不明白为什么 glmm 随机效应方差为零。看过类似的问题还是不明白

机器算法验证 r 方差 lme4-nlme 咕噜咕噜 生态
2022-04-01 18:42:49

我研究一种在殖民地筑巢的鸟类。我正在尝试为巢址选择研究执行 GLMM 的 AICc 评估。我在巢穴和配对的随机地点收集数据。我想评估栖息地特征（固定效应），同时考虑菌落地点（随机效应）来预测巢址（响应变量）。巢址是二元的（巢址/成对的随机位置）。所有数据都经过标准化，以便具有正态分布 ((value-mean)/sd)。

我有 5 个菌落，其中有 4-9 个巢穴，我在分析中使用了 38 个巢穴和 38 个随机对。当我使用 lme4（glmer 函数）在 r 中运行 GLMM 时，我得到一个结果，说我对 60 个模型中的每一个模型都有 0 方差。这是有问题的，因为这不会发生在我移除随机效应的情况下，但我想保留它以解释菌落中的空间自相关。

我已经阅读了这个网站上的一些文章和其他关于低组数如何抑制我的方差的文章，但是，当我运行模型时略有不同（结构不正确，但基本上将所有配对的站点集中在一起形成一个“群体“有 38 个观察值并将其与 4-9 个巢的其他小群体进行比较），它确实给了我一个差异。所以我不认为群体数量本身是罪魁祸首。也许是群体中的小样本和小群体的结合？尽管如此..任何帮助都会很棒。如果我无法弄清楚，我想我不得不求助于没有菌落的 glm 固定效应模型作为随机效应。这是不正确的，但这是一个开始？如果你不能让它发挥作用，你会怎么做。谢谢你。

我想强调一点：我不是统计学家或程序员。我抓鸟。我上过统计课，我对一些事情有一定的把握，但让我们面对现实吧，有时它转瞬即逝。如果答案不能太深奥，请更多地关注实用性，例如“你应该这样做。你应该这样做”，我将非常感激。我正在尝试尽快完成这项分析。再次感谢你。

我咨询过的一些链接：

广义线性混合模型中的随机效应等于 0

为什么我的混合模型中随机效应的方差为零，尽管数据有一些变化？

https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#singular-models-random-effect-variances-estimated-as-zero-or-correlations-estimated-as---1

http://rpubs.com/bbolker/4187

下面是我的代码示例。

helpmeobiwan <-list(NestPlot = c(1, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0),NumDeadJun = c( 0.1409216, -0.1932639,-0.5274494,-0.5274494, 0.1409216, -0.5274494, -0.5274494 , 0.4751071, -0.5274494 , 2.1460347 ,-0.5274494, -0.1932639, 0.8092926, -0.5274494, -0.5274494 ,-0.5274494 ,-0.1932639, 0.1409216, -0.5274494, -0.5274494 ,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216, -0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.1932639, -0.5274494,  0.4751071 , 0.1409216 ,-0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494 ,-0.5274494 ,-0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494, -0.5274494,  0.1409216, -0.5274494, -0.5274494  ,3.1485912 , 2.4802202,  1.4776637, -0.5274494 , 2.8144057, -0.5274494, -0.5274494,  1.1434781,  3.8169623,  3.8169623 ,-0.1932639, -0.5274494  ,1.4776637 , 1.8118492, -0.5274494),RandomPair = c(  "Madera2" ,  "Starfire1", "Madera2" ,  "Madera3" ,  "Starfire1" ,"Starfire1", "Starfire2", "Madera1" , "Madera3"   ,"Starfire2" ,"Starfire2", "Madera1",   "Madera2",  "Starfire1", "Starfire1" ,"Starfire1", "Madera1",   "Madera2" ,  "Starfire1", "Starfire1", "Starfire1", "Madera1" ,  "Starfire1", "Starfire1", "Madera1",   "Madera1" , "Starfire1", "Madera2" ,  "Madera1",   "Madera2" ,  "Madera1" ,  "Madera1"  , "Starfire1" ,"Starfire1", "Starfire1" ,"Starfire1" ,"Madera2"  , "Madera2",   "Starfire2" ,"Starfire2", "Starfire2" ,"Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Starfire2", "Starfire2", "Starfire2", "Starfire2" ,"Starfire2", "Madera3",  "Madera3" ,  "Starfire2", "Madera3" ,  "Madera1"  , "Starfire2" ,"Starfire1", "Madera2" ,  "Madera3" ,  "Madera3"  , "Madera2"  , "Madera3"   ,"Starfire2", "Madera3",   "Starfire1", "Madera3"  , "Starfire2", "Starfire1", "Madera3",   "Starfire1", "Starfire2" ,"Madera1" ,  "Starfire2", "Starfire2", "Madera1"  ))




