I have a dataset of patients who all underwent one of three types of treatment. The disease they have can be graded on a scale of 1-4. The question I would like to answer is:

Is there a difference between recurrence/survival outcomes between the three different treatments when adjusted for disease grade?

I've looked into ggadjustedcurves in the 'survminer' package and have been able to make adjusted survival curves, but I was wondering if there's a statistical test to run to see if the adjusted curves are different. Or is there a statistical test to run to see if 5-year adjusted survival outcomes are different or not.



1 Answer 1


I would recommend staying away from the survminer package, at least until you know a lot more about survival analysis. It has serious software glitches that weren't fixed when last I looked, and I don't think it provides a very good introduction to the basic principles of survival analysis.

Start instead with the basic survival package that comes with R. Although it doesn't fit so nicely into the "tidyverse" as survminer, it has a superb set of vignettes that explain the basis of how to accomplish many types of tasks in survival analysis. The main package vignette is a particularly useful overview.

When you do that, you will find that a Cox model with treatments and grade as predictors will probably accomplish what you want. For example, the model

mod1 <-  coxph(Surv(time, status) ~ treatment + factor(grade))

will allow you to "adjust for" grade in a way that's the same for all treatments. The baseline survival curve will be for the baseline levels of treatment and grade; the 2 regression coefficients for treatment will indicate the differences in log-hazard from baseline associated each of those treatment types after "adjusting for" grade. As a bonus, you get estimates of the associations of other levels of grade with outcome, "adjusted for" treatment.

You can evaluate whether there are any differences associated with treatment by an anova() comparison of the above model against

mod0 <- coxph(Surv(time, status) ~ factor(grade))

provided that you do build both models on the same individuals. Post-modeling software like that of the emmeans package can evaluate pairwise differences among treatments.

  • $\begingroup$ Appreciate your response. I like the idea of using an ANOVA test to compare the two models to look at the impact of treatment type on outcomes. Two questions. (1) What is the downside of using an ANOVA test to do this? (2) How would you go about showing this graphically? $\endgroup$ Feb 5, 2023 at 0:08
  • $\begingroup$ @user2930701 there is no downside, provided that you make sure both models are built on the exact same data set. That's then a simple, reliable likelihood-ratio test between nested models. (The potential problem is that if some individuals have missing data on treatment they might be included in mod0 but silently excluded from mod1.) Graphical display would be predicted survival curves (with confidence intervals) for all 3 treatment types at some particular value of grade. With the additive mod1 it won't matter which grade level you choose. $\endgroup$
    – EdM
    Feb 5, 2023 at 3:35

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.