\documentclass{article}
\usepackage{amsmath,amssymb,graphics}
\usepackage{color}
\usepackage{booktabs}
\newcommand{\bs}{\boldsymbol}
% \VignetteIndexEntry{Population Stochastic Modelling: Model definition, description and examples}
% \VignetteDepends{PSM}
% \usepackage{boxedminipage}

% Sweave hack by inserting %\usepackage{Sweave} then Sweave does not insert it's own line
%\usepackage{Sweave}
<<Sweavepath,include=FALSE,echo=FALSE>>=
sweavepath <- gsub("[\\]",'/',R.home("share/texmf/"))
@
%\usepackage{\Sexpr{sweavepath}Sweave} % not used - inserted below instead.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN Sweave.sty
\NeedsTeXFormat{LaTeX2e}
\ProvidesPackage{Sweave}{}

\RequirePackage{ifthen}
\newboolean{Sweave@gin}
\setboolean{Sweave@gin}{true}
\newboolean{Sweave@ae}
\setboolean{Sweave@ae}{true}

\DeclareOption{nogin}{\setboolean{Sweave@gin}{false}}
\DeclareOption{noae}{\setboolean{Sweave@ae}{false}}
\ProcessOptions

\RequirePackage{graphicx,fancyvrb}

\ifthenelse{\boolean{Sweave@gin}}{\setkeys{Gin}{width=0.8\textwidth}}{}%
\ifthenelse{\boolean{Sweave@ae}}{%
  \RequirePackage[T1]{fontenc}
  \RequirePackage{ae}
}{}%

\DefineVerbatimEnvironment{Sinput}{Verbatim}{fontshape=sl}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{}
\DefineVerbatimEnvironment{Scode}{Verbatim}{fontshape=sl}

\newenvironment{Schunk}{}{}
% END Sweave.sty
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%
% look and feel of the R code
%%%
\DefineVerbatimEnvironment{Sinput}{Verbatim}{fontsize=\small,fontshape=sl}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\small}
\DefineVerbatimEnvironment{Scode}{Verbatim}{fontsize=\small,fontshape=sl}
%%%
% put R code and R-code output in a boxed minipage:
%%%
\renewenvironment{Schunk}
{\begin{center} \begin{minipage}{0.95\textwidth} }
{\end{minipage}\end{center}}


\begin{document}

<<PackageDescription,include=FALSE,echo=FALSE>>=
options(keep.source = TRUE, width = 60)
foo <- packageDescription("PSM")
library("PSM")
foo[['Built']] <- gsub("[#$%^&_{}~]",replacement="-", x=foo[['Built']])
options(digits=4)
@


\title{Population Stochastic Modelling (PSM): \\
Model definition, description and examples}

\date{\Sexpr{format(as.Date(foo[['Date']]),"%B %e, %Y")}}

\author{Stig Mortensen and S\o ren Klim}
\maketitle

\begin{center}
%\begin{tabular}{lp{6cm}}
\begin{tabular}{ll}
Package:& \Sexpr{foo[['Package']]}, version \Sexpr{foo[['Version']]}\\
%Built:  & \Sexpr{foo[['Built']]} \\
URL:    & \Sexpr{foo[['URL']]}
\end{tabular}
\end{center}

%\setcounter{tocdepth}{1}
\tableofcontents



\section{Introduction}

\Sexpr{foo[['Description']]}

Some of the most essential parts of the implentation, namely
the Kalman filter, is for linear models run using compiled code written in Fortran,
which gives significant improvements in the parameter estimation times
in R. However, otherwise this version is almost entirely created in
\emph{R}, and estimation times are thus in no way comparable to
state-of-the-art software for similar types of models based on
ordinary differential equations.



\section{Model definition}

A mixed-effects model is used to describe data with the following
general structure

\begin{equation}
\bs y_{ij}, \quad i = 1, ..., N, \quad j = 1, ..., n_i
\end{equation}

\noindent where $\bs y_{ij}$ is a vector of measurements at time $t_{ij}$ for
individual $i$, $N$ is the number of individuals and $n_i$ is the
number of measurements for individual $i$. In a mixed-effects model
the variation is split into intra-individual variation and
inter-individual variation, which is modelled by a first
and second stage model. For further detail regarding the model
definition please refer to \cite{mortensen07, overgaard05, kristensen03}.

\subsubsection*{First stage model}

The first stage model for a mixed effects model can be
written in the form of a state space model. A state space model
consists of two parts, namely a set of continuous state equations
defining the dynamics of the system and a set of discrete measurement
equations, which defines a functional relationship between the states
of the system and the obtained measurements. In the linear form the
state space equations are written as

\begin{eqnarray}
d\bs x_t & = & ( \bs A(\bs\phi_i) \bs x_t + \bs B(\bs\phi_i) \bs u_t)dt +\bs\sigma_\omega(\bs\phi_i)d\bs\omega_t \label{sssm1}\\
\bs y_{ij} & = & \bs C(\bs\phi_i) \bs x_{ij} + \bs D(\bs\phi_i) \bs u_{ij}+\bs e_{ij}
\label{sssm2}
\end{eqnarray}

\noindent and for a general non-linear model as

\begin{eqnarray}
  d\bs x_t &=& \bs f(\bs x_t,\bs u_t,t,\bs\phi_i)dt + \bs\sigma(u_t,t,\phi_i) d\bs\omega_t \\
  \bs y_{ij} &=& \bs g(\bs x_{ij},\bs u_{ij},t_{ij},\bs\phi_i) + \bs e_{ij}
\end{eqnarray}

