如何在R的ordinal包中获取累积混合模型随机截距标准差的标准误？

数据与模型拟合 #

数据集 #

# A tibble: 7,500 × 4
   cluster cluster_size     x ordinal_y
     <int>        <int> <dbl> <fct>    
 1       1            1  2.57 4        
 2       1            2  3.59 4        
 3       1            3  2.51 1        
 4       1            4  3.80 3        
 5       1            5  2.84 4        
 6       1            6  2.59 1 

拟合的模型 #

mod1 <- ordinal::clmm(ordinal_y ~ x + (1 | cluster),
                      data = dat1,
                      link = "logit")

模型结果 #

mod1

Cumulative Link Mixed Model fitted with the Laplace approximation

formula: ordinal_y ~ x + (1 | cluster)
data:    dat1

 link  threshold nobs logLik   AIC      niter     max.grad
 logit flexible  7500 -7634.28 15278.56 364(1095) 1.45e-02

Random effects:
 Groups  Name        Variance Std.Dev.
 cluster (Intercept) 0.2913   0.5398  
Number of groups:  cluster 150 

Coefficients:
     x 
0.7586 

Thresholds:
   1|2    2|3    3|4 
0.5955 1.1092 1.6108 

通过mod1$ST获取随机截距标准差：

$cluster
            (Intercept)
(Intercept)   0.5397603

问题 #

想要获取该随机截距标准差的标准误，使用vcov(mod1)得到如下结果：

vcov(mod1)

             1|2          2|3          3|4            x          ST1
1|2 0.0038323677 0.0037161142 0.0036710729 0.0006564625 0.0001370114
2|3 0.0037161142 0.0039511707 0.0038837447 0.0007111263 0.0002270493
3|4 0.0036710729 0.0038837447 0.0041400738 0.0007645963 0.0003178104
x   0.0006564625 0.0007111263 0.0007645963 0.0003864496 0.0001702783
ST1 0.0001370114 0.0002270493 0.0003178104 0.0001702783 0.0058518536

请问ST1*ST1是我想要的结果吗？

不是，ST1*ST1不是你要的随机截距标准差的标准误，具体说明如下：

1. vcov矩阵中ST1的含义
   在ordinal::clmm的方差协方差矩阵里，ST1对应的是随机效应标准差的对数尺度估计值的方差——clmm默认采用对数变换来估计随机效应的标准差，以此保证估计值非负。

2. 正确计算标准误的方法
   要得到原始尺度下随机截距标准差的标准误，需要用delta方法做尺度转换：

   • 先提取对数标准差的方差：从vcov(mod1)中取ST1行ST1列的元素，即0.0058518536；
   • 再结合已得到的标准差估计值0.5397603，通过公式计算：

     sigma_hat <- mod1$ST$cluster[1,1]
     var_log_sigma <- vcov(mod1)["ST1", "ST1"]
     se_sigma <- sigma_hat * sqrt(var_log_sigma)
     

     计算结果约为0.0413。

3. 更简便的验证方式
   直接运行ordinal::summary(mod1)，输出结果中会直接给出随机效应标准差的标准误，无需手动计算。

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