Hi all,

I'd like to plot shaded confidence bands for a parametric survival curve fitted with survreg.

It's easy to plot dashed confidence intervals:

library(survival)
library(flexsurv)
# Generate data
dat = data.frame(s = rgamma(n=50, shape=9, rate=.5), event=c(rep(1,40),rep(0,10)))
# Fit a Kaplan-Meier
km <- survfit(Surv(s, event=event) ~ 1, data = dat)
plot(km, conf.int=F, main="", ylab="Probability of Survival", xlab="Overall Survival")
# Fit the parametric curve
f_OS = flexsurvreg(Surv(s, event=event) ~ 1, data=dat, dist="lnorm");
lines(f_OS, t = seq(0,40, by=0.1), ci=T)

However, I'd like the confidence bands to be shaded instead. The code below is the closest I've gotten; It plots 95% confidence intervals based on quantiles and their SEs at each time point (quantile +/- 1.96*SE). This might arguably be an acceptable approximation, however, I'd like the upper and lower confidence limits to be equidistant from each point survival proportion (like those obtained when using lines.flexsurvreg), instead of around each quantile.

f_OS = survreg(Surv(s, event=event) ~ 1, data=dat, dist="lognormal");
lines(predict(f_OS, type="quantile",p=seq(.999,.001,by=-.001))[1,],1-seq(.999,.001,by=-.001),lwd=2, col=c("darkblue"));
newx <- seq(0.0001, 0.999, length.out=100)
preds <- predict(f_OS, type="quantile", p = newx, se.fit=T)
q=preds[[1]][1,]
SEs=preds[[2]][1,]
polygon(c(rev(q-1.96*SEs), q+1.96*SEs), 1-c(rev(newx), newx), col=rgb(1, 0, 0,0.5), border = NA)

Any help is much appreciated! Thank you