$$ \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}}}$$

In this first example, we use the Monolix project iov_project.mlxtran.

iov.project <- 'monolixRuns/iov_project.mlxtran'

d <- readDatamlx(project = iov.project)
print(names(d))
##  [1] "populationParameters" "treatment"            "y1"                  
##  [4] "occasion"             "covariate.iov"        "id"                  
##  [7] "N"                    "catNames.iov"         "covariate.iiv"       
## [10] "catNames.iiv"

1, 2 or 3 occasions of length 12 hours are defined for each individual in this project:

head(d$occasion, 13)
##     id time occ
## 1    1    0   1
## 15   1   24   2
## 29   2    0   1
## 42   2   24   2
## 55   3    0   1
## 64   3   24   2
## 73   4    0   1
## 83   5    0   1
## 93   5   24   2
## 103  6    0   1
## 113  6   24   2
## 123  6   48   3
## 133  7    0   1

Covariates C1 and C3 do not change over time

head(d$covariate.iiv)
##   id C1 C3
## 1  1  5  A
## 2  2 10  A
## 3  3 15  A
## 4  4 20  A
## 5  5 25  A
## 6  6 30  A

while C2 and C4 change over time. Note that the column OCC which defines the occasions is also defined as a time varying covariate

head(d$covariate.iov,12)
##     id time C2 C4 OCC
## 1    1    0  7  A   1
## 15   1   24 12  B   2
## 29   2    0  9  A   1
## 42   2   24 14  B   2
## 55   3    0 11  A   1
## 64   3   24 16  B   2
## 73   4    0 13  A   1
## 83   5    0 15  A   1
## 93   5   24 20  B   2
## 103  6    0 17  A   1
## 113  6   24 22  B   2
## 123  6   48 27  B   3

The model used in this project assumes that individual parameters depend on constant and/or time varying covariates. Furthermore, IOV is assumed for parameters ka and V:

[COVARIATE]
input = {C1, C2, C3, C4, sC4, OCC, sOCC}
C3 = {type=categorical, categories={A, B}}
C4 = {type=categorical, categories={A, B}}
sC4 = {type=categorical, categories={A, A_B, A_B_B}}
OCC = {type=categorical, categories={1, 2, 3}}
sOCC = {type=categorical, categories={1, 1_2, 1_2_3}}

[INDIVIDUAL]
input = {ka_pop, beta_ka_C1, beta_ka_C2, beta_ka_C4_B, C1, C2, C4, omega_ka, gamma_ka, V_pop, beta_V_C2, beta_V_OCC_2, beta_V_OCC_3, OCC, omega_V, gamma_V, Cl_pop, beta_Cl_C1, beta_Cl_C3_B, C3, omega_Cl}
C4 = {type=categorical, categories={A, B}}
OCC = {type=categorical, categories={1, 2, 3}}
C3 = {type=categorical, categories={A, B}}

DEFINITION:
ka = {distribution=lognormal, typical=ka_pop, covariate={C1, C2, C4}, coefficient={beta_ka_C1, beta_ka_C2, {0, beta_ka_C4_B}}, varlevel={id, id*occ}, sd={omega_ka, gamma_ka}}
V = {distribution=lognormal, typical=V_pop, covariate={C2, OCC}, coefficient={beta_V_C2, {0, beta_V_OCC_2, beta_V_OCC_3}}, varlevel={id, id*occ}, sd={omega_V, gamma_V}}
Cl = {distribution=lognormal, typical=Cl_pop, covariate={C1, C3}, coefficient={beta_Cl_C1, {0, beta_Cl_C3_B}}, sd=omega_Cl}

[LONGITUDINAL]
input = {a, b}

file = 'lib:oral1_1cpt_kaVCl.txt'

DEFINITION:
y1 = {distribution=normal, prediction=Cc, errorModel=combined1(a, b)}

New PK data can be simulated using this model with the Monolix project:

res <- simulx(project=iov.project, output= list(name=c("ka","V","Cl","OCC")))
names(res)
## [1] "y1"            "occasion"      "treatment"     "originalId"   
## [5] "parameter.iiv" "parameter.iov" "population"
print(head(res$parameter.iiv))
##   id         Cl       ka0        V0
## 1  1 0.01999307 1.1096880 0.1765979
## 2  2 0.02469007 0.8044610 0.1001704
## 3  3 0.03802401 0.6276444 0.2095770
## 4  4 0.02122040 1.1802836 0.1248067
## 5  5 0.03699083 0.4156609 0.1440094
## 6  6 0.01890576 0.9660440 0.3776085
print(head(res$parameter.iov))
##    id time occ OCC        ka          V
## 1   1    0   1   2 1.1176357 0.20196827
## 2   1   24   2   3 1.2612656 0.18501740
## 27  2    0   1   2 0.7212083 0.09624956
## 28  2   24   2   3 0.7975725 0.09591072
## 52  3    0   1   2 0.6493514 0.19988634
## 53  3   24   2   3 0.8474597 0.26648521

We can then plot both the original and the simulated PK data:

y1 <- subset(d$y1, id %in% 1:10)
y1$out <- "Original data"
y2 <- subset(res$y1, id %in% 1:10)
y2$out <- "Simulated data"
y <- rbind(y1, y2)

pl <- ggplot(data=y, aes(x=time, y=y1)) + geom_point( aes(colour=id)) + geom_line( aes(colour=id)) +
  facet_wrap(~out) +  theme(legend.position="none") 

print(pl)