glmmTMB 是否返回随机效应方差分量(如 glmmADMB)的标准误差?

Question

注意：以下分析仅适用于可重现的对象 - 并不是分析来自MASSbacteria数据的合法方法。

library(glmmADMB)
library(glmmTMB)
data(bacteria,package="MASS")
bacteria$present <- as.numeric(bacteria$y)-1
bacteria$early <- factor(as.numeric(bacteria$week > 3) + 1)
bfit2 <-  glmmadmb(present ~ trt , 
random = ~ (1 | ID) + (1 | early),
family = "binomial", data = bacteria)
bfit2$S
bfit2$sd_S
bfit3 <- glmmTMB(present ~ trt + (1 | ID) + (1 | early),
family = "binomial", data = bacteria)
summary(bfit3)$varcor
confint(bfit3)

在我的理解中，

• bfit2$S包含对glmmADMB随机效应方差的估计
值

• bfit3$varcor包含对glmmTMB随机效应的标准差的估计值(即bfit2$S元素的sqrt())

• bfit2$sd_S包含glmmADMB的bfit2$S估计值的标准误差(如本SO帖子中所述)

存储在glmmTMB对象中的bfit3$varcor的标准错误在哪里？更新：confint是为glmmTMB对象实现的，所以如果计算95%CI是最终目标，那么这是可用的(感谢kaskr)。

> bfit2$S
$ID
(Intercept)
(Intercept)      1.4047
$early
(Intercept)
(Intercept)     0.51465
> bfit2$sd_S
$ID
(Intercept)
(Intercept)     0.94743
$early
(Intercept)
(Intercept)     0.61533
> summary(bfit3)$varcor
Conditional model:
Groups Name        Std.Dev.
ID     (Intercept) 1.18513 
early  (Intercept) 0.71733 
> confint(bfit3)
2.5 %     97.5 %   Estimate
cond.(Intercept)                1.1934429  4.1351114  2.6642771
cond.trtdrug                   -2.6375284 -0.0503639 -1.3439462
cond.trtdrug+                  -2.0819454  0.5325683 -0.7746885
cond.Std.Dev.ID.(Intercept)     0.6118984  2.2953603  1.1851276
cond.Std.Dev.early.(Intercept)  0.2222685  2.3150437  0.7173293

如我们所见，sqrt(1.4047) = 1.18513和sqrt(0.51465) = 0.71733表示bfit2$S给出方差的估计值，summary(bfit3)$varcor给出标准差的估计值。

第二次更新：

经过一番挖掘，我意识到bfit3$sdr返回了对数-sd-scale上的方差分量，以及标准误差。因此，一个想法是避免confint并通过在对数 sd 尺度上计算 95%CI，然后转换为所需的尺度，然后将 CI 的宽度除以 2*1.96 来回计算 SE。

## to get the standard errors from glmmTMB:
bfit3$sdr
## note that theta is just log(sd)
exp(summary(bfit3$sdr, "fixed")[4,1])
exp(summary(bfit3$sdr, "fixed")[5,1])
## calculate the (wald) lower and upper on the log(sd) scale:
low.log.sd.id    <- summary(bfit3$sdr, "fixed")[4,1] - 1.96*summary(bfit3$sdr, "fixed")[4,2]
low.log.sd.early <- summary(bfit3$sdr, "fixed")[5,1] - 1.96*summary(bfit3$sdr, "fixed")[5,2]
upp.log.sd.id    <- summary(bfit3$sdr, "fixed")[4,1] + 1.96*summary(bfit3$sdr, "fixed")[4,2]
upp.log.sd.early <- summary(bfit3$sdr, "fixed")[5,1] + 1.96*summary(bfit3$sdr, "fixed")[5,2]
## convert to variance scale by taking exp and then squaring
low.var.id <- exp(low.log.sd.id)^2
upp.var.id <- exp(upp.log.sd.id)^2
low.var.early <- exp(low.log.sd.early)^2
upp.var.early <- exp(upp.log.sd.early)^2
## back calculate SEs
(upp.var.id - low.var.id) / (2*1.96)
(upp.var.early - low.var.early) / (2*1.96)
## see how they compare to the confint answers for sd:
sqrt(c(low.var.id, upp.var.id))
sqrt(c(low.var.early, upp.var.early))

运行它：

