I'm working with cox proportional hazard (cox-PH) in R. I'm using the smoothed scaled Schoenfeld residual plot to assess the PH assumption.

Suppose I have this model in R, disease as outcome and genetics as exposure:

#Residual plot for genetics:
cox_model <- coxph(Surv(start, end, disease)~1+factor(genetics)+strata(gender), robust =T, data=mydata)

#Residual plot: 
test.ph <- cox.zph(cox_model)
plot(test.ph, resid = FALSE) 

The covariate genetics has levels in it. Meaning I have made a loop for the cox regression so it can run (make the analysis) for every single level in the covariate genetics. The levels are named like this geneA and geneB. So the covariate genetics has two levels.

When I run a PH test in R, and I look at the residual plot, then it says that the covariate on the plot is genetics. It does not show me if it is geneA or geneB. Then I tried to specify the model, so I can visualise the residual plot seperatly for geneA and geneB, like this:

#Residual plot for geneA:
cox_model <- coxph(Surv(start, end, disease)~1+factor(genetics==geneA)+strata(gender), robust =T, data=mydata)

#Residual plot
test.ph <- cox.zph(cox_model)
plot(test.ph, resid = FALSE) 

#Residual plot for geneB:
cox_model <- coxph(Surv(start, end, disease)~1+factor(genetics==geneB)+strata(gender), robust =T, data=mydata)

#Residual plot
test.ph <- cox.zph(cox_model)
plot(test.ph, resid = FALSE) 

I get two different plots (which is not a surprise). But what I'm confused about is, that the first plot I get with genetics as the covariate, looks different from the plots of geneA and geneB. So the conclusion is I get three different residual plots, one for geneA, second for geneB and a third for the covariate genetics with both levels?

My question is, how should I interpret the smooth scaled Schoenfeld residual plot that shows the covariate genetics? Is it a combination of both geneA and geneB? And what does that combination mean?

Many thanks in advance !!! :)

Update: I'm using a 3 level co-factor, meaning genetics has 3 levels saying geneA, geneB and none. EdM has answered my Q below

  • $\begingroup$ Please edit the question to add a sample model summary and the corresponding residual plots. I’m not sure that I understand the situation just from the written description. It might have to do with the choice of reference level for the categorical predictor. $\endgroup$
    – EdM
    Commented Jan 6 at 10:12
  • $\begingroup$ If a categorical predictor has more than 2 levels, the cox.zph() function (in recent versions) can provide either an overall plot for the predictor or plots for individual levels versus the reference. When you edit the question, please specify the version of the software that you are using. $\endgroup$
    – EdM
    Commented Jan 6 at 15:16
  • $\begingroup$ Thanks for answering @EdM. I'm working in a closed system that makes me unable to copy paste codes. I would have to manually write more than 30 lines of coding. I hope my following comment can clarify my Q. $\endgroup$
    – Devi Sita
    Commented Jan 7 at 12:01
  • $\begingroup$ It's the newest version of cox.zph() I'm using. It is exactly what it's providing me - an overall plot for the predictor, genetics (which has two levels, geneA and geneB). So it's showing me one plot with one curve for both geneA and geneB. My Q is, how should I interpret the overall plot? Because when I look at the individual levels plot (geneA and geneB separately), they look different from each other. So what is the overall plot showing? The average HR over time for both geneA and geneB? Or how should I think about it? $\endgroup$
    – Devi Sita
    Commented Jan 7 at 12:03
  • $\begingroup$ I have updated the question with adding a short example of what I have written in R. $\endgroup$
    – Devi Sita
    Commented Jan 7 at 12:10

1 Answer 1


Your attempt to isolate individual levels of genetics isn't doing what you think. The question that you raise about interpretation of Schoenfeld plots is the same as you might have raised about interpretation of model coefficients.

For a 2-level factor you are evaluating the difference in outcome between the 2 levels. You are always thus including both levels in the calculations and plots. The only issue is the direction of the difference, set by the choice of reference level. All 3 models you show in your question are the same except for their choice of reference level of the genetics factor.

If you don't specify otherwise, the reference level for a factor is the lowest in alphanumeric order of the level names, presumably geneA in your data. The regression coefficient is for the log-hazard difference between the second level and the reference level, or for geneB-geneA in this case.

