Understanding the output from coxph when used with pspline

Question

When I run:

res<-coxph(Surv(time_to,event)~pspline(covariate,df=3)
           ,data=a)

I get:

                              coef se(coef)      se2    Chisq   DF       p
pspline(covariate, df = 3), lin 5.42e-02 3.47e-03 3.47e-03 2.44e+02 1.00 < 2e-16
pspline(covariate, df = 3), non                            1.57e+01 2.04 0.00042

What exactly does the p-value on the second line (i.e. 0.00042) refer to (i.e. the name of the statistical test performed) and what are the implications of it being "significant"?

Answer 1

The vignettes of the R survival package are your friends here. With pspline() you are asking for a particular type of smoothing, with a potentially large number of regression coefficients whose magnitudes are reduced from what they would be otherwise (penalized) to provide the smoothing. Alternatives are natural splines or restricted cubic splines, which use fewer but unpenalized coefficients.

As the main survival vignette demonstrates near the end of Section 3.1, the 2 lines in the report about pspline provide a "simple check for linearity" with respect to the continuous predictor you are modeling with the spline. The first line is the linear ("lin")* part of the association in the spline model. That's clearly significant in your case.

The second line ("non") reports the significance of the nonlinear aspects of the fit. This provides a compact report on the multiple coefficients associated with the penalized spline fit. Technically, this is a Wald-type test, using the variance-covariance matrix of the coefficient estimates to determine if any of the nonlinear coefficient estimates differs significantly from 0. In your case, the low p-value indicates significant nonlinearity in the relationship between your covariate and log-hazard.

As a comment on the question indicates, the splines vignette in the survival package provides additional guidance for using this type of fit.

The details of how this is done are a bit tricky as, for example, you have non-integer degrees of freedom for the test and there are different estimates available for the variance-covariance matrix, leading to the two sets of coefficient standard errors reported, se(coef) and se2 (apparently identical in your case, but different in general). If you want to get into those details, see Therneau and Grambsch, which explains the basis for the survival package (then written for the S-Plus language, an older relative of R). Section 5.5 introduces penalized splines, and Section 5.8 shows how penalized models are fit in general, with details on penalized splines in Section 5.8.3.

---

• If the name of the predictor had been shorter, you would see more of the full words "linear" and "nonlinear" displayed in the printout.