Computes a likelihood function for a scalar parameter of interest using the modified profile likelihood (MPL) of Barndorff-Nielsen, as implemented by likelihoodAsy::logMPL(). The ordinary profile likelihood - which curve_lik() constructs via the ProfileLikelihood package - treats nuisance parameter estimates as known, and can be noticeably biased when nuisance parameters are numerous relative to the sample size (e.g., stratified models, variance components). The MPL adjusts for nuisance parameter estimation and is generally preferable in those settings; for variance parameters in mixed models it closely tracks REML.

curve_mpl(data, mle, floglik, datagen, indpsi, lo, hi, steps = 50,
  fscore = NULL, R = 100, seed = NULL, table = TRUE, ...)

Arguments

data

A list containing the data, in the format expected by the user-supplied floglik and datagen functions (see the likelihoodAsy documentation).

mle

Numeric vector with the full maximum likelihood estimate of the model parameter.

floglik

A function of (theta, data) returning the log likelihood at theta.

datagen

A function of (theta, data) returning a copy of data with the response simulated from the model at theta.

indpsi

Integer index of the parameter of interest within theta.

lo

Lower limit of the grid of values for the parameter of interest over which the modified profile likelihood is evaluated.

hi

Upper limit of the grid of values for the parameter of interest.

steps

Number of grid points at which the modified profile likelihood is evaluated. The default is 50; more points give a smoother function at proportionally greater computational cost.

fscore

Optional function of (theta, data) returning the gradient of the log likelihood; supplying it speeds up computation.

R

Number of Monte Carlo draws used for the Skovgaard-type covariance approximations inside likelihoodAsy::logMPL(). The default is 100, which suffices for (curved) exponential family models; increase for stability in other models.

seed

Optional seed for the Monte Carlo computation, for reproducibility.

table

Indicates whether or not a table output with some relevant statistics should be generated. The default is TRUE and generates a table which is included in the list object.

...

Further arguments passed to likelihoodAsy::logMPL().

Value

A list with 2 items where the dataframe of values is in the first object, and the table for the values in the second if table = TRUE.

Details

The user must supply the same ingredients required by likelihoodAsy: a function evaluating the log likelihood, a function that simulates a data set from the fitted model, the full maximum likelihood estimate, and the index of the parameter of interest within the parameter vector.

Because the MPL, like any likelihood, is defined only up to a multiplicative constant, the returned likelihood column is normalized so that its maximum is 1 and is therefore identical to the support column.

References

Barndorff-Nielsen OE. On a formula for the distribution of the maximum likelihood estimator. Biometrika. 1983;70:343-365.

Pierce DA, Bellio R. Modern likelihood-frequentist inference. International Statistical Review. 2017;85:519-541.

Examples

if (FALSE) { # \dontrun{
# Logistic regression on the crying babies data (cond package),
# parameter of interest = coefficient of lull (index 19)
lik <- curve_mpl(
  data = data.obj, mle = coef(mod.glm),
  floglik = loglik.logit, datagen = gendat.logit,
  indpsi = 19, lo = -0.3, hi = 3.7, seed = 2020
)
ggcurve(lik[[1]], type = "l1")
} # }