\noindent where $t$ is the
continuous time variable, the states of the model and the optional
inputs  at time $t$ are denoted $\bs x_t$ and $\bs u_t$
respectively and $\bs \omega_t$ is a standard Wiener process such that
$\bs\omega_{t_2} - \bs\omega_{t_1} \in  N(\bs 0,|t_2-t_1|\bs I)$. Both
the state, measurement and input can be
multi-dimensional, and are in such cases thus represented by a vector
at time $t_{ij}$. The input is assumed constant between sample times
(zero-order hold). The individual model parameters are denoted
$\bs\phi_i$. Measurements are assumed observed with a Gaussian white
noise  measurement error,  that is $\bs e_{ij}\in N(\bs 0,\bs
S(\bs\phi_i))$. For a non-linear model the covariance matrix may also
depend on input and time, that is $S(\bs u_t, t, \bs\phi_i))$.

In the evaluation of the non-linear model it is necessary to specify a
Jacobian matrix function with first-order partial derivatives for $\bs
f$ and $\bs g$. These functions are defined as

\begin{equation}
\frac{\partial \bs f}{\partial \bs x_t} ~~, ~~ \frac{\partial \bs g}{\partial \bs x_t}
\end{equation}

\noindent and must be given with the model specification. PSM will
check the user defined Jacobian functions with numerical evaluations of
the Jacobians of $\bs f$ and $\bs g$ to ensure that they are correct.
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 help file in \verb+R+ for
\verb+PSM.estimate+ for details regarding this.

The initial state of the model is given as a function of $t_1$,
$\bs\phi_i$, and $u_{i1}$ and defines the model state at time $t_1$
before update based on the first observation. The initial state can thus
be included in the parameter estimation as necessary. The covariance
matrix of the initial state  is
set to the integral of the Wiener process and the dynamics of the system over
the first sample interval $t_2-t_1$ as also done in \cite{kristensen03}.

The concept of states is essential to the understanding of the
model setup. The state vector describes the state of the entire system
and is only observable through measurement noise. The actual relation
between measurements and states is defined in the measurement equation
\eqref{sssm2}. A state can represent many different aspects
of the system of interest, e.g. concentrations or amounts in
compartments, a volume, a parameter with unknown time varying behavior
or an input to the system, that we wish to estimate.

\subsubsection*{Second stage model}

The second stage model describes the variation of the individual
parameters $\bs\phi_i$ between individuals and is defined as

\begin{equation}
\bs \phi_i = \bs h(\bs \theta, \bs \eta_i, \bs Z_i)
\label{eqInd}
\end{equation}

\noindent where $\bs \eta_i$ is the multivariate random effect
parameter for the
$i$th individual, which is assumed Gaussian distributed with mean
zero and covariance $\bs \Omega$, i.e.  $\bs \eta_i \in N(\bs 0,\bs
\Omega)$.  The fixed effect parameter of the model is $\bs
\theta$ and $\bs Z_i$ is a vector of co-variates for the $i$th individual.



\section{Estimation}

Parameter estimation is done using maximum likelihood. The likelihood
function will only be outlined briefly here, so please refer to
\cite{mortensen07, overgaard05, kristensen03} for a detailed
description.

The full set of model parameters to be estimated for the final mixed
effects model based on SDEs are the matrices $\bs\Sigma$,
$\bs\sigma_\omega$, $\bs\Omega$ and the fixed effect parameters in the
vector $\bs\theta$. The three matrices are
usually fixed to some degree so that only the diagonals or other
partial structure remains to be estimated. In PSM the parameters in $\bs\Sigma$ and
$\bs\sigma_\omega$ are included in $\bs\theta$.

In PSM the function \verb|ModelPar| defines which part of the model
parameters should estimated. These parameters are denoted $\bs\Theta$ (in
PSM: \verb|THETA|) such that

\begin{equation}
\texttt{ModelPar}:  ~~ \bs\Theta \rightarrow (\bs\theta, \bs\Omega) .
\end{equation}

The exact
population likelihood function cannot be evaluated analytically and
thus a second-order Taylor expansion is made of the individual a
posteriori log-likelihood function around the value of $\hat{\bs
\eta}_i$ that maximizes it. The objective function for PSM is thereby
given as

\begin{eqnarray}
  - \log L(\bs \Theta)
  &\approx& \sum_{i=1}^N \left (    \frac{1}{2} \log \left|
  \frac{-\bs\Delta  l_i}{2\pi}\right | - l_i
  \right ) \label{PSMobj}
\end{eqnarray}

\noindent where $l_i$ is the a posteriori log-likelihood function for
the $i$th individual. This likelihood function is evaluated using the
Kalman Filter which gives an exact solution for linear models. For
non-linear models the Extended Kalman filter (EKF) is used which is only an
approximation. The 2nd derivative $\bs\Delta l_i$
is approximated using the First-Order Conditional Estimation (FOCE) method, in
the same way as it is normally done in mixed effects models based on
ordinary differential equations (ODEs).


\subsubsection*{Uncertainty of parameters}

PSM estimates the uncertainty for the parameter estimates based on the
observed Fisher information.  The parameters to be
estimated are denoted $\bs \Theta$ and the observed information is
then defined as

\begin{equation}
  \bs j ( \bs\Theta ) = -
  \frac{\partial^2}{\partial\bs\Theta\partial\bs\Theta^T} \log L
  (\bs\Theta) = - \bs\nabla^2 \log L(\bs\Theta)
\end{equation}

\noindent which is equal to the Hessian matrix of the negative
log-likelihood function. If the parameters maximizing the
likelihood function are called $\hat{\bs\Theta}$ they will
asymptotically have the distribution

\begin{equation}
  \hat{\bs\Theta} \sim N(\bs\Theta,\bs j ( \hat{\bs\Theta} )^{-1} ). \label{thetadist}
\end{equation}

This is used in PSM to provide a Wald 95\% confidence interval,
standard error and correlation matrix for the estimates. The Hessian
is evaluated using \verb+hessian+ in the \verb+numDeriv+ package.

\subsubsection*{State estimates}

