hoa (version 2.1.4.1)

tem: Tangent exponential model: Higher Order Likelihood Approximation

Description

This function computes the tangent exponential model approximation for higher order likelihood inference about a scalar interest parameter of a parametric model. The function creates an object of class fr.

Usage

tem(psi = NULL, nlogL, phi, make.V, th.init, data, 
              tol = 10^(-5), n.psi = 50)

Arguments

psi

A scalar value for the interest parameter. If NULL (the default) the vector is taken as a grid of values centred on the maximum likelihood estimate (MLE).

nlogL

A function to compute the negative log likelihood for the model of interest. It is a function of three quantities: psi is the scalar interest parameter, lam is the scalar or vector nuisance parameter, and data is the data. See below for an example.

phi

A function to compute the canonical parameter of a local exponential family approximation to the density of interest. It requires values of the parameter theta = (psi,lam), the quantities V, and the data data. The output is a vector of the same length as theta.

make.V

A function to compute V, using as input the parameter theta = (psi,lam) and the data data. The output is a matrix whose rows correspond to individual observations, and whose columns correspond to the elements of theta.

th.init

Initial value(s) of the parameter theta.

data

Data frame or other object containing the data.

tol

Tolerance used for numerical differentiation of phi.

n.psi

Number of values of psi at which the likelihood is computed, if psi is NULL. Avoid odd values, which may give numerical instabilities near the MLE.

Value

normal

The MLE of the interest parameter, and its standard error

th.hat

MLEs of parameters (psi,lam)

th.hat.se

Standard errors of MLEs, based on observed information

th.rest

Restricted MLEs (psi,lam) on grid of values of psi

r

Values of likelihood root corresponding to psi

psi

Values of interest parameter psi

q

Values of likelihood modification

rstar

Values of modified likelihood root

Details

The function computes quantities used for higher order likelihood approximations, which are intended to provide highly accurate inferences on scalar parameters in parametric statistical models. The key aspect is maximisation of the likelihood over a grid of values of the interest parameter psi, and computation of likelihood modifications based on local exponential family approximation to the density. If \(n\) is the sample size, then the resulting inferences should be accurate to order \(n^{-3/2}\) in continuous models and to order \(n^{-1}\) in discrete models, and in many cases they are very close to exact results. The approximations rely on numerical computation of observed information matrices and of derivatives, and may fail in certain cases. The confidence intervals themselves and useful plots are produced using the functions summary and plot. For technical background and further details, see Sections 2.4 and 8.4.2 of the book cited below, which has many further references.

References

Brazzale, A. R., Davison, A. C. and Reid, N. (2007). Applied Asymptotics: Case Studies in Small-Sample Statistics. Cambridge University Press, Cambridge.

See also http://statwww.epfl.ch/AA.

See Also

plot.fr, summary.fr, lik.ci

Examples

# NOT RUN {
# Cost data example from Section 3.5 of "Applied Asymptotics"	
cost <- data.frame(
          f = factor( c(rep(1,13), rep(2,18)) ),
          y = c( 30,172,210,212,335,489,651,1263,1294,1875,2213,2998, 
                 4935,121,172,201,214,228,261,278,279,351,561,622,
                 694,848,853,1086,1110,1243,2543 ) )  
nlogL <- function(psi, lam, data)  { 
  s1 <- exp(lam[2])
  m2 <- lam[1]
  s2 <- exp(lam[3])
  m1 <- psi + m2 + s2^2/2 - s1^2/2
  -sum( dnorm(log(data$y), mean=ifelse(data$f==1, m1, m2),
  	                sd=ifelse(data$f==1, s1, s2), log=TRUE) )
}
phi <- function(th, V, data)  {
  psi <- th[1]
  lam <- th[-1]
  s1 <- exp(lam[2])
  m2 <- lam[1]
  s2 <- exp(lam[3])
  m1 <- psi + m2 + s2^2/2 - s1^2/2
  c( m1/s1^2, 1/s1^2, m2/s2^2, 1/s2^2 ) 
}
make.V <- function(th, data) NULL
cost.lnorm.rat <- tem(psi = NULL, nlogL = nlogL, phi = phi, 
                    make.V = make.V, th.init = c(0, 5, 2, 5), data = cost)
plot(cost.lnorm.rat, psi = 0, all = TRUE)
summary(cost.lnorm.rat)
# }