我的问题可能看起来很简单,但我没有找到令人满意的解决方案。我已经在这个问题上停留了几天。如何在使用库中的glmer函数时获得固定和随机效应估计的协方差矩阵lme4。
我试过vcov(.., full = TRUE)没有成功。
有没有一种函数或方法来计算这个方差-协方差矩阵?
编辑
我需要的协方差矩阵是。对于回归参数估计,。归一化的部分似然信息矩阵的极限值。简而言之,我需要在处评估的观察到的逆信息矩阵。
我的问题可能看起来很简单,但我没有找到令人满意的解决方案。我已经在这个问题上停留了几天。如何在使用库中的glmer函数时获得固定和随机效应估计的协方差矩阵lme4。
我试过vcov(.., full = TRUE)没有成功。
有没有一种函数或方法来计算这个方差-协方差矩阵?
编辑
我需要的协方差矩阵是。对于回归参数估计,。归一化的部分似然信息矩阵的极限值。简而言之,我需要在处评估的观察到的逆信息矩阵。
在GLMMadaptive包中,该vcov()方法返回固定效应系数的最大似然估计的协方差矩阵和随机效应的方差-协方差矩阵的参数(在 log-Cholesky 因子尺度中)。
例如,请查看此处。
除非您特意不计算 Hessian,否则它会隐藏在输出模型结构中。您可以查看lme4:::vcov.merMod这些计算的来源(有什么更复杂的,因为它处理一堆边缘情况;它还仅提取与固定效应相关的协方差矩阵的一部分......)
例子:
library(lme4)
object <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp,
family = binomial)
这会提取 Hessian 矩阵,将其反转,然后将其加倍(因为 Hessian 矩阵是在(-2 对数似然)尺度上计算的。这h+t(h)是在加倍时提高对称性的巧妙方法(如果我没记错的话......)
h <- object@optinfo$derivs$Hessian
h <- solve(h)
v <- forceSymmetric(h + t(h))
检查固定效应部分是否一致(随机效应参数在前):
all.equal(unname(as.matrix(vcov(object))),
unname(as.matrix(v)[-1,-1])) ## TRUE
警告:随机效应在 Cholesky 尺度上进行参数化(即,参数是随机效应协方差矩阵的 Cholesky 因子的下三角形,按列优先顺序排列)......如果您在方差协方差参数化中需要这个,或者在标准偏差相关参数化中,这将需要更多的工作。(如果你只有一个标量随机效应,那么参数就是标准差。)