Extract the position of knots from penalized spline in cox regression in R

31 Views Asked by At

I am calculating a cox model with two predictor variables (var1 and var2). var1 is a numeric variable and enters the cox regression without any modification. var2 on the other hand enters the cox regression after transformation into a penalized spline. See the following code:

library(splines)
fit <- coxph(Surv(time2event, Mortality)~ var1 + pspline(var2,2), data = data)
summary(fit)

The output of the summary function is the following:

                          coef    se(coef) se2      Chisq DF   p      
var1                      0.05215 0.005283 0.005266 97.42 1.00 5.6e-23
pspline(var2, 2), linea   0.43643 0.065597 0.065597 44.26 1.00 2.9e-11
pspline(var2, 2), nonli                              2.90 1.07 9.7e-02

            exp(coef) exp(-coef) lower .95 upper .95
var1            1.054   0.949191    1.0427     1.064
ps(var2)3       3.179   0.314526    0.7537    13.413
ps(var2)4      10.098   0.099031    1.0029   101.670
ps(var2)5      29.274   0.034160    2.2205   385.929
ps(var2)6      61.328   0.016306    4.8178   780.665
ps(var2)7     101.100   0.009891    8.0312  1272.687
ps(var2)8     136.941   0.007302    9.7441  1924.514
ps(var2)9     182.140   0.005490    8.0620  4115.002

Can someone explain how to interpret the output? And is it possible to extract the position of knots chosen in the penalized spline?

Thanks!

1

There are 1 best solutions below

0
PBulls On

The interpretation of such effects isn't really on topic here (a link to Cross Validated was raised in the comments), though I can say that the individual basis functions aren't too useful. You'd rather check the joint test for the non-linear parameters (lines 2-3 in the printout) and then have methods like predict reconstruct the actual spline for you.

Obtaining the knots directly doesn't seem possible, although some of the code in the non-exported but still documented survival:::psplineinverse seems to suggest that at one point this was an attribute of the spline object. (Through the magic of version control we can see that this attribute was added in commit f5a8c35b and removed again in commit 0ce58fc6, both dating back well over a decade)

The next best thing is to look at the relevant parts of the source code:

pspline_knots <- function(x, df=4, nterm=2*df, degree=3,
                          Boundary.knots=range(x), ...) {
  if (df == 0) nterm <- 15 else nterm <- round(nterm)
  dx <- (Boundary.knots[2] - Boundary.knots[1])/nterm
  knots <- c(Boundary.knots[1] + dx * ((-degree):(nterm - 1)),
             Boundary.knots[2] + dx * (0:degree))
  return(knots)
}

Note that this doesn't implement all of the input validation of the original function and only defines the relevant arguments.