I am trying to get hazard ratio predictions from a coxph model at different levels of a given predictor. I know this is possible in SAS using the HAZARDRATIO function, but I do not know how to do this in R. I am using a large biomedical dataset, and in order to properly assess the relationship between my predictor of interest (inflammation) and outcome (mortality), I need to adjust for categorical values such as education level, ethnicity, sex, location (there are a limited number of assessment centers), and smoking status.

I have fit my hazard model as follows: coxph(Surv(Follow_up, Mortality) ~ Leukocyte_Count + Platelet_Count + Age + Sex + Smoking_Status + Education + Ethnicity + Location + Inflammation_Score, dataframe, na.action=na.exclude)

Essentially, I want to generate risk predictions for a few given "individuals" with different inflammation scores. To do this, I am using the predict.coxph function in the following way:

predict(inflammation_model, newdata = mean_DF, se.fit = TRUE, type = “risk”)

Where "mean_DF" is:

mean_DF <- with(dataframe,
                data.frame(Leukocyte_Count = rep(mean(Leukocyte_Count),5),
                           Platelet_Count = rep(mean(Platelet_Count),5),
                           Age = rep(mean(Age),5),
                           Sex = c("Male","Male","Male","Male","Male"),
                           Location = c("Boston","Boston","Boston","Boston","Boston"),
                           Smoking_Status = c("Never","Never","Never","Never","Never"),
                           Education = c("High School","High School","High School","High School","High School"),
                           Ethnicity = c("White","White","White","White","White"),
                           Inflammation_Score = c(-10,-5,0,5,10)

I am trying to ascertain the impact of only different inflammation scores on the outcome, so I have set all continuous covariates to the mean for each "individual" and all the factors to the same level. However, in the documentation for this function it states that the reference value for the predict function is the mean within each strata. Since I cannot (or do not know how to) calculate the mean for categorical variables, I am stuck using just the same "level" for each. I have three questions, although the third may be a different topic.

  1. How does predict.coxph calculate the "mean" for categorical variables to serve as the reference when predicting a hazard ratio?

  2. How can I set my categorical variables to this "mean" value? This is mostly for data visualization purposes - I do not want to make a misleading figure by selecting, e.g. "High School" as the level for education and having that be different from the mean value and therefore impact the predicted risk relative to the "mean"

  3. How can I plot a predicted survival curve based on the output of predict.coxph?


2 Answers 2

  1. The reference level is given by the $means component of the fitted model
  2. I don't think you can: a categorical variable only has the categories that it has.
  3. You can't. That's what survfit.coxph is for. The predicted curve is absolute, not relative, so you don't need to worry about means; just set the covariates to the values you want.

Essentially, what you are looking for is a plot of the marginal effect of the inflammation score on the survival probability, adjusted for further variables. Using the mean of variables does not really do the trick here, even if you only had continuous variables. I developed an R-package specifically for this purpose called contsurvplot.

It uses a specified model (e.g. your cox model) to perform g-computation to estimate the marginal survival probabilities for the population and plots it using different methods. For example, you could create the plot you are looking for using code like this:


# using data from the survival package
data(nafld, package="survival")

# fit cox-model with age
model <- coxph(Surv(futime, status) ~ age + male + bmi, data=nafld1, x=TRUE)

# plot effect of age on survival using defaults
                horizon=c(15, 30, 40))

lines This example shows the marginal survival probability over time for different values of the body-mass-index. Instead of picking some arbitrary values you could use a survival area plot:



The methodology and further plots are explained in detail in the associated publication: https://arxiv.org/pdf/2208.04644.pdf (soon to be published in "Epidemiology").


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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