I'm implementing a Weibull survival analysis fitter, and have successfully estimated the parameters and their standard errors. I can also produce the fitted survival curve. My question is how can I use the fitted parameters and the standard errors to compute confidence bounds of the survival curve?
$\begingroup$
$\endgroup$
asked Jan 31, 2015 at 4:43
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
I just worked this out as I had the same question myself. My answer is largely thanks to Mara Tableman (p65 of Survival Analysis Using S/R).
Let's say you have the survival function $S(t)$ using the Weibull distribution, the lower and upper bounds of the 95% confidence interval can be calculated for every $t$ of the follow-up:
$$CI_{lo} = exp\left\{ log\left(S(t)\right)\times e^{z^*/\sqrt{n_u}}\right\}$$ $$CI_{hi} = exp\left\{ log\left(S(t)\right)\times e^{-z^*/\sqrt{n_u}}\right\}$$
where $z^*$ is the relevant quantile of the normal distribution ($z^*\approx1.96$ for a 95% CI) and $n_u$ is the number of events observed for the follow-up time (i.e., not the number of patients/subjects).
The below is R code to calculate the confidence bounds for a vector, named S_t, of survival proportions (a vector of $S(t)$ at specified follow-up times, $t$).
nu # number of events, needs to be set
z_st<-qnorm(0.975)
ci_lo <- exp(log(S_t)*exp(z_st/sqrt(nu)))
ci_hi <- exp(log(S_t)*exp(-z_st/sqrt(nu)))
$\begingroup$
$\endgroup$
I found that @tystanza's answer (Tableman's equations) give a too conservative confidence interval -- that is, they are too wide. When I did a Monte Carlo simulation to estimate the coverage of the Tableman confidence interval, I found that nominal 95% CI's gave 100% coverage.
I found another option from the package flexsurv. I don't think they directly expose the confidence intervals, but they can be found from "summary". Here is the hack workaround I found to get CI's. I did verify through MC sampling that the coverage of these CI's is just about right (I found the CI's contained the true value in 930/1000 cases for 95% confidence intervals -- pretty close).
require(flexsurv)
require(survival)
sWei = flexsurvreg(Surv(samples, rep(1,N)) ~ 1, dist='weibull', cl = 1-alpha)
s = unlist(summary(sWei))
lowerCI = s[(2*N+1):(3*N)]
upperCI = s[(3*N+1):(4*N)]
