glmer过度离散的模型检验和测试

Question

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

我的响应变量(n)是每株植物每个柱头的pollengrains(log10 + 1)的数量,平均因为我每株植物收集3个柱头.
数据不适合泊松分布,因为(i)不是整数,(ii)方差远高于平均值(比率= 911.0756).所以,我适合作为negative.binomial.

模型选择后,我有:

m4a<-glmer(n ~ habitat*stigmatype + 

(1 | stigmaspecies /代码),家庭= 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号图形的趋势.这是我的最后一个问题:

3. 我可以用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)

希望这有所帮助....

干杯

Answer 2

对于那些可能通过谷歌搜索发现这一点的人来说，这只是 Q1 的补充：该blmco dispersion_glmer功能似乎已经过时。为此目的最好使用@Ben_Bolker 的函数：

overdisp_fun <- function(model) {
rdf <- df.residual(model)
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

来源：https : //bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion。

突出显示的概念：

请注意这里提到的通常和额外的警告：这是过度分散参数的近似估计。

附注。为什么过时了？

该lme4软件包包括residuals这些天的函数，据说 Pearson 残差对于这种类型的计算比偏差残差更稳健。将blmeco::dispersion_glmer偏差残差与u立方相加，除以残差自由度并取值（函数）的平方根：

 dispersion_glmer <- function (modelglmer) 
{
    n <- length(resid(modelglmer))
    return(sqrt(sum(c(resid(modelglmer), modelglmer@u)^2)/n))
}

该blmeco解决方案提供了比 Bolker 函数高得多的偏差/df 比率。由于 Ben 是该lme4软件包的作者之一，尽管我没有资格将统计原因合理化，但我更相信他的解决方案。

x <- InsectSprays
x$id <- rownames(x)

mod <- lme4::glmer(count ~ spray + (1|id), data = x, family = poisson)

blmeco::dispersion_glmer(mod)
# [1] 1.012649
overdisp_fun(mod)
#     chisq      ratio        rdf          p 
# 55.7160734  0.8571704 65.0000000  0.7873823