如何在不加载coxme包时计算coxme模型的std.error

背景 #

我正尝试将他人编写的函数整合到自己的R包中，但遇到了问题。

问题 #

计算std.error的代码会报错，但加载coxme包后错误就消失了。编写R包函数时无法使用library(coxme)，需要找到成功计算std.error的方法。

可复现代码 #

可复现示例如下：

library(survival)
# 拟合coxme模型
fit <- coxme::coxme(Surv(time, status) ~ ph.ecog + age + (1|inst), lung)

# estimate、nvar和nfrail用于计算std.error
estimate <- coxme::fixef(fit)
nvar <- base::length(estimate)
nfrail <- base::nrow(fit$var) - nvar

运行以下代码会报错：

std.error <- base::sqrt(diag(fit$var)[nfrail + 1:nvar])
# 错误信息：
# Error in as.integer(x) : 
#  cannot coerce type 'S4' to vector of type 'integer'

加载coxme包后，代码可成功运行：

library(coxme)
std.error <- base::sqrt(diag(fit$var)[nfrail + 1:nvar]) # 运行成功

已尝试的操作 #

• 曾以为diag()函数来自coxme包，但确认后发现coxme中使用的diag()实际来自base包。
• 计算std.error的代码参考自Stack Overflow的相关回答。

解决方案 #

问题核心是fit$var是coxme包定义的S4对象，未加载coxme包时，base包的diag()无法识别该类的处理方法。以下是三种可行解决办法：

方法1：显式调用coxme的diag方法 #

无需加载coxme包，直接调用其内部的diag.coxmeVar方法处理fit$var：

std.error <- base::sqrt(diag(coxme:::diag.coxmeVar(fit$var))[nfrail + 1:nvar])

方法2：添加coxme为包依赖 #

在R包的DESCRIPTION文件中添加：

Imports: coxme

之后在代码中使用coxme::前缀调用相关方法，确保S4对象的方法能被正确识别。

方法3：手动提取协方差矩阵 #

直接访问S4对象的底层数据，fit$var的协方差矩阵存储在@x槽中，可手动构造矩阵后计算：

# 构造完整协方差矩阵
cov_matrix <- matrix(fit$var@x, nrow = fit$var@Dim[1], ncol = fit$var@Dim[2])
# 计算标准误
std.error <- base::sqrt(diag(cov_matrix)[nfrail + 1:nvar])

内容的提问来源于stack exchange，提问作者zhiwei li