tem: Tangent exponential model: Higher Order Likelihood...

In hoa: Higher Order Likelihood Inference

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.

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.

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

Author(s)

Anthony Davison <[e-mail]> Alex-Antoine Fortin <[e-mail]>

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

# 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)

hoa documentation built on May 2, 2019, 8:56 a.m.