A key feature of the SDE approach to population modelling is the
ability to give improved estimates of the system states given the
individual parameters and also to provide confidence bands for the
states. Confidence bands at a time point $t$ are directly given by the
estimated state  covariance matrix $\hat{\bs P}_{i(t|...)}$ from the
EKF, where $t$ can be both at or between measurements.

There are four types of state and state covariance estimates available
when using the EKF, each of which differs in the way data is used. The four
types are:

\begin{itemize}
\item Simulation estimate: $\hat{\bs x}_{i(j|0)}$, $\hat{\bs P}_{i(j|0)}$ \\
Provides an estimate of the state evolution for a repeated experiment,
without updating based on measurements. This is an ODE-like estimate,
but it also yields a confidence band for the state evolution.
\item Prediction estimate: $\hat{\bs x}_{i(j|j-1)}$, $\hat{\bs P}_{i(j|j-1)}$ \\
The prediction is used here to give the conditional density for
the next observation at time $t_{ij}$ given the observations up to
$t_{i(j|j-1)}$.
\item Filtering estimate: $\hat{\bs x}_{i(j|j)}$, $\hat{\bs P}_{i(j|j)}$ \\
Best estimate at time $t_{ij}$ given the observations up to
time $t_{ij}$.
\item Smoothing estimate: $\hat{\bs x}_{i(j|N)}$, $\hat{\bs P}_{i(j|N)}$\\
Optimal estimate at time $t_{ij}$ utilizing all observations both prior
to and after time $t_{ij}$.

\end{itemize}

For a conventional ODE model the state is found by the simulation
estimate, which is entirely given  by the (possibly ML-estimated)
initial state of the system. The covariance matrix for the states is
$\bs 0$ since no system noise is estimated.

With SDEs three new types of estimates, apart from the simulation
estimate, also become available. In the present setup the
prediction estimate is used to give conditional Gaussian
densities to form the likelihood function. The filter estimate is
the best obtainable state estimate during the experiment, where
the subsequent observations are not present. The third type of state
estimate is the smoothed estimate. This provides the optimal state and
state covariance estimate ($\hat{\bs x}_{i(j|N)}$ and $\hat{\bs P}_{i(j|N)}$) based on all obtained observations, both prior and
subsequent to the time of interest. The smoothed estimate is therefore
often the natural estimate of choice when studying the behavior of
the system in \emph{post hoc} analysis \cite{mortensen07}.

\subsubsection*{Comparison to NONMEM}

The NONMEM software is a widely used tool for mixed effects modelling
based on ODEs \cite{nonmem}. NONMEM also performs maximum likelihood
estimation using the FOCE approximation and for many it might thus be
of interest to know how the two objective functions are related.

The objective function in NONMEM ($l_{NM}$) is advertised as $-2\log L$ but
in fact it lacks a constant equal to the likelihood of the data. The PSM
objective function ($l_{PSM}$) is  $-\log L$ as seen in
Eq. \eqref{PSMobj} and the relation thereby becomes $l_{NM} = 2\cdot
l_{PSM} - \log(2\pi) \cdot \sum n_i $. The relation has been tested
for a number of models \cite{mortensen07}.

\subsubsection*{Models with no random effects}

A special case arises in models where no random effects are specified.
This may be defined in PSM by setting \verb|OMEGA = NULL| in the
\verb|ModelPar| function. This greatly simplifies the likelihood
function, as it is no longer necessary to integrate out the random
effects. The likelihood function is thus reduced to a product of
conditional probabilities of the observations. In the form of a
log-likelihood function this may be written as

\begin{equation}
  - \log L(\bs \Theta) = - \sum_{i=1}^N \sum_{j=1}^{n_i} \log
  p(y_{ij}|\cdot) \label{CTSMobj}
\end{equation}

\noindent where the $\cdot$ indicates conditioning on parameters and past
observations for individual $i$. The latter conditioning is necessary
due to the inclusion of SDEs in the model.

The likelihood function in Eq. \eqref{CTSMobj} is identical to the one
used in CTSM \cite{kristensen03}, and is sometimes referred to as a
pooled likelihood. The inner sum of the likelihood function is also
equal to the a priori individual log-likelihood function in the mixed
effects framework.

\subsubsection*{Implementation issues}
The estimation algorithm in PSM for linear models is based on the
ordinary Kalman filter, which
has been written in Fortran for faster execution times. However, the
Fortran code does not support a singular $A$-matrix, and will in these
cases fall back on an R version of the Kalman filter. This may be
circumvented by adding a very small value to the diagonal of $A$,
at least in order to find the first rough parameter estimates.

