Good morning, I have a question about how to calculate the standard error and confidence intervals for predicting the probability of survival at a point in time, for a multivariate prognostic model. I have recently fitted an AFT model (lognormal), which allows to calculate survival probabilities.

The E code would be something like:

fit <- psm(Surv(d.time,death) ~ X1,X2,X3, dist=“lognormal”, x=TRUE, y=TRUE)

Thus, applying the equation of the AFT model it is possible to estimate a survival probability, in subjects with variables X1, X2 and X3.

The equations for the survival function of the AFT model are easily obtained in the rms package, for example with Function(fit) or latex(fit).

To estimate confidence intervals, I have used the percentiles of 1000 replications per bootstrap (boot package). This method works very well, but is somewhat slow to perform to get individual predictions. Since the website would need to access an excel document with 1000 different models, and perform 3x1000 calculations, that calculator would take some time.

Besides, it doesn’t allow me to draw a survival function with a confidence band, according to the patient’s covariates, because it takes too long.

To solve it I would like to obtain an equation for the confidence intervals as a function of time by the delta method. This can easily be done using the R functions (e., predict with conf.int=T), but I need to do it by hand, to get it out of R, to a web-based calculator.

Professor Harrell explains how to apply the delta-method to an AFT model in Regression Modelling Strategies, chapter 18, page 439. But it is still complicated to implement by hand the formula 18.40: Var(f) = F V F '.

Since estimating model errors is an important issue, can someone help me get the formula in R by hand to get the standard error and prediction confidence interval?

Thanks.

Umm… not a very popular question… plof !

You might tweet a brief statement of the question with a link to the URL for this topic

1 Like

I have some doubts that twitter is a good source of mathematical help for this topic, but I find it interesting as an experiment. I’m going to try to see what happens but I bet you that nobody responds.

I’m going to try to provide an example R code:

```
library(survival)
library(rms)
data(lung)
head(lung)
ddist <- datadist(lung)
options(datadist='ddist')
lung$SurvObj <- with(lung, Surv(time, status == 2))
fit <- psm(SurvObj ~ rcs(age,3)+sex+ph.karno,
dist="lognormal", x=T , y=T , data=lung)
# then I would like to obtain the 365-day probability of survival for patients aged 70, males, with Karnosfy 60%, so....
newdata <- expand.grid(age=70, sex=1, ph.karno=60) # I try to predict survival for a specific profile
sur<-survest(fit,newdata,times = 365)
sur # so the SE is 0.161, then I would like to replicate that by hand...
xb<-sur$linear.predictors
survival365 <- 1-pnorm((log(365)-xb)/fit$scale) # I obtain the same survival rate 'by hand'
pdf <- (dnorm((log(365)-xb)/fit$scale))/(365* fit$scale)
x <- c(70,70,70,1,60) # 5 elements to codify restricted cubic splines for age
se <- sqrt (pdf * t(x) %*% vcov(fit) %*% x * pdf )
se # 0.059 not similar to the result with survest ??
```

Anyone know how to obtain 0.16 by hand?