通过带有 DHARMa 包的 glmmTMB 检查 beta 回归模型

机器算法验证 r 残差 随机效应模型 咕噜咕噜 glmmtmb
2022-03-30 04:21:42

我想澄清一下我的模型是否明确指定（因为我对 Beta 回归模型没有太多经验）。

我的变量是假牙中污垢面积的百分比。对于每位患者，牙医都会在假牙的左侧或右侧使用特殊产品（将另一侧作为安慰剂）以去除污垢区域。

之后，他计算了假牙每一侧的总面积，以及每一侧的总污垢面积。

我需要测试产品是否能有效去除污垢。

我的初始模型（prop.bio 是污垢区域的比例）：

library(glmmTMB)    
m1 <- glmmTMB(prop.bio ~ Product*Side + (1|Pacients), data, family=list(family="beta",link="logit"))

更新：

我通过 TRV 测试手动向后选择后的最终模型（这也是研究人员的主要问题）：

m1.f <- glmmTMB(prop.bio ~ Product + (1|Pacients), data, family=list(family="beta",link="logit"))

我的残留诊断使用DHARMa：

library(DHARMa)
res = simulateResiduals(m1.f)
plot(res, rank = T)

在此处输入图像描述

根据我对DHARMa小插图的阅读，基于正确的情节，我的模型可能是错误的。那我该怎么办？（我的型号规格错了吗？）

提前致谢！

数据：

structure(list(Pacients = structure(c(5L, 6L, 2L, 11L, 26L, 29L, 
20L, 24L, 8L, 14L, 19L, 7L, 13L, 4L, 3L, 5L, 6L, 2L, 11L, 26L, 
29L, 20L, 24L, 8L, 14L, 19L, 7L, 13L, 4L, 3L, 23L, 25L, 12L, 
21L, 10L, 22L, 18L, 27L, 15L, 9L, 17L, 28L, 1L, 16L, 23L, 25L, 
12L, 21L, 10L, 22L, 18L, 27L, 15L, 9L, 17L, 28L, 1L, 16L), .Label = c("Adlf", 
"Alda", "ClrW", "ClsB", "CrCl", "ElnL", "Gema", "Héli", "Inác", 
"Inlv", "InsS", "Ircm", "Ivnr", "Lnld", "Lrds", "LusB", "Mart", 
"Mrnz", "Murl", "NGc1", "NGc2", "Nlcd", "Norc", "Oliv", "Ramr", 
"Slng", "Svrs", "Vldm", "Vlsn"), class = "factor"), Area = c(3942, 
3912, 4270, 4583, 2406, 2652, 2371, 4885, 3704, 3500, 4269, 3743, 
3414, 4231, 3089, 4214, 3612, 4459, 4678, 2810, 2490, 2577, 4264, 
4287, 3487, 4547, 3663, 3199, 3836, 3237, 3846, 4116, 3514, 3616, 
3609, 4053, 3810, 4532, 4380, 4103, 4552, 3745, 3590, 3386, 3998, 
4449, 3367, 3698, 3840, 4457, 3906, 4384, 4000, 4156, 3594, 3258, 
4094, 2796), Side = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L), .Label = c("Right", "Left"), class = "factor"), Biofilme = c(1747, 
1770, 328, 716, 1447, 540, 759, 1328, 2320, 1718, 1226, 977, 
1193, 2038, 1685, 2018, 1682, 416, 679, 2076, 947, 1423, 1661, 
1618, 1916, 1601, 1833, 1050, 1780, 1643, 1130, 2010, 2152, 812, 
2550, 1058, 826, 1526, 2905, 1299, 2289, 1262, 1965, 3016, 1630, 
1823, 1889, 1319, 2678, 1205, 472, 1694, 2161, 1444, 1062, 819, 
2531, 2310), Product = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L), .Label = c("No", "Yes"), class = "factor"), prop.bio = c(0.443176052765094, 
0.452453987730061, 0.0768149882903981, 0.156229543966834, 0.601413133832086, 
0.203619909502262, 0.320118093631379, 0.271852610030706, 0.626349892008639, 
0.490857142857143, 0.287186694776294, 0.261020571733903, 0.349443468072642, 
0.481682817300874, 0.545483975396568, 0.478879924062648, 0.465669988925803, 
0.0932944606413994, 0.145147498931167, 0.738790035587189, 0.380321285140562, 
0.552192471866511, 0.389540337711069, 0.377420107301143, 0.549469457986808, 
0.352100285902793, 0.5004095004095, 0.328227571115974, 0.464025026068822, 
0.507568736484399, 0.293811752470099, 0.488338192419825, 0.612407512805919, 
0.224557522123894, 0.706566916043225, 0.261041204046385, 0.216797900262467, 
0.336716681376876, 0.66324200913242, 0.316597611503778, 0.502855887521968, 
0.3369826435247, 0.547353760445682, 0.890726520968695, 0.407703851925963, 
0.409755001123848, 0.561033561033561, 0.356679286100595, 0.697395833333333, 
0.270361229526587, 0.12083973374296, 0.386405109489051, 0.54025, 
0.347449470644851, 0.295492487479132, 0.251381215469613, 0.618221787982413, 
0.82618025751073)), row.names = c(NA, -58L), class = "data.frame")
3个回答

