Adding pspline term to coxph model changes summary table in R

Question

I'm trying to use knitr/xtable to produce tables from coxph() objects in a report. When I don't include a pspline term in the model everything works as expected. In a single chunk:

<<results = 'asis'>>=
require(survival, quietly = T)
require(xtable, quietly = T)
data(cancer,  package = "survival")
fit0 <- coxph(Surv(time, status) ~ meal.cal + ph.ecog + age, cancer)
# construct data frame for tables - no spline
fit0table <- data.frame(Variable = c("Calories Consumed", "ECOG Performance Score","Age"), RiskRatio = summary(fit0)$conf.int[,1], Lower = summary(fit0)$conf.int[,3],   Upper = summary(fit0)$conf.int[,4], Pval = summary(fit0)$coeff[,5])
# print latex table
print(xtable(fit0table, digits = 3), include.rownames = F)
@

But when I include a penalized spline term, the structure of the summary() object changes and the $conf.int and $coeff slots are no longer available.

> fit1 <- coxph(Surv(time, status) ~ meal.cal + ph.ecog + pspline(age, 3), cancer)
> str(summary(fit0))
List of 14
$ call        : language coxph(formula = Surv(time, status) ~ meal.cal + ph.ecog + age, data = cancer)
$ fail        : NULL 
$ na.action   :Class 'omit'  Named int [1:48] 3 5 12 13 14 16 23 25 33 44 ...
  .. ..- attr(*, "names")= chr [1:48] "3" "5" "12" "13" ...
$ n           : int 180
$ loglik      : num [1:2] -574 -567
$ nevent      : num 133
$ coefficients: num [1:3, 1:5] 3.84e-05 4.00e-01 1.10e-02 1.00 1.49 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "meal.cal" "ph.ecog" "age"
  .. ..$ : chr [1:5] "coef" "exp(coef)" "se(coef)" "z" ...
$ conf.int    : num [1:3, 1:4] 1 1.491 1.011 1 0.671 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "meal.cal" "ph.ecog" "age"
.. ..$ : chr [1:4] "exp(coef)" "exp(-coef)" "lower .95" "upper .95"
$ logtest     : Named num [1:3] 13.2142 3 0.0042
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ sctest      : Named num [1:3] 13.46468 3 0.00373
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ rsq         : Named num [1:2] 0.0708 0.9983
..- attr(*, "names")= chr [1:2] "rsq" "maxrsq"
$ waldtest    : Named num [1:3] 13.28 3 0.00407  
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ used.robust : logi FALSE
$ concordance : Named num [1:2] 0.6061 0.0291
..- attr(*, "names")= chr [1:2] "concordance.concordant" "se.std(c-d)"
- attr(*, "class")= chr "summary.coxph"


> str(summary(fit1))
Call: 
coxph(formula = Surv(time, status) ~ meal.cal + ph.ecog + pspline(age, 
3), data = cancer)

n= 180, number of events= 133 
(48 observations deleted due to missingness)

                        coef     se(coef) se2      Chisq DF   p     
meal.cal                3.65e-05 0.000228 0.000228 0.03  1.00 0.8700
ph.ecog                 3.98e-01 0.131938 0.131738 9.10  1.00 0.0026
pspline(age, 3), linear 1.07e-02 0.010694 0.010694 1.00  1.00 0.3200
pspline(age, 3), nonlin                            2.90  2.07 0.2500

          exp(coef) exp(-coef) lower .95 upper .95
meal.cal       1.00     1.0000     1.000      1.00
ph.ecog        1.49     0.6717     1.150      1.93
ps(age)3       1.75     0.5717     0.473      6.47
ps(age)4       3.03     0.3302     0.365     25.14
ps(age)5       4.49     0.2228     0.395     50.96
ps(age)6       4.65     0.2150     0.405     53.43
ps(age)7       3.96     0.2526     0.363     43.12
ps(age)8       3.84     0.2604     0.360     41.01
ps(age)9       4.44     0.2250     0.413     47.84
ps(age)10      5.39     0.1855     0.486     59.82
ps(age)11      7.94     0.1260     0.599    105.23
ps(age)12     12.25     0.0816     0.537    279.91

Iterations: 4 outer, 12 Newton-Raphson
     Theta= 0.836 
Degrees of freedom for terms= 1.0 1.0 3.1 
Concordance= 0.616  (se = 0.029 )
Rsquare= 0.092   (max possible= 0.998 )
Likelihood ratio test= 17.5  on 5.06 df,   p=0.00389
Wald test            = 15.9  on 5.06 df,   p=0.0073
 NULL


> coefficients(fit1) # doesn't give p-values
meal.cal      ph.ecog     ps(age)3     ps(age)4     ps(age)5     ps(age)6         ps(age)7     ps(age)8 
3.647054e-05 3.980039e-01 5.590767e-01 1.108052e+00 1.501557e+00 1.537249e+00     1.375833e+00 1.345564e+00 
ps(age)9    ps(age)10    ps(age)11    ps(age)12 
1.491454e+00 1.684622e+00 2.071641e+00 2.505932e+00 

> confint(fit1) # getting closer
                  2.5 %       97.5 %
meal.cal  -0.0004104346 0.0004833757
ph.ecog    0.1394097867 0.6565980826
ps(age)3  -0.7493022459 1.8674555679
ps(age)4  -1.0084545140 3.2245588375
ps(age)5  -0.9278798219 3.9309933396
ps(age)6  -0.9038092211 3.9783071434
ps(age)7  -1.0122388810 3.7639051908
ps(age)8  -1.0226368192 3.7137644105
ps(age)9  -0.8849251510 3.8678337954
ps(age)10 -0.7221442743 4.0913878825
ps(age)11 -0.5129062130 4.6561876883
ps(age)12 -0.6226068259 5.6344701023

Answer 1

I do not think there is a single number (or even two or three) that would be meaningful to describe the confidence interval for a penalized spline term fit suitable for inclusion in a table, and I certainly do (Edit: meant to say not) think that long list of intervals produced by confint is meaningful. (There is no confint.coxph.penal function.) When a similar question (albeit one asking for a graphical display) was posed 7 years ago on R-help, Terry Therneau posted this code for displaying what he thought would be a meaningful, which I have modified to fit your names and display the fit and CI for 'age':

fit1 <- coxph(Surv(time, status) ~ meal.cal + ph.ecog + pspline(age, 3), na.omit(cancer) )
temp <- predict(fit1, type='terms', se=TRUE) 
matplot(na.omit(cancer)$age, exp(cbind( temp$fit[, 3], 
                                  temp$fit[,3] - 2* temp$se.fit[,3], 
                                  temp$fit[,3] + 2* temp$se.fit[,3])), 
        log='y', xlab="Age", ylab="Estimated Relative Risk", col=c('red',"blue","blue") )

BTW: There is nothing returned by summary(fit0) except invisible(), so all you are seeing from str(summary(fit1)) is the output sent to the console by the cat calls followed by that lonely little NULL. If you doubt my veracty just review the code with getAnywhere(summary.coxph.penal).