PSM.estimate: Estimate population parameters

Description

Estimates population parameters in a linear or non-linear mixed effects model based on stochastic differential equations by use of maximum likelihood and the Kalman filter.

Usage

PSM.estimate(Model, Data, Par, CI = FALSE, trace = 0, control=NULL, fast=TRUE)

Arguments

Model

A list containing the following elements:

Matrices = function(phi)

Only in linear models.

Defines the matrices $A$, $B$, $C$ and $D$ in the model equation. Must return a list of matrices named matA, matB, matC and matD. If there is no input, matB and matD may be omitted by setting them to NULL. Note, if the matrix $A$ is singular the option fast is set to FALSE, as this is not supported in the compiled Fortran code.

Functions

Only in non-linear models.

A list containing the functions f(x,u,time,phi), g(x,u,time,phi), df(x,u,time,phi) and dg(x,u,time,phi).

The functions f and g defines the system and df and dg are the Jacobian matrices with first-order partial derivatives for f(x) and g(x) which is needed to evaluate the model. A warning is issued if df or dg appear to be incorrect based on a numerical evaluation of the Jacobians of f(x) and g(x).

It is possible to avoid specifying the Jacobian functions in the model and use numerical approximations instead, but this will increase estimation time at least ten-fold. See the section ‘Numerical Jacobians of f and g’ below for more information.

X0 = function(Time, phi, U)

Defines the model state at Time[1] before update. Time[1] and U[,1] can be used in the evaluation of X0. Must return a column matrix.

SIG = function(phi)

in linear models and SIG = function(u,time,phi) in non-linear models. It defines the matrix $\sigma$ for the diffusion term. Returns a square matrix.

S = function(phi)

in linear models and S = function(u,time,phi) in non-linear models. It defines a covariance matrix for the observation noise. Returns a square matrix.

h = function(eta,theta,covar)

Second stage model. Defines how random effects (eta) and covariates (covar) affects the fixed effects parameters (theta). In models where OMEGA=NULL (no random-effects) h must still be defined with the same argument list to allow for covariates to affect theta, but the function h is evaluated with eta=NULL. Must return a list (or vector) phi of individual parameters which is used as input argument in the other user-defined functions.

ModelPar = function(THETA)

Defines the population parameters to be optimized. Returns a list containing 2 elements, named:

theta: A list of fixed effects parameters $\theta$ which are used as input to the function h listed above.
OMEGA: A square covariance matrix $\Omega$ for the random effects. If OMEGA is missing or NULL then no 2nd stage model is used. However, the function h must still be defined, see above.

Data

An unnamed list where each element contains data for one individual. Each element in Data is a list containing:

Time

A vector of timepoints for measurements

Y

A matrix of multivariate observations for each timepoint, where each column is a multivariate measurement. Y may contain NA for missing observations and a column may consist of both some or only NAs. The latter is useful if a dose is given when no measurement is taken.

U

A matrix of multivariate input to the model for each timepoint. U is assumed constant between measurements and may not contain any NA. If U is ommitted, the model is assumed to have no input and matB and matD need no to be specified.

Dose

A list containing the 3 elements listed below. If the element Dose is missing or NULL, no dose is assumed.

Time: A vector of timepoints for the dosing. Each must coinside with a measurement time. Remember to insert a missing measurement in Y if a corresponding timepoint is not present. Dose is considered added to the system just after the measurement.
State: A vector with indexes of the state for dosing.
Amount: A vector of amounts to be added.

Par

A list containing the following elements:

Init: A vector with initial estimates for THETA, vector of population parameters to be optimized.
LB, UB: : Two vectors with lower and upper bounds for parameters. If ommitted, the program performs unconstrained optimization. It is highly recommended to specify bounds to ensure robust optimization.

Boolean. If true, the program estimates 95% confidence intervals, standard deviation and correlation matrix for the parameter estimates based on the Hessian of the likelihood function. The Hessian is estimated by hessian in the numDeriv package.

trace

Non-negative integer. If positive, tracing information on the progress of the optimization is produced. Higher values produces more tracing information.

control

A list of control parameters for the optimization of the likelihood function. The list has one required component, namely:

optimizer: A string value equal to either 'optim' or 'ucminf'. This gives the choise of optimizer. Default is optimizer = 'optim'.

The remaining components in the list are given as the control argument for the chosen optimizer. See corresponding help file for further detail.

fast

Boolean. Use compiled Fortran code for faster estimation.

Value

A list containing the following elements:

NegLogL

Value of the negative log-likelihood function at optimum.

THETA

Population parameters at optimum

95% confidence interval for the estimated parameters

Standard deviation for the estimated parameters

COR

Correlation matrix for the estimated parameters

sec

Time for the estimation in seconds

opt

Raw output from optim

Numerical Jacobians of f and g

