Limit on number of features in Cox Proportional Hazards model with penalized spline

Question

I am trying to fit a Cox proportional hazards model with penalized splines using the survival R package with the following, but it seems like I am hitting a limit on the number of coefficients that can be estimated given the total number of observations in the dataset.

For example,

library(survival)

data(lung)
lung <- na.omit(lung)
lung <- data.frame(apply(lung, 2, as.numeric))
nrow(lung) # 167 rows
model <- coxph(Surv(time, status) ~ pspline(inst) + pspline(age) 
               + pspline(ph.ecog) + pspline(ph.karno) + pspline(pat.karno)
               + pspline(meal.cal) + pspline(wt.loss) + sex, data=lung,
               method='breslow')

This yields Error in coxpenal.fit(X, Y, istrat, offset, init = init, control, weights = weights, : NA/NaN/Inf in foreign function call (arg 5).

But when I double the number of observations:

lung.2x <- as.data.frame(lapply(lung, rep, 2))
nrow(lung.2x) # 334 rows
model.2x <- coxph(Surv(time, status) ~ pspline(inst) + pspline(age)
                  + pspline(ph.ecog) + pspline(ph.karno) + pspline(pat.karno)
                  + pspline(meal.cal) + pspline(wt.loss) + sex, data=lung.2x,
                  method='breslow')

This runs.

Could anyone provide some insight as to why the model won't fit on the smaller dataset? Looking at model.2x$coefficients, I see that there are 85 coefficients estimated which is still fewer than the number of observations in the smaller dataset.

(If it helps - the default values for pspline() are 4 degrees of freedom and degree = 3, i.e. a cubic spline.)

Answer 1

You are trying to apply splines to predictors that aren't continuous variables, while not recognizing that one of those variables had only a limited number of distinct values.

In the lung data "inst" is an institution code that should be treated as an unordered 17-level factor or, better, as a random effect; "ph.ecog" is a potentially 5-level ordered performance score that only takes on 4 levels in your complete-cases subset.

I have no problem with fitting the model when those 2 predictors are handled appropriately: frailty(inst) to treat institutions as a random effect, and ordered(ph.ecog) to use ECOG score as an ordered categorical predictor. The model:

coxph(formula = Surv(time, status) ~ frailty(inst) + pspline(age) + 
    ordered(ph.ecog) + pspline(ph.karno) + pspline(pat.karno) + 
    pspline(meal.cal) + pspline(wt.loss) + sex, data = lung, 
    method = "breslow")

is probably terribly overfit with only 120 events, but there is no convergence problem.

Your convergence problem seems to have stemmed from the "ph.ecog" predictor. It only has 4 levels in your complete-cases subset, and the default for pspline()is thus too many at 4 df. Taking your original model and just changing the corresponding term to pspline(ph.ecog,2) led to no technical difficulty.