Profile-Likelihood Inference for Fitted Parametric Model Objects

Profile-likelihood inference for fitted model objects.

Usage

# S3 method for default
profLik(object,
        fun,
        level = 0.95,
        deriv = TRUE,
        trace = 0,
        diagno = TRUE,
        ftol_abs = 1e-12, ftol_rel = 1e-8,
        ...)

Arguments

object: A fitted model object. See Details.
fun: A function of the parameter vector for which the profile-likelihood will be carried over. This function must have the arguments: theta for the vector of parameters and object for the model object; so the function can use the velue of some of slots of object see Details. The function must return a list with two elements with names "objective" and "gradient" as required by nloptr see Examples. If needed, a wrapper function can be used to use more arguments.
level: Level of confidence. Can be of length > 1.
deriv: Logical. If TRUE, the function fun is assumed to provide a gradient vector as an attribute named "gradient" of the result. For now deriv can only be TRUE, which implies that fun must compute the gradient.
trace: Level of verbosity; the value 0 prints nothing.
diagno: Logical. When TRUE two diagnostics are returned as two attributes "diagno" and "theta" of the returned array. The array in attribute "diagno" contains the return codes of nloptr (named "status"), the values of the objective and of the constraint functions as well as a value named "gradDist" whcih measures the distance between the two directions of the two gradients: objective and constraint. The array in attribute "theta" contains the value of the Poisson-GP parameter that was found to maxi/minimise the function fun under the constraint of a high log-likelihood. See section Note.
ftol_abs, ftol_rel: Absolute and relative tolerance to stop the constrained optimisation nloptr. These apply to the objective of the constrained optimisation that is to the value of fun. Remind that ftol_abs is thus given with the same unit as fun.
...: Not used yet.

Value

An array with the value of the function and the corresponding Lower and Upper end-points for the given confidence levels. This array has two attributes with names "diagno"

and "psi" which both are arrays. The attributes provide information about the numerical optimisation and the values of the vector of parameter that maximised or minimised the function fun.

Details

Compute the lower and upper end-points of a profile-likelihood based confidence interval. The (apparently new) method used here relies on maximising and minimising the function of interest, say $\eta(\boldsymbol(\theta)$, under the constraint that the log-likelihood is greater than the maximal log-likelihood minus a positive quantity $\delta$ depending on the confidence level. This differs from the usual method which relies on an univariate zero-finding for the profile-likelihood function (minus a constant). Remind that each evaluation of the profile requires a $p-1$ dimensional optimisation. As a major advantage, the new method does not require a re-parameterisation of the model.

The requirement for the method to work is that object has suitable "slots", i.e. is a list with suitable elements. These are: the number p of parameters, the vector parNames of the parameter names, the vector estimate of ML estimates, the closure negLogLikFun computing the negative log-likelihood function and the value negLogLik of the minimised negative log-likelihood. These requirements are fulfilled for poisGP objects of the package. The result returned by object$negLogLikFun should have a "gradient" attribute to allow the use of derivatives.

Note

For each confidence limit the numerical optimisation may fail, in which case the limit will be NA. Using trace = 1 can be useful to further check the optimisation. We encourage the user to inspect the set of important diagnostics returned when diagno = TRUE, especially those in the attribute of the result named "diagno". The optimisation status value in status should be between 1 and 4 to indicate that small changes on the parameter or on the objective were eventually obtained. On the other hand, a value of 5 for status indicates that the maximal number of iterations was reached, which is considered here as a failure. Both the constraint and gradDist values should be small because the constraint must be active at the optimum and the two gradients for the objective and the constraint must be colinear at the optimum (Lagrange conditions). In other words the optimum parameter vector must lie on the boundary of the high-likelihood region corresponding to the chosen confidence level.

References

Deville Y. (2017) "Profile-Likelihood Using Constrained Optimisation". Unpublished Tech. Report.

Johnson S.G. The NLopt Nonlinear-Optimization Package. https://github.com/stevengj/nlopt.

Section Return values in the manual NLOPT Reference.

Author

Yves Deville

Examples

object <- poisGP(Garonne)
#> Warning: 'threshold' is missing and set just below the smallest obsevation
#> Warning: 'threshold' is smaller than the smallest observation
#> Warning: 'threshold' is smaller than the smallest observation

## =========================================================================
## Define a function of the parameter vector: here the first component of
## the "PP" parameter vector.
## This is for code illustration only, since the the result can be obtained
## using the 'confint' method with 'method = "proflik"', which
## gives the profile-likelihood confidence intervals for each of the
## three parameters.
## =========================================================================

numPP <- 2
myfun <- function(theta, object) {
    thetaStar <- nieve::poisGP2PP(lambda = theta[1],
                                  scale = theta[2],
                                  shape = theta[3],
                                  loc = object$threshold, deriv = TRUE)
    res <- thetaStar[numPP]
    grad <- attr(thetaStar, "gradient")[1, numPP, c(1, 3, 4)]
    attr(res, "gradient") <- grad
    res
}

pl <- profLik(object = object, fun = myfun, level = c(0.95, 0.70))