smooth scaled Schoenfeld residual plot - covariates with levels

Question

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

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.

library(survival)
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().