Solved – P-values for spline terms in coxph in R

chunk-testcox-modelsplinessurvival

I have a Cox PH model including a smoothing spline as follows, where x and z are covariates:

fit<-coxph(Surv(start,end,exit) ~ x + pspline(z))

While I understand that getting a single coefficient to summarise a spline term is not meaningful and, accordingly, I intend to plot these non-linear relationships (see this previous question), I do want to get p-values for the spline term.

I have seen displayed in e.g., Keele, L., 2010. Proportionally Difficult: Testing for Nonproportional Hazards in Cox Models. Political Analysis, 18 (2), 189–205. Having replicated the coxph models in that article, it is not obvious to me how the significance values for the spline terms are derived.

I'd appreciate any advice you can give.

Best Answer

I suggest you use the ANOVA for nested models.

fit<-coxph(Surv(start,end,exit) ~ x + z)
nonlinear_fit <- update(fit, .~.-z+pspline(z))
anova(fit, nonlinear_fit, tst="Chisq")

This is at least how I deal with it in my addNonlinearity function in the Greg package.