dMod (version 0.1)

profile: Profile-likelihood (PL) computation

Description

Profile-likelihood (PL) computation

Usage

profile(obj, pars, whichPar, alpha = 0.05, limits = c(lower = -Inf, upper =
  Inf), stepControl = NULL, algoControl = NULL, optControl = NULL,
  verbose = FALSE, ...)

Arguments

obj

Objective function obj(pars, fixed, ...) returning a list with "value", "gradient" and "hessian".

pars

Parameter vector corresponding to the log-liklihood optimum.

whichPar

Numeric or character. The parameter 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]

stepControl

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

algoControl

List of arguments controlling the fast PL algorithm. defaults to list(gamma = 1, W = c("hessian", "identity"), reoptimize = FALSE, correction = 1, reg = 1e-6)

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 mor

verbose

Logical, print verbose messages.

...

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), 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 matrix reg*Id is added to the Hessian.

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.

Examples

## ----------------------
## Example 1
## ----------------------
trafo <- c(a = "exp(loga)", b = "exp(logb)",c = "exp(loga)*exp(logb)*exp(logc)")
p <- P(trafo)
obj <- function(pOuter, fixed = NULL)
   constraintL2(p(pOuter, fixed), c(a =.1, b = 1, c = 10), 1)

ini <- c(loga = 1, logb = 1, logc = 1)
myfit <- trust(obj, ini, rinit=1, rmax=10)
profiles <- sapply(1:3, function(i)
   profile(obj, myfit$argument, whichPar = i, limits = c(-5, 5),
                 algoControl=list(gamma=1, reoptimize=FALSE), verbose=TRUE))
plotProfile(profiles)
plotPaths(profiles)

## ----------------------------
## Example 2
## ----------------------------
trafo <- c(a = "exp(loga)", b = "exp(logb)",c = "exp(loga)*exp(logb)*exp(logc)")
p <- P(trafo)
obj <- function(pOuter, fixed = NULL, sigma)
  constraintL2(p(pOuter, fixed), c(a =.1, b = 1, c = 10), 1) +
  constraintL2(pOuter, mu = c(loga = 0, logb = 0), sigma = sigma, fixed = fixed)


ini <- c(loga = 1, logb = 1, logc = 1)
myfit <- trust(obj, ini[-1], rinit=1, rmax=10, fixed = ini[1], sigma = 10)
profiles.approx <- sapply(1:2, function(i)
  profile(obj, myfit$argument, whichPar = i, limits = c(-10, 10),
                algoControl=list(gamma=1, reoptimize=FALSE),
                verbose=TRUE, fixed = ini[1], sigma = 10))
profiles.exact  <- sapply(1:2, function(i)
  profile(obj, myfit$argument, whichPar = i, limits = c(-10, 10),
                algoControl=list(gamma=0, reoptimize=TRUE),
                verbose=TRUE, fixed = ini[1], sigma = 10))

plotProfile(profiles.approx, profiles.exact)