R script: analytical.R

Mlxtran code: model/analytical.txt

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

1 Introduction

Functions of time are implemented in a block EQUATION of the section [LONGITUDINAL] of a script Mlxtran.

We consider parametric functions of time of the form \(f(t;\phi)\), where \(\phi\) is a vector of structural parameters.

The [LONGITUDINAL] section starts with the definition of the vector of input parameters \(\phi\) in a list input.

Then, the simulx script should provide

Results are returned as a list of data frames.


2 Example

We consider in this example a model where two functions \(f_1\) and \(f_2\) depending on a vector of parameters \(\phi=(ka,V,k)\) are defined:

\[ \begin{aligned} f_1(t) &= \frac{D}{V}e^{-k\, t} \\ f_2(t) &= \frac{D\, k_a}{V(k_a - k)}(e^{-k\, t} - e^{-k_a \, t}) \end{aligned} \] where \(D=100\) is given.

This model is implemented in the model file model/analytical.txt:

[LONGITUDINAL]
input = {ka, V, k}

EQUATION:
D  = 100
f1 = D/V*exp(-k*t)
f2 = D*ka/(V*(ka-k))*(exp(-k*t) - exp(-ka*t))

We will start evaluating only \(f_1\) every 0.1h between time 0 and time 25, with \(k_a=0.5\), \(V=10\) and \(k=0.2\):

f <- list(name = 'f1',  time = seq(0, 25, by=0.1))

p <- c(ka = 0.5, V = 10, k = 0.2)

res <- simulx(model     = 'model/analytical.txt', 
              parameter = p, 
              output    = f)

res is a list with element f1 which is a data frame with 251 rows and 2 columns:

res$f1[1:5,]
##   time        f1
## 1  0.0 10.000000
## 2  0.1  9.801987
## 3  0.2  9.607894
## 4  0.3  9.417645
## 5  0.4  9.231163
res$f1[248:251,]
##     time         f1
## 248 24.7 0.07154598
## 249 24.8 0.07012928
## 250 24.9 0.06874063
## 251 25.0 0.06737947

we can plot this data using ggplot

print(ggplot(data=res$f1) + geom_line(aes(x=time, y=f1)))

We way now want evaluate both \(f_1\) and \(f_2\) at the same time points. We can still define the output in a single list, but the field name has now two elements f1 and f2:

f <- list(name = c('f1','f2'), time = seq(0, 25, by=0.1))

res <- simulx(model     = 'model/analytical.txt', 
              parameter = p, 
              output    = f)

print(ggplot() + geom_line(data=res$f1, aes(x=time, y=f1), color="black") +
                geom_line(data=res$f2, aes(x=time, y=f2), color="red") +
                ylab('concentration'))

It can be useful to reshape the data before plotting it with ggplot. We merge the two outputs into a unique data frame r with three columns: time, f (which takes the values f1 or f2) and value.

library(reshape2)
r <- merge(res$f1,res$f2)
r <- melt(r ,  id = 'time', variable.name = 'f', value.name="concentration")
r[c(1:4),]
##   time  f concentration
## 1  0.0 f1     10.000000
## 2  0.1 f1      9.801987
## 3  0.2 f1      9.607894
## 4  0.3 f1      9.417645
r[c(250:253),]
##     time  f concentration
## 250 24.9 f1    0.06874063
## 251 25.0 f1    0.06737947
## 252  0.0 f2    0.00000000
## 253  0.1 f2    0.48282081

This format is now suitable for ggplot:

print(ggplot(r, aes(time,concentration)) + geom_line(aes(colour = f),size=1) +
        guides(colour=guide_legend(title=NULL)) + theme(legend.position=c(.9, .8)))

It is also possible to evaluate \(f_1\) and \(f_2\) at different time points. The two outputs should be then defined in two separate lists:

f1 <- list(name = 'f1', time = seq(0, 15, by=0.1))
f2 <- list(name = 'f2', time = seq(10, 25, by=1))

res <- simulx(model     = 'model/analytical.txt', 
              parameter = p, 
              output    = list(f1,f2))

plot(ggplot() + geom_line(data=res$f1, aes(x=time, y=f1), color="black") +
                geom_line(data=res$f2, aes(x=time, y=f2), color="red") +
                ylab('concentration'))

We may again prefer to reshape the data before plotting it with ggplot

r <- merge(res$f1,res$f2,all=TRUE)
r <- melt(r ,  id = 'time', variable.name = 'f', value.name="concentration")
r=r[(!is.na(r$concentration)),]
print(ggplot(r, aes(time,concentration)) + geom_line(aes(colour = f),size=1) +
       guides(colour=guide_legend(title=NULL)) + theme(legend.position=c(.9, .8)))