m1 <- glmer(NestPair ~ NumDeadJun+ (1|RandomPair), family=binomial, data=helpmeobiwan)
2个回答

函数glmer()默认使用拉普拉斯近似，这对于二分数据不是最优的。一个更好的选择是自适应高斯正交。您可以通过将参数设置nAGQ为glmer()更高的数字（例如，11 或 15）或使用GLMMadaptive包来使用此方法。在您的示例中，它给出：

library("GLMMadaptive")
helpmeobiwan <- list(NestPlot = c(1, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0),NumDeadJun = c( 0.1409216, -0.1932639,-0.5274494,-0.5274494, 0.1409216, -0.5274494, -0.5274494 , 0.4751071, -0.5274494 , 2.1460347 ,-0.5274494, -0.1932639, 0.8092926, -0.5274494, -0.5274494 ,-0.5274494 ,-0.1932639, 0.1409216, -0.5274494, -0.5274494 ,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216, -0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.1932639, -0.5274494,  0.4751071 , 0.1409216 ,-0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494 ,-0.5274494 ,-0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494, -0.5274494,  0.1409216, -0.5274494, -0.5274494  ,3.1485912 , 2.4802202,  1.4776637, -0.5274494 , 2.8144057, -0.5274494, -0.5274494,  1.1434781,  3.8169623,  3.8169623 ,-0.1932639, -0.5274494  ,1.4776637 , 1.8118492, -0.5274494),RandomPair = c(  "Madera2" ,  "Starfire1", "Madera2" ,  "Madera3" ,  "Starfire1" ,"Starfire1", "Starfire2", "Madera1" , "Madera3"   ,"Starfire2" ,"Starfire2", "Madera1",   "Madera2",  "Starfire1", "Starfire1" ,"Starfire1", "Madera1",   "Madera2" ,  "Starfire1", "Starfire1", "Starfire1", "Madera1" ,  "Starfire1", "Starfire1", "Madera1",   "Madera1" , "Starfire1", "Madera2" ,  "Madera1",   "Madera2" ,  "Madera1" ,  "Madera1"  , "Starfire1" ,"Starfire1", "Starfire1" ,"Starfire1" ,"Madera2"  , "Madera2",   "Starfire2" ,"Starfire2", "Starfire2" ,"Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Starfire2", "Starfire2", "Starfire2", "Starfire2" ,"Starfire2", "Madera3",  "Madera3" ,  "Starfire2", "Madera3" ,  "Madera1"  , "Starfire2" ,"Starfire1", "Madera2" ,  "Madera3" ,  "Madera3"  , "Madera2"  , "Madera3"   ,"Starfire2", "Madera3",   "Starfire1", "Madera3"  , "Starfire2", "Starfire1", "Madera3",   "Starfire1", "Starfire2" ,"Madera1" ,  "Starfire2", "Starfire2", "Madera1"  ))
helpmeobiwan <- as.data.frame(helpmeobiwan)

fm <- mixed_model(NestPlot ~ NumDeadJun, random = ~ 1 | RandomPair, 
                  family = binomial(), data = helpmeobiwan)

