Inter occasion variability (mlxR >= 3.2.2)
Marc Lavielle August 31st, 2017 Â
R script: iov.R
$$ \newcommand{\esp}[1]{\mathbb{E}\left(#1\right)} \newcommand{\var}[1]{\mbox{Var}\left(#1\right)} \newcommand{\deriv}[1]{\dot{#1}(t)} \newcommand{\prob}[1]{ \mathbb{P}\!(#1)} \newcommand{\eqdef}{\mathop{=}\limits^{\mathrm{def}}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bc}{\boldsymbol{c}} \newcommand{\bpsi}{\boldsymbol{\psi}} \def\pmacro{\texttt{p}} \def\like{{\cal L}} \def\llike{{\cal LL}} \def\logit{{\rm logit}} \def\probit{{\rm probit}} \def\one{{\rm 1\!I}} \def\iid{\mathop{\sim}_{\rm i.i.d.}} \def\simh0{\mathop{\sim}_{H_0}} \def\df{\texttt{df}} \def\res{e} \def\xomega{x} \newcommand{\argmin}[1]{{\rm arg}\min_{#1}} \newcommand{\argmax}[1]{{\rm arg}\max_{#1}} \newcommand{\Rset}{\mbox{$\mathbb{R}$}} \def\param{\theta} \def\setparam{\Theta} \def\xnew{x_{\rm new}} \def\fnew{f_{\rm new}} \def\ynew{y_{\rm new}} \def\nnew{n_{\rm new}} \def\enew{e_{\rm new}} \def\Xnew{X_{\rm new}} \def\hfnew{\widehat{\fnew}} \def\degree{m} \def\nbeta{d} \newcommand{\limite}[1]{\mathop{\longrightarrow}\limits_{#1}} \def\ka{k{\scriptstyle a}} \def\ska{k{\scriptscriptstyle a}} \def\kel{k{\scriptstyle e}} \def\skel{k{\scriptscriptstyle e}} \def\cl{C{\small l}} \def\Tlag{T\hspace{-0.1em}{\scriptstyle lag}} \def\sTlag{T\hspace{-0.07em}{\scriptscriptstyle lag}} \def\Tk{T\hspace{-0.1em}{\scriptstyle k0}} \def\sTk{T\hspace{-0.07em}{\scriptscriptstyle k0}} \def\thalf{t{\scriptstyle 1/2}} \newcommand{\Dphi}[1]{\partial_\pphi #1} \def\asigma{a} \def\pphi{\psi} \newcommand{\stheta}{{\theta^\star}} \newcommand{\htheta}{{\widehat{\theta}}}$$
1 Introduction
Assumption that an individual parameter remains constant over time can be weakened by introducing intra-individual variability, or inter-occasion variability (IOV), of individual parameters into the model.
The model then consists of splitting the study into \(K\) time periods or occasions and assuming that individual parameters can vary from occasion to occasion but remain constant within occasions.
In other words, we assume that there exists \(\tau_0=0 < \tau_1 < \tau_2 < \ldots\) such that the individual parameter for individual \(i\) is \(\psi_{ik}\) within occasion \(k\), i.e. between time \(\tau_{k-1}\) and time \(\tau_k\).
We will assume that the (possibly transformed) individual parameters are normally distributed. Then, there exists for each individual \(i\),
a vector \(\eta_i^{(0)}\), of random effects which describes the random inter-individual variability of the individual parameters,
a vector \(\eta_{ik}^{(1)}\), of random effects which describes the random intra-individual variability of the individual parameters in occasion \(k\), for each \(1\leq k \leq K\).
Here and in the following, the superscript \({(0)}\) is used to represent inter-individual variability, i.e., variability at the individual level, while superscript \({(1)}\) represents inter-occasion variability, i.e., variability at the occasion level for each individual.
Then, model for parameter \(\psi_{ik}\) assumes that there exists a transformation \(h\) (log, logit, probit, identity, …) such that \[ h(\psi_{ik}) = h(\tilde{\psi}_{ik}) + \eta_i^{(0)} + \eta_{ik}^{(1)} \] where \(\tilde{\psi}_{ik}\) is the prediction of \(\psi_{ik}\), i.e. of parameter \(\psi\) for individual \(i\) within occasion \(k\).
Assume, for example, that the predicted individual parameter \(\tilde{\psi}_{ik}\) has been defined in the EQUATION block, then the definition of \(\psi_{ik}\) will include the variability levels (id for individual and id*occ for occasion, for example), and the standard deviations of the random effects for each variability level:
DEFINITION: psi = {distribution = lognormal, prediction = psi_pred, varlevel = {id, id*occ}, sd = {oV0, oV1}}
Remark 1: mlxR 3.2.2 supports inter-occasion variability for the individual parameters defined in section [INDIVIDUAL] of the Mlxtrancode and for the individual covariates defined in section [COVARIATE]
Remark 2: mlxR 3.2.2 only supports nested levels of variability. When there is only one level of variability, as described in this article, the variability levels should be defined by varlevel={id, id*occ} or just varlevel={id*occ} if there is no IIV in the model.
2 A basic example
2.1 Defining the same occasions for all the individuals
We consider in this example a one compartment PK model with IIV for absorption rate constant \(ka\) and volume \(V\), and IIV and IOV for clearance \(Cl\).
mymodel <- inlineModel("
[LONGITUDINAL]
input = {ka, V, Cl}
EQUATION:
Cc = pkmodel(ka, V, Cl)
[INDIVIDUAL]
input = {ka_pop, omega_ka, V_pop, omega_V, Cl_pop, omega_Cl, gamma_Cl}
DEFINITION:
ka = {distribution=lognormal, reference=ka_pop, sd=omega_ka}
V = {distribution=lognormal, reference=V_pop, sd=omega_V}
Cl = {distribution=lognormal, reference=Cl_pop, varlevel={id, id*occ}, sd={omega_Cl,gamma_Cl}}
")We consider in this example three occasions of 12 hours. The beginnings of these 3 occasions should then be defined in a vector:
occ <- list(name="occ", time=c(0,12,24))Parameters, outputs and design are defined as usual
param <- c(ka_pop=1, omega_ka=0.05, V_pop=10, omega_V=0.05, Cl_pop=5, omega_Cl=0.3, gamma_Cl=0.2)
trt <- list(time=c(0,12,24), amount=100)
out.cc <- list(name= 'Cc', time=seq(0,36,by=0.1))
out.param <- list(name=c( 'ka', 'V', 'Cl'))simulx now uses
res <-simulx(model = mymodel,
parameter = param,
treatment = trt,
varlevel = occ,
output = list(out.cc, out.param),
group = list(size=3, level='individual'),
settings = list(seed=123123))
print(names(res))## [1] "Cc" "group" "occasion" "treatment"
## [5] "parameter.iiv" "parameter.iov"An additional output
print(res$occasion)## id time occ
## 1 1 0 1
## 2 1 12 2
## 3 1 24 3
## 4 2 0 1
## 5 2 12 2
## 6 2 24 3
## 7 3 0 1
## 8 3 12 2
## 9 3 24 3Parameters defined as outputs are now splitted into constant individual parameters (with inter individual variability) and time varying individual parameters (with inter occasion variability).
Note that the three typical individual clearances, denoted
print(res$parameter.iiv)## id ka V Cl0
## 1 1 0.9766671 9.491637 4.997328
## 2 2 1.0314359 9.847188 10.638905
## 3 3 0.9731704 10.241550 5.201438print(res$parameter.iov)## id time occ Cl
## 1 1 0 1 4.596113
## 2 1 12 2 5.297767
## 3 1 24 3 7.694320
## 4 2 0 1 12.190594
## 5 2 12 2 12.572875
## 6 2 24 3 9.927218
## 7 3 0 1 5.160954
## 8 3 12 2 5.683008
## 9 3 24 3 4.118505Plotting the predicted concentrations for the three subjects shows how the clearance varies from one occasion to another.
print(ggplot(res$Cc) +geom_line(aes(x=time,y=Cc,color=id)))
2.2 Defining individual occasions
Individual occasions can also be defined in a data frame:
occ <- inlineDataFrame("
id time occ
1 0 1
1 12 2
1 24 3
2 0 1
2 24 2
3 0 1
")
res <-simulx(model = mymodel,
parameter = param,
treatment = trt,
varlevel = occ,
output = list(out.cc, out.param),
settings = list(seed=123))
print(res$parameter.iov)## id time occ Cl
## 1 1 0 1 4.996669
## 2 1 12 2 6.848626
## 3 1 24 3 4.190739
## 4 2 0 1 4.165440
## 5 2 24 2 3.157346
## 6 3 0 1 9.325054print(ggplot(res$Cc) +geom_line(aes(x=time,y=Cc,color=id)))
3 Combining IOV and time varying covariates
3.1 Simulating individual covariates
In this example, the weight wt is defined as a piecewise constant individual covariate: it changes from one occaion to another but remains constant within each occasion.
Then, any combination beetween time-varying/constant covariates and time-varying/constant parameters is possible
mymodel <- inlineModel("
[LONGITUDINAL]
input = {F, ka, V, Cl}
EQUATION:
Cc = pkmodel(p=F, ka, V, Cl)
[INDIVIDUAL]
input = {wt, age, wt_pop, age_pop}
input = {F_pop, omega_F, gamma_F, beta_F_age, ka_pop, omega_ka, beta_ka_wt,
V_pop, omega_V, Cl_pop, omega_Cl, gamma_Cl, beta_Cl_wt, beta_Cl_age}
EQUATION:
lwt = log(wt/wt_pop)
lage = log(age/age_pop)
DEFINITION:
F = {distribution=logitnormal, reference=F_pop, covariate=lage, coefficient=beta_F_age,
varlevel={id, id*occ}, sd={omega_F,gamma_F}}
ka = {distribution=lognormal, reference=ka_pop, covariate=lwt, coefficient=beta_ka_wt, sd=omega_ka}
V = {distribution=lognormal, reference=V_pop, sd=omega_V}
Cl = {distribution=lognormal, reference=Cl_pop,
covariate={lwt, lage}, coefficient={beta_Cl_wt,beta_Cl_age},
varlevel={id, id*occ}, sd={omega_Cl,gamma_Cl}}
[COVARIATE]
input = {wt_pop, omega_wt, gamma_wt, age_pop, omega_age}
DEFINITION:
wt = {distribution=normal, reference=wt_pop, varlevel={id, id*occ}, sd={omega_wt,gamma_wt}}
age = {distribution=normal, reference=age_pop, sd=omega_age}
")Three occasions are defined in this example:
occ <- list(name="occ", time=c(0,12,24))3 individuals are simulated with the following parameters and design:
param <- c(F_pop=0.8, omega_F=0.5, gamma_F=0.4, ka_pop=1, omega_ka=0.2,
V_pop=10, omega_V=0.2, Cl_pop=5, omega_Cl=0.2, gamma_Cl=0.2,
wt_pop=70, omega_wt=10, gamma_wt=0.1, age_pop=45, omega_age=5,
beta_F_age=0.5, beta_Cl_wt=1, beta_ka_wt=1, beta_Cl_age=0.4)
trt <- list(time=c(0,12,24), amount=100)
out.cc <- list(name= 'Cc', time=seq(0,36,by=0.1))
out.param <- list(name=c( 'ka', 'V', 'Cl'))
out.cov <- list(name=c('wt', 'age'))
res <-simulx(model = mymodel,
parameter = param,
treatment = trt,
varlevel = occ,
output = list(out.cc, out.param, out.cov),
group = list(size=3, level='covariate'),
settings = list(seed=12345))
print(names(res))## [1] "Cc" "group" "occasion" "treatment"
## [5] "parameter.iiv" "parameter.iov"Constant and time-varying parameters now include the covariates defined as outputs
print(res$parameter.iiv)## id V age ka0 Cl0 wt0
## 1 1 14.431164 45.41041 0.7339646 5.927781 61.35785
## 2 2 8.282027 49.25530 0.8538765 4.969817 57.54195
## 3 3 12.695709 44.84863 0.7449855 4.557900 72.29923print(res$parameter.iov)## id time occ Cl ka wt
## 1 1 0 1 4.055155 0.6418336 61.21323
## 2 1 12 2 6.264927 0.6437035 61.39157
## 3 1 24 3 4.531852 0.6425998 61.28631
## 4 2 0 1 4.442053 0.7030259 57.63341
## 5 2 12 2 4.142937 0.7010613 57.47235
## 6 2 24 3 4.420663 0.7023575 57.57862
## 7 3 0 1 4.312119 0.7693441 72.28877
## 8 3 12 2 4.898974 0.7696795 72.32029
## 9 3 24 3 2.853995 0.7692118 72.27634print(ggplot(res$Cc) +geom_line(aes(x=time,y=Cc,color=id)))
3.2 Using existing individual covariates
If, instead of simulating the individual covariates, existing individual covariates are used in the model, then the same model can be used by removing the section [COVARIATE]
mymodel <- inlineModel("
[LONGITUDINAL]
input = {F, ka, V, Cl}
EQUATION:
Cc = pkmodel(p=F, ka, V, Cl)
[INDIVIDUAL]
input = {wt, age, wt_pop, age_pop}
input = {F_pop, omega_F, gamma_F, beta_F_age, ka_pop, omega_ka, beta_ka_wt,
V_pop, omega_V, Cl_pop, omega_Cl, gamma_Cl, beta_Cl_wt, beta_Cl_age}
EQUATION:
lwt = log(wt/wt_pop)
lage = log(age/age_pop)
DEFINITION:
F = {distribution=logitnormal, reference=F_pop, covariate=lage, coefficient=beta_F_age,
varlevel={id, id*occ}, sd={omega_F,gamma_F}}
ka = {distribution=lognormal, reference=ka_pop, covariate=lwt, coefficient=beta_ka_wt, sd=omega_ka}
V = {distribution=lognormal, reference=V_pop, sd=omega_V}
Cl = {distribution=lognormal, reference=Cl_pop,
covariate={lwt, lage}, coefficient={beta_Cl_wt,beta_Cl_age},
varlevel={id, id*occ}, sd={omega_Cl,gamma_Cl}}
")Both time varying covariates and constant covariates should be defined in a data frame
covariate <- inlineDataFrame("
id time wt age
1 0 60 50
1 12 65 50
1 24 70 50
2 0 75 60
2 18 75 60
3 0 80 55
3 24 90 55
")Note that it is mandatory to also define occasions. These occasions should be consistent with those defined by the time varying covariates
occ <- inlineDataFrame("
id time occ
1 0 1
1 12 2
1 24 3
2 0 1
2 18 2
3 0 1
3 24 3
")res <-simulx(model = mymodel,
parameter = list(param,covariate),
treatment = trt,
varlevel = occ,
output = list(out.cc, out.param, out.cov),
settings = list(seed=12345))
print(res$parameter.iiv)## id V age ka0 Cl0
## 1 1 7.022771 50 1.2974183 5.082759
## 2 2 14.431164 60 0.7339646 5.927781
## 3 3 8.282027 55 0.8538765 4.969817print(res$parameter.iov)## id time occ Cl ka wt
## 1 1 0 1 3.822883 1.1120729 60
## 2 1 12 2 3.686424 1.2047456 65
## 3 1 24 3 5.671433 1.2974183 70
## 4 2 0 1 5.554202 0.7863907 75
## 5 2 18 2 8.555922 0.7863907 75
## 6 3 0 1 6.444115 0.9758589 80
## 7 3 24 3 7.221586 1.0978412 90