1
$\begingroup$

Hello I am performing survival analysis I really like SAS's probability plots from lifereg as a way to visualize a distributions' fit and was curious if there is a way to do it in R? Here is the documentation I found about it but I got lost in the symbiology of it all https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.4/statug/statug_lifereg_details19.htm

enter image description here

Could someone explain in laymen's terms what's being plotted in the plot so I can make the plot by hand in R or direct me to a package or function that can help me here? Ideally I would like to know how to do it by hand so my knowledge can grow. I imagine it is similar to Q-Q plot but I may be wrong.

I am working in survival in R and here is some sample code that will be my starting point For making these probability plot. I imagine the estimated parameters from survreg are used in the Probability plot.

library(survival)
set.seed(123)
size <-  1000
deathtime <- rweibull(size, shape = 3, scale = 2)
(wei.fit <- survreg(Surv(deathtime) ~ 1, dist = "weibull"))

enter image description here

$\endgroup$
1
$\begingroup$

This is a plot of a type of residual against $\log T$, for accelerated failure time (AFT) models of the form:

$$\log T = \alpha + \sigma W ,$$

where $\alpha$ is a location parameter, $\sigma$ is a scale parameter, and $W$ has a specified probability distribution. For a Weibull survival model, $W$ has a minimum extreme-value distribution. In that case, residuals calculated as

$$ \frac{\log T - \hat \alpha}{\hat \sigma},$$

where $\hat \alpha$ and $\hat \sigma$ are the modeled estimates of those parameter values, should fit a minimum extreme-value distribution. A plot of those residuals against $\log T$ should have a slope of $1/\hat \sigma$ and an intercept of $-\hat \alpha/\hat \sigma$. Starting with your definition of deathtime, calculate and plot the residuals

residsWei <- (log(deathtime)-0.6963)/.3308
plot(log(deathtime),residsWei,pch=19,bty="n")
abline(-0.6963/.3308,1/.3308,col="red",lwd=3)

residuals against log T

This is the same type of plot as you show, except that the x-axis is labeled with values of $\log T$ instead of $T$, the y-axis is labeled with actual residual values instead of their percentiles, and it doesn't show the confidence interval for the linear fit.

Your example doesn't include censored survival times. With censoring you need to do extra work in defining the x-axis values, with the SAS default the "modified Kaplan-Meier method" described in the documentation to which you link. Censoring times are noted at the bottom of the plot you display. The documentation isn't very clear about how to define the x-axis values when covariates are involved and there is censoring, where $\hat \alpha$ in the above is replaced with $X' \hat \beta$, for covariate values $X$ and coefficient estimates $\hat \beta$.

I'm not a great fan of this type of plot for survival data, as it visually over-emphasizes the lowest survival times, it squishes together most of the observations, and it doesn't handle censoring very intuitively. I'd recommend a different way to evaluate parametric survival model fits that nicely incorporates censoring.

The trick is to examine the survival curve of the residual distribution, which is simply 1 minus its CDF, and compare that against the survival curve expected for the distribution of $W$. For your data and the minimum extreme value distribution of $W$ that corresponds to a Weibull survival model:

plot(survfit(Surv(residsWei)~1),xlab="Residual",ylab="1 - extreme value CDF",bty="n")
curve(exp(-exp(x)),from=-6,to=2,add=TRUE,col="red",lwd=2)

survival curve of residuals

This directly handles residuals from models with covariates and spreads out the observations more evenly along the y-axis than the type of residual plot you show from SAS. Right-censored cases contribute to the survival curve of residuals similarly to how they contribute to the survival curve of event times, and can be plotted directly on that curve, so they fit very nicely into this type of display. Sections 18.3.6 (see figure 18.11) and 19.2 of Harrell's Regression Modeling Strategies, second edition show how to use this approach as part of evaluating AFT models.

$\endgroup$
4
  • $\begingroup$ It is okay to have negative values in the Kaplan Meir curve? As seen here survfit(Surv(residsWei)~1) $\endgroup$
    – Vefeagins
    Nov 22 at 19:15
  • $\begingroup$ Also that book Regression Modeling Strategies is very nice thank you for sharing $\endgroup$
    – Vefeagins
    Nov 22 at 19:43
  • 1
    $\begingroup$ @Vefeagins if you think of a survival curve $S(x)$ as being the complement of a cumulative probability distribution function $F(x)$, that is $S(x) = 1 - F(x)$, then the values of $x$ can take on any values in the domain of that probability distribution. In survival analysis such distributions often are limited to non-negative values of $x$, but there's no reason in general why $x$ can't cover the entire real line. In this case, with $x$ representing values of residuals, the x-axis necessarily includes negative values. $\endgroup$
    – EdM
    Nov 23 at 15:10
  • $\begingroup$ That makes sense thank you! $\endgroup$
    – Vefeagins
    2 days ago

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.