tl; dr您担心是有道理的，但是在查看了各种不同的图形诊断后，我认为一切看起来都不是很好。我的回答将说明一系列其他查看合身的方法glmmTMB- 比 DHARMa 更复杂/更不方便，但尽可能多地以不同的方式看待合身是很好的。

首先让我们看一下原始数据（我称之为dd）：

library(ggplot2); theme_set(theme_bw())
ggplot(dd,aes(Product,prop.bio,colour=Side))+
    geom_line(colour="gray",aes(group=Pacients))+
    geom_point(aes(shape=Side))+
    scale_colour_brewer(palette="Dark2")

在此处输入图像描述

我的第一点是由（通常是所有预测与残差图）制作的右手图DHARMa正在寻找模型中的偏差，即残差相对于均值具有系统模式的模式。对于只有分类预测变量的模型（假设它包含预测变量的所有可能交互作用），这绝不应该发生，因为该模型对每个可能的拟合值都有一个参数......我们将在下面看到，如果我们在总体水平而不是个体水平上查看拟合与残差...

获得拟合与残差图的最快方法（例如，类似于 base-R 的plot.lm()方法或lme4's plot.merMod()）是通过broom.mixed::augment()+ ggplot：

library(broom.mixed)
aa <- augment(m1.f, data=dd)
gg2 <- (ggplot(aa, aes(.fitted,.resid))
    + geom_line(aes(group=Pacients),colour="gray")
    + geom_point(aes(colour=Side,shape=Product))
    + geom_smooth()
)

在此处输入图像描述

这些拟合值和残差值处于个体患者水平。它们确实显示出温和的趋势（我承认目前我不理解），但相对于数据的分散性而言，总体趋势似乎并不大。

为了检查这种现象确实是由患者而不是人口水平的预测引起的，并检验上面的论点，即人口水平效应在拟合与残差图中的趋势应该恰好为零，我们可以破解glmmTMB预测来构建人口水平的预测和残差（下一个版本glmmTMB应该更容易）：

aa$.fitted0 <- predict(m1.f, newdata=transform(dd,Pacients=NA),type="response")
aa$.resid0 <- dd$prop.bio-aa$.fitted0
gg3 <- (ggplot(aa, aes(.fitted0,.resid0))
    + geom_line(aes(group=Pacients),colour="gray")
    + geom_point(aes(colour=Side,shape=Product))
    + geom_smooth()
)

（请注意，如果您运行此代码，您会收到很多警告geom_smooth()，当预测变量 [即拟合值] 仅具有两个唯一级别时运行时，您会感到不高兴）

在此处输入图像描述

现在残差的平均值（几乎？）对于两个级别（Product=="No"和Product=="Yes"）都为零。

只要我们还在，让我们检查一下随机效应的诊断：

lme4:::dotplot.ranef.mer(ranef(m1.f)$cond)

在此处输入图像描述

这看起来不错：没有不连续跳跃的迹象（表明随机效应中可能存在多模态）或​​异常患者。

其他的建议

  • 我不赞成根据哪些术语似乎很重要来减少模型的一般原则（例如Side，在运行后从模型中删除anova()）：通常，数据驱动的模型减少会混淆推理。

看看DHARMa的小插图中关于glmmTMB的部分。考虑到随机效应，如何计算预测似乎是一个问题。

作为替代方案，您可以尝试使用GLMMadaptive包。您可以在此处找到使用DHARMa 的示例。

我是 DHARMa 的开发者。Dimitris 和 Ben 是正确的，该模式源于 glmmTMB （目前）不允许仅基于固定效应进行预测的已知问题，这有时会产生这种模式。我希望我们可以在下一个版本的 glmmTMB 中解决这个问题，它应该允许固定效果预测。

[编辑 11 月 21 日：这个问题大约 1 年前在 glmmTMB 中得到修复。]

在您的情况下，很明显模型中的预测变量基于固定效应和随机效应，因为您的固定效应只有一个分类预测变量，因此您的 x 轴上应该只有 2 个值。我们可以轻松地手动生成仅使用固定效应作为预测变量的图：

plotResiduals(data$Product, res$scaledResiduals)

这导致一个看起来不错的情节

在此处输入图像描述

顺便说一句，同意 Ben 的观点，我不会根据重要性进行模型选择，这本质上是 p-hacking。如果您从 Product*Side 开始，请报告此模型，除非您认为推理存在严重问题。

其它你可能感兴趣的问题
上一篇使用 scikit-learn 库的 StandardScaler 函数 下一篇对最大操作的期望