Plotting survival plots after brm

Interfaces brms
1

Hello, is there any code available to plot survival curves from a model with brms, based on covariates?

2

I am not aware of code for this. There is a longstanding issue in brms to implement survival curves though: https://github.com/paul-buerkner/brms/issues/232

3

I thought the code could be implemented relatively easily and would be a nice improvement for medical statistics. So I wanted to ask if this is correct, particularly the credibility interval for survival rates. I wish someone could help, as we could leave the code written for other users to plot the curves. So, if possible, I would ask for some feedback.
For example, using a multivariable lognormal model with brm, the fitted function could be used to estimate the linear predictor. For example:

   fit<-brm_multiple(SGm | cens(1-Die)  ~ treat+X1+X2+X3, prior=prior, data = imp,family = lognormal(), 
    iter=2000)

newdata = as.data.frame( expand.grid(
      treat=c(0,1), 
      X1=c(1),
      X2=c(1),
      X3=TRUE )) 

    fi <- data.frame(fitted(fit, newdata=newdata, scale="linear",summary = FALSE))
    LETTERS2<-c(LETTERS [1:2])
    colnames(fi) <- LETTERS2
    fi$Effect <- exp( fi$B - fi$A)
    fi$y<-"Effect"

After that, survival could be estimated by applying the equation of the lognormal model. My doubt is whether or not for the survival credible interval it would be correct to use the upper and lower limits of the linear predictor credible interval, as follows.

tt<-seq(1,36,1)
st<-1 - pnorm((log(tt) - median_qi(fi$B)$y) / sigma)
sti<-1 - pnorm((log(tt) - median_qi(fi$B)$ymin) / sigma)
stu<-1 - pnorm((log(tt) - median_qi(fi$B)$ymax) / sigma)
sd<-1 - pnorm((log(tt) - median_qi(fi$A)$y) / sigma)
sdi<-1 - pnorm((log(tt) - median_qi(fi$A)$ymin) / sigma)
sdu<-1 - pnorm((log(tt) - median_qi(fi$A)$ymax) / sigma)

After that it should be trivial plotting the survival curve with something like this. The geom_ribbon function of ggplot2 would paint the confidence band, but I’m not sure if it’s correct. Visually it looks a bit wide to me.

datst<-data.frame(sur=st, sui=sti, suu=stu, tt=tt, reg="Treatment 1")
datsd<-data.frame(sur=sd, sui=sdi, suu=sdu, tt=tt, reg="Treatment 2")
dat<-rbind(datst,datsd)
dat1<-dat [dat$reg=="Treatment 1",]
dat2<-dat [dat$reg=="Treatment 2",]
mycolors<-c("steelblue3", "black")

s<-ggplot(dat,aes(x=tt,y=sur,group=reg))+
  geom_line(aes(colour=factor(reg),linetype = factor(reg)), size=1.25, 
  show.legend=FALSE) +
  xlab("Months")+
  ylab("Survival probability")+
  geom_ribbon(aes(ymin=sui, ymax=suu, fill = factor(reg)), alpha = 0.1)+    
  theme_minimal_grid() +theme(text = element_text(size=18))+
  scale_fill_manual("",values = mycolors)+
  scale_color_manual("",values = mycolors) +
  scale_x_continuous(breaks =seq(0,36,6))

Median OS to add in the text could be:
text<-paste(“Median OS with DPF:”,round(as.numeric(median_qi(exp(fi$B))[1]),1),“months”)
text2<-paste(“Median OS with PF:”,round(as.numeric(median_qi(exp(fi$A))[1]),1),“months”)

I added an eye plot and obtained this:

dibujo prueba
dibujo prueba500×750 97.8 KB

But not sure if the credible bands are accurate.

4

I’ve been working through the survival analysis sections of Singer and Willett’s (2003) text with brms here. Chapters 9 and above are the most relevant.

5

Has this been implemented or is @Solomon’s great resource on Survival analysis still the way to go. A Bayesian Kaplan-Meier curve is implemented in 13.3 but I cannot find anything about whether this is easier now the cox family has been added to brms for a couple of years?

Could you use bayesplot::ppc_km_overlay() (PPC censoring — PPC-censoring • bayesplot) and then calculate the credible intervals at each time point to estimate the Kaplain meier curve? It does not look like it would do exactly as the KM curve from using survminer::ggsurvplot() for example.