profile: Profile-likelihood (PL) computation
In dMod: Dynamic Modeling and Parameter Estimation in ODE Models

Description Usage Arguments Details Value Examples

Profile-likelihood (PL) computation

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,
  ...
)

`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()

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.

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).

## 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")


## End(Not run)

Loading required package: cOde

Attaching package: 'dMod'

The following object is masked from 'package:stats':

    profile

The following local files were dynamically loaded: 
The following local files were dynamically loaded: 
The following local files were dynamically loaded: 
The following local files were dynamically loaded: 
The following local files were dynamically loaded: 
The following local files were dynamically loaded: 
The following local files were dynamically loaded: 
The following local files were dynamically loaded: 
     name         value    lower     upper
loga loga -1.917574e+00     -Inf 0.2485995
logb logb  5.069280e-17     -Inf 0.7769882
logc logc  4.220159e+00 1.716189       Inf