• [主要编辑以斜体显示]*

使用R I从正态分布生成数据，但是，当我将随机截距混合效应模型拟合到它们时，* 正态性的Kolmogorov-Smirnoff检验和分位数-分位数图 * 都由R包DHARMa执行，表明偏离正态性。为什么会这样？

测试1

• 我生成了一组随机截距数据，其中重复次数和组内标准差随随机效应的每个水平而变化：*

rm(list=ls())
library(DHARMa)
library(glmmTMB)
library(ggplot2)

# "regression" with an 8-level RE with variable slopes:
set.seed(666)
ii=8 # parameter set size; no. groups
reps <- sample(x=8:12, size=ii, replace=T)
slope <- 3
intercept <- runif(ii, min = 70, max = 140)
group.sd <- runif(ii, min = 5, max = 15)
group.ID <- LETTERS[1:ii]

xx <- 1:15

out.list <- list() # empty list to store simulated data

for (i in 1:ii){
        df00 <- data.frame(Group=rep(group.ID[i], reps[i]*length(xx)),
                                    x=rep(xx, reps[i]))
        df00$y = intercept[i] + df00$x * slope + rnorm(length(df00[,1]), 0, group.sd[i])
        out.list[[i]] = df00
        }

# to turn out.list into a data.frame:
df <- do.call(rbind.data.frame, out.list)

# data visualisation
ggplot(df, aes(x = x, y = y, colour = Group)) + 
    geom_point() + geom_smooth(method="lm") + theme_classic()

如果我将随机截距模型拟合到我的数据集：

glmmTMB_ri <- glmmTMB(y~x + (1|Group), data=df)

以下是DHARMa生成的诊断图和测试：

mem1Resids <- simulateResiduals(glmmTMB_ri, n = length(df[, 1]) * 3)
par(mfrow = c(1, 2))
plotQQunif(mem1Resids)
mtext(text = "(a)", side = 3, adj = 0, line = 2)
plotResiduals(mem1Resids, quantreg = T)
mtext(text = "(b)", side = 3, adj = 0, line = 2)

测试2

• 我认为问题可能是模拟残差的分辨率太低。* 通过simulateResiduals()将模拟观测值的数量从数据集大小的三倍增加到数据集大小的10倍并没有改善诊断图：

测试3

我认为问题可能不是DHARMa的设置，而是我的数据集 （可能太小或太可变）。 我生成了一个新的数据集，每个 * 随机 * 组具有***更大，统一的复制水平***（每组每个x值有40个观察结果），并且具有***更低，统一的组间变异性***（所有 * 随机 * 组的SD = 5）：

reps <- rep(40, length.out=ii)
# let's also reduce and uniformise maximum within-group variability:
group.sd <- rep(5, ii)

out.list <- list() # empty list to store simulated data

for (i in 1:ii){
        df00 <- data.frame(Group=rep(group.ID[i], reps[i]*length(xx)),
                                    x=rep(xx, reps[i]))
        df00$y = intercept[i] + df00$x * slope + rnorm(length(df00[,1]), 0, group.sd[i])
        out.list[[i]] = df00
        }

# to turn out.list into a data.frame:
df2 <- do.call(rbind.data.frame, out.list)

# data visualisation
ggplot(df2, aes(x = x, y = y, colour = Group)) +
    geom_point() + geom_smooth(method="lm") + theme_classic()

如果我将随机截距模型拟合到我的新数据集：

glmmTMB_ri2 <- glmmTMB(y~x + (1|Group), data=df2)

以下是DHARMa生成的诊断图和测试：

mem2Resids <- simulateResiduals(glmmTMB_ri2, n = length(df2[, 1]) * 3)
par(mfrow = c(1, 2))
plotQQunif(mem2Resids)
mtext(text = "(a)", side = 3, adj = 0, line = 2)
plotResiduals(mem2Resids, quantreg = T)
mtext(text = "(b)", side = 3, adj = 0, line = 2)

如果有什么不同的话，诊断图显示的情况比以前更糟。
为什么DHARMa表明我的模型的残差分布的非正态性，即使我生成的数据是高斯分布？如果问题不在于样本量或DHARMa，也许我模拟的数据不正确，但似乎也不是这样。