R: Estimate population parameters

PSM.estimate {PSM}

R Documentation

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 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 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 `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`	The list is passed to `"optim"` to control the settings for the optimization. As default the max iteration limit is set to 100 and the remaining options are the same as in `"optim"`.
`fast`	Boolean. Use compiled Fortran code for faster estimation.

Details

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

Value

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`

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.

Note

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

Author(s)

Stig B. Mortensen and Søren Klim

References

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

Examples

cat("\nExamples are included in the package vignette.\n")

[Package PSM version 0.8-3 Index]