使用 glmmadmb 解释 mcmc 输出

问题描述

我正在尝试使用 glmmadmb 评估负二项式混合模型的输出.为了总结输出，我将汇总函数与 mcmc 选项的输出进行比较.我已经运行了这个模型:

I'm trying to evaluate the output from a negative binomial mixed model using glmmadmb. To summarize the output I'm comparing the summary function with output forom the mcmc option. I have run this model:

         pre1 <- glmmadmb(walleye~(1|year.center) + (1|Site) ,data=pre,
             family="nbinom2",link="log",
             mcmc=TRUE,mcmc.opts=mcmcControl(mcmc=1000))

我有两个随机截取:年份和站点.年份有 33 个级别，站点有 15 个.

I have two random intercepts: year and site. Year has 33 levels and site has 15.

来自 summary(pre1) 的站点和年份的随机效应参数估计似乎与来自 mcmc 输出的后验分布不一致.我使用 50% 置信区间作为估计值，该估计值应该与汇总函数的参数估计值一致.那不正确吗?有没有办法使用汇总函数来获得随机效应参数周围的误差来衡量这是否是方差问题?我尝试将 postvar=T 与 ranef 一起使用，但这不起作用.另外，有没有办法用信息丰富的行名称格式化 mcmc 输出以确保我使用正确的估计?

The random effect parameter estimate for site and year from summary(pre1) do not seem to agree with the posterior distribution from the mcmc output. I am using the 50% confidence interval as the estimate that should coincide with the parameter estimate from the summary function. Is that incorrect? Is there a way to obtain an error around the random effect parameter using the summary function to gauge whether this is variance issue? I tried using postvar=T with ranef but that did not work. Also, Is there a way to format the mcmc output with informative row names to ensure I'm using the proper estimates?

glmmabmb 的摘要输出:总结(pre1)

summary output from glmmabmb:summary(pre1)

Call:
glmmadmb(formula = walleye ~ (1 | year.center) + (1 | Site),
data = pre, family = "nbinom2", link = "log", mcmc = TRUE,
mcmc.opts = mcmcControl(mcmc = 1000))

AIC: 4199.8

Coefficients:
        Estimate Std. Error z value Pr(>|z|)
 (Intercept)    3.226      0.154      21   <2e-16 ***

 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of observations: total=495, year.center=33, Site=15
Random effect variance(s):
Group=year.center
             Variance StdDev
 (Intercept)   0.1085 0.3295
 Group=Site
             Variance StdDev
 (Intercept)   0.2891 0.5377

 Negative binomial dispersion parameter: 2.0553 (std. err.: 0.14419)

 Log-likelihood: -2095.88

mcmc 输出:m <- as.mcmc(pre1$mcmc)CI

mcmc output:m <- as.mcmc(pre1$mcmc)CI <- t(apply(m,2,quantile,c(0.025,0.5,0.975)))

                    2.5%          50%         97.5%
(Intercept)  2.911667943  3.211775843  3.5537371345
tmpL.1       0.226614903  0.342206509  0.4600328729
tmpL.2       0.395353518  0.554211483  0.8619127547
alpha        1.789687691  2.050871824  2.3175742167
u.01         0.676758365  0.896844797  1.0726750539
u.02         0.424938481  0.588191585  0.7364795440

这些估计继续到 u.48 以包括年份和地点特定系数.

these estimates continue to u.48 to include year and site specific coefficients.

预先感谢您对此问题的任何想法.蒂芙尼

Thank you in advance for any thoughts on this issue.Tiffany

推荐答案

这不是 50% 的置信区间，而是 50% 的分位数(即中位数).年间和站点间标准差的拉普拉斯近似值分别为{0.3295,0.5377}，这似乎非常接近MCMC中值估计{0.342206509,0.554211483} ... 如下所述，MCMC tmpL 参数是随机效应标准偏差，而不是方差——这可能是导致您困惑?

It's not the 50% confidence interval, it's the 50% quantile (i.e. the median). The point estimates from the Laplace approximation of the among-year and among-site standard deviations respectively are {0.3295,0.5377}, which seem quite close to the MCMC median estimates {0.342206509,0.554211483} ... as discussed below, the MCMC tmpL parameters are the random-effects standard deviations, not the variances -- this might be the main cause of your confusion?

有没有办法使用汇总函数来获得随机效应参数周围的误差来判断这是否是方差问题?我尝试将 postvar=T 与 ranef 一起使用，但这不起作用.

lme4 包(不是 glmmadmb 包)允许估计条件模式的方差(即与特定级别)通过 ranef(...,condVar=TRUE)(postVar=TRUE 现在已弃用).关于条件模式不确定性的等效信息可通过 ranef(model,sd=TRUE) 获得(参见 ?ranef.glmmadmb).

The lme4 package (not the glmmadmb package) allows estimates of the variances of the conditional modes (i.e. the random effects associated with particular levels) via ranef(...,condVar=TRUE) (postVar=TRUE is now deprecated). The equivalent information on the uncertainty of the conditional modes is available via ranef(model,sd=TRUE) (see ?ranef.glmmadmb).

但是，我认为您可能正在寻找 $S(方差-协方差矩阵)和 $sd_S(方差-协方差估计的 Wald 标准误差)(虽然如上所述，我认为这真的没有问题).

However, I think you might be looking for the $S (variance-covariance matrices) and $sd_S (Wald standard errors of the variance-covariance estimates) instead (although as stated above, I don't think there's really a problem).

另外，有没有办法用信息丰富的行名称格式化 mcmc 输出以确保我使用正确的估计?

见第vignette("glmmADMB",package="glmmADMB") 的 15 个:

glmmADMB 中的 MCMC 输出未完全翻译.它按顺序包括:

• pz 零通胀参数(原始)
• fixed Effect parameters 以与coef() 的结果相同的方式命名或fixef().
• tmpL 方差(标准偏差量表)
• tmpL1 方差-协方差矩阵的相关/off?-对角元素(相关矩阵的 Cholesky 因子的 off?-对角元素).(如果您需要将这些转换为相关性，则需要在对角线上构建相关矩阵，并计算叉积 CC^T(请参阅 tcrossprod)；如果这没有意义给你，联系维护者)
• alpha 过度分散/尺度参数
• u 随机效应(未缩放:这些可以使用来自 VarCorr() 的估计随机效应标准偏差进行缩放)
• pz zero-inflation parameter (raw)
• fixed effect parameters Named in the same way as the results of coef() orfixef().
• tmpL variances (standard-deviation scale)
• tmpL1 correlation/off-diagonal elements of variance-covariance matrices (off-diagonal elements of the Cholesky factor of the correlation matrix'). (If you need to transform these to correlations, you will need to construct the relevant matrices with 1 on the diagonal and compute the cross-product, CC^T (see tcrossprod); if this makes no sense to you, contact the maintainers)
• alpha overdispersion/scale parameter
• u random effects (unscaled: these can be scaled using the estimated random-effects standard deviations from VarCorr())

这篇关于使用 glmmadmb 解释 mcmc 输出的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！