R scripts: joint1.R ; joint2.R


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

A joint model is one that allows us to simultaneously describe the distribution of different types of observations made on the same individual.

Suppose that we have \(L\) different types of observation: \(y^{(1)}=(y_{j}^{(1)},1\leq j \leq n_{1})\), \(y^{(2)}=(y_{j}^{(2)},1\leq j \leq n_{2})\), …, \(y^{(L)}=(y_{j}^{(L)},1\leq j \leq n_{L})\), where \(n_{\ell}\) is the number of observations of type \(\ell\).

Note that the numbers of observations \((n_{\ell})\) and the observation times \((t_{j}^{(\ell)})\) may be different.

The joint model is the joint distribution of \(y = (y^{(1)}, y^{(2)}, \ldots, y^{(L)}).\)

This joint distribution is defined in the block DEFINITIONof the Section [LONGITUDINAL].


2 Examples

2.1 Joint continuous PKPD model

We consider in this example a joint model for continuous PK and PD data. We assume here that the predicted effect \(E\) is function of the predicted concentration \(C\):

\[ \begin{align} E(t) &= k_{\rm in}/k_{\rm out} \quad \text{for} \ t\leq 0 \\ \deriv{E} &= k_{\rm in}\left(1 - \frac{C(t)}{IC_{50}+C(t)}\right) - k_{\rm out}E(t) \quad \text{for} \ t\geq 0 \end{align} \]

We then assume that the measured concentrations \((y_{j}^{(1)})\) are log-normally distributed and the measured effects \((y_{j}^{(2)})\) normally distributed:

\[ \begin{align} \log(y_{j}^{(1)}) &= \log(C(t_{j}^{(1)})) + a_1 \varepsilon_{j}^{(1)}, \quad ; \quad 1 \leq j \leq n_1 \\ y_{j}^{(2)} &= f_2(t_{j}^{(2)}) + a_2 \varepsilon_{j}^{(2)}, \quad ; \quad 1 \leq j \leq n_2 \end{align} \]

This model is implemented in joint1.R:

joint.model1 <- inlineModel("
[LONGITUDINAL]
input = {ka, V, Cl, IC50, kin, kout, a1, a2}  

EQUATION:
C     = pkmodel(ka, V, Cl)
t0    = 0
E_0   = kin/kout
ddt_E = kin*(1 - C/(IC50+C)) - kout*E

DEFINITION:
Concentration = {distribution = lognormal, 
                 prediction   = C, 
                 sd           = a1}

Effect = {distribution = normal, 
          prediction   = E, 
          sd           = a2}
")

We use simulx for computing \(C\) and \(E\) every hour between 0h and 100h and for generating measured concentrations every 12 hours between time 4h and 100h and measured effects every 10 hours between time 5h and 100h.

A dose of 10 mg is administrated orally every 12 hours, starting at time 2h.

p <- c(ka=0.5, V=8, Cl=1.5, IC50=0.5, kin=10, kout=0.1, a1=0.1, a2=5)

a <- list(amount = 10, time = seq(2,120,by=12))

f <- list(name = c('C', 'E'),     time = seq(0,100,by=1))
c <- list(name = 'Concentration', time = seq(4,100,by=12))
e <- list(name = 'Effect',        time = seq(5,100,by=10))

res <- simulx(model     = 'model/joint1.txt',
              treatment = a,
              parameter = p,
              output    = list(f, c, e),
              settings  = list(seed = 1234))

Here, res is a list with four elements:

names(res)
## [1] "C"             "E"             "Concentration" "Effect"       
## [5] "treatment"

Let us plot this data:

library("gridExtra")
plot1 = ggplot() + geom_line(data=res$C, aes(x=time, y=C)) +
        geom_point(data=res$Concentration, aes(x=time, y=Concentration),colour="red") 
plot2 = ggplot() + geom_line(data=res$E, aes(x=time, y=E)) +
        geom_point(data=res$Effect, aes(x=time, y=Effect),colour="red") 
grid.arrange(plot1, plot2, ncol=2)


2.2 Joint PK and time-to-event data model

We consider now a joint model for continuous PK data and repeated events. The hazard \(h\) is function of the predicted concentration \(C\): \[ h(t) = u\, e^{v\,C(t)} \]

We also assume that the events are exactly observed until time 100h (right censoring time) This model is implemented in joint2.R:

joint.model2 <- inlineModel("
[LONGITUDINAL]
input = {ka, V, Cl, u, v, a1}  

EQUATION:
C = pkmodel(ka, V, Cl)
h = u*exp(v*C)

DEFINITION:
Concentration = {distribution = lognormal, 
                 prediction   = C, 
                 sd           = a1}

Hemorrhaging  = {type               = event, 
                 rightCensoringTime = 100,  
                 hazard             = h}
")

We use Simulx for computing \(C\) and \(h\) every hour between 0h and 100h and for generating measured concentrations every 12 hours between time 4h and 100h and repeated events, starting at time 0.

p <- c(ka=0.5, V=8, Cl=1.5, u=0.003, v=3, a1=0.05)

a <- list(amount = 10, time = seq(2,120,by=12))

f <- list(name = c('C', 'h'),     time = seq(0,100,by=1))
c <- list(name = 'Concentration', time = seq(4,100,by=12))
e <- list(name = 'Hemorrhaging',  time = 0)

res1 <- simulx(model     = joint.model2,
               treatment = a,
               parameter = p,
               output    = list(f, c, e),
               settings  = list(seed = 12345))

Let us see the simulated events:

print(res1$Hemorrhaging)
##        time Hemorrhaging
## 1   0.00000            0
## 2  86.95924            1
## 3  99.83272            1
## 4 100.00000            0

Indeed, in this example,

  • the experiment starts at time 0, which means that the first event will occur after time 0
  • there is one observed event at time 86.96h
  • there is one observed event at time 99.83h
  • the next event will occur after time 100

We can also plot the predicted and observed concentrations and the hazard function

plot1 = ggplot() + geom_line(data=res1$C, aes(x=time, y=C)) +
        geom_point(data=res1$Concentration, aes(x=time, y=Concentration),colour="red") 
plot2 = ggplot() + geom_line(data=res1$h, aes(x=time, y=h)) +
        ylab("hazard")  + theme(legend.position="none")
grid.arrange(plot1, plot2, ncol=2)

Instead of simulating data for a single individual, we can simulate data for \(N=20\) individuals with the same model, using the input argument group

N   <- 20
res2 <- simulx(model     = joint.model2,
               treatment = a,
               parameter = p,
               output    = list(f, c, e),
               group     = list(size = N, level='longitudinal'),
               settings  = list(seed = 121212))

We can plot the predicted and observed concentrations

plot3 = ggplot() + geom_line(data=res2$C, aes(x=time, y=C, group=id)) +
  geom_point(data=res2$Concentration, aes(x=time, y=Concentration),colour="red") 
plot(plot3)

We can also plot the observed events for the 20 individuals:

h1    = res2$Hemorrhaging[res2$Hemorrhaging[,3]==1,]
plot4 = ggplot()+geom_point(data=h1, aes(x=time,y=id), size=3) + 
        xlab("event time") + ylab("replicate")  + theme(legend.position="none")
plot(plot4)