# smooth scaled Schoenfeld residual plot - how to interpret the plot which is based on a stratified cox model?

I understand what a Smoothed scaled Schoenfeld residual plot is, and what it shows. Before I go into details with my question, I will provide an example.

I use R to code.

Here is a cox model:

cox <- coxph(Surv(start, end, status) ~ 1 + genetics + strata(gender), robust=T, data = my_data)


Then I check if the proportional hazard assumption is hold

test_ph <- cox.zph(cox)
cox


From this I get results for the covariate genetics. The PH assumptions is not violated as you can see

          chisq df        p
genetics   6.75  1   0.0602
GLOBAL     6.75  1   0.0602


Then I make a smoothed scaled Schoenfeld residual plot (which is also a way of inspecting the PH assumption), like this

plot(test_ph, resid = FALSE, main = "residual plot")


I get one plot showing a line which is the HR over time for the covariate genetics. The line is approximately around 0 and almost straight. From this we can also say that the PH assumption is not violated.

Now to the question: When doing the PH test, shouldn't I get two results for the covariate genetics? One for females and one for males, as I choose gender as my strata? Also, shouldn't I also get two residual plots showing how HR varies over time for each gender? But instead I get one plot showing the HR over time. How do I interpret that in relation to stratification in the cox model?

Sidenote: If anyone sees this question and wants to understand more of what a residual plot is showing, then you can see my previous question about the same model. I asked for help to understand how to interpret the residual plot when a covariate has 3 levels. Here is the link smooth scaled Schoenfeld residual plot - covariates with levels

When doing the PH test, shouldn't I get two results for the [3-level] covariate genetics?

For a multi-level categorical predictor, the terms argument to cox.zph() determines whether the proportional hazards test is on its individual coefficients (terms=FALSE) or is combined over all coefficients (terms=TRUE, the default).

Here's a reproducible example inspired by a model in Section 6.3 of Therneau and Grambsch. The cell type factor has 4 levels, with squamous as the reference level.

library(survival)
fit.vet.CT <- coxph(Surv(time, status) ~ celltype, data=veteran)


With a 4-level categorical predictor there are 3 associated regression coefficients. Each event thus has 3 associated Schoenfeld residuals. With terms=FALSE the test is performed on each coefficient separately.

zphNoTerms <- cox.zph(fit.vet.CT,terms=FALSE)
zphNoTerms
#                   chisq df     p
# celltypesmallcell  1.44  1 0.229
# celltypelarge      4.70  1 0.030
# GLOBAL             8.87  3 0.031


That's what you seem to expect. The smoothed plot for each of the 3 coefficients displays the smoothed scaled Schoenfeld residuals added to the coefficient estimate, thus providing a visual representation of the change in the coefficient as a function of time. The plot.cox.zph() function labels the y-axis with the name of the individual coefficient. See this page for details.

With the default terms=TRUE instead (as you implicitly chose), those 3 residuals are combined into a single value.

zphTerms <- cox.zph(fit.vet.CT)
zphTerms
#cox.zph(fit.vet.CT)
#         chisq df     p
# celltype  8.87  3 0.031
# GLOBAL    8.87  3 0.031


The code for cox.zph() shows that this starts by evaluating the linear predictor associated with all levels of the factor at the values of the residuals. That is: weight each of the (unscaled) Schoenfeld residuals by the corresponding regression coefficient and add up all the weighted residuals. Then scale that value by a similarly combined variance estimate. Finally add a value of 1, the coefficient for a linear predictor.

You can think about this as an estimate of the overall residual associated with the categorical predictor, as each individual residual is weighted by its "importance" in the sense of its regression coefficient. The plot shows the overall deviation from proportional hazards (PH) over time; under PH, the plot should be flat in time with a y-axis value of 1. The plot.cox.zph() function labels the y-axis with the name of the categorical predictor itself.

shouldn't I also get two residual plots showing how HR varies over time for each gender?

No, unless you have covariate by strata interactions.

With an additive strata() term like you show, a Cox regression coefficient is assumed to be the same for all strata. You thus want an estimate of the PH violation associated with that single regression coefficient. The scaling of the Schoenfeld residuals is done within strata, however, and the object produced by cox.zph() has a "strata" value that you can use to identify which rows of the matrix of residuals come from each stratum and generate your own plots.

If you specify a model having an interaction of a multi-level categorical predictor with strata, however, then you can get what you expect if you also specify terms=FALSE to cox.zph(): one row in the output and one plot provided for each combination of coefficient and stratum. Play with the following to see how this works for additive versus interaction terms for strata:

fit.vet.CTstrat <- coxph(Surv(time, status) ~ celltype + strata(factor(trt)), data=veteran)
fit.vet.CTstratInt <- coxph(Surv(time, status) ~ celltype * strata(factor(trt)), data=veteran)

• A BIG thank you @EdM!! You help me understanding things much better! Commented Jan 12 at 13:24
• Hi again @EdM ! So looking back at this, and I was thinking, that in the calculation for each Schoenfeld residual for a covariate, does the calculation also (somehow), taking a second covariate into account? in other words, a model which estimates the risk of covariateA, and adjusted for covariateB, does the schoenfeld residuals for covariateA, also adjusts for covariateB? I mean, I do understand that the adjusted estimate for covariateA is part of the calculations for each schoenfeld residual, so I'm assuming the answer is yes, that the schoenfeld residuals for covariateA is adjusted? Commented Jun 11 at 12:39
• @DeviSita although the residual calculation is done for one covariate at a time, it uses the risk-weighted values of that covariate, where the risk-weighting depends on all the covariates. See this page. For details of the risk-weighting, which is part of getting the solution to the Cox model, see this page.
– EdM
Commented Jun 11 at 12:47
• In the last sentence I meant that, I'm assuming that the Schoenfeld residuals for covariateA are considered as adjusted residuals, as the estimate for covariateA is adjusted. I hope my question makes sense :) Commented Jun 11 at 12:54
• my last comment was sent before I saw your answer Commented Jun 11 at 12:54