\newpage
\section{User's guide to PSM}

PSM is built around two key objects. These are

\begin{itemize}
\item a data object and
\item a model object.
\end{itemize}

The data object contains sample times, observations and possible
input, covariates and dosing regimen for all individuals and the model
object contains everything related to the model.

\subsubsection*{Model object}

Before setting up a model in PSM it is a good idea to write it down on
paper and note the dimensions of the state, observations and possible
input and random effects. When this is done the function
\verb+PSM.template()+ can be used for both linear and non-linear
models as shown below to print a template for the model.

<<PSMtemplateLinear>>=
PSM.template(Linear=TRUE,dimX=3,dimY=1,dimU=0,dimEta=3)
@

The structure of the model object and what each function should return
can be derived from the template
shown above. It is important to keep the input arguments of all
functions unchanged even though a particular model may not need use
every argument in a function.

The input arguments to the functions are defined as follow:


\begin{center}
\begin{tabular}{cp{7cm}}
\toprule
Input arg. & Type \\
\hline
\verb+THETA+ &  Vector (possibly named) containing the parameters to be
maximum likelihood estimated.\\
%\hline
\verb+theta+ &  Vector of all parameters in the model.\\
%\hline
\verb+phi+ & As standard a list of individual parameters, but can also
  be a vector as defined in \verb+$h+. \\
%\hline
\verb+eta+ &  An unnamed vector of random effects.\\
%\hline
\verb+Time, time+ &  A scalar value.\\
%\hline
\verb+covar+ &  As specified in the user defined data object.\\
%\hline
\verb+x,u,U+ &  Column matrices with state or input at a sample time.\\
\bottomrule
\end{tabular}
\end{center}

%\vspace{6mm}

\subsubsection*{Data object}

The data object is a list with one element for each individual in the
data set. Each element in the list must contain:

\begin{itemize}
  \item \verb+Y+ - Matrix with observations. Each columns holds one
    (possibly multi-dimensional) observation. Missing observations
    should be marked as NA.
  \item \verb+Time+ - Vector of sample times. The length must correspond
    to the number of columns in \verb+Y+ for the individual.
\end{itemize}

\noindent Each element in the list can optionally contain:

\begin{itemize}
  \item \verb+U+ - Matrix with input at sample times. The input cannot contain
    missing values and is assumed constant between sample times
    (zero-order hold). It must have the same number of columns as \verb+Y+.
  \item \verb+covar+ - A vector/list with covariates to be used in the
    function \verb+$h+.
  \item \verb+Dose+ - A list containing three vectors: Time, State and
    Amount. See help(PSM.estimate) for more detail.
\end{itemize}

The data object is illustrated with a small example. The object shown
contains 4 and 5 observations for two individuals sampled at different
times and it also has 'BMI' as a covariate for each individual.

<<DataObj>>=
MyData <- list()
MyData[[1]] <- list(Time=1:4,Y=matrix(c(2.1,3.2,3.4,3.7),nrow=1),covar=c(BMI=20.1))
MyData[[2]] <- list(Time=3:7,Y=matrix(c(1.9,2.1,2.0,2.9,3.5),nrow=1),covar=c(BMI=23.4))
@

\subsubsection*{Main functions}

The PSM program is accessed in R through five main functions:

\begin{itemize}
\item \verb+PSM.simulate(Model, Data, THETA, deltaTime)+ \\
  Simulates data for multiple individuals. The number of individuals is
  determined by length(Data). The simulation is based on the Euler method
  to be able to simulate SDEs and thus a short time step should be chosen.
\item \verb+PSM.estimate(Model, Data, Par, CI)+\\
  Estimates population parameters for any linear or non-linear
  model. The Par argument is a list containing initial estimates and
  bounds for the parameter search.
\item \verb+PSM.smooth(Model, Data, THETA, subsample)+\\
  Optimal estimates of model states based on estimated parameters. It
  returns both the predicted, filtered and smoothed state
  estimates and for models with random effects an estimate of these are
  also returned.
\item \verb+PSM.plot(Data, Smooth, indiv, type)+\\
  Creates a matrix plot with a column for each individual. The rows
  can show observations, inputs, simulated and estimated states, residuals and
  auto-correlation functions. The x and/or y axis can be on log-scale and
  it is possible to list simulated or estimated random effects on the plot.
\item \verb+PSM.template(Linear,dimX,dimY,dimU,dimEta,file)+\\
  Creates a template with R-syntax to help setup a model in PSM. It
  works for both linear and non-linear models and it can output the
  resulting template to the screen or a file.

\end{itemize}

For detailed information please refer to the help files in R.



\newpage
\section{Examples}

\begin{tabular}{llr}
\ref{dose2comp} & Dosing in two-compartment model (Linear)  . . . .& \pageref{dose2comp} \\
\ref{isr} & Extraction of insulin secretion rate (Linear)  . . . .& \pageref{isr} \\
\end{tabular}

\vspace{5mm}

\noindent The following provides two examples to illustrate the use of PSM. It
may also be useful when trying to set up a new model, by using the
code shown as model templates.
The document is written using
Sweave\footnote{http://www.ci.tuwien.ac.at/$\sim$leisch/Sweave/}, and thus
all R-code in the examples shown here can be extracted into an R-script
file by calling \verb+Stangle()+ on the Sweave file, which is found in
the PSM package installation folder.

To save time in the processing of the document all computer intensive
operations are skipped by setting a flag \verb|Redo = FALSE| and
instead the outcome is loaded from a saved .Rdata file. Changing the
flag to \verb|TRUE| in PSM.Rnw and writing \verb|Sweave("PSM.Rnw")| in
R will generate a new version of the document with analysis and plots
based on the new simulated data sets.

<<Redo.flag>>=
Redo = FALSE
@

\newpage

\subsection{Dosing in two-compartment model (Linear)}
\label{dose2comp}

This example illustrates how a standard two-compartment model with a
random diffusion between the compartments can be set up. An overview
of the model is shown in Figure \ref{fig:modeldose}.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.8\textwidth]{images/comp2dose}
  \caption{Model layout}
  \label{fig:modeldose}
\end{figure}

The model is used to simulate data based on two doses of $1500mg$
given after 30 and 180 minutes. In state space formulation the model
is described as

\begin{eqnarray}
  dA_1 &=& \left(-\frac{CL}{V_1^i}A_1 - \frac{CL_d}{V_1^i}A_1 +
    \frac{CL_d}{V_2}A_2 \right )dt + \sigma_1d\omega \\
  dA_2 &=& \left( \frac{CL_d}{V_1^i}A_1 - \frac{CL_d}{V_2}A_2
\right) dt - \sigma_1d\omega\\
  Y &=& A_1/V_1^i + e
\end{eqnarray}

\noindent or in matrix notation

\begin{eqnarray}
  d\bs A_t &=& \left[ \begin{array}{cc}
      -(CL/V_1^i + CL_d/V_1^i) & CL_d/V_2  \\
    CL_d/V_1^i & - CL_d/V_2 \end{array} \right]\bs A_tdt + \left[
  \begin{array}{cc} \sigma_1 & 0 \\ -\sigma_1 & 0 \end{array}
\right]d\bs\omega\\
  Y_{ij} &=& [ \begin{array}{cc} 1/V_1^i & 0 \end{array}]\bs A_{ij} + e_{ij}
\end{eqnarray}

\noindent where $\bs A_t = [ A_1 ~ A_2]^T$ is the amount in each compartment and thus
$A_1/V_1^i$ is the measured concentration. It is seen that the
elimination will follow a normal two-compartment model, but with a
small random diffusion between the compartments. The mass is preserved
since the diffusion terms are equal with opposite signs.

For simplicity the individual variation is modelled as

\begin{equation}
V_1^i =  V_1\exp(\eta_1)
\end{equation}

The parameters of the model is $V_1= 5L$, $V_2 = 10L$, $CL_d =  0.005
L/min$, $CL = 0.002 L/min$, $\sigma_1 = 10$, $S = 20 mg^2/L^2$ and
$\Omega = 0.5$.

The model can be defined in R as shown below.

<<Dose.ModelDefinition>>=
Model.SimDose = list()
Model.SimDose$Matrices = function(phi) {
  V1i <- phi$V1i; V2=phi$V2; CL = phi$CL; CLd = phi$CLd;
  matA <- matrix(c(-(CL+CLd)/V1i ,  CLd/V2 ,
                   CLd/V1i , -CLd/V2 ) ,nrow=2,byrow=T)
  matC <- matrix(c(1/V1i,0),nrow=1)
  list(matA=matA,matC=matC)
}
Model.SimDose$X0 = function(Time=Na,phi,U=Na) {
  matrix(0,nrow=2)
}
Model.SimDose$SIG = function(phi) {
  sig1 <- phi[["sig1"]]
  matrix(c( sig1,0,
           -sig1,0), nrow=2, byrow=T)
}
Model.SimDose$S = function(phi) {
  matrix(phi[["S"]])
}
Model.SimDose$h = function(eta,theta,covar) {
  phi <- theta
  phi$V1i <- theta$V1*exp(eta[1])
  phi
}
Model.SimDose$ModelPar = function(THETA){
  V2 <- 10
  CLd <- 0.1
  list(theta=list(V1 = THETA['V1'],V2=V2,CLd=CLd,CL=THETA['CL'], sig1=THETA['sig1'], S=THETA['S']),
       OMEGA=matrix(THETA['OMEGA1']) )
}
SimDose.THETA <-  c(CL=0.05,V1 = 5, sig1 = 10 , S = 20 , OMEGA1 = .2)
@

Five parameters in the model will be estimated, as it can be
seen from the \verb|ModelPar| function above. The parameters to be estimated
are $\bs\Theta$ = (\verb|CL|, \verb|V1|, \verb|sig1|, \verb|S|, \verb|OMEGA1|).

For this example 5 individuals will be simulated. They will
all be sampled every $10min$ for $400min$ which is described as below.

<<SamplingScheme>>=
N = 5
SimDose.Data <- vector(mode="list",length=N)
for (i in 1:N) {
  SimDose.Data[[i]]$Time <- seq(from=10,by=10,to=400)
  SimDose.Data[[i]]$Dose <-list(
                                Time = c(30,180),
                                State = c(1, 1),
                                Amount = c(1500,1500)
                                )
}
@

Everything is now setup and the simulation can be performed.

<<Dose.SimulateData>>=
if(Redo) {
  SimDose.Data <- PSM.simulate(Model.SimDose, SimDose.Data, SimDose.THETA, deltaTime=.1)
} else
  load("simdose.RData")
@

The simulated data are shown in Figure \ref{fig:fig1} using the
\verb+PSM.plot+ function. The first row
shows the observations for individuals 1 and 2, the next two show
state 1 and 2 which we wish to estimate and the simulated
values of $\eta_1$ is shown.

\begin{figure}[htb]
\begin{center}
%<<label=figure1,fig=TRUE,echo=FALSE>>=
<<label=figure1,fig=TRUE,echo=TRUE,height=5.5,width=6>>=
PSM.plot(SimDose.Data,indiv=1:2,type=c('Y','longX','eta'))
#par(mfcol=c(3,2),mar = c(2, 4, 2, 2)+.1)
#for(id in 1:2) {
#  plot(SimDose.Data[[id]]$Time , SimDose.Data[[id]]$Y,
#         ylab="Observations", main=paste('Individual ',id,', eta= ',
#                                round(SimDose.Data[[id]]$eta,3),sep=""))
#  for(i in 1:2) {
#    plot(SimDose.Data[[id]]$longTime , SimDose.Data[[id]]$longX[i,],type="l",
#         ylab=paste('State',i))
#    rug(SimDose.Data[[id]]$Time)
#  }
#}
@
\end{center}
\caption{Simulated data and states.}
\label{fig:fig1}
\end{figure}

As initial guess for the parameters in $\bs\Theta$ the true parameters
are used and the bounds are $\pm$ a factor away.

<<Dose.EstimateParameters>>=
parA <- list(LB=SimDose.THETA*.5, Init=SimDose.THETA , UB=SimDose.THETA*1.5 )
if(Redo) fitA <- PSM.estimate(Model.SimDose,SimDose.Data,parA,CI=T)
fitA[1:5]
SimDose.THETA
@

It is seen that the estimated parameters are reasonably close to the
true values in \verb+SimDose.THETA+ and the 95\% confidence intervals
include the true values. In particular the first parameter $\sigma_1$
is significantly different from zero which shows that deviations from
a normal ODE two-compartment model are significant.

Based on the estimated parameters it is possible to obtain an estimate
of the model states by using \verb|PSM.smooth|. The
estimates are shown in Figure \ref{fig:fig2}. The structure of the
output from the smoothing function is also shown using the
\verb|names| command. Please refer to \verb|help(PSM.smooth)| for a
more detailed description of the output.

<<Dose.EstimateStates>>=
if(Redo)
  out <- PSM.smooth(Model.SimDose, SimDose.Data, fitA$THETA, subsample = 20)
# View the data structure
names(out[[1]])
@

By comparing the smoothed estimates of the states to the
true simulated states, it can be seen that they are very close. This
shows that the system noise and observation noise has been separated
properly in the reconstruction.

\begin{figure}[hptb]
\begin{center}
<<label=figure2,fig=TRUE,echo=FALSE>>=
par(mfcol=c(3,2),mar = c(2, 4, 2, 2)+.1)
pos <- c('topright','bottomright')
for(id in 1:2) {
  tmp <- out[[id]]$Ys
  plot(SimDose.Data[[id]]$Time , SimDose.Data[[id]]$Y,
         ylab="Observations", main=paste('Individual',id),ylim=range(tmp))
  legend(x=pos[1], legend=c('Smooth est.'),
         lty=c('71A1'),col='blue', xjust=1,yjust=1,bty="n")
  lines(out[[id]]$Time , out[[id]]$Ys,lty='71A1',col='blue')
  for(i in 1:2) {
    plot(SimDose.Data[[id]]$longTime , SimDose.Data[[id]]$longX[i,],type="l",
         ylab=paste('State',i))
    lines(out[[id]]$Time , out[[id]]$Xs[i,],lty='71A1',col="blue")
    rug(SimDose.Data[[id]]$Time)
    legend(x=pos[i], legend=c('Simulation','Smooth est.'),
           lty=c('solid','71A1'),col=c("black","blue"),xjust=1,
           yjust=1,bty="n",y.intersp=.8)
  }
}
@
\end{center}
\caption{Smoothed estimate of states.}
\label{fig:fig2}
\end{figure}

\newpage
\subsection{Extraction of insulin secretion rate (Linear)}
\label{isr}

Insulin secretion rates (ISR) can be estimated based on measurements
of the concentration of C-peptide in the blood, since insulin and
C-peptide are secreted in equi-molar amounts. This example will first
illustrate a way to simulate C-peptide data based on a model for ISR,
and then how ISR can be estimated again using a more simple model. The
models used in the example are described in further detail in
\cite{klim06}.

The simulated measurements of C-peptide spans over 24H, during
which the patients recieves three meals at 8 a.m., 12 a.m. and 6
p.m. These meals will give rise to an increase in insulin secretion,
which will be modelled and estimated.

The simulation model for the C-peptide measurements is based on the
commonly used two compartment model for C-peptide as shown in Figure
\ref{fig:modelisr}. The kinetic parameters are set equal to the Van
Cauter estimates \cite{vancauter92}.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.65\textwidth]{images/isrmodel}
  \caption{Model layout}
  \label{fig:modelisr}
