添加观察水平随机项会混淆残差与拟合图。为什么?
谢谢你更新你的帖子,查理。我使用了一些过度分散的泊松数据,以了解在 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(尤其是左侧面板)的诊断图清楚地表明过度分散是一个问题。

有关这些类型的绘图的更多详细信息,请参阅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。
此特定示例表明,添加观察水平随机效应可能会恶化残差与拟合值的图的外观,同时生成其他看起来不错的诊断图。我想知道这个站点上的其他人是否能够进一步了解在这种情况下应该如何进行,而不是报告使用过度分散校正时每个诊断图会发生什么。