summary(fm)
#> 
#> Call:
#> mixed_model(fixed = NestPlot ~ NumDeadJun, random = ~1 | RandomPair, 
#>     data = helpmeobiwan, family = binomial())
#> 
#> Data Descriptives:
#> Number of Observations: 76
#> Number of Groups: 5 
#> 
#> Model:
#>  family: binomial
#>  link: logit 
#> 
#> Fit statistics:
#>   log.Lik      AIC      BIC
#>  -46.2248 98.44959 97.27791
#> 
#> Random effects covariance matrix:
#>                StdDev
#> (Intercept) 0.0477673
#> 
#> Fixed effects:
#>             Estimate Std.Err z-value  p-value
#> (Intercept)  -0.1568  0.2829 -0.5544 0.579304
#> NumDeadJun   -1.2274  0.4917 -2.4961 0.012558
#> 
#> Integration:
#> method: adaptive Gauss-Hermite quadrature rule
#> quadrature points: 11
#> 
#> Optimization:
#> method: hybrid EM and quasi-Newton
#> converged: TRUE

我尝试使用 glmmADMB 包，它是用于线性混合建模的 lme4 的替代方案。您可以使用以下代码安装此软件包：

install.packages("R2admb")
install.packages("glmmADMB", 
repos=c("http://glmmadmb.r-forge.r-project.org/repos",
        getOption("repos")),
type="source")

然后你去：

library(glmmADMB)
helpmeobiwan <-list(NestPlot = c(1, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0),NumDeadJun = c( 0.1409216, -0.1932639,-0.5274494,-0.5274494, 0.1409216, -0.5274494, -0.5274494 , 0.4751071, -0.5274494 , 2.1460347 ,-0.5274494, -0.1932639, 0.8092926, -0.5274494, -0.5274494 ,-0.5274494 ,-0.1932639, 0.1409216, -0.5274494, -0.5274494 ,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216, -0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.1932639, -0.5274494,  0.4751071 , 0.1409216 ,-0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494 ,-0.5274494 ,-0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494, -0.5274494,  0.1409216, -0.5274494, -0.5274494  ,3.1485912 , 2.4802202,  1.4776637, -0.5274494 , 2.8144057, -0.5274494, -0.5274494,  1.1434781,  3.8169623,  3.8169623 ,-0.1932639, -0.5274494  ,1.4776637 , 1.8118492, -0.5274494),RandomPair = c(  "Madera2" ,  "Starfire1", "Madera2" ,  "Madera3" ,  "Starfire1" ,"Starfire1", "Starfire2", "Madera1" , "Madera3"   ,"Starfire2" ,"Starfire2", "Madera1",   "Madera2",  "Starfire1", "Starfire1" ,"Starfire1", "Madera1",   "Madera2" ,  "Starfire1", "Starfire1", "Starfire1", "Madera1" ,  "Starfire1", "Starfire1", "Madera1",   "Madera1" , "Starfire1", "Madera2" ,  "Madera1",   "Madera2" ,  "Madera1" ,  "Madera1"  , "Starfire1" ,"Starfire1", "Starfire1" ,"Starfire1" ,"Madera2"  , "Madera2",   "Starfire2" ,"Starfire2", "Starfire2" ,"Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Starfire2", "Starfire2", "Starfire2", "Starfire2" ,"Starfire2", "Madera3",  "Madera3" ,  "Starfire2", "Madera3" ,  "Madera1"  , "Starfire2" ,"Starfire1", "Madera2" ,  "Madera3" ,  "Madera3"  , "Madera2"  , "Madera3"   ,"Starfire2", "Madera3",   "Starfire1", "Madera3"  , "Starfire2", "Starfire1", "Madera3",   "Starfire1", "Starfire2" ,"Madera1" ,  "Starfire2", "Starfire2", "Madera1"  ))
dontworryLeia <- helpmeobiwan
dontworryLeia$RandomPair <- as.factor(dontworryLeia$RandomPair)
attach(dontworryLeia)
mod <- glmmadmb(NestPlot ~ NumDeadJun + (1|RandomPair), family='binomial', data=dontworryLeia)
mod
summary(mod)
drop1(mod)

First RandomPair不被视为一个因素，这解释了这种转变。我猜你在你的 m1 模型中指的是NestPlot而不是NestPair。无论如何，这应该工作！

其它你可能感兴趣的问题
上一篇方差XX和方差日志( X)log⁡(X). 如何将它们联系起来？ 下一篇使用无法区分的主题作为预测变量/随机效应