$$ \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: statout.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=100, level='individual')
s <- list(seed=123456)
out.y <- list(name='y', time=seq(0,to=25,by=5))
out.e <- list(name='e', time=0)

We can compute the mean and standard deviation of the individual parameters \(V\) and \(Cl\), compute over the \(N=100\) individuals:

out.p <- list(name=c("V", "Cl"), FUN= c("mean", "sd"))

out   <- list(out.p, out.y, out.e)
res1a <- simulx(model=modelPK, treatment=adm, parameter=p, output=out, group=g, settings=s)

res1a$parameter
##     V.mean  Cl.mean     V.sd     Cl.sd
## 1 10.30865 1.051291 2.046022 0.2146551

Available functions are mean, sd, var, median, quantile.

Let us compute now the percentiles of order 5%, 50% and 95%:

out.p <- list(name=c("V", "Cl"), FUN="quantile", probs=c(0.05, 0.5, 0.95)) 
out   <- list(out.p, out.y, out.e)
res1b <- simulx(model=modelPK, treatment=adm, parameter=p, output=out, group=g, settings=s)
res1b$parameter
##       V.p5    V.p50    V.p95     Cl.p5   Cl.p50  Cl.p95
## 1 7.541696 9.845827 13.98242 0.7404126 1.041992 1.37721

It is possible to compute several functions with a single call:

out.p <- list(name=c("V", "Cl"),  FUN=c("mean","quantile"), probs=c(0.05,0.95)) 
out   <- list(out.p, out.y, out.e)
res1c <- simulx(model=modelPK, treatment=adm, parameter=p, output=out, group=g, settings=s)
res1c$parameter
##     V.mean  Cl.mean     V.p5    V.p95     Cl.p5  Cl.p95
## 1 10.30865 1.051291 7.541696 13.98242 0.7404126 1.37721

Remark: we would get the same results by running first Simulx with \(V\) and \(Cl\) as outputs, and computing then the summary statistics using statmlx:

out.p <- c("V", "Cl")
out   <- list(out.p, out.y, out.e)
res1d <- simulx(model=modelPK, treatment=adm, parameter=p, output=out, group=g, settings=s)
statmlx(res1d$parameter, FUN=c("mean","quantile"), probs=c(0.05,0.95))
##     V.mean  Cl.mean     V.p5    V.p95     Cl.p5  Cl.p95
## 1 10.30865 1.051291 7.541696 13.98242 0.7404126 1.37721

Summary statistics can also be computed for the observations. In this example, we compute

out.p <- list(name=c("V", "Cl"),  FUN=c("mean","quantile"), probs=c(0.05,0.95)) 
out.y <- list(name='y', time=seq(0,to=25,by=5), FUN = c("mean", "sd", "quantile"), probs = c(0.05, 0.95))
out.e <- list(name='e', time=0, type="event", surv.time=c(10,20))
out   <- list(out.p, out.y, out.e)
res2 <- simulx(model=modelPK, treatment=adm, parameter=p, output=out, group=g, settings=s)
res2$y
##   time     y.mean      y.sd      y.p5     y.p95
## 1    0 10.0438921 2.1153454 6.5011208 13.452909
## 2    5  5.8266127 1.0655046 4.2174192  7.640414
## 3   10  3.4799493 0.7009569 2.3488558  4.544024
## 4   15  2.1444653 0.7108819 1.1351004  3.266180
## 5   20  1.3189945 0.5670839 0.4787293  2.278551
## 6   25  0.8440965 0.4477673 0.2153766  1.737770
res2$e
##   nbEv.mean S10.mean S20.mean
## 1      0.63     0.63     0.45

nbEv.mean is the mean number of events per individual during the whole study (i.e. between time 0 and time 30).

When replicates and/or groups are defined, the summary statistics are computed over each group of each replicate.

g1 <- list(size=100, level="individual", treatment=list(amount=100, time=0))
g2 <- list(size=100, level="individual", treatment=list(amount=50,  time=0))
g  <- list(g1, g2)
res3a <- simulx(model=modelPK, parameter=p, output=out, group=g, nrep=3, settings=s)
res3a$parameter
##   rep group    V.mean  Cl.mean     V.p5    V.p95     Cl.p5   Cl.p95
## 1   1     1 10.308647 1.051291 7.541696 13.98242 0.7404126 1.377210
## 2   1     2 10.463900 1.039790 7.385356 14.15422 0.7002062 1.445878
## 3   2     1 10.383206 1.003870 7.654567 13.22531 0.7272242 1.329500
## 4   2     2 10.256424 1.028617 7.446694 13.12068 0.6726915 1.438391
## 5   3     1  9.933524 1.067755 6.782763 13.18811 0.7077879 1.499508
## 6   3     2  9.981329 1.059125 6.902989 13.00620 0.7276924 1.491159
res3a$e
##   rep group nbEv.mean S10.mean S20.mean
## 1   1     1      0.63     0.63     0.45
## 2   1     2      0.60     0.74     0.53
## 3   2     1      0.50     0.77     0.64
## 4   2     2      0.48     0.76     0.61
## 5   3     1      0.53     0.66     0.53
## 6   3     2      0.52     0.72     0.60

When the simulated data is saved in a data file, the complete data is saved and the summary statistics are returned by Simulx:

res3b <- simulx(model=modelPK, parameter=p, output=out, group=g, nrep=3, settings=s,
                result.file = "res3b.csv")
head(read.table("res3b.csv", header=T, sep=","))
##   rep id group time       y ytype amount     V   Cl
## 1   1  1     1    0       .     .    100 8.817 1.64
## 2   1  1     1    0  9.8693     1      . 8.817 1.64
## 3   1  1     1    0       0     2      . 8.817 1.64
## 4   1  1     1    5  5.0853     1      . 8.817 1.64
## 5   1  1     1   10  1.9722     1      . 8.817 1.64
## 6   1  1     1   15 0.71898     1      . 8.817 1.64
res3b$parameter
##   rep group    V.mean  Cl.mean     V.p5    V.p95     Cl.p5   Cl.p95
## 1   1     1 10.308647 1.051291 7.541696 13.98242 0.7404126 1.377210
## 2   1     2 10.463900 1.039790 7.385356 14.15422 0.7002062 1.445878
## 3   2     1 10.383206 1.003870 7.654567 13.22531 0.7272242 1.329500
## 4   2     2 10.256424 1.028617 7.446694 13.12068 0.6726915 1.438391
## 5   3     1  9.933524 1.067755 6.782763 13.18811 0.7077879 1.499508
## 6   3     2  9.981329 1.059125 6.902989 13.00620 0.7276924 1.491159
res3b$e
##   rep group nbEv.mean S10.mean S20.mean
## 1   1     1      0.63     0.63     0.45
## 2   1     2      0.60     0.74     0.53
## 3   2     1      0.50     0.77     0.64
## 4   2     2      0.48     0.76     0.61
## 5   3     1      0.53     0.66     0.53
## 6   3     2      0.52     0.72     0.60