# How to interpret p values of a non-linear covariate using pspline in a coxph model

I tested the assumptions for Cox proportional hazards model on my time-to-event data. I found that the assumption of linearity between independent variables and model residuals is violated.
After some reading I realized that I could use pspline with 4 degrees of freedom to handle non-linearity. My Cox model has 2 explanatory continuous variables both of which are non-linear.
Here is the output of the coxph with pspline of both variables:

Call: recsurve <- Surv(timetoevent,convert)
coxph(formula = recsurve ~ pspline(Occipital_lGI, df = 4) +
pspline(Prefrontal_lGI, df = 4), data = C_NC_pref_occ_lGI)

n= 72, number of events= 24

coef  se(coef) se2   Chisq DF   p
pspline(Occipital_lGI, df 1.794 1.594    1.542 1.27  1.00 0.2600
pspline(Occipital_lGI, df                      1.48  3.00 0.6900
pspline(Prefrontal_lGI, d 5.724 2.153    2.096 7.07  1.00 0.0078
pspline(Prefrontal_lGI, d                      5.80  2.99 0.1200

exp(coef) exp(-coef) lower .95 upper .95
ps(Occipital_lGI)3     0.39212     2.5502 4.405e-03 3.491e+01
ps(Occipital_lGI)4     0.20268     4.9339 5.591e-04 7.347e+01
ps(Occipital_lGI)5     0.21215     4.7135 7.094e-04 6.345e+01
ps(Occipital_lGI)6     0.28838     3.4676 1.146e-03 7.256e+01
ps(Occipital_lGI)7     0.19500     5.1281 7.183e-04 5.294e+01
ps(Occipital_lGI)8     0.31445     3.1802 1.099e-03 8.997e+01
ps(Occipital_lGI)9     0.65732     1.5213 2.354e-03 1.836e+02
ps(Occipital_lGI)10    0.69503     1.4388 2.392e-03 2.019e+02
ps(Occipital_lGI)11    0.62975     1.5879 1.540e-03 2.576e+02
ps(Occipital_lGI)12    0.69336     1.4423 3.600e-04 1.335e+03
ps(Occipital_lGI)13    0.85750     1.1662 2.450e-05 3.001e+04
ps(Occipital_lGI)14    1.08123     0.9249 3.405e-07 3.433e+06
ps(Prefrontal_lGI)3    0.61636     1.6224 1.416e-03 2.683e+02
ps(Prefrontal_lGI)4    0.38646     2.5876 1.447e-05 1.032e+04
ps(Prefrontal_lGI)5    0.26852     3.7242 8.383e-07 8.601e+04
ps(Prefrontal_lGI)6    0.25250     3.9603 2.823e-07 2.258e+05
ps(Prefrontal_lGI)7    0.19687     5.0794 2.096e-07 1.849e+05
ps(Prefrontal_lGI)8    0.08341    11.9884 1.155e-07 6.023e+04
ps(Prefrontal_lGI)9    0.17187     5.8183 2.685e-07 1.100e+05
ps(Prefrontal_lGI)10   0.46969     2.1291 7.953e-07 2.774e+05
ps(Prefrontal_lGI)11   0.83136     1.2029 1.313e-06 5.264e+05
ps(Prefrontal_lGI)12   1.34341     0.7444 1.943e-06 9.288e+05
ps(Prefrontal_lGI)13   2.05412     0.4868 2.699e-06 1.563e+06
ps(Prefrontal_lGI)14   3.16502     0.3160 1.929e-06 5.193e+06

Iterations: 4 outer, 13 Newton-Raphson
Theta= 0.1809375
Theta= 0.1955302
Degrees of freedom for terms= 4 4
Concordance= 0.727  (se = 0.061 )
Likelihood ratio test= 18.14  on 7.99 df,   p=0.02


Can anyone explain which p-values I have to look at?
Whether any of the 2 variables is significant??
Why there are 2 p values for each variable marked as df and d?

Any guidance would be appreciated.

• Please note that you have too few events to fit this complicated a model. With only 24 events you can't really fit more than about 2 predictors without overfitting (rule of thumb, 10-20 events per predictor), yet each degree of freedom for each spline counts as a predictor (for 8 predictors total). Perhaps you can find a simple transformation (e.g., log) of one or both predictors that will improve linearity and use only 1 df per predictor. At best, you need to cut back to much less flexible splines. Also, consider using restricted cubic splines (rcs()) in the R rms package instead.
– EdM
Nov 23, 2020 at 20:25

The rest of the display is about the coefficients associated with each of the ps() polynomial terms from the basis functions used to do the penalized spline fitting: with 4 df specified, that's 10 third-degree polynomials (although constraints in putting together the spline lower the number of independent coefficients considerably) for each of your 2 continuous predictors. Those latter values aren't needed for the critical inference. They are shown because you used summary(recsurve) to get this output rather than print(recsurve) (or just typing recsurve at the prompt), which is what is shown for pspline() examples in the R survival package documentation.
Unfortunately, with long variable names the beginning of each of those critical top 4 lines about inference on coefficient values is truncated. pspline(Occipital_lGI, df should have been followed by  = 4, lin in the first line and by  = 4, non in the second. Truncation is even more severe for the third and fourth lines. For each of those pairs of lines, the first shows a test on the linear trend and the second a test on the non-linear terms. So you show no significant association of Occipital_lGI with outcome and a strong linear association of Prefrontal_lGI with outcome. There's also a possible non-linear association of Prefrontal_lGI with outcome (p = 0.12), which you will have to evaluate based on your understanding of the subject matter.
That said, this highly flexible model is almost certainly overfitting with only 24 events. You need to use up many fewer degrees of freedom, down to 2 or maybe 3 total including both predictors. A simple linear fit with your 2 continuous predictors is about as far as one might go without significant overfitting in a non-penalized model with so few events. If you feel compelled to do spline fits, cut back to the minimum number of 2 df per spline and perhaps restrict the spline fit to Prefrontal_lGI.