如何提取JointModel竞争风险联合模型随机效应的标准误

1.1 线性混合模型 #

mod_mml=lme(MDRD~DELAI_MONTH, data=Mixed, random=~1+DELAI_MONTH|N_ID,na.action=na.exclude)

1.2 包含死亡、移植两类竞争风险事件的数据转换 #

SurvieCR=crLong(Survie,"Statut1_60m","Vivant 60m")

1.3 竞争风险生存子模型 #

生存子模型需采用长格式数据，将竞争风险标识作为分层因子纳入：

Surv_comp=coxph(Surv(DELAI_M,status2) ~ 1*strata + strata(strata),
                data=SurvieCR, x=TRUE,method="efron",na.action=na.exclude)

1.4 联合模型拟合 #

将上述两个子模型作为核心传入参数，设置CompRisk=TRUE即可拟合对应联合模型：

JM_comp=jointModel(mod_mml, Surv_comp, timeVar="DELAI_MONTH",
                   method="spline-PH-aGH",
                   interFact=list(value= ~ strata, data=SurvieCR),
                   CompRisk=TRUE)

已尝试调用如下函数提取结果，但未获得目标值：

summary(JM_comp)
random.effects(JM_comp)
coef(JM_comp) 
coef(JM_comp, process = "Longitudinal")

调用random.effects(JM_comp)仅返回各受试者对应截距、DELAI_MONTH斜率的随机效应点估计值，输出示例如下：

> random.effects(JM_comp)
     (Intercept)   DELAI_MONTH
1    3.727740214 -0.2020818077
2    1.070965087  0.0370529800
3    3.440276448  0.1582811172
4   -4.551210225 -0.1568204062
5   -2.943917109 -0.0389127841

需要获取每个个体随机效应中截距、斜率对应的标准差（标准误），寻求可行的实现方案。

JM包返回的jointModel对象中，个体随机效应的后验协方差数组存储在posterior.var.b字段，结构为「个体数 × 随机效应个数 × 随机效应个数」的三维数组，每个二维切片对应单个受试者的随机效应后验协方差矩阵，取矩阵对角线元素开平方即可得到对应随机效应的后验标准差（即需要的标准误）。
直接运行以下代码即可批量提取并合并结果：

# 提取已有的随机效应点估计
re_est <- random.effects(JM_comp)
# 遍历所有个体，计算每个随机效应的标准误
re_se <- t(
  apply(JM_comp$posterior.var.b, 1, function(cov_mat){
    sqrt(diag(cov_mat))
  })
)
# 列名对齐
colnames(re_se) <- c("(Intercept)_se", "DELAI_MONTH_se")
# 合并点估计和标准误为最终结果表
re_final <- cbind(re_est, re_se)

说明：如果需要的是总体层面随机效应分布的截距、斜率标准差，而非个体层面估计值的标准误，直接提取JM_comp$coefficients$D随机效应总体协方差矩阵，对对角线元素开平方即可。

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