lmer() 的事后测试出错：multcomp() 和 emmeans()

Posted 2023-02-18

【中文标题】lmer() 的事后测试出错：multcomp() 和 emmeans()

【英文标题】：Error in post hoc test for lmer(): both multcomp() and emmeans()

【发布时间】：2021-11-20 21:47:34

【问题描述】：

我有一个不同位置的“Y”测量数据集，我试图通过运行lmer() 模型并分析结果来确定变量 Y 如何受变量 A、B 和 D 的影响。但是，当我到达 post hoc 步骤时，我在尝试分析时收到错误消息。

这是我的数据示例：

table <- "   ID location A         B      C     D       Y
1   1       AA 0 0.6181587 -29.67 14.14 168.041
2   2       AA 1 0.5816176 -29.42 14.21 200.991
3   3       AA 2 0.4289670 -28.57 13.55 200.343
4   4       AA 3 0.4158891 -28.59 12.68 215.638
5   5       AA 4 0.3172721 -28.74 12.28 173.299
6   6       AA 5 0.1540603 -27.86 14.01 104.246
7   7       AA 6 0.1219355 -27.18 14.43 128.141
8   8       AA 7 0.1016643 -26.86 13.75 179.330
9   9       BB 0 0.6831649 -28.93 17.03 210.066
10 10       BB 1 0.6796935 -28.54 18.31 280.249
11 11       BB 2 0.5497743 -27.88 17.33 134.023
12 12       BB 3 0.3631052 -27.48 16.79 142.383
13 13       BB 4 0.3875498 -26.98 17.81 136.647
14 14       BB 5 0.3883785 -26.71 17.56 142.179
15 15       BB 6 0.4058061 -26.72 17.71 109.826
16 16       CC 0 0.8647298 -28.53 11.93 220.464
17 17       CC 1 0.8664036 -28.39 11.59 326.868
18 18       CC 2 0.7480748 -27.61 11.75 322.745
19 19       CC 3 0.5959143 -26.81 13.27 170.064
20 20       CC 4 0.4849077 -26.77 14.68 118.092
21 21       CC 5 0.3584687 -26.65 15.65  95.512
22 22       CC 6 0.3018285 -26.33 16.11  71.717
23 23       CC 7 0.2629121 -26.39 16.16  60.052
24 24       DD 0 0.8673077 -27.93 12.09 234.244
25 25       DD 1 0.8226558 -27.96 12.13 244.903
26 26       DD 2 0.7826429 -27.44 12.38 252.485
27 27       DD 3 0.6620447 -27.23 13.84 150.886
28 28       DD 4 0.4453213 -27.03 15.73 102.787
29 29       DD 5 0.3720257 -27.13 16.27 109.201
30 30       DD 6 0.6040217 -27.79 16.41 101.509
31 31       EE 0 0.8770987 -28.62 12.72 239.036
32 32       EE 1 0.8504547 -28.47 12.92 220.600
33 33       EE 2 0.8329484 -28.45 12.94 174.979
34 34       EE 3 0.8181102 -28.37 13.17 138.412
35 35       EE 4 0.7942685 -28.32 13.69 121.330
36 36       EE 5 0.7319724 -28.22 14.62 111.851
37 37       EE 6 0.7014828 -28.24 15.04 110.447
38 38       EE 7 0.7286984 -28.15 15.18 121.831"

#Create a dataframe with the above table
df <- read.table(text=table, header = TRUE)
df

# Make sure location is a factor
df$location<-as.factor(df$location)

这是我的模型：

# Load libraries
library(ggplot2)
library(pscl)
library(lmtest)
library(lme4)
library(car)

mod = lmer(Y ~ A * B * poly(D, 2) * (1|location), data = df)
summary(mod)
plot(mod)

我现在需要确定哪些变量对 Y 有显着影响。因此我从包 car 中运行了 Anova()（输出粘贴在这里）。

Anova(mod)
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Y
# Chisq Df Pr(>Chisq)    
# A                 8.2754  1   0.004019 ** 
# B                 0.0053  1   0.941974    
# poly(D, 2)        40.4618  2  1.636e-09 ***
# A:B               0.1709  1   0.679348    
# A:poly(D, 2)      1.6460  2   0.439117    
# B:poly(D, 2)      5.2601  2   0.072076 .  
# A:B:poly(D, 2)    0.6372  2   0.727175    
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这表明：

A 显着影响 Y

B 对 Y 没有显着影响

D 显着影响 Y

所以接下来我将对这些变量中的每一个进行事后测试，但这是我遇到问题的地方。我已经尝试使用下面的emmeans 和multcomp 包：

library(emmeans)
emmeans(mod, list(pairwise ~ A), adjust = "tukey")
# NOTE: Results may be misleading due to involvement in interactions
# Error in if ((misc$estType == "pairs") && (paste(c("", by), collapse = ",") !=  : 
#  missing value where TRUE/FALSE needed

pairs(emmeans(mod, "A"))
# NOTE: Results may be misleading due to involvement in interactions
# Error in if ((misc$estType == "pairs") && (paste(c("", by), collapse = ",") !=  : 
#  missing value where TRUE/FALSE needed

library(multcomp)
summary(glht(mod, linfct = mcp(A = "Tukey")), test = adjusted("fdr"))
# Error in h(simpleError(msg, call)) : 
#  error in evaluating the argument 'object' in selecting a method for function 'summary': Variable(s) ‘depth’ of class ‘integer’ is/are not contained as a factor in ‘model’.

这是我第一次在lmer() 模型上运行方差分析/事后测试，虽然我已经阅读了一些关于这个模型的介绍性网站，但我不确定我是否正确地测试了它。任何帮助将不胜感激。

【参考方案1】：

如果我正确地查看数据，A 是值为 0、1、...、7 的变量。现在查看您的方差分析表，您会发现 A 只有 1 df，而不是 7它应该是一个具有 8 个级别的因子。这意味着您的模型将 A 作为数字预测器——这毫无意义。将 A 变成一个因子并重新拟合他的模型。你会有更好的运气。

我还认为您的意思是在模型公式的末尾加上 + (1|location)，而不是让随机效应与某些多项式效应相互作用。

【讨论】：

嗨，我尝试进行这些调整，但在尝试运行事后测试时仍然收到错误。

您需要重新拟合模型。然后再次运行 Anova() 因为测试会改变，并确认 A 仍然没有 1 d.f.