如何检查具有负二项分布的 GAM 中的过度分散？答案

【问题标题】：How to check for overdispersion in a GAM with negative binomial distribution?如何检查具有负二项分布的 GAM 中的过度分散？
【发布时间】：2019-12-15 09:20:27
【问题描述】：

我使用 mgcv 包中的 gam 在负二项式系列中拟合了广义加法模型。我有一个数据框，其中包含我的因变量y、一个自变量x、一个因子fac 和一个随机变量ran。我适合以下模型

gam1 <- gam(y ~ fac + s(x) + s(ran, bs = 're'), data = dt, family = "nb"

我在 Negative Binomial Regression 一书中读到，模型仍然可能过度分散。我在glm 中找到了检查过度分散的代码，但我找不到gam 的代码。我也遇到过仅检查 QQ 图和标准化残差与预测残差的建议，但我无法从我的图中确定数据是否仍然过度分散。因此，我正在寻找一个可以解决我的问题的方程式。

【问题讨论】：

标签： r gam mgcv

【解决方案1】：

检查模型与观察数据的比较情况（并因此检查数据相对于模型隐含的条件分布是否过度分散）的好方法是通过根图。

我有一个 blog post 显示如何使用 countreg 包为 glm() 模型执行此操作，但这也适用于 GAM。

适用于模型的 GAM 版本的帖子的主要部分是：

library("coenocliner")
library('mgcv')

## parameters for simulating
set.seed(1)
locs <- runif(100, min = 1, max = 10)     # environmental locations
A0 <- 90                                  # maximal abundance
mu <- 3                                   # position on gradient of optima
alpha <- 1.5                              # parameter of beta response
gamma <- 4                                # parameter of beta response
r <- 6                                    # range on gradient species is present
pars <- list(m = mu, r = r, alpha = alpha, gamma = gamma, A0 = A0)
nb.alpha <- 1.5                           # overdispersion parameter 1/theta
zprobs <- 0.3                             # prob(y == 0) in binomial model

## simulate some negative binomial data from this response model
nb   <- coenocline(locs, responseModel = "beta", params = pars,
                   countModel = "negbin",
                   countParams = list(alpha = nb.alpha))
df <- setNames(cbind.data.frame(locs, nb), c("x", "yNegBin"))

好的，所以我们有一个从负二项式抽样分布中抽取的数据样本，现在我们将用两个模型拟合这些数据：

泊松 GAM

m_pois <- gam(yNegBin ~ s(x), data = df, family = poisson())

负二项式 GAM

m_nb   <- gam(yNegBin ~ s(x), data = df, family = nb())

countreg 包尚未在 CRAN 上，但可以从 R-Forge 安装：

install.packages("countreg", repos="http://R-Forge.R-project.org")

然后加载包并绘制根图：

library("countreg")
library("ggplot2")

root_pois <- rootogram(m_pois, style = "hanging", plot = FALSE)
root_nb   <- rootogram(m_nb, style = "hanging", plot = FALSE)

现在绘制每个模型的根图：

autoplot(root_pois)
autoplot(root_nb)

这就是我们得到的（在使用cowplot::plot_grid() 将两个根图排列在同一个图上之后）

我们可以看到，对于这些数据，负二项式模型在这里比 Poisson GAM 做得更好——在观察到的计数范围内，条形的底部更接近于零。

countreg 包详细介绍了如何在零线周围添加不确定带作为拟合优度测试的一种形式。

您还可以使用每个模型的 Pearson 残差计算分散参数的 Pearson 估计值：

r$>  sum(residuals(m_pois, type = "pearson")^2) / df.residual(m_pois)
[1] 28.61546
r$>  sum(residuals(m_nb, type = "pearson")^2) / df.residual(m_nb)
[1] 0.5918471

在这两种情况下，这些都应该是 1；我们看到 Poisson GAM 中存在显着的过度分散，而负二项式 GAM 中存在一些欠分散。

【讨论】：