我正在尝试使用 knitr/xtable 从报告中的 coxph() 对象生成表格。当我在模型中不包含 pspline 项时,一切都按预期工作。在单个块中:
<<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)
@
但是当我包含一个惩罚样条项时,summary() 对象的结构会发生变化,并且$conf.int
and$coeff
插槽不再可用。
> 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