> ## to get the standard errors from glmmTMB:
> bfit3$sdr
sdreport(.) result
Estimate Std. Error
beta   2.6642771  0.7504394
beta  -1.3439462  0.6600031
beta  -0.7746885  0.6669800
theta  0.1698504  0.3372712
theta -0.3322203  0.5977910
Maximum gradient component: 4.83237e-06 
> ## note that theta is just log(sd)
> exp(summary(bfit3$sdr, "fixed")[4,1])
[1] 1.185128
> exp(summary(bfit3$sdr, "fixed")[5,1])
[1] 0.7173293
> ## calculate the (wald) lower and upper on the log(sd) scale:
> low.log.sd.id    <- summary(bfit3$sdr, "fixed")[4,1] - 1.96*summary(bfit3$sdr, "fixed")[4,2]
> low.log.sd.early <- summary(bfit3$sdr, "fixed")[5,1] - 1.96*summary(bfit3$sdr, "fixed")[5,2]
> upp.log.sd.id    <- summary(bfit3$sdr, "fixed")[4,1] + 1.96*summary(bfit3$sdr, "fixed")[4,2]
> upp.log.sd.early <- summary(bfit3$sdr, "fixed")[5,1] + 1.96*summary(bfit3$sdr, "fixed")[5,2]
> ## convert to variance scale by taking exp and then squaring
> low.var.id <- exp(low.log.sd.id)^2
> upp.var.id <- exp(upp.log.sd.id)^2
> low.var.early <- exp(low.log.sd.early)^2
> upp.var.early <- exp(upp.log.sd.early)^2
> ## back calculate SEs
> (upp.var.id - low.var.id) / (2*1.96)
[1] 1.24857
> (upp.var.early - low.var.early) / (2*1.96)
[1] 1.354657
> ## see how they compare to the confint answers for sd:
> sqrt(c(low.var.id, upp.var.id))
[1] 0.611891 2.295388
> sqrt(c(low.var.early, upp.var.early))
[1] 0.2222637 2.3150935

上面的最后两行与 confint(bfit3) 输出的最后两行非常匹配。 现在我想我只是想知道为什么 glmmADMB 的 SE 是 0.94743 和 0.61533，而 glmmTMB 的反向计算分别为 1.24857 和 1.354657...(?)

Answer 1

不是 100% 确定您的分析，但这是我所做的检查(包括挖掘glmmADMB的内脏并使用glmmTMB稍微模糊的方面)：

• 运行glmmADMB并深入了解 ADMB.std输出文件以检查结果：

n par 估计 sd 4 tmpL
1.6985e-01 3.3727e-01 5
tmpL -3.2901e-01 5.9780e-01
...
62 S 1.4045 9.474e-01 63 S 5.178e-01 6.191e-01

这些线分别是内部参数(tmpL：对数标准差)和变换参数(方差)。

重新执行转换并检查：

tmpL <-    c(0.16985,-0.32901)
cbind(admb.raw=exp(2*tmpL),
admb=unlist(bfit2$S),
tmb =unlist(VarCorr(bfit3)))
##        admb.raw    admb       tmb
## ID    1.4045262 1.40450 1.4045274
## early 0.5178757 0.51788 0.5145613

我们从随机效应的 log-std dev 的标准差乘以变换的导数得到方差的标准差(V = exp(2*logsd)，所以 dV/d(logsd) = 2*exp(2*logsd))

## sd of log-sd from glmmTMB
tmb_sd <- sqrt(diag(vcov(bfit3,full=TRUE)))[4:5]
tmb_logsd <- bfit3$sdr$par.fixed[4:5]
tmpL_sd <- c(0.33727,0.5978)
cbind(admb.raw=tmpL_sd*2*exp(2*tmpL),
admb=unlist(bfit2$sd_S),
tmb=tmb_sd*2*exp(2*tmb_logsd))
##        admb.raw    admb       tmb
## ID    0.9474091 0.94741 0.9474132
## early 0.6191722 0.61918 0.6152003

所以这些似乎都匹配得很好。

Answer 2

$varcor对象是一个列表。标准差存储为属性：

str( summary(bfit3)$varcor )
List of 2
$ cond:List of 2
..$ ID   : num [1, 1] 1.4
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr "(Intercept)"
.. .. ..$ : chr "(Intercept)"
.. ..- attr(*, "stddev")= Named num 1.19
.. .. ..- attr(*, "names")= chr "(Intercept)"
.. ..- attr(*, "correlation")= num [1, 1] 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr "(Intercept)"
.. .. .. ..$ : chr "(Intercept)"
.. ..- attr(*, "blockCode")= Named num 1
.. .. ..- attr(*, "names")= chr "us"
..$ early: num [1, 1] 0.515
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr "(Intercept)"
.. .. ..$ : chr "(Intercept)"
.. ..- attr(*, "stddev")= Named num 0.717
.. .. ..- attr(*, "names")= chr "(Intercept)"
.. ..- attr(*, "correlation")= num [1, 1] 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr "(Intercept)"
.. .. .. ..$ : chr "(Intercept)"
.. ..- attr(*, "blockCode")= Named num 1
.. .. ..- attr(*, "names")= chr "us"
..- attr(*, "sc")= num 1
..- attr(*, "useSc")= logi FALSE
$ zi  : NULL
- attr(*, "sc")= logi FALSE
- attr(*, "class")= chr "VarCorr.glmmTMB"

这将遍历该对象的.$cond元素：

sapply( summary(bfit3)$varcor$cond, function(x) attr( x, "stddev") )
ID.(Intercept) early.(Intercept) 
1.1851276         0.7173293