1

我正在尝试使用 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.intand$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
4

1 回答 1

2

我认为没有一个数字(甚至两个或三个)可以有意义地描述适合包含在表格中的惩罚样条项的置信区间,我当然会(编辑:意思是说不是)认为 confint 产生的一长串区间是有意义的。(没有confint.coxph.penal功能。)当 7 年前在 R-help 上提出类似的问题(尽管是要求图形显示)时,Terry Therneau 发布了这段代码,用于显示他认为有意义的内容,我已将其修改为适合您的姓名并显示“年龄”的适合度和 CI:

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") )

在此处输入图像描述

summary(fit0)顺便说一句: except没有返回任何内容invisible(),因此您所看到的只是通过调用str(summary(fit1))发送到控制台的输出,然后是那个孤独的 little 。如果您怀疑我的真实性,只需使用.catNULLgetAnywhere(summary.coxph.penal)

于 2013-04-18T17:13:47.593 回答