我使用包riskRegression和函数selectCox来执行向后选择模型。
https://rdrr.io/github/tagteam/pec/man/selectCox.html
library(pec)
library(prodlim)
data(GBSG2)
library(survival)
library(riskRegression)
f <- selectCox(Surv(time,cens)~horTh+age+menostat+tsize+tgrade+pnodes+progrec+estrec ,
data=GBSG2)
是否可以显示 95% 置信区间？
我遇到了同样的问题，但 ChatGPT 和我通过解决方案互相帮助，因为 chatGPT 不知道如何从对象中提取相关信息： 首先，您必须从对象中提取系数，然后提取标准误差。对象只有一个方差-协方差矩阵。要从此方差-协方差矩阵计算系数的标准误差 (SE)，您需要提取对角线单元(表示方差)并取其平方根。对角线元素表示每个系数的平方标准误差。然后将系数以及 95% 置信区的上限和下限相增。 系数的名称保存在"命名数字"向量中
coef_estimates <- f[["fit"]][["coefficients"]]
var_cov_matrix <- f[["fit"]][["var"]]
# Calculate standard errors (S.E.) by taking the square root of the diagonal elements
se <- round(sqrt(diag(var_cov_matrix)),2)
z_score <- qnorm(0.975) # Z-score for 95% CI
# Calculate confidence intervals
lower_ci <- round(exp(coef_estimates - z_score * se),2)
upper_ci <- round(exp(coef_estimates + z_score * se),2)
# Print the results (chatGPT)
for (i in seq_along(coef_estimates)) {
coef_name <- names(coef_estimates)[i]
coef_value <- round(exp(coef_estimates),2)[i]
se_value <- se[i]
lower_ci_value <- lower_ci[i]
upper_ci_value <- upper_ci[i]
cat(paste0("Coefficient ", coef_name, ": ", coef_value, "n"))
cat(paste0("Standard Error: ", se_value, "n"))
cat(paste0("95% CI (Exponentiated): [", lower_ci_value, ", ", upper_ci_value, "]nn"))
}
它对我的数据非常有效！
我在 SO 上真诚地搜索了一个好的重复提名，但我没有找到。有趣的是，在selectCox对象中正确组件的先前答案中使用了chatGPT编码，因为它基本上是在重新发明一个名为confint的现有函数默认函数。让我们首先展示应该如何处理这个问题：在本例中，首先要查看名为f的返回对象的结构。
f
#----screen outputs---
$fit
Cox Proportional Hazards Model
rms::cph(formula = newform, data = data, surv = TRUE)
Model Tests Discrimination
Indexes
Obs 686 LR chi2 99.16 R2 0.135
Events 299 d.f. 5 R2(5,686)0.128
Center 0.526 Pr(> chi2) 0.0000 R2(5,299)0.270
Score chi2 115.95 Dxy 0.375
Pr(> chi2) 0.0000
Coef S.E. Wald Z Pr(>|Z|)
horTh=yes -0.3203 0.1257 -2.55 0.0108
tgrade=II 0.6524 0.2488 2.62 0.0087
tgrade=III 0.8084 0.2678 3.02 0.0025
pnodes 0.0542 0.0068 8.01 <0.0001
progrec -0.0022 0.0006 -3.97 <0.0001
$In
[1] "horTh" "tgrade" "pnodes" "progrec"
$call
selectCox(formula = Surv(time, cens) ~ horTh + age + menostat +
tsize + tgrade + pnodes + progrec + estrec, data = GBSG2)
attr(,"class")
[1] "selectCox"
####---- end of output ---
因此，显然有一个名为print.selectCox的函数，它将计算结果捆绑在第一个 copmponet 上，并将它们显示在标题"$fit"下。现在看看实际结构：
str(f$fit)
####--- somewhat more informative ---
List of 30
$ coefficients : Named num [1:5] -0.32028 0.65239 0.80839 0.05425 -0.00221
..- attr(*, "names")= chr [1:5] "horTh=yes" "tgrade=II" "tgrade=III" "pnodes" ...
$ var : num [1:5, 1:5] 1.58e-02 -3.57e-05 1.55e-03 2.82e-05 -3.77e-07 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:5] "horTh=yes" "tgrade=II" "tgrade=III" "pnodes" ...
.. ..$ : chr [1:5] "horTh=yes" "tgrade=II" "tgrade=III" "pnodes" ...
$ loglik : num [1:2] -1788 -1739
$ score : num 116
$ iter : int 5
$ linear.predictors: Named num [1:686] 0.1831 0.0511 0.1795 -0.1094 0.1232 ...
..- attr(*, "names")= chr [1:686] "1" "2" "3" "4" ...
$ residuals : Named num [1:686] 0.18 0.169 0.705 0.399 0.683 ...
..- attr(*, "names")= chr [1:686] "1" "2" "3" "4" ...
$ means : num [1:5] 0.359 0.647 0.235 5.01 109.996
### snipped the remaining items in the list
因此，第一项名为"fit"，它本身就是一个列表，具有名为coefficients的子组件，第二项名为var。
名为confint.default的基本 R 函数可以从任何名为"系数"和"var"的 R 对象返回置信区间，只要其中第二个是维度与"系数"长度对应的方阵。 如果您尝试在整个selectCox类对象上使用它，它会失败，如果您只在fit部分使用它就成功了：
> confint(f)
Error in vcov.default(object) :
object does not have variance-covariance matrix
> confint(f$fit)
2.5 % 97.5 %
horTh=yes -0.566648468 -0.073905052
tgrade=II 0.164830108 1.139941870
tgrade=III 0.283558082 1.333218011
pnodes 0.040965920 0.067528158
progrec -0.003298602 -0.001117193
但是您需要将这些值相加，以便在风险比率尺度上：
> exp( confint(f$fit) )
2.5 % 97.5 %
horTh=yes 0.5674240 0.9287599
tgrade=II 1.1791928 3.1265866
tgrade=III 1.3278460 3.7932304
pnodes 1.0418166 1.0698604
progrec 0.9967068 0.9988834