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.
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
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.
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.
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), θ
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:
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") in R.
cat("\nExamples are included in the package vignette.\n")
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(PSM)
Loading required package: MASS
Loading required package: numDeriv
Loading required package: deSolve
Attaching package: 'deSolve'
The following object is masked from 'package:graphics':
matplot
Loading required package: ucminf
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/PSM/PSM.estimate.Rd_%03d_medium.png", width=480, height=480)
> ### Name: PSM.estimate
> ### Title: Estimate population parameters
> ### Aliases: PSM.estimate
> ### Keywords: htest models multivariate ts
>
> ### ** Examples
>
> cat("\nExamples are included in the package vignette.\n")
Examples are included in the package vignette.
>
>
>
>
>
> dev.off()
null device
1
>