R scripts: groupI.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

It is possible to define groups for several purposes.

  1. We can define several groups (or arms) with different values of the parameters, different values of the regression variables, different treatments and/or different outputs (see the example below).

  2. We can create groups that consist of several replicates of the longitudinal data, several individuals sharing or not the same covariates, and/or several populations.

  3. We can combine these features and create groups of different sizes, with different levels of randomization and different parameters, treatments, outputs, … (see the examples of the next article)


2 Example

We will use in this example a very simple PK model for iv administration:

pk.model <- inlineModel("
[LONGITUDINAL]
input = {V, k}
EQUATION:
C = pkmodel(V,k)
")

We create here three groups with the same PK parameters and the same output, but with three different dosage regimens: 50mg every 6 hours, 100mg every 12 hours, and 150mg every 18 hours.

treatment is not anymore an input argument of simulx since it is defined in each element of the list group:

adm1 <- list(time=seq(0,to=66,by=6), amount=50)
adm2 <- list(time=seq(0,to=66,by=12), amount=100)
adm3 <- list(time=seq(0,to=66,by=18), amount=150)
g1 <- list(treatment=adm1);
g2 <- list(treatment=adm2);
g3 <- list(treatment=adm3);

C <- list(name='C', time=seq(0, 100, by=1))
p <- c(V=10, k=0.2)


res1 <- simulx(model     = pk.model, 
               parameter = p, 
               output    = C, 
               group     = list(g1,g2,g3))

Here, res1 has a unique element C: res1$C is a data frame with three columns, id, time and C where id takes the values 1, 2 and 3:

print(res1$C[1:4,])
##   id group time        C
## 1  1     1    0 5.000000
## 2  1     1    1 4.093654
## 3  1     1    2 3.351600
## 4  1     1    3 2.744058
print(res1$C[102:105,])
##     id group time         C
## 102  2     2    0 10.000000
## 103  2     2    1  8.187308
## 104  2     2    2  6.703200
## 105  2     2    3  5.488116
print(res1$C[203:206,])
##     id group time         C
## 203  3     3    0 15.000000
## 204  3     3    1 12.280961
## 205  3     3    2 10.054801
## 206  3     3    3  8.232175

We can plot the concentration \(C\) using one colour per group (i.e. per individual):

print(ggplot(data=res1$C, aes(x=time, y=C, colour=id)) + geom_line(size=1))

We can instead define the same treatment, but different PK parameter values for the three groups. In this case, parameter is not anymore an input argument of simulx since it is defined in each element of the list group:

adm <- list(time=seq(0,to=66,by=12), amount=100)

g1 <- list(parameter=c(V=10, k=0.2))
g2 <- list(parameter=c(V=15, k=0.1))
g3 <- list(parameter=c(V=20, k=0.05))

res2 <- simulx(model     = pk.model, 
               treatment = adm, 
               output    = C, 
               group     = list(g1,g2,g3))

print(ggplot(data=res2$C, aes(x=time, y=C, colour=id)) + geom_line(size=1))

We can instead define different outputs for the three groups:

adm <- list(time=seq(0,to=66,by=12), amount=100)
p <- c(V=10, k=0.2)
C1 <- list(name='C', time=seq(0,  25, by=1))
C2 <- list(name='C', time=seq(30, 55, by=0.5))
C3 <- list(name='C', time=seq(60,100, by=0.25))
g1 <- list(output=C1)
g2 <- list(output=C2)
g3 <- list(output=C3)

res3 <- simulx(model     = pk.model, 
               treatment = adm, 
               parameter = p, 
               group     = list(g1,g2,g3))

print(ggplot(data=res3$C, aes(x=time, y=C, colour=id)) + geom_line(size=1))

It is also possible to combine these different features by defining different treatments, different outputs and different parameter values in the three groups:

adm1 <- list(time=seq(0,to=66,by=6), amount=50)
adm2 <- list(time=seq(0,to=66,by=12), amount=100)
adm3 <- list(time=seq(0,to=66,by=18), amount=150)
p1 <- c(V=10, k=0.2)
p2 <- c(V=15, k=0.1)
p3 <- c(V=20, k=0.05)
C1 <- list(name='C', time=seq(0,  75, by=1))
C2 <- list(name='C', time=seq(20, 90, by=0.5))
C3 <- list(name='C', time=seq(25,100, by=0.25))
g1 <- list(output=C1, treatment=adm1, parameter=p1)
g2 <- list(output=C2, treatment=adm2, parameter=p2)
g3 <- list(output=C3, treatment=adm3, parameter=p3)

res4 <- simulx(model     = pk.model, 
               group     = list(g1,g2,g3))

print(ggplot(data=res4$C, aes(x=time, y=C, colour=id)) + geom_line(size=1))