\end{figure}

The first two states of the simulation model is concentration in
compartment 1 and 2, $C1$, $C2$. The third state is $ISR$ and is the
secretion which is modelled as a structural part based on the three
meal time plus a constant baseline $B$ and random noise through a
Wiener process. The fourth state $Q$ is used to model ISR, where $Q$
is controlled by an input $u_2$ which is equal to 1 for $30min$ after meal
times. The model can be defined as

\begin{eqnarray}
  dC_1 &=& [-(k_1+k_e)C_1 +k_2C_2 + ISR]dt \\
  dC_2 &=& [k_1C_1 - k_2C_2]dt \\
  dISR  &=& [-a_1ISR + a_1Q + B^i]dt + \sigma_{ISR}d\omega \label{eqisr3}\\
  dQ    &=& [-a_2Q + a_2K^iu_2]dt \label{eqisr4}
\end{eqnarray}

\noindent which is linear and can thus again be written on matrix
form.

The model is initialized in steady state just prior to the
first meal time. The individual variation in the model is included in
the initial concentration $C_1(0)^i$, baseline $B^i$ and height of the
peaks $K^i$ such that

\begin{eqnarray}
  B^i &=& B\exp(\eta_1) \\
  K^i &=& K\exp(\eta_2) \\
  C_1(0)^i &=& C_1(0)\exp(\eta_3)
