2013-04-18 1 views
1

knitr/xtable을 사용하여 보고서의 coxph() 개체에서 테이블을 생성하려고합니다. 모델에 pspline 용어를 포함하지 않으면 모든 것이 예상대로 작동합니다. 하나의 덩어리에서 : 나는 범 스플라인 기간을 포함 할 때coxph 모델 변경 요약 테이블에 pspline 용어 추가 R

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

는하지만, 요약() 객체 변경의 구조와 $conf.int$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 

답변

2

(또는 두 개 또는 세) (테이블에 포함하기에 적합 맞는 벌칙을 스플라인 용어에 대한 신뢰 구간을 설명하는 의미가, 내가 확실히 할 것 편집 : 이 아닌이 아니라고 말하면서 confint에 의해 생성 된 간격의 긴 목록은 의미가 있다고 생각합니다. 그래픽 디스플레이를 요구하는 유사한 질문이 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")) 
BTW

enter image description here

: invisible() 제외 summary(fit0)에 의해 반환 아무것도 그래서 당신은 str(summary(fit1))에서보고있는 모든 전송 출력은 없다 콘솔은 cat 콜에 의해 그 외로운 작은 NULL이 뒤 따른다. 내 증거를 의심한다면 getAnywhere(summary.coxph.penal) 코드를 검토하십시오.

+0

이것은 스플라인 용어에 분명히 도움이되지만 요약을 사용하여 다른 용어의 값에 액세스 할 수는 없습니다. 아마도 나는 그것을 "손으로"계산해야 할 것인가? –

+0

'predict.coxph'를 사용하지 않는 이유가 있습니까? 그것은'se.fit' 인자를 가지며'newdata'를 받아들입니다. 아마도 당신은 관찰 할 수있는 관점에서 테이블에 무엇을보고 싶은지를 지정해야합니다. 예를 들어, 45,55,65,75 세의 예측 세트를 말할까요? 테 르노가 이미 너를 위해 열심히 노력했다고 의심한다. –