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

Question

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.

Answer 1

The statistical issues here are whether there are any associations of the predictors with outcome and whether the penalized splines show any evidence of a non-linear associations. In this display, that fundamental information is in the top 4 lines of coefficient values.

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.