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