\end{eqnarray}

In order to write a model containing a constant in the differential
equations (here $B^i$) on the linear form as defined in Eq.
\eqref{sssm1} and \eqref{sssm2} it is necessary to include a constant
input $u_1=1$ and multiply this with $B^i$. The matrix description of
the model can be found in \cite{klim06} p. 69. Using this, the model
can be defined in PSM as follows.


<<ISR.ModelDefinition>>=
k1 = 0.053; k2 = 0.051; ke = 0.062;
Model.SimISR <- list()
Model.SimISR$Matrices = function(phi) {
  a1  <- phi[["a1"]]
  a2  <- phi[["a2"]]
  B  <- phi[["B"]]
  K  <- phi[["K"]]
  matA <- matrix( c(-(k1+ke) ,  k2 ,   1 ,   0,
                    k1 , -k2 ,   0 ,   0,
                    0 ,   0 , -a1 ,  a1,
                    0 ,   0 ,   0 , -a2),nrow=4,byrow=T)
  matB <- matrix( c(0 , 0   ,
                    0 , 0   ,
                    B , 0   ,
                    0 , a2*K),byrow=T,nrow=4)
  matC <- matrix(c(1,0,0,0),nrow=1)
  matD <- matrix(c(0,0),nrow=1)
  list(matA=matA,matB=matB,matC=matC,matD=matD)
}
Model.SimISR$X0 = function(Time=NA,phi,U=NA) {
  C0 <- phi[["C0"]]
  tmp    <- C0
  tmp[2] <- C0*k1/k2
  tmp[3] <- C0*ke
  tmp[4] <- 0
  matrix(tmp,ncol=1)
}
Model.SimISR$SIG = function(phi) {
  diag( c(0,0,phi[["SIG33"]],0))
}
Model.SimISR$S = function(phi) {
  return( matrix(phi[["S"]]))
}
Model.SimISR$h = function(eta,theta,covar) {
  phi <- theta
  phi[["B"]] <- theta[["B"]]*exp(eta[1])
  phi[["K"]] <- theta[["K"]]*exp(eta[2])
  phi[["C0"]] <- theta[["C0"]]*exp(eta[3])
  return(phi)
}
Model.SimISR$ModelPar = function(THETA){
  list(theta=list(C0=900,S=8500,
                a1=THETA['a1'],a2=THETA['a2'],
                SIG33=THETA['SIG33'],
                K = THETA['K'], B = THETA['B']),
              OMEGA=diag(c(.2,.2,.2))
       )
}

