添加观察水平随机项会混淆残差与拟合图。为什么?

机器算法验证 咕噜咕噜 过度分散
2022-03-26 09:51:15

我为比例数据(响应变量)运行混合效应广义模型。我使用了二项式族和 logit 链接函数。我遭受了过度分散的困扰,因此我添加了一个观察级别的随机效应来吸收部分无法解释的可变性。

问题是:现在我的残差图(残差与拟合)有一个清晰的模式。如果我删除那个新的“观察水平随机项”,情节看起来又好了。

为什么?这是一个常见问题还是取决于我的特定数据?有什么办法可以解决吗?

编辑:我附上了“前后”的情节。 之前,不是 OLRE 之后,使用 OLRE

1个回答

谢谢你更新你的帖子,查理。我使用了一些过度分散的泊松数据,以了解在 glmer 模型中添加观察水平效应对残差与拟合值图的影响。这是R代码:

# generate data like here: https://rpubs.com/INBOstats/OLRE

set.seed(324)
n.i <- 10
n.j <- 10
n.k <- 10
beta.0 <- 1
beta.1 <- 0.3
sigma.b <- 0.5
theta <- 5
dataset <- expand.grid(
X = seq_len(n.i),
b = seq_len(n.j),
Replicate = seq_len(n.k)
)
rf.b <- rnorm(n.j, mean = 0, sd = sigma.b)
dataset$eta <- beta.0 + beta.1 * dataset$X + rf.b[dataset$b]
dataset$mu <- exp(dataset$eta)
dataset$Y <- rnbinom(nrow(dataset), mu = dataset$mu, size = theta)
dataset$OLRE <- seq_len(nrow(dataset))


require(lme4)

m.1 <- glmer(Y ~ X + (1 | b), family=poisson(link="log"), data=dataset)

m.2 <- glmer(Y ~ X + (1 | b) + (1 | OLRE), family=poisson(link="log"),  
             data=dataset)

请注意,模型 m.2 包括一个观察水平随机效应来解释过度分散。

要诊断模型 m.1 中是否存在过度分散,我们可以使用以下命令:

# check for over-dispersion:  
# values greater than 1.4 indicate over-dispersion

require(blmeco)

dispersion_glmer(m.1)

dispersion_glmer 返回的值是 2.204209,它大于 1.4 的截止值,我们开始怀疑存在过度分散。

将dispersion_glmer 应用到模型m.2 时,我们得到的值为1.023656:

dispersion_glmer(m.2)

这是残差(皮尔逊或偏差)与拟合值图的 R 代码:

par(mfrow=c(1,2))
plot(residuals(m.1, type="pearson") ~ fitted(m.1), col="darkgrey")
     abline(h=0, col="red")
plot(residuals(m.2, type="pearson") ~ fitted(m.2), col="darkgrey")
     abline(h=0, col="red")

par(mfrow=c(1,2))
plot(residuals(m.1, type="deviance") ~ fitted(m.1), col="darkgrey")
abline(h=0, col="red")
plot(residuals(m.2, type="deviance") ~ fitted(m.2), col="darkgrey")
abline(h=0, col="red")

如您所见,与模型 m.1 的图相比,模型 m.2 的 Pearson 残差图(包括观察水平随机效应)看起来很可怕。 残差与拟合

我没有显示 m.2 的偏差残差图,因为它看起来大致相同(即可怕)。

以下是模型 m.1 和 m.2 的拟合值与观测响应值的关系图:

par(mfrow=c(1,2))
plot(fitted(m.1) ~ dataset$Y, col="darkgrey", 
     xlim=c(0, 250), ylim=c(0, 250), 
     xlab="Y (response)", ylab="Fitted Values")
abline(a=0, b=1, col="red")
plot(fitted(m.2) ~ dataset$Y, col="darkgrey", 
     xlim=c(0, 250), ylim=c(0, 250), 
     xlab="Y (response)", ylab="Fitted Values")
abline(a=0, b=1, col="red")

在此处输入图像描述

对于模型 m.2,拟合值与实际响应值的图看起来更好。

我们应该检查两个模型对应的摘要:

summary(m.1)

summary(m.2)

正如https://rpubs.com/INBOstats/OLRE中所述,固定效应系数之间的巨大差异,尤其是 b 的随机效应方差之间的差异表明可能有问题。(初始模型中存在的过度分散程度将推动这些差异的程度。)

让我们看一下使用 Dharma 包获得的两个模型的一些诊断图:

require(DHARMa)

fittedModel <- m.1
simulationOutput <- simulateResiduals(fittedModel = fittedModel)
plot(simulationOutput)

#----

fittedModel <- m.2
simulationOutput <- simulateResiduals(fittedModel = fittedModel)
plot(simulationOutput)

模型 m.1(尤其是左侧面板)的诊断图清楚地表明过度分散是一个问题。 m.1 的佛法

模型 m.2 的诊断图显示过度分散不再是问题。 在此处输入图像描述

有关这些类型的绘图的更多详细信息,请参阅https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html

最后,让我们对这两个模型进行后验预测检查(即,在真实响应值 Y 的直方图上绘制从每个模型构建的模拟数据集获得的拟合值),如http://www.flutterbys 中所述。 com.au/stats/ws/ws12.html

range(dataset$Y)  # Actual response values Y range from 0 to 247

set.seed(1234567)
glmer.sim1 <- simulate(m.1, nsim = 1000)
glmer.sim2 <- simulate(m.2, nsim = 1000)

out <- matrix(NA, ncol = 2, nrow = 251)
cnt <- 0:250
for (i in 1:length(cnt)) {
 for (j in 1:2) {
     eval(parse(text = paste("out[i,", j, "] <- 
          mean(sapply(glmer.sim", j,",\nFUN = function(x) {\nsum(x == cnt  
      [i]) }))", sep = "")))
 }
}


plot(table(dataset$Y), ylab = "Frequency", xlab = "Y", lwd = 2, 
     col="darkgrey")
lines(x = 0:250, y = out[, 1], lwd = 2, lty = 2, col = "red")    
lines(x = 0:250, y = out[, 2], lwd = 2, lty = 2, col = "blue")    

结果图表明,这两个模型在逼近 Y 的分布方面都做得很好。

在此处输入图像描述

当然,还有其他可以查看的预测检查,包括蜈蚣图,它将显示具有观察水平随机效应的模型会在哪里失败(例如,模型倾向于低估 Y 的低值):http: //newprairiepress.org/cgi/viewcontent.cgi?article=1005&context=agstatconference

此特定示例表明,添加观察水平随机效应可能会恶化残差与拟合值的图的外观,同时生成其他看起来不错的诊断图。我想知道这个站点上的其他人是否能够进一步了解在这种情况下应该如何进行,而不是报告使用过度分散校正时每个诊断图会发生什么。