为何DHARMa对高斯生成数据的混合模型检测出正态性偏离？

我使用R生成了正态分布数据，但对其拟合随机截距混合效应模型后，R包DHARMa执行的Kolmogorov-Smirnoff正态性检验和分位数-分位数图均显示数据偏离正态分布。这是为什么？

测试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()生成的模拟观测数从数据集大小的3倍增加到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，难道是我模拟数据的方式有误？但看起来并非如此。

内容的提问来源于stack exchange，提问作者Marco Plebani