@

For this example two individuals will be simulated. They will both be
sampled at predefined time points during 24H. This defined as below
together with the input data for each individual.

<<ISR.SimulateData>>=
Sim.Data <- vector(mode="list",length=2)
for (i in 1:2) {
  Sim.Data[[i]]$Time <- c( 0,15,30,45,60,75,90,120,150,180,210,240,270,300,330,
                          360,420,480,600,615,630,645,660,675,690,720,750,780,810,
                          840,960,1140,1320,1410,1440)
  Sim.Data[[i]]$U <- matrix(c( rep(1,35) ,
                              as.numeric( Sim.Data[[i]]$Time %in% c(0,15,240,600,615)) )
                            ,byrow=T,nrow=2)
}
@

Both the model, sample times and input is now prepared and the
simulation can be performed. The parameter estimates are taken from
\cite{klim06} p. 70. The simulated data are shown in Figure
\ref{fig:fig3}.

<<ISR.SimulateData>>=
Sim.THETA <-  c(a1=0.02798, a2=0.01048, SIG33=4 , K=427.63 , B=1.7434)
if(Redo) {
  Sim.Data <- PSM.simulate(Model.SimISR, Sim.Data, Sim.THETA, deltaTime=.1 )
} else
  load("simisr.RData")
@


\begin{figure}[htb]
\begin{center}
<<label=figure3,fig=TRUE,echo=FALSE>>=
par(mfcol=c(4,2),mar = c(2, 4, 2, 2)+.1)
for(id in 1:2) {
  for(i in 1:4) {
    plot(Sim.Data[[id]]$longTime , Sim.Data[[id]]$longX[i,],type="l", ylab=paste('State',i))
    rug(Sim.Data[[id]]$Time)
    if(i==1) {
      title(main=paste('Individual ',id,', eta=(',paste(round(Sim.Data[[id]]$eta,3),
              collapse=","),')',sep=""))
      points(Sim.Data[[id]]$Time, Sim.Data[[id]]$Y)
      legend(1440, y = max(Sim.Data[[id]]$X[i,]), legend=c('Obs.'),
           pch=21,xjust=1, yjust=1)
    }
  }
}
@
\end{center}
\caption{Simulated data and states.}
\label{fig:fig3}
\end{figure}

The next step is to generate a model for estimation of ISR. It is
again based on the model illustrated in Figure \ref{fig:modelisr} only
this time the ISR is simply modelled as a random walk. Thus no
information about the meal times is used in the estimation of ISR.

The model simplification is done by replacing Eq. \eqref{eqisr3} and
\eqref{eqisr4} by Eq. \eqref{eq:isrny} below. The model for estimation
can thus be seen as estimating the outcome of the random walk for
$ISR$ based on the observed (simulated) data for $C1$.

\begin{equation}
  dISR = \sigma_{ISR} \label{eq:isrny}
\end{equation}

<<ISR.EstimationModel>>=
Model.Est <- list(
            Matrices=function(phi) { list(
              matA=matrix(c(-(k1+ke),  k2, 1,
                             k1     , -k2, 0,
                             0      ,   0, 0  ),ncol=3,byrow=T),
              matB=NA,
              matC=matrix(c(1,0,0),nrow=1),
              matD=NA )  },
            X0 = function(Time=NA,phi=NA,U=NA) {
              C0 <- phi[["C0"]]
              tmp    <- C0
              tmp[2] <- C0*k1/k2
              tmp[3] <- C0*ke
                    return(matrix(tmp,ncol=1) )} ,
            SIG = function(phi) {
                      return( diag( c(1e-3,1e-3,phi[["SIG33"]])) ) } ,
            S = function(phi) {
                      return( matrix(phi[["S"]])) } ,
            h = function(eta,theta,covar) {
              phi <- theta
              phi[["C0"]] <- theta[["C0"]]*exp(eta[1])
              return(phi) } ,
            ModelPar = function(THETA){
              return(list(theta=list(C0=THETA['C0'],S=THETA['S'],SIG33=THETA['SIG33']),
                          OMEGA=matrix(THETA['OMEGA'])))}
            )
@

Looking at the \verb|ModelPar|-function it is seen that it is chosen
to include the average initial concentration $C_1(0)$, measurement
variation $S$, the coefficient of the random walk for $ISR$
$\sigma_{ISR}$ and the variance of the random effect on $C_1(0)$
denoted $\Omega_{C_1(0)}$ in the likelihood estimation.

Since the model now does not use any input, this must be removed from
the simulated data before estimation.

<<ISR.RemoveInputFromData>>=
Pop.Data <- Sim.Data
for (i in 1:2)
  Pop.Data[[i]]$U <- NULL
@

The data and model for estimation is now prepared, and the model can
be estimated by calling \verb|PSM.estimate|. This is done below and
the output in \verb|obj1[1:3]| containing the log-likelihood value,
parameter estimates and confidence intervals is shown as output.

<<ISR.ParameterEstimation>>=
par1 <- list(LB   = c(C0= 200, S= 50^2, SIG33= 0, OMEGA=.0 ),
             Init = c(C0=1000, S=100^2, SIG33=10, OMEGA=.25),
             UB   = c(C0=3000, S=150^2, SIG33=15, OMEGA=.50))
if(Redo) obj1 <- PSM.estimate(Model.Est, Pop.Data, par1,CI=T,trace=1)
obj1[1:5]
@



