Hi all, I noticed that the scaled Schoenfeld residuals produced by residuals.coxph(fit, type="scaledsch") were different from those returned by cox.zph for a model where robust standard errors have been estimated. Looking at the source code for both functions suggests this is because residuals.coxph uses the naive variance to scale the Schoenfeld residuals whereas cox.zph uses the robust version when it is available. Lines 20-21 of the version of residuals.coxph currently on github: vv <- drop(object$naive.var) if (is.null(vv)) vv <- drop(object$var) i.e. the naive variance is used even when a robust version is available. Why is this the case? Have I missed something? Am I right in thinking that using the robust variance is the better choice if the intention is to check the proportional hazards assumption? Here is a reproducible example using the heart data: data(heart) fit <- coxph(Surv(start, stop, event) ~ year + age + surgery + cluster(id), data=jasa1) # Should return True since both produce the scaled Schoenfeld residuals all(residuals(fit, type='scaledsch') == cox.zph(fit)$y) Thanks for your help. [[alternative HTML version deleted]]