Survival Analysis – Test for PH Assumption vs Visual Inspection of Schoenfeld Residuals

cox-modelproportional-hazardsschoenfeld-residualssurvival

I am inspecting the PH assumption of a Cox model, both by testing it (with cox.zph in R) and visualizing the Schoenfeld residuals. I found many references with cases where the test is significant, but looking at the residuals, it looks there is no violation. This is the case with large sample size.

I have the other way around. I ran a Cox model adjusting for multiple variables with quite a large sample size, 6000 observations, 2000 events. One of them shows no violation of PH according to the test, but looking at the plot, I am wondering if there could be something. On the the left with the data, on the right without. To me, it looks like the variability of beta(t) is quite important at the end, after 25-28 days. So in favor of a time-varying coefficient. because a log-hazard ratio of zero around 11-15 days and 1 at the end is quite substantial.

Does it make sense to have such a behaviour? And no accounting for the result of the test? Scientifcally speaking, it would also make sense to observe this change.

Best Answer

This might be an artifact of the way that plot.cox.zph() generates its smoothed curves. The p-values for cox.zph() do not depend on the displayed smoothed curves, so an apparent discrepancy is possible.

By default, the smoothed curve is a natural spline with 4 degrees of freedom (df = 4). That's a reasonable choice for typical survival studies with several dozen to a few hundred events. It might be too small for a large data set like yours.

Natural splines enforce linearity outside the outer knots. I don't know how plot.cox.zph() chooses the knot locations, but the behavior of the smoothed curve you show suggests that the linear part of the curve at the far right is a linear extrapolation beyond an outer knot that might not well represent the data you have.

Try specifying more knots first (by changing df in the call to plot.cox.zph()), or playing with the knot locations (if that's possible). If you still see similar behavior, then what might be going on is that the number of events at those late time points is too small to overwhelm the wavy behavior around the coefficient estimate from many earlier time points, in the score test performed and reported by cox.zph().

You could then proceed with some handling of time-dependent behavior at late times if you think it's important, but that might then affect the validity of the PH assumption for other predictors.

Related Solutions

Solved – Test Cox proportional hazard assumption (Bad Schoenfeld residuals)

It is likely that the large sample size is responsible for the seemingly strong evidence against the PH assumption. P-values are a function of sample size, and their usefulness declines when sample size grows very large as the null hypothesis is never exactly true. They don't help too much with your question here, which is not "is the PH assumption satisfied" but "is the deviation from the PH assumption so large that inference is impaired".

One way to assess this for categorical variables, which your model seems to mostly contain, is by a log-minus-log plot. This is explained e.g. in this book and easily implemented in R using the rms library

library(rms)

myfit = survfit(Surv(time, status) ~ catvar)
survplot(myfit, loglog=T)

When the PH assumption holds, the lines are parallel, and their vertical distance is the log hazard ratio.

Convergent curves are seen when the difference between the groups decreases with time, and divergent curves when it increases, which indicates some deviation from the PH assumption. Crossing of the curves indicates a more severe deviation, with the effect of group membership changing signs.

Solved – Schoenfeld residual test for model with time varying coefficients

I believe I have answered my own question, and it just required thinking through what Schoenfeld residuals represent. Schoenfeld residuals can be thought of as observed minus expected values of the covariates at each failure time. In a standard Cox model, these residuals can be inspected for temporal trends to determine if any of the covariates have a time varying effect.

However, when you incorporate a time-varying coefficient, the time-varying coefficient is a function of time. For example, in the model I've been working with, the time-varying coefficient is described with two parameters: a slope and intercept defining a linear relationship between the coefficient and time. It seems that Schoenfeld residuals can technically be computed with time-varying coefficients (and cox.zph allows you to do it), however I don't see how they can be interpreted sensibly. You will get one set of Schoenfeld residuals for the intercept and another set for the slope, but it makes no sense to ask how the observed/expected values of the slope and intercept change through time, because these parameters themselves describe how a coefficient changes through time: the intercept corresponds to the expected value of the time varying coefficient when t = 0, and together with the slope, the two parameters define the expected value of the coefficient at every time t.

If I'm off base here or if there is a better way to think through my question, I would still really appreciate answers from people more experienced with Cox models. For a few weeks now, I've been struggling to develop intuition for the various types of residuals defined for Cox models and best practices for assessing model fit with time varying coefficients (particularly as implemented with the time-transform feature in coxph).

Best Answer

Related Solutions

Solved – Test Cox proportional hazard assumption (Bad Schoenfeld residuals)

Solved – Schoenfeld residual test for model with time varying coefficients

Related Question