1
$\begingroup$

What is the relationship between the df in survival::pspline and the number of knots? Say I want to fit a curve made up of cubic polynomials and has N internal knots. What would should I set for df? I don't really understand this from the R documentation.

Secondly, when I fit a model, say fit <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,3), lung), and look at fit$coefficients, there terms ps(age)3 until ps(age)12. What does this mean?

I have read this post but I am not sure I fully understand how it translates to my case.

$\endgroup$
2
  • 2
    $\begingroup$ Close-voters: I believe the first question is primarily about its statistical content, and this should therefore stay open. See here. $\endgroup$ Commented Feb 4, 2023 at 6:21
  • 1
    $\begingroup$ Although the second part is software-specific to R, it's more about how to interpret the output of statistical software, which should be on-topic too. $\endgroup$
    – Silverfish
    Commented Feb 4, 2023 at 12:24

1 Answer 1

1
$\begingroup$

You seem to be thinking about pspline() similarly to a cubic regression spline. That's not how it works. It's a simple form of a smoothing spline. This page compares different smoothing methods. This page compares pspline smoothing and the regression splines implemented by the rcs() function in the rms package, in the context of a Cox model.

If you want to specify knot positions between which unpenalized cubic functions are fit, then you are talking about regression splines. The approach of pspline() is different. It sets up a series of basis functions spaced evenly across the limits of the predictor values and produces a penalized regression coefficient for each of the basis functions.

If you specify a value for df and don't specify boundary knot locations, then the default nterm argument is round(2.5*df) (8 in your example) and knots to join the basis functions are placed according to the following code:

dx <- (Boundary.knots[2] - Boundary.knots[1])/nterm
knots <- c(Boundary.knots[1] + dx * ((-degree):(nterm - 1)), 
        Boundary.knots[2] + dx * (0:degree))

where the default Boundary.knots[2] and Boundary.knots[1] are the limits of the data range and degree is 3 for a cubic fit.

The multiple coefficients that you see reported are for each of the corresponding basis functions between the knots. Each of those coefficient values, however, is penalized or reduced in magnitude from what it would be in a fully fit model. The result is a less wiggly curve than an unpenalized fit would produce.

The df argument to pspline() effectively determines how wiggly the final fit is. That's useful if you have decided to spend a certain number of degrees of freedom on modeling a continuous predictor and prefer the smoothing spline to the regression spline approach. The penalization then results in a number of degrees of freedom that is substantially less than the number of knots. Alternatively, if you specify df = 0 the modeling will find optimal smoothing based on the Akaike Information Criterion.

This page goes into more detail about what a Cox model fit with a pspline() term reports.

$\endgroup$
4
  • $\begingroup$ Thanks @EdM -- love your answers! I don't understand exactly what pspline's "df" means then? If df=3 gives 8 coefficients, it's not really interpretable as degrees of freedom? If we use rms::rcs(x,1), it says the number of knots must be >=3. Does this mean this number includes the boundary knots? Also if the location is not specified, Is this fitted by AIC? $\endgroup$
    – Snoot
    Commented Feb 6, 2023 at 13:55
  • $\begingroup$ @Snoot the df value is the effective degrees of freedom devoted to that predictor in the model. With penalized coefficients it's the effective df that matters for statistical analysis and estimates of overfitting, not the total number of associated coefficients. Without penalization, those are the same. The more the penalization, the lower the effective df. Frank Harrell give a brief summary here, with a link to more information. $\endgroup$
    – EdM
    Commented Feb 6, 2023 at 14:24
  • $\begingroup$ @Snoot parameterization differs among implementations of splines. The rcs() function uses unpenalized regression splines, not the penalized splines of pspline(). For rcs() you specify the total number of knots including "boundary" knots, but its default is to place the "boundary" knots somewhere within the range of data values instead of at the extremes. That's unlike the ns() function in R, which by default places them at the extremes. With a regression spline, the fit is forced to be linear outside the boundary knots. $\endgroup$
    – EdM
    Commented Feb 6, 2023 at 14:31
  • $\begingroup$ Thanks a million @EdM, you are a wealth of knowledge! $\endgroup$
    – Snoot
    Commented Feb 6, 2023 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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