R script: optimcode.R

$$\newcommand{\esp}{\mathbb{E}\left(#1\right)} \newcommand{\var}{\mbox{Var}\left(#1\right)} \newcommand{\deriv}{\dot{#1}(t)} \newcommand{\prob}{ \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}{{\rm arg}\min_{#1}} \newcommand{\argmax}{{\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}{\mathop{\longrightarrow}\limits_{#1}} \def\ka{k{a}} \def\ska{k{\scriptscriptstyle a}} \def\kel{k{e}} \def\skel{k{\scriptscriptstyle e}} \def\cl{C{\small l}} \def\Tlag{T\hspace{-0.1em}{lag}} \def\sTlag{T\hspace{-0.07em}{\scriptscriptstyle lag}} \def\Tk{T\hspace{-0.1em}{k0}} \def\sTk{T\hspace{-0.07em}{\scriptscriptstyle k0}} \def\thalf{t{1/2}} \newcommand{\Dphi}{\partial_\pphi #1} \def\asigma{a} \def\pphi{\psi} \newcommand{\stheta}{{\theta^\star}} \newcommand{\htheta}{{\widehat{\theta}}}$$

# 1 Introduction

1. When simulx is called for the first time with a given Mlxtran model, a C++ code is automatically generated from this model and automatically compiled. This process takes some time (a few seconds). If simulx is called again with the same model, this process is not performed anymore. The run is then much faster than the first one.

2. For each run of simulx, the design of the simulation must be created and loaded in memory. This process is also time consuming. When several calls of simulx are done successively, a test is automatically performed to check if the design has changed and should be loaded again, or if the design already in memory can be used. The user can also force simulx to create and load the design with the field load.design of settings that takes the values TRUE/FALSE.

3. simulx converts the inputs of the function defined in several lists into a unique list which can be processed by the C++ mlxLibrary. It is possible to bypass the execution of this R code and use directly the API (Application Programming Interface) of mlxLibrary instead of the simulx standard one.

All these techniques for reducing the computation time of simulx are shown to be extremly efficient for simulating a very large number of clinical trials.

# 2 Example

We will use in this example a quite complex PK model which requires to solve a system of ODEs (a two compartement model with transit compartments and non linear elimination).

Furthermore, each call of simulx will simulate PK data for 100 individuals.

myModel = inlineModel("
[LONGITUDINAL]
input = {Mtt, Ktr, ka, V, Vm, Km, k12, k21}
EQUATION:
C = pkmodel(Mtt, Ktr, ka, V, Vm, Km, k12, k21)
[INDIVIDUAL]
input = {V_pop, omega_V}
DEFINITION:
V = {distribution=lognormal, prediction=V_pop, sd=omega_V}
")

adm <- list(time=seq(0, 200, by=12), amount=100)
p   <- c(V_pop=10, omega_V=0.2, Mtt=2, Ktr=0.5, ka=1, Vm=10, Km=1, k12=0.5, k21=0.3)
C   <- list(name='C',time=seq(100, 200, by=10))
g   <- list(size = 100, level='individual')

Let us run simulx and display the computation time for this first run which requires the creation and the compilation of a C++ code:

ptm <- proc.time()
res <- simulx(model = myModel,parameter = p, output = C, treatment = adm, group = g)
print(proc.time() - ptm)
##    user  system elapsed
##    0.14    0.02    1.33

A second run of simulx with the same model will be much faster since the model is already compiled:

ptm <- proc.time()
res <- simulx(model = myModel,parameter = p, output = C, treatment = adm, group = g)
print(proc.time() - ptm)
##    user  system elapsed
##    0.11    0.00    0.03

Computation time becomes critical if we need to run simulx a large number of times (with different parameters values, or if we want to simulate several replicates of the same experience for instance).

As an illustration, let us run 500 times Simulx and display the total running time:

M=500

ptm <- proc.time()
for(i in seq(1,M)){
res <- simulx(model = myModel,
parameter = p,
output = C,
group = g)
}
print(proc.time() - ptm)
##    user  system elapsed
##   16.33    1.19   10.61

Simulx is quite fast beacause it uses the same compiled model for each run, but also because a test, automatically performed, avoids to load the same design at each run.

Forcing Simulx to load the design at each run significantly increases the computation time:

ptm <- proc.time()
for(i in seq(1,M)){
res <- simulx(model = myModel,
parameter = p,
output = C,
group = g,
}
print(proc.time() - ptm)
##    user  system elapsed
##   23.25    4.56   18.07

Simulxwill return the input of mlxLibrary as a list by setting the field data.in to TRUE:

dataIn <- simulx(model   = myModel,
parameter = p,
output    = C,
settings  = list(data.in=TRUE))
We can then bypass the R code by using Simulx with the input argument data:
ptm <- proc.time()
print(proc.time() - ptm)
##    user  system elapsed
##   14.68    0.87    7.40