Accelerated Failure Time Regression Performance (Survival Analysis)

Question

Issue is around high intercept with AFT Regression. Let me explain below:

Suppose you are modelling the time to an event via an Accelerated Failure Time Regression i.e. given survival time $T$, suppose we have observed values of covariates $x_{i1}, ..., x_{ip}$ and possibly censored survival time $t_i$, then:
$$ log(t_i) = \beta_0 + \beta_1 x_{i1} + ... + \beta_p x_{ip}+ \sigma \epsilon_i $$

Suppose we are looking at a Weibull AFT i.e. where $\epsilon_i $ are IID according to a Gumbel Distribution (Extreme Value Type 1).

You are looking at the case of time varying covariates (assume just one for now) e.g. you have a dataset like the following example with a single time dependent covariate (TDC_1). Where Start is the enter time (period start) and End is the period end (exit time) and UNIT_ID is the ID for the entity in the study:

START END EVENT UNIT_ID TDC_1 
0     1   0     1       0.1 
1     2   0     1       0.2
2     3   0     1       0.3
...
19    20  1     1       1.9
0     1   0     2       0.1
1     2   0     2       0.2
2     3   0     2       0.3
...
19    20  1     2       1.9

With the aftreg function from the eha library in R you can construct a Weibull AFT e.g.

model <- aftreg(Surv(START, END, EVENT) ~ TDC_1, dist="weibull", data=df, id=UNIT_ID, param='lifeExp')

Calling model.coefficients gives:

             model.coefficients
TDC_1        -0.905
log(scale)    9.393
log(shape)    0.046

The expected time to event when $T$ follows a Weibull distrubtion is given by: $$E(T|X_i) = exp \left( \beta_0 + x_i \beta_1 \right)\Gamma(1 + \sigma) = exp \left( 9.393 - 0.905*TDC_1 \right)*0.98 $$

As $\beta_0 = log(scale)$ and $\sigma = \frac{1}{exp(log(shape))}$

