logPL: Profile likelihood computation
In likelihoodAsy: Functions for Likelihood Asymptotics

Description Usage Arguments Details Value References Examples

View source: R/likelihoodAsy.R

This function evaluates the profile likelihood for a subset of the model parameter. The result is optionally returned with a minus sign, so the function can be used directly as input to a general-purpose optimizer.

1 2	logPL(psival, data, thetainit, floglik, fscore=NULL, indpsi, minus=FALSE, onestep=FALSE, jhat=NULL, trace=FALSE)

`psival`	A numerical vector containing the value of the parameter of interest.
`data`	The data as a list. All the elements required to compute the likelihood function at a given parameter value should be included in this list. The required format of such list will be determined by the user-provided function `floglik`.
`thetainit`	A numerical vector with the size of the entire model parameter, that will be used as starting point in the constrained optimization performed to obtain the maximum likelihood estimate under the null. The specific meaning of `thetainit` is determined by the specification of `floglik`.
`floglik`	A function which returns the log likelihood function at a given parameter value. In particular, for a certain parameter value contained in a numerical vector `theta`, a call `floglik(theta, data)` should return a scalar numerical value, the log likelihood function at `theta`. Note that the parameter of interest should be a subset of the coordinates of `theta`.
`fscore`	An optional function which returns the score function at a given parameter value. It must return a numerical vector of the same length as `thetainit`. For a certain parameter value contained in a numerical vector `theta`, a call `fscore(theta, data)` should return the gradient of the log likelihood function at `theta`. Default is `NULL`, implying that numerical differentiation will be employed.
`indpsi`	A vector of integers in the range `1:length(theta)` containing the indexes of the parameter of interest, so that the parameter of interest will be given by `theta[indpsi]`.
`minus`	Logical. Should the profile likelihood be multiplied by -1? This may be useful for usage with optimizers. Default is `FALSE`.
`onestep`	Logical. If set to `TRUE` the constrained estimate of the nuisance parameter is replaced by a one-step approximation around the maximum likelihood estimate. Default is `FALSE`.
`jhat`	A squared matrix with dimension equal to `length(mle)` containing the observed information matrix evaluated at `mle`. It is employed only when `onestep=TRUE`. Default is `NULL`.
`trace`	Logical. When set to `TRUE` will cause the printing of the MPL value, which can be useful to monitor optimization. Default is `FALSE`.

This function is designed to be used with external functions, such as optimizers and evaluators over a grid of points.

A scalar value, minus the profile likelihood at psival.

Severini, T.A. (2000). Likelihood Methods in Statistics. Oxford University Press.

# A negative binomial example, taken from Venables and Ripley (2002, MASS4 book)
library(MASS)
# The quine data are analysed in Section 7.4
data(quine)
# We fit a model with just the main effects
quine.nb1 <- glm.nb(Days ~ Eth + Sex + Age + Lrn, data = quine) 
# The data list includes the design matrix and the response vector
quinedata<-list(X=model.matrix(quine.nb1), y=quine$Days)      
# Let us define the various functions
# Log likelihood, log link
logLikNbin <- function(theta,data) 
{
  y <- data$y
  X <- data$X
  eta <- X %*% theta[1:ncol(X)] 
  mu <- exp(eta)
  alpha <- theta[ncol(X)+1]
  l <- sum(lgamma(y + alpha) + y * log(mu) - (alpha + y) * log(alpha + mu) 
            - lgamma(alpha) + alpha * log(alpha))
  return(l)
}


# Score function
gradLikNbin <- function(theta,data) 
{
  y <- data$y
  X <- data$X
  eta <- X %*% theta[1:ncol(X)] 
  mu <- exp(eta)
  alpha <- theta[ncol(X)+1]
  g <-rep(0,ncol(X)+1)
  g[1:ncol(X)] <- t(y - (alpha+y)*mu / (alpha+mu)) %*% X
  g[ncol(X)+1] <- sum(digamma( y + alpha) - log(alpha + mu) - (alpha + y) / (alpha + mu) 
                  - digamma(alpha) + 1 + log(alpha))
  return(g)
}
# Data generator
genDataNbin<- function(theta,data)
{
  out <- data
  X <- data$X
  eta<- X %*% theta[1:ncol(X)] 
  mu <- exp(eta)
  out$y <- rnegbin(length(data$y), mu=mu, theta=theta[ncol(X)+1])
  return(out)
}		
# First we refine the maximum likelihood estimates
mleFull <- optim( c(coef(quine.nb1),quine.nb1$theta), logLikNbin, gr=gradLikNbin,
           method="BFGS", data=quinedata, control=list(fnscale=-1), hessian=TRUE) 
 # Then we can plot the profile likelihood
list.psi <- seq(0.90, 1.70, l=30)
list.prof <- sapply(list.psi, logPL, data=quinedata, thetainit=mleFull$par, floglik=logLikNbin, 
                    fscore=gradLikNbin, indpsi=8, trace=FALSE) 
plot(list.psi, list.prof-max(list.prof), type="l", xlab=expression(psi), ylab="Log likelihood")