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

R scripts: replicate.R


We will use the following joint model for continuous PK and time-to-event data to illustrate this feature of Simulx:

modelPK <- inlineModel("
[LONGITUDINAL] 
input={V,Cl,alpha, beta,b}

EQUATION:
C = pkmodel(V, Cl)
h = alpha*exp(beta*C)
g = b*C

DEFINITION:
y = {distribution=normal, prediction=C, sd=g}
e = {type=event, maxEventNumber=1, rightCensoringTime=30, hazard=h}

;-----------------------------------------
[INDIVIDUAL]
input={V_pop,Cl_pop,omega_V,omega_Cl}

DEFINITION:
V     = {distribution=lognormal,   prediction=V_pop,    sd=omega_V}
Cl    = {distribution=lognormal,   prediction=Cl_pop,   sd=omega_Cl}
")

We can define the design and the parameters values as usual

adm  <- list(amount=100, time=0)
p <- c(V_pop=10, Cl_pop=1, omega_V=0.2, omega_Cl=0.2, alpha=0.02, beta=0.1, b=0.1)
g <- list(size=5, level='individual')

out.y <- list(name='y', time=seq(0,to=25,by=5))
out.e <- list(name='e', time=0)
out.p <- c("V", "Cl")
out   <- list(out.p, out.y, out.e)

The number of replicates is defined with the input argument nrep:

res1 <- simulx(model=modelPK, treatment=adm, parameter=p, output=out, group=g, nrep=3)
print(head(res1$parameter))
##   rep id         V        Cl
## 1   1  1 17.644951 0.8920649
## 2   1  2 10.829684 0.8960151
## 3   1  3  9.931543 0.7888289
## 4   1  4 10.583175 1.2900281
## 5   1  5 12.511587 0.9921788
## 6   2  1  8.777186 0.9528762
print(summary(res1$parameter))
##  rep   id          V                Cl        
##  1:5   1:3   Min.   : 6.051   Min.   :0.6459  
##  2:5   2:3   1st Qu.: 8.071   1st Qu.:0.8940  
##  3:5   3:3   Median : 9.901   Median :0.9882  
##        4:3   Mean   :10.047   Mean   :0.9830  
##        5:3   3rd Qu.:10.746   3rd Qu.:1.1094  
##              Max.   :17.645   Max.   :1.2900

When summary statistics are defined for some output, they are computed over each replicate

out.p2 <- list(name=c("V", "Cl"), FUN="quantile", probs=c(0.05, 0.5, 0.95)) 
out2   <- list(out.p, out.y, out.e)
res2 <- simulx(model=modelPK, treatment=adm, parameter=p, output=out2, group=g, nrep=3)
print(res2$parameter)
##    rep id         V        Cl
## 1    1  1 10.878496 1.0420109
## 2    1  2 12.341423 0.9757487
## 3    1  3 11.102384 1.0929885
## 4    1  4  8.181268 0.7755158
## 5    1  5  8.155356 0.9148135
## 6    2  1 10.772539 1.4129120
## 7    2  2  6.995855 0.9762706
## 8    2  3 14.091847 1.0377984
## 9    2  4 11.621393 1.1055963
## 10   2  5  7.291591 0.8617778
## 11   3  1 11.432198 0.9642706
## 12   3  2  8.423529 1.0038978
## 13   3  3  9.676487 1.0423061
## 14   3  4  8.219045 0.7509642
## 15   3  5  9.229826 0.9421577

Groups can be defined for each replicate:

g1 <- list(size=3, level="individual", treatment=list(amount=100, time=0))
g2 <- list(size=3, level="individual", treatment=list(amount=50,  time=0))
g  <- list(g1, g2)
res3 <- simulx(model=modelPK, parameter=p, output=out, group=g, nrep=2)
print(res3$parameter, digits=3)
##    rep id group     V    Cl
## 1    1  1     1 10.59 1.262
## 2    1  2     1 13.87 1.154
## 3    1  3     1 11.49 0.812
## 4    1  4     2  8.15 1.432
## 5    1  5     2  8.31 0.802
## 6    1  6     2  9.78 0.999
## 7    2  1     1  9.90 1.292
## 8    2  2     1  8.62 0.974
## 9    2  3     1 10.16 1.196
## 10   2  4     2 17.76 0.945
## 11   2  5     2 10.03 1.131
## 12   2  6     2 12.70 0.808
print(head(res3$e))
##   rep id group        time e
## 1   1  1     1  0.00000000 0
## 2   1  1     1  4.33002566 1
## 3   1  2     1  0.00000000 0
## 4   1  2     1 30.00000000 0
## 5   1  3     1  0.00000000 0
## 6   1  3     1  0.07264378 1