My question is around these parameter estimates (in particular the the intercept term $\beta_0 = log(scale)$. No matter how I change the error term parameterisation e.g. if $\epsilon_i$ are distributed normally (then $T$ lognormal) or if $\epsilon_i$ ~ Logistic etc, the intercept is exceptionally high and appears not to be optimal in terms of minimising error on time to event.

For example, if I manually subtract 2 from the intercept (9.393 - 2) I can reduce the root mean squared error on the time to event on the dataset fit:.

Intercept TIME_TO_EVENT_RMSE
9.393     776 days
7.393     97 days

Here TIME_TO_EVENT_RMSE is calculated as (with a dataset that only contains non-censored events):

$$ RMSE = \sqrt{\sum_{i}^{n} \frac{(exp \left( \beta_0 + x_i \beta_1 \right)\Gamma(1 + \sigma) - t_i)^2}{n}} $$

For further illustration, suppose you model directly using exponential regression (i.e. linear regression and logging the target variable) with exactly the same dataset (only using non-censored events so the two are comparable). I know they are minimising different loss functions and aren't directly comparable, but just for illustration purposes:

TIME_TO_EVENT UNIT_ID TDC_1 
19             1      0.1 
18             1      0.2
17             1      0.3
...

Here we have:

$$E(T|X_i) = exp \left( \beta_0 + x_i \beta_1 \right) = exp \left( 8.03 - 0.5*x_i \right) $$

I know that AFT Regression is not directly minimising RMSE, and that with the AFT regression the TDC_1 coefficient magnitude is larger in addition to a larger intercept, however with the intercept as high as it is, the model isn't particularly useful (significantly over-predicting the time to event).

Questions:

1. Has anyone experienced this before and have any advice on how to improve the AFT model?
2. Is there anyway to fix the scale with time varying covariates in AFTRegression?

Answer 1

This question shows why use of RMSE as proposed is not a reliable measure of survival model performance in real-world data.

This type of result can be seen in the mort data, analyzed by the aftreg() function in the R eha package. Load the package and the mort data frame.

> library(survival)
> library(eha)
> data(mort)
> head(mort)
  id  enter   exit event birthdate   ses
1  1  0.000 20.000     0  1800.010 upper
2  2  3.478 17.562     1  1800.015 lower
3  3  0.000 13.463     0  1800.031 upper
4  3 13.463 20.000     0  1800.031 lower
5  4  0.000 20.000     0  1800.064 lower
6  5  0.000  0.089     0  1800.084 lower

So ses is a time-varying categorical covariate, as seen for id=3. Now use aftreg to model survival as a function of ses with the default Weibull model:

> aft1.id <- aftreg(Surv(enter, exit, event) ~ ses, param = "lifeExp", data = mort,id=id)
> aft1.id
Call:
aftreg(formula = Surv(enter, exit, event) ~ ses, data = mort, 
    id = id, param = "lifeExp")

Covariate          W.mean      Coef Life-Expn  se(Coef)    Wald p
ses 
           lower    0.416     0         1           (reference)
           upper    0.584     0.365     1.440     0.089     0.000 

Baseline parameters:
log(scale)                    3.587               0.071     0.000 
log(shape)                    0.339               0.058     0.000 
Baseline life expectancy:  

Events                    276 
Total time at risk         17038 
Max. log. likelihood      -1390.6 
LR test statistic         17.3 
Degrees of freedom        1 
Overall p-value           3.14473e-05

For comparison, do a corresponding Cox model:

> cox1 <- coxph(Surv(enter, exit, event) ~ ses, data = mort)

Now, calculate (predicted - observed) event times and RMSE from the aftreg() fit as in the question, limiting analysis to rows with events:

> mortEvents <- subset(mort,subset=event==1)
> mortEvents[,"diff1"] <- exp(3.587+ 0.365*(mortEvents$ses=="upper")*gamma(1+(1/exp(0.339))))- mortEvents$exit
> mean(mortEvents$diff1^2)
[1] 1051.286

Following the approach in the question, reduce the magnitude of log(scale), say by a factor of 2, and repeat the RMSE calculation:

> mortEvents[,"diff2"] <- exp((3.587/2)+ 0.365*(mortEvents$ses=="upper")*gamma(1+(1/exp(0.339))))- mortEvents$exit
> mean(mortEvents$diff2^2)
[1] 51.90771

That sure did reduce the RMSE, calculated as proposed in the question. But now look at the associated survival curves, plotted for an ses value of "lower" out to the maximum exit time of 20. Start with the survival function from the aftreg() fit in black:

> weibull_S1 <- function(x) 1 - pweibull(x,exp(0.339),exp(3.587))
> curve(weibull_S1, from = 0, to = 20,xlab="Time (years)",ylab="Fraction surviving",bty="n")

Then superimpose the corresponding Cox fit with 95% CI for ses of "lower", in blue:

> lines(survfit(cox1,newdata=data.frame(enter=0,exit=20,event=0,ses="lower")),col="blue")

and finally add the corresponding Weibull survival curve with the reduced intercept value of scale in red, with annotation:

> weibull_S2 <- function(x) 1 - pweibull(x,exp(0.339),exp(3.587/2))
> curve(weibull_S2, from = 0, to = 20, add = TRUE, col="red")
> legend(4.5,.75,legend="blue, Cox\nblack, AFT Weibull, RMSE = 1051\nred, lower Weibull Scale, RMSE = 52",bty="n")

You get the following:

There's little question that the original aftreg() model is a superior fit to the data, despite the much higher RMSE. Here's what's happening.

The aftreg() mean survival for ses of "lower" is:

> exp(3.587*gamma(1+(1/exp(0.339))))
[1] 26.25835

which makes a lot of sense given the y-axis limits on the survival curve. With the reduced value of scale, you get a clearly erroneous mean survival:

> exp((3.587/2)*gamma(1+(1/exp(0.339))))
[1] 5.12429

The predicted mean survival is for a cohort followed until survival hits 0. As with most real-world data sets, the cases in this data set (particularly those with events) are not representative of such a complete cohort. With a maximum time of 20, the entire mort data set doesn't even extend out in time to the predicted mean survival of 26.3 for the most at-risk subset; all event times will be less than that mean value.

So if you choose a scale value that provides a mean survival time closer to the observed event times, you can certainly improve the RMSE for this particular data sample. But the RMSE doesn't represent how well the model fits over the entire range of survival times.