dMod (version 1.0.0)

profile: Profile-likelihood (PL) computation

Description

Profile-likelihood (PL) computation

Usage

profile(obj, pars, whichPar, alpha = 0.05, limits = c(lower = -Inf,
  upper = Inf), method = c("integrate", "optimize"),
  stepControl = NULL, algoControl = NULL, optControl = NULL,
  verbose = FALSE, cores = 1, ...)

Arguments

obj

Objective function obj(pars, fixed, ...) returning a list with "value", "gradient" and "hessian". If attribute "valueData" and/or "valuePrior are returned they are attached to the return value.

pars

Parameter vector corresponding to the log-liklihood optimum.

whichPar

Numeric or character vector. The parameters for which the profile is computed.

alpha

Numeric, the significance level based on the chisquare distribution with df=1

limits

Numeric vector of length 2, the lower and upper deviance from the original value of pars[whichPar]

method

Character, either "integrate" or "optimize". This is a short-cut for setting stepControl, algoControl and optControl by hand.

stepControl

List of arguments controlling the step adaption. Defaults to integration set-up, i.e. list(stepsize = 1e-4, min = 1e-4, max = Inf, atol = 1e-2, rtol = 1e-2, limit = 100)

algoControl

List of arguments controlling the fast PL algorithm. defaults to list(gamma = 1, W = "hessian", reoptimize = FALSE, correction = 1, reg = .Machine$double.eps)

optControl

List of arguments controlling the trust() optimizer. Defaults to list(rinit = .1, rmax = 10, iterlim = 10, fterm = sqrt(.Machine$double.eps), mterm = sqrt(.Machine$double.eps)). See trust for more details.

verbose

Logical, print verbose messages.

cores

number of cores used when computing profiles for several parameters.

...

Arguments going to obj()

Value

Named list of length one. The name is the parameter name. The list enty is a matrix with columns "value" (the objective value), "constraint" (deviation of the profiled paramter from the original value), "stepsize" (the stepsize take for the iteration), "gamma" (the gamma value employed for the iteration), "valueData" and "valuePrior" (if specified in obj), one column per parameter (the profile paths).

Details

Computation of the profile likelihood is based on the method of Lagrangian multipliers and Euler integration of the corresponding differential equation of the profile likelihood paths.

algoControl: Since the Hessian which is needed for the differential equation is frequently misspecified, the error in integration needs to be compensated by a correction factor gamma. Instead of the Hessian, an identity matrix can be used. To guarantee that the profile likelihood path stays on the true path, each point proposed by the differential equation can be used as starting point for an optimization run when reoptimize = TRUE. The correction factor gamma is adapted based on the amount of actual correction. If this exceeds the value correction, gamma is reduced. In some cases, the Hessian becomes singular. This leads to problems when inverting the Hessian. To avoid this problem, the pseudoinverse is computed by removing all singular values lower than reg.

stepControl: The Euler integration starts with stepsize. In each step the predicted change of the objective function is compared with the actual change. If this is larger than atol, the stepsize is reduced. For small deviations, either compared the abolute tolerance atol or the relative tolerance rtol, the stepsize may be increased. max and min are upper and lower bounds for stepsize. limit is the maximum number of steps that are take for the profile computation. stop is a character, usually "value" or "data", for which the significance level alpha is evaluated.

Examples

# NOT RUN {
## Parameter transformation
trafo <- eqnvec(a = "exp(loga)", 
                b = "exp(logb)", 
                c = "exp(loga)*exp(logb)*exp(logc)")
p <- P(trafo)

## Objective function
obj1 <- constraintL2(mu = c(a = .1, b = 1, c = 10), sigma = .6)
obj2 <- constraintL2(mu = c(loga = 0, logb = 0), sigma = 10)
obj <- obj1*p + obj2

## Initialize parameters and obtain fit
pars <- c(loga = 1, logb = 1, logc = 1)
myfit <- trust(obj, pars, rinit = 1, rmax = 10)
myfit.fixed <- trust(obj, pars[-1], rinit = 1, rmax = 10, fixed = pars[1])

## Compute profiles by integration method
profiles.approx <- do.call(
  rbind, 
  lapply(1:3, function(i) {
    profile(obj, myfit$argument, whichPar = i, limits = c(-3, 3),
            method = "integrate")
  })
)

## Compute profiles by repeated optimization 
profiles.exact <- do.call(
  rbind, 
  lapply(1:3, function(i) {
    profile(obj, myfit$argument, whichPar = i, limits = c(-3, 3),
            method = "optimize")
  })
)

## Compute profiles for fit with fixed element by integration method
profiles.approx.fixed <- do.call(
  rbind, 
  lapply(1:2, function(i) {
    profile(obj, myfit.fixed$argument, whichPar = i, limits = c(-3, 3),
            method = "integrate",
            fixed = pars[1])
  })
)

## Plotting
plotProfile(profiles.approx)
plotProfile(list(profiles.approx, profiles.exact))
plotProfile(list(profiles.approx, profiles.approx.fixed))

plotPaths(profiles.approx, sort = TRUE)
plotPaths(profiles.approx, whichPar = "logc")
plotPaths(list(profiles.approx, profiles.approx.fixed), whichPar = "logc")

## Confidence Intervals
confint(profiles.approx, val.column = "value")

# }