r - glmer 过度离散的模型检验和检验

Question

我正在测试不同生境和柱头类型的植物柱头上花粉粒数量的差异。
我的样本设计包括两个栖息地，每个栖息地有 10 个站点。
在每个地点，我有多达 3 种柱头类型(湿、干和半干)，对于每种柱头类型，我有不同数量的植物物种，每种植物物种的个体数量不同(代码)。
所以，我最终得到了如下的嵌套设计:habitat/site/stigmatype/stigmaspecies/code由于这是一项描述性研究，柱头类型、柱头种类和代码因地点而异。

我的响应变量 (n) 是每株植物每个柱头的花粉粒数 (log10+1)，平均值，因为我从每株植物收集了 3 个柱头。
数据不符合泊松分布，因为 (i) 不是整数，并且 (ii) 方差远高于平均值(比率 = 911.0756)。所以，我拟合为负二项式。

模型选择后，我有:

m4a <- glmer(n ~ habitat*stigmatype + (1|stigmaspecies/code),
             family=negative.binomial(2))

    > summary(m4a)
    Generalized linear mixed model fit by maximum likelihood ['glmerMod']
    Family: Negative Binomial(2) ( log )
    Formula: n ~ habitat * stigmatype + (1 | stigmaspecies/code) 
    
          AIC       BIC    logLik  deviance
     993.9713 1030.6079 -487.9856  975.9713 

    Random effects:
    Groups             Name        Variance  Std.Dev. 
    code:stigmaspecies (Intercept) 1.034e-12 1.017e-06
    stigmaspecies      (Intercept) 4.144e-02 2.036e-01
    Residual                       2.515e-01 5.015e-01
    Number of obs: 433, groups: code:stigmaspecies, 433; stigmaspecies, 41

    Fixed effects:  
                                         Estimate Std. Error t value Pr(>|z|)    
    (Intercept)                          -0.31641    0.08896  -3.557 0.000375 ***
    habitatnon-invaded                   -0.67714    0.10060  -6.731 1.68e-11 ***
    stigmatypesemidry                    -0.24193    0.15975  -1.514 0.129905    
    stigmatypewet                        -0.07195    0.18665  -0.385 0.699885    
    habitatnon-invaded:stigmatypesemidry  0.60479    0.22310   2.711 0.006712 ** 
    habitatnon-invaded:stigmatypewet      0.16653    0.34119   0.488 0.625491    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

    Correlation of Fixed Effects:
                     (Intr) hbttn- stgmtyps stgmtypw hbttnn-nvdd:stgmtyps
hbttnn-nvdd          -0.335                                            
stgmtypsmdr          -0.557  0.186                           
stigmatypwt          -0.477  0.160  0.265                                
hbttnn-nvdd:stgmtyps  0.151 -0.451 -0.458   -0.072                       
hbttnn-nvdd:stgmtypw  0.099 -0.295 -0.055   -0.403    0.133  

两个问题:

1. 如何检查此输出是否过度分散？
2. 在这里进行模型验证的最佳方法是什么？

我一直在使用:

qqnorm(resid(m4a))
hist(resid(m4a))
plot(fitted(m4a),resid(m4a))

虽然qqnorm() 和hist() 看起来没问题，但在第3 个图上有异方差的倾向。这是我的最后一个问题:

1. 我可以在 glmer 中使用此图进行模型验证吗？或者有更好的方法吗？如果不是，我应该为第三张图担心多少？

Answer 1

在 glmer 中检查过度分散的一种简单方法是:

> library("blmeco") 
> dispersion_glmer(your_model) #it shouldn't be over
> 1.4

为了解决过度分散，我通常会添加一个观察级别的随机因子

对于模型验证，我通常从这些图开始……但随后取决于您的特定模型……

par(mfrow=c(2,2))
qqnorm(resid(your_model), main="normal qq-plot, residuals")
qqline(resid(your_model))

qqnorm(ranef(your_model)$id[,1])
qqline(ranef(your_model)$id[,1])

plot(fitted(your_model), resid(your_model)) #residuals vs fitted
abline(h=0)


dat_kackle$fitted <- fitted(your_model)    #fitted vs observed
plot(your_data$fitted, jitter(your_data$total,0.1))
abline(0,1)

希望对您有所帮助....

干杯