# Predicting absolute risk using cox regression

I am trying to use R to predict the absolute risk of developing adverse events in a cohort, and to compare that with the observed outcome. Should I use survreg or coxph to do this? Anyone kind enough to explain how to do this with R code?

The mean follow up period of my cohort is only up to 6 years, so am I able to predict the absolute risk for each individual at the end of the follow up period (including both censored and non censored data)?

• By absolute risk do you mean the probability at some specific point in time? Nov 7, 2013 at 17:43
• @ Michael Mayer Thanks for your comment. Yes, it is the probability of having an adverse event at the time point where the follow up ends (for both censored and non-censored data). I wonder if we could calculate the probability for every single individual in the survival analysis? – Nov 9, 2013 at 9:47
• Would such lower bound of the absolute risk be useful? The longer the follow up, the higher the probability. Nov 9, 2013 at 12:21
• @ michael mayer I don't quite understand the meaning of lower bound of absolute risk. But I need to fit models to predict the absolute risk of developing a clinical event in order to compare the strength of different models. Nov 10, 2013 at 5:10
• The absolute risk will increase with increasing follow-up. So you end ip with a lower bound. Nov 10, 2013 at 7:07

So by "absolute risk" I'm going to assume you mean either the cumulative probability of an event at time t, or the hazard at time t.

In short, no, a Cox proportional hazards model doesn't really give you back that information - the model itself doesn't calculate the underlying hazard, just the relative difference in the hazard between covariate values. This rather nicely frees you from having to specify the underlying hazard of your outcome, which in many cases is unknown, not particularly of interest in the first place, or difficult to specify using a parametric model.

If you do want to estimate the underlying hazard function, you need to use parametric survival models, such as those used in survreg. There are a large number of tutorials online, including code, for survreg and parametric models.

This is one of my favorites, as it includes some theoretical treatment and a good bit of code. The 'Survival' package documentation is also a good place to start.

• I wonder if we could get the predicted probability for the censored cases as well, i.e. those who did not develop any event at the end of follow up? Nov 8, 2013 at 12:03
• @user32454 Parametric survival models can handle censoring - they can actually do it slightly better than Cox models, IIRC. Nov 8, 2013 at 18:26
• @ EpiGrad Yes, as I am trying to predict the probability of developing an adverse event for each individual in the survival analysis, i.e. the cumulative probability of an event at the time point where the follow up period ends. Nov 9, 2013 at 9:51
• @Fpmite would it not be possible doing predictions using the Nelson-Aalen estimator of the baseline hazard-function (described here stats.stackexchange.com/questions/36015/… by @ocram) ? Mar 21, 2015 at 22:34
• @ErosRam Honestly, I'm not familiar enough at the moment to definitively assert that's true. Mar 22, 2015 at 9:20

I agree with the above answers that Cox was not primarily meant for absolute risk calculations, but I believe it was more for historical reasons (people focused on relative risks and not absolute), and don't see a legitimate reason why not to use it for predictions. There are works in prediction modelling that successfully use Cox already, but I think they prefer glmnet package for Cox family, with flexibility of regularization terms and hdnom package for survival probabilities (hdnom:::glmnet_survcurve).

Back to the question, going in cycles for survfit basic Cox model, I came up with this code to back out estimated individual probabilities of an event:

bh=basehaz(cox_model_fit, newdata = data, centered = FALSE) #baseline cumulative hazard function for all covariates set to 0 (centered = false
lps <- predict(cox_model_fit, data, type = "lp") #linear predictors, betas x covariates for each observation

time_of_interest = 6 #if you are interested in survival at t=6

i_time_of_interest = match(1, round(bh$time,1) == time_of_interest, nomatch = -100) #annoyingly, no time argument in basehaz function, so find that time in the list (it builds at all censored and event times, so either choose a point like that or need to interpolate between available points) event_prob_t = 1-exp(-bh[i_time_of_interest,1])^exp(lps) #calculate risk of event from the baseline hazard and linear predictors. returns a vector of probabilities for each person in the "data"  You can then check event_prob_t separately for cases and controls / censored and test the difference etc. Or, use this function to calculate concordance statistic for survival data- this gives probability that two random observations from cases and controls have probabilities in the right order (lower chances of event for controls than for the cases). For binary outcome c-index is the same as ROC_AUC, so very common measure of how good a model is for classification/binary models.  survConcordance(Surv(data$survtimed, data$outcome) ~ lps)$concordance