Looking at the estimated confidence intervals, it is seen that the
values used in the simulation $\bs\Theta$ = (900,8500,4,0.2) are
nicely contained within the limits.

The estimation time including the confidence interval is about 3
minutes on a 2GHz computer. Since the matrix $A$ in the estimation
model is singular, the estimation cannot make use of the compiled
Fortran code. As mentioned, this may be circumvented by adding
e.g. $10^{-6}$ to the diagonal. This reduces the estimation time to 12
sec. It changes the maximum log-likelihood value to 499.4781, and thus
yields virtually no difference in parameter estimates.

Using the estimated model parameters it is possible to give smoothed
estimates of the three model states $C1$, $C2$ and $ISR$. This is done
below and the result is plottet in Figure \ref{fig:fig4}. In Figure \ref{fig:fig5}
the smoothed $ISR$ state is plotted together with the estimated
uncertainty. For both figures the true simulated states is also
plotted for reference.

<<ISR.EstimateStates>>=
if(Redo)
  Data.Sm <- PSM.smooth( Model.Est , Pop.Data, obj1$THETA, subsample=10)
@


\begin{figure}[hptb]
\begin{center}
<<label=figure4,fig=TRUE,echo=FALSE>>=
par(mfcol=c(4,2),mar = c(2, 4, 2, 2)+.1)
for(id in 1:2) {
  for(i in 1:4) {
    plot(Sim.Data[[id]]$longTime , Sim.Data[[id]]$longX[i,],type="l",
         ylab=paste('State',i), xlab=paste('Individual',id))
    if(i<=3) lines(Data.Sm[[id]]$Time, Data.Sm[[id]]$Xs[i,],col="blue",lty='71A1')
    rug(Sim.Data[[id]]$Time)
    legend(1440, y = max(Sim.Data[[id]]$X[i,]), legend=c('Simul.','Sm. est.'),
       lty=c('solid','71A1'),col=c("black","blue"),xjust=1,
       yjust=1,bty="n",y.intersp=.8)
    if(i==1)  {
      title(main=paste('Individual',id))
      points(Sim.Data[[id]]$Time, Sim.Data[[id]]$Y)
    }
  }
}
@
\end{center}
\caption{Smoothed estimate of states.}
\label{fig:fig4}
\end{figure}

\begin{figure}[hptb]
\begin{center}
<<label=figure5,fig=TRUE,echo=FALSE,height=4.5,width=6>>=
par(mfcol=c(2,1),mar = c(3,3, 0, 2)+.1)
for(ID in 1:2) {
  CMT <- 3
  Data <- Data.Sm[[ID]]
  polyCI = c(Data$Xs[CMT,]+sqrt(abs(Data$Ps[CMT,CMT,])),
    rev(Data$Xs[CMT,]-sqrt(abs(Data$Ps[CMT,CMT,]))))
  ymax=max(polyCI)
  plot.new()
  plot.window(xlim=range(Data$Time),ylim=c(0,ymax))
  axis(1);axis(2);box()
  title(xlab="minutes",ylab="pmol/min",line=2)
  polygon( c(Data$Time,rev(Data$Time)) ,polyCI,col=rgb(.9,.9,1),border=NA)
  polygon( c(Data$Time,rev(Data$Time)) ,polyCI,col=rgb(.6,.6,1),density=0)
  lines( Data$Time, Data$Xs[CMT,], type="l",lwd=2,col="blue")
  lines( Sim.Data[[ID]]$longTime, Sim.Data[[ID]]$longX[3,],lwd=.5)
  #points( Sim.Data[[ID]]$Time, Sim.Data[[ID]]$X[3,],pch=20,col=1,cex=.7)
  legend(1440, y = max(polyCI), legend=c('Simulation','Smooth est.'),
         col=c("black","blue"),xjust=1,
         yjust=1,bty="n",y.intersp=.8,lwd=c(1,2),pch=c(NA,NA))
  box()
}
@
\end{center}
\caption{Smoothed estimate of insulin secretion rate $\pm$1SD for
  individual 1 and 2 compared with the true simulated ISR.}
\label{fig:fig5}
\end{figure}

\newpage

\small
\begin{thebibliography}{99}
  \addcontentsline{toc}{section}{\refname}
  \bibitem{mortensen07}
    Mortensen SB, Klim S, Dammann B, Kristensen NR, Madsen H,
    Overgaard RV (2007) A Matlab framework for estimation of nlme
    models using stochastic differential equations: applications for
    estimation of insulin secretion rates. J of Parmacokinet
    Pharmacodyn 34:623-642

  \bibitem{overgaard05}
    Overgaard RV, Jonsson N, Torn\o e CW, Madsen H (2005) Non-linear
    mixed-effects models with stochastic differential equations:
    implementation of an estimation algorithm. J of Parmacokinet
    Pharmacodyn 32(1):85-107

  \bibitem{kristensen03}
    Kristensen NR, Madsen H (2003) Continous time stochastic
    modelling: CTSM 2.3 mathematics guide, Technical University of
    Denmark \\
    http://www2.imm.dtu.dk/ctsm/MathGuide.pdf

  \bibitem{nonmem}
    Beal SL, Sheiner LB (2004) NONMEM\textregistered Users Guide. University of
    California, NONMEM Project Group.

  \bibitem{klim06}
    Klim S, Mortensen SB (2006) Stochastic PK/PD Modelling. M.Sc. thesis,
    Informatics and Mathematical Modelling, Technical University of
    Denmark. \\
    http://www2.imm.dtu.dk/pubdb/views/edoc\_download.php/4533/pdf/imm4533.pdf

  \bibitem{vancauter92}
    Van Cauter E, Mestrez F, Sturis J, Polonsky KS (1992) Estimation of
    insulin secretion rates from C-peptide levels. Comparison of
    individual and standard kinetic parameters for C-peptide
    clearance. Diabetes, 41(3), pp. 368-77.


\end{thebibliography}

\end{document}
