如何在多变量 Coxph 中获取特定变量的 wald 检验？
How to get the wald test of a specific variable in a multivariate Coxph?
我使用 R 生存包拟合了一个多变量 Cox 模型，如下所示：
library(survival)
data(lung)
res.cox1 <- coxph(Surv(time, status) ~ sex + ph.karno + wt.loss, data = lung)
res.cox1
Call:
coxph(formula = Surv(time, status) ~ sex + ph.karno + wt.loss,
data = lung)
coef exp(coef) se(coef) z p
sex -0.521839 0.593428 0.174454 -2.991 0.00278
ph.karno -0.015243 0.984873 0.005988 -2.546 0.01091
wt.loss -0.002523 0.997480 0.006233 -0.405 0.68558
Likelihood ratio test=16.42 on 3 df, p=0.0009298
n= 214, number of events= 152
(14 observations deleted due to missingness)
如何在多变量 Cox 模型 (sex + ph.karno + wt.loss )?
我试图查看 coxph 的结构和 coxph 对象的摘要，但我发现 wald 测试只有一个值 $wald.test : num 16.5、$ waldtest : Named num [1:3] 1.65e+01 3.00 8.81e-04 ..- attr(*, "names")= chr [1:3] "test" "df" "pvalue" !
这个测试值对应的是什么？
如何获得 Wald 性别检验的 3 个值，ph.karno 和 wt.loss ?
str(res.cox1)
List of 20
$ coefficients : Named num [1:3] -0.52184 -0.01524 -0.00252
..- attr(*, "names")= chr [1:3] "sex" "ph.karno" "wt.loss"
$ var : num [1:3, 1:3] 3.04e-02 -6.78e-05 2.77e-05 -6.78e-05 3.59e-05 ...
$ loglik : num [1:2] -680 -672
$ score : num 16.9
$ iter : int 4
$ linear.predictors: num [1:214] 0.0756 0.0756 0.0857 -0.039 0.7232 ...
$ residuals : Named num [1:214] -0.147 -2.93 0.58 -1.613 -5.599 ...
..- attr(*, "names")= chr [1:214] "2" "3" "4" "5" ...
$ means : Named num [1:3] 1.4 82.06 9.83
..- attr(*, "names")= chr [1:3] "sex" "ph.karno" "wt.loss"
$ method : chr "efron"
$ n : int 214
$ nevent : num 152
$ terms :Classes 'terms', 'formula' language Surv(time, status) ~ sex + ph.karno + wt.loss
.. ..- attr(*, "variables")= language list(Surv(time, status), sex, ph.karno, wt.loss)
.. ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "Surv(time, status)" "sex" "ph.karno" "wt.loss"
.. .. .. ..$ : chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..- attr(*, "term.labels")= chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..- attr(*, "specials")=Dotted pair list of 2
.. .. ..$ strata: NULL
.. .. ..$ tt : NULL
.. ..- attr(*, "order")= int [1:3] 1 1 1
.. ..- attr(*, "intercept")= num 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(Surv(time, status), sex, ph.karno, wt.loss)
.. ..- attr(*, "dataClasses")= Named chr [1:4] "nmatrix.2" "numeric" "numeric" "numeric"
.. .. ..- attr(*, "names")= chr [1:4] "Surv(time, status)" "sex" "ph.karno" "wt.loss"
$ assign :List of 3
..$ sex : int 1
..$ ph.karno: int 2
..$ wt.loss : int 3
$ wald.test : num 16.5
$ concordance : Named num [1:7] 11071 6046 96 22 0 ...
..- attr(*, "names")= chr [1:7] "concordant" "discordant" "tied.x" "tied.y" ...
$ na.action : 'omit' Named int [1:14] 1 20 36 44 56 63 108 138 178 183 ...
..- attr(*, "names")= chr [1:14] "1" "20" "36" "44" ...
$ y : 'Surv' num [1:214, 1:2] 455 1010+ 210 883 1022+ 310 361 218 166 170 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:214] "2" "3" "4" "5" ...
.. ..$ : chr [1:2] "time" "status"
..- attr(*, "type")= chr "right"
$ timefix : logi TRUE
$ formula :Class 'formula' language Surv(time, status) ~ sex + ph.karno + wt.loss
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
$ call : language coxph(formula = Surv(time, status) ~ sex + ph.karno + wt.loss, data = lung)
- attr(*, "class")= chr "coxph"
str(summary(res.cox1))
List of 14
$ call : language coxph(formula = Surv(time, status) ~ sex + ph.karno + wt.loss, data = lung)
$ fail : NULL
$ na.action : 'omit' Named int [1:14] 1 20 36 44 56 63 108 138 178 183 ...
..- attr(*, "names")= chr [1:14] "1" "20" "36" "44" ...
$ n : int 214
$ loglik : num [1:2] -680 -672
$ nevent : num 152
$ coefficients: num [1:3, 1:5] -0.52184 -0.01524 -0.00252 0.59343 0.98487 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..$ : chr [1:5] "coef" "exp(coef)" "se(coef)" "z" ...
$ conf.int : num [1:3, 1:4] 0.593 0.985 0.997 1.685 1.015 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..$ : chr [1:4] "exp(coef)" "exp(-coef)" "lower .95" "upper .95"
$ logtest : Named num [1:3] 16.42029 3 0.00093
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ sctest : Named num [1:3] 1.69e+01 3.00 7.52e-04
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ rsq : Named num [1:2] 0.0739 0.9983
..- attr(*, "names")= chr [1:2] "rsq" "maxrsq"
$ waldtest : Named num [1:3] 1.65e+01 3.00 8.81e-04
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ used.robust : logi FALSE
$ concordance : Named num [1:2] 0.646 0.0274
..- attr(*, "names")= chr [1:2] "C" "se(C)"
- attr(*, "class")= chr "summary.coxph"
谢谢！
http://www.sthda.com/english/wiki/cox-proportional-hazards-model
"z" 列与多变量 Cox 模型中每个协变量的 Wald 检验统计量相同。
也可以这样调用Cox模型统计：
summary(res.cox1)
A "Wald test" 基于回归过程的参数值呈正态分布的假设。您检查系数估计值 ("coef") 除以估计值的标准误差 ("coef(se)") 的比率，并查看该比率的 95% 置信区间是否包含零值。操作说明：采用 coef +/- 1.96*se(coef) 并查看间隔是否包括零。或者等价地，您可以取比率：coef/se(coef)，看看它的绝对值是否大于 1.96。当我说 "test" 是 yes/no 结果，回答问题 "does the ratio value lie in a critical interval or not"，而 "test statistic" 就像 z 值一样时，也许我是在迂腐纯数.
您构建的摘要中实际上报告了 4 个 Wald 检验。其中三个用于单个系数，其中一个用于整个模型，即名为 "wald" 的那个。但是您不想对整体模型进行 Wald 检验。您想要 summary() 处理结果的 "coefficient" 矩阵的结果（而不是 coxph() 结果的 "coefficient" 值。）当您采用此类比率时，它被分析为z 检验，因此您不对统计数据进行平方（当然，除非您想使用卡方 table，此时 Z^2 将用于评估。）
summ.coef <- summary(res.cox1)$coefficients
( Wald.ratios <- summ.coef[,"coef"]/summ.coef[,"se(coef)"] )
sex ph.karno wt.loss
-2.9912645 -2.5456273 -0.4048609
identical(Wald.ratios, summ.coef[, "z"])
#[1] TRUE
如果您想按名称关注单个变量：
summ.coef["sex", "coef"]/summ.coef["sex", "se(coef)"]
我使用 R 生存包拟合了一个多变量 Cox 模型，如下所示：
library(survival)
data(lung)
res.cox1 <- coxph(Surv(time, status) ~ sex + ph.karno + wt.loss, data = lung)
res.cox1
Call:
coxph(formula = Surv(time, status) ~ sex + ph.karno + wt.loss,
data = lung)
coef exp(coef) se(coef) z p
sex -0.521839 0.593428 0.174454 -2.991 0.00278
ph.karno -0.015243 0.984873 0.005988 -2.546 0.01091
wt.loss -0.002523 0.997480 0.006233 -0.405 0.68558
Likelihood ratio test=16.42 on 3 df, p=0.0009298
n= 214, number of events= 152
(14 observations deleted due to missingness)
如何在多变量 Cox 模型 (sex + ph.karno + wt.loss )?
我试图查看 coxph 的结构和 coxph 对象的摘要，但我发现 wald 测试只有一个值 $wald.test : num 16.5、$ waldtest : Named num [1:3] 1.65e+01 3.00 8.81e-04 ..- attr(*, "names")= chr [1:3] "test" "df" "pvalue" !
这个测试值对应的是什么？ 如何获得 Wald 性别检验的 3 个值，ph.karno 和 wt.loss ?
str(res.cox1)
List of 20
$ coefficients : Named num [1:3] -0.52184 -0.01524 -0.00252
..- attr(*, "names")= chr [1:3] "sex" "ph.karno" "wt.loss"
$ var : num [1:3, 1:3] 3.04e-02 -6.78e-05 2.77e-05 -6.78e-05 3.59e-05 ...
$ loglik : num [1:2] -680 -672
$ score : num 16.9
$ iter : int 4
$ linear.predictors: num [1:214] 0.0756 0.0756 0.0857 -0.039 0.7232 ...
$ residuals : Named num [1:214] -0.147 -2.93 0.58 -1.613 -5.599 ...
..- attr(*, "names")= chr [1:214] "2" "3" "4" "5" ...
$ means : Named num [1:3] 1.4 82.06 9.83
..- attr(*, "names")= chr [1:3] "sex" "ph.karno" "wt.loss"
$ method : chr "efron"
$ n : int 214
$ nevent : num 152
$ terms :Classes 'terms', 'formula' language Surv(time, status) ~ sex + ph.karno + wt.loss
.. ..- attr(*, "variables")= language list(Surv(time, status), sex, ph.karno, wt.loss)
.. ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "Surv(time, status)" "sex" "ph.karno" "wt.loss"
.. .. .. ..$ : chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..- attr(*, "term.labels")= chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..- attr(*, "specials")=Dotted pair list of 2
.. .. ..$ strata: NULL
.. .. ..$ tt : NULL
.. ..- attr(*, "order")= int [1:3] 1 1 1
.. ..- attr(*, "intercept")= num 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(Surv(time, status), sex, ph.karno, wt.loss)
.. ..- attr(*, "dataClasses")= Named chr [1:4] "nmatrix.2" "numeric" "numeric" "numeric"
.. .. ..- attr(*, "names")= chr [1:4] "Surv(time, status)" "sex" "ph.karno" "wt.loss"
$ assign :List of 3
..$ sex : int 1
..$ ph.karno: int 2
..$ wt.loss : int 3
$ wald.test : num 16.5
$ concordance : Named num [1:7] 11071 6046 96 22 0 ...
..- attr(*, "names")= chr [1:7] "concordant" "discordant" "tied.x" "tied.y" ...
$ na.action : 'omit' Named int [1:14] 1 20 36 44 56 63 108 138 178 183 ...
..- attr(*, "names")= chr [1:14] "1" "20" "36" "44" ...
$ y : 'Surv' num [1:214, 1:2] 455 1010+ 210 883 1022+ 310 361 218 166 170 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:214] "2" "3" "4" "5" ...
.. ..$ : chr [1:2] "time" "status"
..- attr(*, "type")= chr "right"
$ timefix : logi TRUE
$ formula :Class 'formula' language Surv(time, status) ~ sex + ph.karno + wt.loss
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
$ call : language coxph(formula = Surv(time, status) ~ sex + ph.karno + wt.loss, data = lung)
- attr(*, "class")= chr "coxph"
str(summary(res.cox1))
List of 14
$ call : language coxph(formula = Surv(time, status) ~ sex + ph.karno + wt.loss, data = lung)
$ fail : NULL
$ na.action : 'omit' Named int [1:14] 1 20 36 44 56 63 108 138 178 183 ...
..- attr(*, "names")= chr [1:14] "1" "20" "36" "44" ...
$ n : int 214
$ loglik : num [1:2] -680 -672
$ nevent : num 152
$ coefficients: num [1:3, 1:5] -0.52184 -0.01524 -0.00252 0.59343 0.98487 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..$ : chr [1:5] "coef" "exp(coef)" "se(coef)" "z" ...
$ conf.int : num [1:3, 1:4] 0.593 0.985 0.997 1.685 1.015 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "sex" "ph.karno" "wt.loss"
.. ..$ : chr [1:4] "exp(coef)" "exp(-coef)" "lower .95" "upper .95"
$ logtest : Named num [1:3] 16.42029 3 0.00093
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ sctest : Named num [1:3] 1.69e+01 3.00 7.52e-04
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ rsq : Named num [1:2] 0.0739 0.9983
..- attr(*, "names")= chr [1:2] "rsq" "maxrsq"
$ waldtest : Named num [1:3] 1.65e+01 3.00 8.81e-04
..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
$ used.robust : logi FALSE
$ concordance : Named num [1:2] 0.646 0.0274
..- attr(*, "names")= chr [1:2] "C" "se(C)"
- attr(*, "class")= chr "summary.coxph"
谢谢！
http://www.sthda.com/english/wiki/cox-proportional-hazards-model
"z" 列与多变量 Cox 模型中每个协变量的 Wald 检验统计量相同。
也可以这样调用Cox模型统计：
summary(res.cox1)
A "Wald test" 基于回归过程的参数值呈正态分布的假设。您检查系数估计值 ("coef") 除以估计值的标准误差 ("coef(se)") 的比率，并查看该比率的 95% 置信区间是否包含零值。操作说明：采用 coef +/- 1.96*se(coef) 并查看间隔是否包括零。或者等价地，您可以取比率：coef/se(coef)，看看它的绝对值是否大于 1.96。当我说 "test" 是 yes/no 结果，回答问题 "does the ratio value lie in a critical interval or not"，而 "test statistic" 就像 z 值一样时，也许我是在迂腐纯数.
您构建的摘要中实际上报告了 4 个 Wald 检验。其中三个用于单个系数，其中一个用于整个模型，即名为 "wald" 的那个。但是您不想对整体模型进行 Wald 检验。您想要 summary() 处理结果的 "coefficient" 矩阵的结果（而不是 coxph() 结果的 "coefficient" 值。）当您采用此类比率时，它被分析为z 检验，因此您不对统计数据进行平方（当然，除非您想使用卡方 table，此时 Z^2 将用于评估。）
summ.coef <- summary(res.cox1)$coefficients
( Wald.ratios <- summ.coef[,"coef"]/summ.coef[,"se(coef)"] )
sex ph.karno wt.loss
-2.9912645 -2.5456273 -0.4048609
identical(Wald.ratios, summ.coef[, "z"])
#[1] TRUE
如果您想按名称关注单个变量：
summ.coef["sex", "coef"]/summ.coef["sex", "se(coef)"]