如何检查具有负二项分布的 GAM 中的过度分散?
How to check for overdispersion in a GAM with negative binomial distribution?
我使用 mgcv
包中的 gam
在负二项式族中拟合广义加性模型。我有一个包含因变量 y
、自变量 x
、因子 fac
和随机变量 ran
的数据框。我适合以下型号
gam1 <- gam(y ~ fac + s(x) + s(ran, bs = 're'), data = dt, family = "nb"
我在负二项式回归书上看到模型仍然有可能过度分散。我在 glm
中找到了检查过度分散的代码,但我没能在 gam
中找到它。我也遇到过只检查 QQ 图和标准化残差与预测残差的建议,但如果数据仍然过度分散,我无法从我的图中决定。因此,我正在寻找可以解决我的问题的方程式。
检查模型与观测数据的比较情况(并因此检查数据是否相对于模型隐含的条件分布过度分散)的好方法是通过根图。
我有一个 blog post 展示了如何使用 countreg 包对 glm()
模型执行此操作,但这也适用于 GAM。
应用于 GAM 版本模型的 post 的显着部分是:
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()
将两个根图排列在同一个图上绘制之后)
我们可以看到,对于这些数据,负二项式模型在这里比泊松 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;我们在泊松 GAM 中看到大量过度分散,在负二项式 GAM 中看到一些分散不足。
我使用 mgcv
包中的 gam
在负二项式族中拟合广义加性模型。我有一个包含因变量 y
、自变量 x
、因子 fac
和随机变量 ran
的数据框。我适合以下型号
gam1 <- gam(y ~ fac + s(x) + s(ran, bs = 're'), data = dt, family = "nb"
我在负二项式回归书上看到模型仍然有可能过度分散。我在 glm
中找到了检查过度分散的代码,但我没能在 gam
中找到它。我也遇到过只检查 QQ 图和标准化残差与预测残差的建议,但如果数据仍然过度分散,我无法从我的图中决定。因此,我正在寻找可以解决我的问题的方程式。
检查模型与观测数据的比较情况(并因此检查数据是否相对于模型隐含的条件分布过度分散)的好方法是通过根图。
我有一个 blog post 展示了如何使用 countreg 包对 glm()
模型执行此操作,但这也适用于 GAM。
应用于 GAM 版本模型的 post 的显着部分是:
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()
将两个根图排列在同一个图上绘制之后)
我们可以看到,对于这些数据,负二项式模型在这里比泊松 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;我们在泊松 GAM 中看到大量过度分散,在负二项式 GAM 中看到一些分散不足。