机器算法验证 - 难以找到适合具有混合效应的计数数据的良好模型 - ZINB 或其他？ - 吾爱随笔录

难以找到适合具有混合效应的计数数据的良好模型 - ZINB 或其他？

机器算法验证计数数据负二项分布混合模式零通胀 lme4-nlme

2022-03-19 10:45:05

我有一个非常小的关于孤独蜜蜂丰度的数据集，我无法对其进行分析。它是计数数据，几乎所有计数都在一种处理中，而另一种处理中的大多数零。还有几个非常高的值（六个站点中的两个站点各一个），因此计数的分布具有极长的尾巴。我在 R 工作。我使用了两个不同的包：lme4 和 glmmADMB。

泊松混合模型不适合：未拟合随机效应时模型过度分散（glm 模型），而拟合随机效应时模型分散不足（glmer 模型）。我不明白为什么会这样。实验设计需要嵌套的随机效应，所以我需要包括它们。泊松对数正态误差分布没有改善拟合。我使用 glmer.nb 尝试了负二项式误差分布，但无法拟合 - 达到了迭代限制，即使使用 glmerControl(tolPwrss=1e-3) 更改了容差。

因为很多零是由于我根本没有看到蜜蜂（它们通常是微小的黑色东西），所以我接下来尝试了一个零膨胀模型。ZIP 不太合适。ZINB 是迄今为止最合适的模型，但我对模型的拟合仍然不太满意。我不知道接下来要尝试什么。我确实尝试了一个跨栏模型，但无法将截断分布拟合到非零结果——我认为因为很多零都在控制处理中（错误消息是“model.frame.default 中的错误（公式 = s.bee ~ tmt + lu + : 可变长度不同（为“治疗”找到）”）。

此外，我认为我所包含的交互对我的数据做了一些奇怪的事情，因为系数小得不切实际——尽管当我使用包 bbmle 中的 AICctab 比较模型时，包含交互的模型是最好的。

我包含了一些 R 脚本，它们几乎可以重现我的数据集。变量如下：

d=儒略日期，df=儒略日期（作为因子），d.sq=df 平方（蜜蜂数量在整个夏天增加然后下降），st=site，s.bee=蜜蜂计数，tmt=治疗，lu=土地利用类型，hab=周围景观中半自然栖息地的百分比，ba=田野边界区域。

任何关于如何获得良好模型拟合的建议（替代误差分布、不同类型的模型等）都将非常感激！

谢谢你。

d <- c(80,  80,  121, 121, 180, 180, 86,  86,  116, 116, 144, 144, 74,  74, 143, 143, 163, 163, 71, 71,106, 106, 135, 135, 162, 162, 185, 185, 83,  83,  111, 111, 133, 133, 175, 175, 85,  85,  112, 112,137, 137, 168, 168, 186, 186, 64,  64,  95,  95,  127, 127, 156, 156, 175, 175, 91,  91, 119, 119,120, 120, 148, 148, 56, 56)
df <- as.factor(d)
d.sq <- d^2
st <- factor(rep(c("A", "B", "C", "D", "E", "F"), c(6,12,18,10,14,6)))
s.bee <- c(1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,4,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,3,0,0,0,0,5,0,0,2,0,50,0,10,0,4,0,47,3)
tmt <- factor(c("AF","C","C","AF","AF","C","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","AF","C","AF","C","C","AF","C","AF","C","AF","AF","C","AF","C",
"AF","C","AF","C","AF","C"))
lu <- factor(rep(c("p","a","p","a","p"), c(6,12,28,14,6)))
hab <- rep(c(13,14,13,14,3,4,3,4,3,4,3,4,3,4,15,35,37,35,37,35,37,35,37,0,2,1,2,1,2,1), 
        c(1,2,2,1,1,1,1,2,2,1,1,1,1,1,18,1,1,1,2,2,1,1,1,14,1,1,1,1,1,1))
ba <-  c(480,6520,6520,480,480,6520,855,1603,855,1603,1603,855,855,12526,855,5100,855,5100,2670,7679,7679,2670,
2670,7679,2670,7679,7679,2670,2670,7679,2670,7679,2670,7679,2670,7679,1595,3000,1595,3000,3000,1595,1595,3000,1595
,3000,4860,5460,4860,5460,5460,4860,5460,4860,5460,4860,4840,5460,4840,5460,3000,1410,3000,1410,3000,1410)
data <- data.frame(st,df,d.sq,tmt,lu,hab,ba,s.bee)
with(data, table(s.bee, tmt) )

