# How to construct probability plots for survival distributions in R? An example plot from SAS

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

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"))


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)


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")

• It is okay to have negative values in the Kaplan Meir curve? As seen here survfit(Surv(residsWei)~1) Nov 22 '21 at 19:15
• @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.