PSM.estimate: Estimate population parameters
In PSM: Non-Linear Mixed-Effects Modelling using Stochastic Differential Equations

Description Usage Arguments Details Value Numerical Jacobians of f and g Note Author(s) References See Also Examples

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.

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

`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 SIG 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 θ 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 `NA`s. 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.
`CI`	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.

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

dx = (A(phi)*x + B(phi)*u)dt + SIG(phi)*dw

y = C(phi)*x + D(phi)*u + e

and for non-linear models as

dx = f(x,u,t,phi)dt + SIG(u,t,phi)dw

y = g(x, u, t, phi) + e

where e ~ N(0,S(x, u, t)) and w is a standard Brownian motion.

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

phi = h(eta,theta,Z)

where eta ~ N(0,OMEGA), θ are the fixed effect parameters and Z are covariates for individual i. In a model without random-effects the function h is only used to include possible covariates in the model.

A list containing the following elements:

`NegLogL`	Value of the negative log-likelihood function at optimum.
`THETA`	Population parameters at optimum
`CI`	95% confidence interval for the estimated parameters
`SD`	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`

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.

For further details please also read the package vignette pdf-document by writing vignette("PSM") in R.

Stig B. Mortensen and Soeren Klim

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

PSM, PSM.smooth, PSM.simulate, PSM.plot, PSM.template

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