$$ \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 script: writeresults.R



1 Writing simulated data

1.1 Writing results in a single file

Let’s start with a PKPD model:

pkpd.model <- inlineModel("
[LONGITUDINAL]
input = {V, Cl, EC50, a1, a2}

EQUATION:
Cc = pkmodel(V, Cl)
E  = 100*Cc/(Cc+EC50)

DEFINITION:
y1 ={distribution=lognormal, prediction=Cc, sd=a1}
y2 ={distribution=normal,    prediction=E,  sd=a2}

[INDIVIDUAL]
input={V_pop,o_V,Cl_pop,o_Cl,EC50_pop,o_EC50}

DEFINITION:
V   ={distribution=lognormal, prediction=V_pop,   sd=o_V}
Cl  ={distribution=lognormal, prediction=Cl_pop,  sd=o_Cl}
EC50={distribution=lognormal, prediction=EC50_pop,sd=o_EC50}
")

p <- c(V_pop=10, o_V=0.1, Cl_pop=1, o_Cl=0.2, EC50_pop=3, o_EC50=0.2, a1=0.1, a2=1)
adm  <- list(amount=100, time=seq(0,50,by=12))
y1 <- list(name=c('y1'), time=seq(5,to=50,by=5))
y2 <- list(name='y2', time=seq(2,to=50,by=6))

We can save the outputs y1 and tt>y2 in a file “res1a.csv” (using the standard NONMEM/Monolix format), thanks to the additional input argument result.file:

res1 <- simulx(model     = pkpd.model,
              treatment = adm,
              parameter = p,
              group     = list(size=5, level="individual"),
              output    = list(y1,y2),
              result.file   ="res1a.csv",
              settings      = list(seed = 32323))

print(head(read.table("res1a.csv", header=T, sep=",")))
##   id time      y ytype amount
## 1  1    0      .     .    100
## 2  1    2   77.2     2      .
## 3  1    5 6.4569     1      .
## 4  1    8 64.444     2      .
## 5  1   10 3.6855     1      .
## 6  1   12      .     .    100

Remark: These outputs y1 and tt>y2 are not anymore outputs of simulx:

names(res1)
## [1] "group"

Results can be saved in various types of files, using different separators. The separator and the number of digits are additional settings:

res1 <- simulx(model     = pkpd.model,
              treatment = adm,
              parameter = p,
              group     = list(size=5, level="individual"),
              output    = list(y1,y2),
              result.file   ="res1a.txt",
              settings      = list(seed = 32323, sep="\t", digits=1))
print(head(read.table("res1a.txt", header=T, sep="\t")))
##   id time    y ytype amount
## 1  1    0    .     .    100
## 2  1    2 77.2     2      .
## 3  1    5  6.5     1      .
## 4  1    8 64.4     2      .
## 5  1   10  3.7     1      .
## 6  1   12    .     .    100

1.2 Writing results in separated files

Results can be saved in separated files instead of a single one:

res1 <- simulx(model     = pkpd.model,
              treatment = adm,
              parameter = p,
              group     = list(size=5, level="individual"),
              output    = list(y1,y2),
              result.folder ="res1a",
              settings      = list(seed = 32323))
print(list.files(path="res1a"))
## [1] "treatment.csv" "y1.csv"        "y2.csv"

Both the output file and the output folder can be defined

res1 <- simulx(model     = pkpd.model,
              treatment = adm,
              parameter = p,
              group     = list(size=5, level="individual"),
              output    = list(y1,y2),
              result.file   ="res1b.csv",
              result.folder ="res1b",
              settings      = list(seed = 32323))
print(head(read.table("res1b.csv", header=T, sep=",")))
##   id time      y ytype amount
## 1  1    0      .     .    100
## 2  1    2   77.2     2      .
## 3  1    5 6.4569     1      .
## 4  1    8 64.444     2      .
## 5  1   10 3.6855     1      .
## 6  1   12      .     .    100
print(list.files(path="res1b"))
## [1] "treatment.csv" "y1.csv"        "y2.csv"

1.3 Adding individual parameters

individual parameters can be saved as well in the data file

res1 <- simulx(model     = pkpd.model,
               treatment = adm,
               parameter = p,
               group     = list(size=5, level="individual"),
               output    = list(y1,y2, c("V" ,"Cl", "EC50")),
               result.file   ="res1c.csv",
               settings      = list(seed = 32323))
print(head(read.table("res1c.csv", header=T, sep=",")))
##   id time      y ytype amount      V      Cl   EC50
## 1  1    0      .     .    100 9.4832 0.90628 2.5527
## 2  1    2   77.2     2      . 9.4832 0.90628 2.5527
## 3  1    5 6.4569     1      . 9.4832 0.90628 2.5527
## 4  1    8 64.444     2      . 9.4832 0.90628 2.5527
## 5  1   10 3.6855     1      . 9.4832 0.90628 2.5527
## 6  1   12      .     .    100 9.4832 0.90628 2.5527

1.4 Writing results using writeDatamlx

Instead of saving the results ‘’within’’ simulx, it is possible to first use simulx as usual and then, write the results using the writeDatamlxfunction.

res2 <- simulx(model     = pkpd.model,
              treatment = adm,
              parameter = p,
              group     = list(size=5, level="individual"),
              output    = list(y1,y2),
              settings  = list(seed = 32323))

writeDatamlx(res2, result.file = "res2.csv")
print(head(read.table("res2.csv", header=T, sep=",")))
##   id time        y ytype amount
## 1  1    0        .     .    100
## 2  1    2  77.1999     2      .
## 3  1    5  6.45691     1      .
## 4  1    8 64.44424     2      .
## 5  1   10  3.68555     1      .
## 6  1   12        .     .    100
writeDatamlx(res2, result.folder = "res2")
print(list.files(path="res2"))
## [1] "treatment.csv" "y1.csv"        "y2.csv"


2 Using a Monolix project

2.1 Using a simplified format

project.file <- 'monolixRuns/theophylline_project.mlxtran'
res3 <- simulx(project=project.file, result.file="theo1.csv", settings=list(seed=12345))
print(head(read.csv("theo1.csv")))
##   id time      y amount
## 1  1 0.00      .    320
## 2  1 0.25 2.2614      .
## 3  1 0.57 4.4562      .
## 4  1 1.12 6.2658      .
## 5  1 2.02 8.8488      .
## 6  1 3.82 6.6862      .

2.2 Using the original format of the data

res4 <- simulx(project=project.file, result.file="theo2.csv", settings=list(format.original=TRUE,seed=12345))
print(head(read.csv("theo2.csv", sep="\t")))
##   ID AMT AMT.KG TIME   CONC WEIGHT SEX
## 1  1 320   4.02 0.00      .   79.6   M
## 2  1   .      . 0.25 2.2614   79.6   M
## 3  1   .      . 0.57 4.4562   79.6   M
## 4  1   .      . 1.12 6.2658   79.6   M
## 5  1   .      . 2.02 8.8488   79.6   M
## 6  1   .      . 3.82 6.6862   79.6   M