Diagnostics of a parametric Survival Regression in R

Question

I am doing a Survival Analysis in R with the "survival" package and I don't know how to do any plots of the results for diagnostic purposes.
Here is my model (I have given the variables self-explanatory names, I believe):
my_fit <- survreg(Surv(time, status) ~ feature_1 + feature_2 + feature_3, data = my_data_frame, dist = "loglogistic")
I use summary(my_fit) to see the coefficients and the corresponding p-values together with a couple of other statistical outputs like the Chi^2 of the model. However, with a linear regression model, for example, one can use the plot function to check out the residuals. That can give one a clue about possible nonlinearities. In this case plot(my_fit) gives an error:
"Error in xy.coords(x, y, xlabel, ylabel, log):
'x' is a list, but does not have components 'x' and 'y'"
I know that the survreg function uses the maximum likelihood method, so one doesn't have residuals, but perhaps there is still some useful function for a visual inspection of the fitted model.
Another question is how to compare two different models (i.e. models that differ only in their predictors). I used the anova function but am not sure how to interpret the output and I believe this makes sense only if the models are self-contained. Hence, if one of the models has predictors feature_1, feature_2 and feature_3, and the other one uses predictors feature_1, feature_2 and feature_4, then the anova function does not give a meaningful result. Is that correct?

Answer 1

Model Diagnostics

With a parametric survival regression model you do have residuals. You just might have to work a bit harder to get them than you do with ordinary least squares.

Chapter 18 of Frank Harrell's course notes and of his book on regression modeling strategies goes into detail about diagnostics for parametric survival models. A parametric model expresses the survival time $T$ as a function of the covariates $X$ (with associated coefficients $\beta$), a scale parameter $\sigma$, and a specific form of an error distribution, $\epsilon$, related to your choice of parametric family:

$$ \log (T) = X \beta + \sigma \epsilon$$

Thus for each event time you can calculate a residual from the parameter values estimated from your model fit:

$$\frac {\log(T)-X \beta} {\sigma} $$

and check whether the residuals have the assumed form. You also can check whether the covariates have the assumed linear relationship with $\log(T)$. Chapter 19 of Harrell's notes and book works through a detailed example with a log-normal model, with principles readily applied to your log-logistic model.

Model comparison

You are correct that anova() comparisons can't be used for non-nested models. I'm reluctant to recommend comparison of non-nested survival models, given the risk of omitted-variable bias when you omit any predictor associated with survival from a model.* It's generally best to include as many predictors associated with survival as you can, without overfitting. You might consider Akaike's An Information Criterion (AIC) for such non-nested comparisons, but please read this discussion before you do. You also could examine which model has the best predictive performance, for example repeating both models on multiple bootstrapped samples of your data and evaluating the competing models against the full data set.

---

*Unlike ordinary linear regression, but similar to logistic regression,you can have such bias even if the omitted variable isn't correlated with the included variables.