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.


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.

  • $\begingroup$ Forgive my ignorance, but it as far as I understand your ANOVA tests the difference between a model with a linear z fit and a model with a spline z fit. Would this p value tell you if the model with the spline z fit was significantly better than the model with the linear z fit? Would it not make more sense to just test the model without the z term at all versus the model with spline z to get an overall p value for the spline term? $\endgroup$ Apr 1 at 19:25
  • $\begingroup$ By default most spline model have more than one additional term and the key thing is that the models are truly different. When you add variables usually the entire solution changes and you want to look at the entire solution and not just 1-3 columns $\endgroup$
    – Max Gordon
    Apr 3 at 8:26

Tests involving spline terms are of two types: chunk tests for nonlinearity and chunk tests for total association (chunk = multiple degree of freedom). You can get both of these by comparing full and reduced models to get the likelihood ratio test (best) or by using general contrasts on Wald tests. For the latter, the R rms package makes it easy:

f <- cph(Surv( ) ~ rcs(age, 5) + sex)
anova(f)  # gets 3 d.f. test of nonlinearity, 4 d.f. association test

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.