Robust standard errors with splines

Question

I realize that large changes in model results between using robust and non-robust standard errors can suggest a misspecified model.

My case refers to using a Cox regression and I have experimented with using both robust and non-robust SE. The difference in SE and p-values in my first run is small. However, to account for nonlinearity and proportional hazard violation, I use penalized splines (psplines()).

If I test this run both with and without standard errors, the results for my continuous variables with splines differ largely, and some pairwise differences become significant after the addition of splines that were not before in my categorical variables.

I am struggling to interpret what this means. On the one hand, some continuous variables have complex and nonlinear relationships with the outcome, so assuming linearity in all continuous variables (as in the first run) could make the model so bad that specifying robust=TRUE or not makes little difference to anything, whereas when I use splines to capture these relationships these differences become more apparent.

On the other hand, the changes only occur in continuous variables with splines; those with linear relationships not. I wonder whether these findings are just an artifact of using splines and not true results that suggest a misspecified model.

Does anyone have any experience with this? How might I test this?

Answer 1

You are dealing with regresssion effects, not with relaxing the proportional hazards assumption by using splines, unless you add time-dependent spline terms to the model. Huber-White robust sandwich covariance estimators are not guaranteed to be more accurate than ordinary Wald covariance estimators. In binary logistic regression William Gould has shown that the variance of robust standard errors can be many times larger than the variance of Wald standard errors. So it's not clearly the case that using a robust method improves things in your situation. It is true that you might get slightly better covariance estimates with sandwich estimates if proportional hazards is violated, but then you're doing what David Freedman wrote about: getting the right standard errors on the wrong quantities.

When looking at prop. hazards for splines don't look at individual coefficients. Use capabilities of the R survival package to get multiple degree of freedom measures of non-PH using cox.zph. And use a smoother to get a plot of the effect of the overall spline term collapsed to one degree of freedom, against time.