Automatic numerical approximations of the Jacobians of f and g can be used in PSM. In the folliwing, the name of the model object is assumed to be MyModel.

First define the functions MyModel$Functions$f and MyModel$Functions$g. When these are defined in MyModel the functions df and dg can be added to the model object by writing as below:

  MyModel$Functions$df = function(x,u,time,phi) {
    jacobian(MyModel$Functions$f,x=x,u=u,time=time,phi=phi)
  }
  MyModel$Functions$dg = function(x,u,time,phi) {
    jacobian(MyModel$Functions$g,x=x,u=u,time=time,phi=phi)
  }

This way of defining df and dg forces a numerical evaluation of the Jacobians using the numDeriv package. It may be usefull in some cases, but it should be stressed that it will probably give at least a ten-fold increase in estimation times.

Details

The first stage model describing intra-individual variations is for linear models defined as

$$ dx_t = (A(\phi_i)x_t + B(\phi_i)u_t)dt + \sigma(\phi_i) d\omega_t$$ $$ y_{ij} = C(\phi_i)x_{ij} + D(\phi_i)u_{ij} + e_{ij}$$

and for non-linear models as

$$ dx_t = f(x_t,u_t,t,\phi_i)dt + \sigma(u_t,t,\phi_i) d\omega_t$$ $$y_{ij} = g(x_{ij},u_{ij},t_{ij},\phi_i) + e_{ij}$$

where $e_{ij} \sim N(0,S(u_{ij},t_{ij},\phi_i))$ and $\omega_t$ is a standard Brownian motion.

The second stage model describing inter-individual variations is defined as:

$$\phi_i = h(\eta_i,\theta,Z_i)$$

where $\eta_i \sim N(0,\Omega)$, $\theta$ are the fixed effect parameters and $Z_i$ are covariates for individual i. In a model without random-effects the function $h$ is only used to include possible covariates in the model.

References

Please visit http://www.imm.dtu.dk/psm or refer to the main help page for PSM.

Examples

Run this code

# NOT RUN {
#detailed examples are provided in the package vignette

#Theophylline data from Boeckmann et al (1994)
#objective: recover the administered doses

library(datasets)
data(Theoph) 


#reshape data to PSM format

TheophPSM = list()
for(i in 1:length(unique(Theoph$Subject))){
  TheophPSM[[i]] = with(
    Theoph[Theoph$Subject == i,],
    list(Y = matrix(conc, nrow=1), Time = Time)
  )
}


#specify a simple pharmacokinetic model comprised of
#2 state equations and 1 observation equation
#initial value of 1 state eq. varies randomly across individuals

mod = vector(mode="list")
mod$Matrices = function(phi) {
  list(
    matA=matrix(c(-phi$ka, 0, phi$ka, -phi$ke), nrow=2, ncol=2, byrow=TRUE),
    matC=matrix(c(0, 1), nrow=1, ncol=2)
  )
}
mod$h = function(eta, theta, covar) {
  phi = theta
  phi$dose = theta$dose * exp(eta[1])
  phi
}
mod$S = function(phi) {
  matrix(c(phi$sigma), nrow=1, ncol=1)
}
mod$SIG = function(phi) {
  matrix(c(0, 0, 0, phi$omega), nrow=2, ncol=2, byrow=TRUE)
}
mod$X0 = function(Time, phi, U) {
  matrix(c(phi$dose, 0), nrow=2, ncol=1)
}
mod$ModelPar = function(THETA) {
  list(theta=list(dose = THETA["dose"], ka = THETA["ka"], ke = THETA["ke"],
                  omega = THETA["omega"], sigma = THETA["sigma"]),
       OMEGA=matrix(c(THETA["BSV_dose"]), nrow=1, ncol=1)
  )
}


#specify the search space of the fitting algorithm

parM = c(ka = 1.5, ke = 0.1, dose = 10, omega = .3, sigma = 1,
         BSV_dose = 0.015)
pars = list(LB=parM*.25, Init=parM, UB=parM*1.75)

#fit model and predict data

fit = PSM.estimate(mod, TheophPSM, pars, trace = 1, fast = TRUE,
  control=list(optimizer="optim", maxit=1))
pred = PSM.smooth(mod, TheophPSM, fit$THETA)

#visualize recovery performance

true_dose = tapply(Theoph$conc, Theoph$Subject, mean)
true_dose = true_dose[order(as.numeric(names(true_dose)))]
est_dose = fit$THETA["dose"] * exp(unlist(lapply(pred, function(x) x$eta)))
plot(true_dose, est_dose,
  xlab="actually administered dose", ylab= "recovered dose")
abline(lm(est_dose ~ true_dose), lty=2)

# }

Run the code above in your browser using DataLab

Description

Usage

Arguments

Value

Numerical Jacobians of f and g

Details

References

See Also

Examples