When you specify factor(genetics==geneA) instead, you are simply reversing the choice of reference level. Both geneA and geneB cases are still included, but the reference level is now for a value of 0 (FALSE) for genetics==geneA; that is, the reference level is now geneB and the regression coefficient is for the geneA-geneB log-hazard difference. Specifying factor(genetics==geneB) just re-reverses the reference-level choice and gives the same model as you get with factor(genetics).

Here's are reproducible examples based on the lung dataset built into the R survival package, similar to your models. Sex is coded Male=1, Female=2, providing a 2-level categorical predictor like your genetics. I included 2 strata based on the inst numbers distinguishing treating institutions. There's no need to include the 1 in the formula, but I kept it to match your code.

mod1 <- coxph(Surv(time, status) ~ 1 + factor(sex) + strata(inst<10), data=lung)
mod2 <- coxph(Surv(time, status) ~ 1 + factor(sex==1)+ strata(inst<10), data=lung)
mod3 <- coxph(Surv(time, status) ~ 1 + factor(sex==2)+ strata(inst<10), data=lung)

The coefficients for mod1 and mod3 are the same; their scaled Schoenfeld residual plots are the same except for the y-axis label. The coefficient for mod2 is just the negative of that for the other models; its scaled Schoenfeld residual plot is a mirror image of the others. Both levels of genetics are involved in all models and plots.

You should know that things are a bit more complicated with a factor having 3 or more levels. Then the terms argument to cox.zph() determines whether you evaluate all levels together or evaluate each other level relative to the reference level.

In response to comments

Comments on this answer indicate that the genetics predictor actually was a 3-level categorical predictor, with levels none (the reference), geneA (for one type of variant), and geneB (for a different variant).

In that case the first residual plot produced by the code in the question, with the default terms=TRUE argument to cox.zph(), combines the scaled Schoenfeld residuals over all levels of genetics. That's analogous to how an anova() function combines a test over all levels of a categorical predictor, and is probably the most useful analysis and plot for evaluating PH for genetics.

With the 3-level genetics predictor, the second and third models in the question still didn't do what might have been intended. The factor(genetics==geneA) model compared geneA cases against the combination of none and geneB cases, and similarly for the factor(genetics==geneB) model. I suspect that neither the models nor their PH assessments would be very useful in practice.

If for some reason you want separate analyses and plots for geneA versus none and for geneB versus none, specify terms=FALSE in cox.zph().

  • $\begingroup$ As always, thanks a lot EdM! You are so good at explaining stuff about cox regression. I finally understand it now. I was confused before because I thought that the reference level was those who did not carry the cofactor genetics (meaning neither geneA or geneB). That's why I was confused about what the plot was showing. what if I wanted to see if the PH assumption is hold true between those carring the genetics and those who does not (the outcome would still be disease). Would that make sense to compare? $\endgroup$
    – Devi Sita
    Commented Jan 7 at 18:46
  • $\begingroup$ @SAphi11 perhaps you are already using a 3-level genetics predictor? That would make sense if the only (and mutually exclusive) options were none, geneA and geneB. Look at summary(factor(my.data$genetics)). If you aren't, I would recommend that. In that case you could set the reference level (via relevel()) to be those without either of geneA or geneB. You would choose, with the terms argument to cox.zph(), whether to get an overall estimate of PH combining all levels (my favorite) or to evaluate each of geneA and geneB separately with respect to the reference level. $\endgroup$
    – EdM
    Commented Jan 7 at 19:01
  • $\begingroup$ I missed that detail. I am using a 3-level genetics as you say. It's a none, geneA and geneB. So wouldn't that change the interpretation of the first plot I got from R, by visualising the PH test by only specifying the cofactor genetics (the plot that provided me results with genetics on the y-axis). Because then I don't understand which of the levels it is using as a reference and what it is compared to. When you say "get an overall estimate of PH combining all levels". Does that mean that level 0, A and B are all combined? Or the combination of geneA or geneB and comparing it to not carrying $\endgroup$
    – Devi Sita
    Commented Jan 7 at 19:32
  • $\begingroup$ I just found out that, the reference level by default is none (as it is assigned a zero as factor). So I'm suspecting that the first residual plot I got, used none as ref, and combined the HR of geneA and B. Does that make sense? and was it that one you pointed out as you favorite? $\endgroup$
    – Devi Sita
    Commented Jan 7 at 20:09
  • $\begingroup$ @SAphi11 I've added a bit to the end of the answer to cover the 3-level situation. $\endgroup$
    – EdM
    Commented Jan 7 at 20:48

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.