# below is a much abbreviated summary of attempted models:

library(MASS)
library(lme4)
library(glmmADMB)
library(coefplot2)

###
### POISSON MIXED MODEL

    m1 <- glmer(s.bee ~ tmt + lu + hab + (1|st/df), family=poisson)
    summary(m1)

    resdev<-sum(resid(m1)^2)
    mdf<-length(fixef(m1)) 
    rdf<-nrow(data)-mdf
    resdev/rdf
# 0.2439303
# underdispersed. ???



###
### NEGATIVE BINOMIAL MIXED MODEL

    m2 <- glmer.nb(s.bee ~ tmt + lu + hab + d.sq + (1|st/df))
# iteration limit reached. Can't make a model work.



###
### ZERO-INFLATED POISSON MIXED MODEL

    fit_zipoiss <- glmmadmb(s.bee~tmt + lu + hab + ba + d.sq +
                    tmt:lu +
                    (1|st/df), data=data,
                    zeroInflation=TRUE,
                    family="poisson")
# has to have lots of variables to fit
# anyway Poisson is not a good fit



###
### ZERO-INFLATED NEGATIVE BINOMIAL MIXED MODELS

## BEST FITTING MODEL SO FAR:

    fit_zinb <- glmmadmb(s.bee~tmt + lu + hab +
                    tmt:lu +
                    (1|st/df),data=data,
                    zeroInflation=TRUE,
                    family="nbinom")
    summary(fit_zinb)
# coefficients are tiny, something odd going on with the interaction term
# but this was best model in AICctab comparison

# model check plots
    qqnorm(resid(fit_zinb))
    qqline(resid(fit_zinb))

    coefplot2(fit_zinb)

    resid_zinb <- resid(fit_zinb , type = "pearson")
    hist(resid_zinb)

    fitted_zinb <- fitted (fit_zinb)
    plot(resid_zinb ~ fitted_zinb)


## MODEL WITHOUT INTERACTION TERM - the coefficients are more realistic:

    fit_zinb2 <- glmmadmb(s.bee~tmt + lu + hab +
                    (1|st/df),data=data,
                    zeroInflation=TRUE,
                    family="nbinom")

# model check plots
    qqnorm(resid(fit_zinb2))
    qqline(resid(fit_zinb2))

    coefplot2(fit_zinb2)

    resid_zinb2 <- resid(fit_zinb2 , type = "pearson")
    hist(resid_zinb2)

    fitted_zinb2 <- fitted (fit_zinb2)
    plot(resid_zinb2 ~ fitted_zinb2)



# ZINB models are best so far
# but I'm not happy with the model check plots

1个回答

这篇文章有四年了，但我想关注 fsociety 在评论中所说的话。GLMM 中残差的诊断并不简单，因为标准残差图可以显示非正态性、异方差性等，即使模型被正确指定。有一个 R 包，DHARMa特别适用于诊断这些类型的模型。

该软件包基于一种模拟方法，可从拟合的广义线性混合模型生成缩放残差，并生成不同的易于解释的诊断图。这是一个小示例，其中包含来自原始帖子和第一个拟合模型 (m1) 的数据：

library(DHARMa)
sim_residuals <- simulateResiduals(m1, 1000)
plotSimulatedResiduals(sim_residuals)

左边的图显示了缩放残差的 QQ 图，以检测与预期分布的偏差，右边的图表示残差与预测值，同时执行分位数回归以检测与均匀性的偏差（红线应为水平线，位于 0.25 , 0.50 和 0.75)。

此外，该软件包还具有测试过度/不足分散和零通货膨胀等特定功能：

testOverdispersionParametric(m1)

Chisq test for overdispersion in GLMMs

data:  poisson
dispersion = 0.18926, pearSS = 11.35600, rdf = 60.00000, p-value = 1
alternative hypothesis: true dispersion greater 1

testZeroInflation(sim_residuals)

DHARMa zero-inflation test via comparison to expected zeros with 
simulation under H0 = fitted model


data:  sim_residuals
ratioObsExp = 0.98894, p-value = 0.502
alternative hypothesis: more

其它你可能感兴趣的问题

上一篇Bootstrap：估计值超出置信区间下一篇人口稳定指数 - 除以零