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

Plotting several prediction distributions on a same plot

Simulate first data from a joint model (PK + TTE) defining 3 treatment arms:

jointModel <- inlineModel("
input={Tk0, V, Cl, alpha, beta, a}

C = pkmodel(Tk0, V, Cl)
h = exp(-alpha - beta*C)
H_0 = 0
ddt_H = h
S = exp(-H)
y = {distribution=lognormal, prediction=C, sd=a}
e = {type=event, maxEventNumber=1, rightCensoringTime=200, hazard=h}

Tk0   = {distribution=lognormal,   prediction=Tk0_pop,  sd=omega_Tk0}
V     = {distribution=lognormal,   prediction=V_pop,    sd=omega_V}
Cl    = {distribution=lognormal,   prediction=Cl_pop,   sd=omega_Cl}
alpha = {distribution=lognormal,   prediction=alpha_pop,sd=omega_alpha}

t.time <- seq(24,120,by=24)
g1 <- list(size=N, level='individual', treatment=list(time=t.time, amount=0))
g2 <- list(size=N, level='individual', treatment=list(time=t.time, amount=25))
g3 <- list(size=N, level='individual', treatment=list(time=t.time, amount=50))

pop.param   <- c(
  Tk0_pop   = 3,    omega_Tk0   = 0.3,
  V_pop     = 10,   omega_V     = 0.3,
  Cl_pop    = 1,    omega_Cl    = 0.3,
  alpha_pop = 6,    omega_alpha = 0.01,
  beta      = 0.4,  a           = 0.1

f <- list(name=c('C','S'), time=seq(0,to=180,by=3))
e <- list(name="e", time=0)

res1 <- simulx(model     = jointModel,
               parameter = pop.param,
               group     = list(g1, g2, g3),
               output    = list(f,e))

prctilemlx displays the 3 prediction intervals on 3 subplots


or on the same plot:

print(prctilemlx(res1$S, facet=FALSE))

Alternatively, the function plot.prediction.intervals can be used if we want to display these prediction intevals with different colors.

Remark: the R function plot.prediction.intervals is not part of the mlxR package. This is just an example of what can be done: this function can be freely used/modified by anyone.

Show the R code
plot.S1 <- plot.prediction.intervals(res1$S, 
                                    labels       = c("placebo","25 mg","50mg"), 
                                    legend.title = "arm")
plot.S1 <- plot.S1  + ylab("Survival prediction interval") + theme(legend.position=c(0.9,0.8))

The median can be removed

plot.S2 <- plot.prediction.intervals(res1$S, 
                                     labels       = c("placebo","25 mg","50mg"), 
                                     legend.title = "arm",
                                     plot.median  = FALSE)
plot.S2 <- plot.S2  + ylab("Survival prediction interval") + theme(legend.position=c(0.9,0.8))

and the colors can be defined by the user

plot.S3 <- plot.prediction.intervals(res1$S, 
                                     labels       = c("placebo","25 mg","50mg"), 
                                     legend.title = "arm",
                                     colors       = c('#01b7a5', '#c17b01', '#a00159'))
plot.S3 <- plot.S3  + ylab("Survival prediction interval") + theme(legend.position=c(0.9,0.8))

Plotting several prediction intervals of a Kaplan Meier plot

We now simulate 100 replicates of a same trial (with only 2 arms in this example)

res2 <- simulx(model     = jointModel,
               parameter = pop.param,
               group     = list(g1,g3),
               nrep      = 100,
               output    = e)

Prediction intervals of the Kaplan Meier estimate can then be computed for each group and displayed using these 100 replicates using the R function plot.survival.intervals

Remark: the R function plot.survival.intervals is not part of the mlxR package.

Show the R code
plot.SPI <- plot.survival.intervals(res2$e, labels=c("placebo", "100 mg"), level=80, legend.title="arm") + 
  ylab("Survival prediction interval") + theme(legend.position=c(0.9,0.8))