# 1 Plotting several prediction distributions on a same plot

A R function for plotting several prediction distributions on a same plot:

plot.prediction.intervals <- function(r, plot.median=TRUE, level=90, labels=NULL,
legend.title=NULL, colors=NULL) {
P <- prctilemlx(r, number=1, level=level, plot=FALSE)
if (is.null(labels))  labels <- levels(r$group) if (is.null(legend.title)) legend.title <- "group" names(P$y)[2:4] <- c("p.min","p50","p.max")
pp <- ggplot(data=P$y)+ylab(NULL)+ geom_ribbon(aes(x=time,ymin=p.min, ymax=p.max,fill=group),alpha=.5) if (plot.median) pp <- pp + geom_line(aes(x=time,y=p50,colour=group)) if (is.null(colors)) { pp <- pp + scale_fill_discrete(name=legend.title, breaks=levels(r$group),
labels=labels)
pp <- pp + scale_colour_discrete(name=legend.title,
breaks=levels(r$group), labels=labels, guide=FALSE) } else { pp <- pp + scale_fill_manual(name=legend.title, breaks=levels(r$group),
labels=labels,
values=colors)
pp <- pp + scale_colour_manual(name=legend.title,
breaks=levels(r$group), labels=labels, guide=FALSE,values=colors) } return(pp) } Simulate first data from a joint model (PK + TTE) defining 3 treatment arms: jointModel <- inlineModel(" [LONGITUDINAL] input={Tk0, V, Cl, alpha, beta, a} EQUATION: C = pkmodel(Tk0, V, Cl) h = exp(-alpha - beta*C) H_0 = 0 ddt_H = h S = exp(-H) DEFINITION: y = {distribution=lognormal, prediction=C, sd=a} e = {type=event, maxEventNumber=1, rightCensoringTime=200, hazard=h} [INDIVIDUAL] input={Tk0_pop,V_pop,Cl_pop,alpha_pop,omega_Tk0,omega_V,omega_Cl,omega_alpha,a} DEFINITION: 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} ") #----------------------------------------------- N=100 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 print(prctilemlx(res1$S))

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 plot.S1 <- plot.prediction.intervals(res1$S,
labels       = c("placebo","25 mg","50mg"),
legend.title = "arm")
plot.S <- plot.S1  + ylab("Survival prediction interval") + theme(legend.position=c(0.9,0.8))
print(plot.S1)

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)) print(plot.S2) 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))
print(plot.S3)

# 2 Plotting several prediction intervals of a Kaplan Meier plot

A R function for plotting several prediction intervals of a Kaplan Meier plot on a same plot:

plot.survival.intervals <- function(r, plot.median=TRUE, level=90, labels=NULL,
legend.title=NULL, colors=NULL) {
n <- 100
surv.time <- seq(min(r$time), max(r$time), length.out=n)
km <-  kmplotmlx(r, time=surv.time, plot=F)$surv attr(km,"name") <- "S" n.groups <- nlevels(km$group)
n.rep    <- nlevels(km$rep) n.id <- n.groups*n.rep km$id    <- factor(rep((1:n.id), each=n))
pp <- plot.prediction.intervals(km, plot.median=plot.median, level=level, labels=labels,
legend.title=legend.title, colors=colors)
return(pp)
}

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

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))
print(plot.SPI)