r - glmmTMB，事后测试和 glht

Question

我正在使用 glmmTMB 分析负二项式广义线性混合模型 (GLMM)，其中因变量是过度分散的计数数据 (CT)。

相关数据框中有 115 个样本（行）。有两个固定效应（F1，F2）和一个随机截距（R），其中嵌套了另一个随机效应（NR）。还有一个偏移量，由每个样本中总计数的自然对数 (LOG_TOT) 组成。

数据框 df 的示例是：

CT  F1  F2  R   NR  LOG_TOT
77  0   0   1   1   12.9
167 0   0   2   6   13.7
289 0   0   3   11  13.9
253 0   0   4   16  13.9
125 0   0   5   21  13.7
109 0   0   6   26  13.6
96  1   0   1   2   13.1
169 1   0   2   7   13.7
190 1   0   3   12  13.8
258 1   0   4   17  13.9
101 1   0   5   22  13.5
94  1   0   6   27  13.5
89  1   25  1   4   13.0
166 1   25  2   9   13.6
175 1   25  3   14  13.7
221 1   25  4   19  13.8
131 1   25  5   24  13.5
118 1   25  6   29  13.6
58  1   75  1   5   12.9
123 1   75  2   10  13.4
197 1   75  3   15  13.7
208 1   75  4   20  13.8
113 1   8   1   3   13.2
125 1   8   2   8   13.7
182 1   8   3   13  13.7
224 1   8   4   18  13.9
104 1   8   5   23  13.5
116 1   8   6   28  13.7
122 2   0   1   2   13.1
115 2   0   2   7   13.6
149 2   0   3   12  13.7
270 2   0   4   17  14.1
116 2   0   5   22  13.5
94  2   0   6   27  13.7
73  2   25  1   4   12.8
61  2   25  2   9   13.0
185 2   25  3   14  13.8
159 2   25  4   19  13.7
125 2   25  5   24  13.6
75  2   25  6   29  13.5
121 2   8   1   3   13.0
143 2   8   2   8   13.8
219 2   8   3   13  13.9
191 2   8   4   18  13.7
98  2   8   5   23  13.5
115 2   8   6   28  13.6
110 3   0   1   2   12.8
123 3   0   2   7   13.6
210 3   0   3   12  13.9
354 3   0   4   17  14.4
160 3   0   5   22  13.7
101 3   0   6   27  13.6
69  3   25  1   4   12.6
112 3   25  2   9   13.5
258 3   25  3   14  13.8
174 3   25  4   19  13.5
171 3   25  5   24  13.9
117 3   25  6   29  13.7
38  3   75  1   5   12.1
222 3   75  2   10  14.1
204 3   75  3   15  13.5
235 3   75  4   20  13.7
241 3   75  5   25  13.8
141 3   75  6   30  13.9
113 3   8   1   3   12.9
90  3   8   2   8   13.5
276 3   8   3   13  14.1
199 3   8   4   18  13.8
111 3   8   5   23  13.6
109 3   8   6   28  13.7
135 4   0   1   2   13.1
144 4   0   2   7   13.6
289 4   0   3   12  14.2
395 4   0   4   17  14.6
154 4   0   5   22  13.7
148 4   0   6   27  13.8
58  4   25  1   4   12.8
136 4   25  2   9   13.8
288 4   25  3   14  14.0
113 4   25  4   19  13.5
162 4   25  5   24  13.7
172 4   25  6   29  14.1
2   4   75  1   5   12.3
246 4   75  3   15  13.7
247 4   75  4   20  13.9
114 4   8   1   3   13.1
107 4   8   2   8   13.6
209 4   8   3   13  14.0
190 4   8   4   18  13.9
127 4   8   5   23  13.5
101 4   8   6   28  13.7
167 6   0   1   2   13.4
131 6   0   2   7   13.5
369 6   0   3   12  14.5
434 6   0   4   17  14.9
172 6   0   5   22  13.8
126 6   0   6   27  13.8
90  6   25  1   4   13.1
172 6   25  2   9   13.7
330 6   25  3   14  14.2
131 6   25  4   19  13.7
151 6   25  5   24  13.9
141 6   25  6   29  14.2
7   6   75  1   5   12.2
194 6   75  2   10  14.2
280 6   75  3   15  13.7
253 6   75  4   20  13.8
45  6   75  5   25  13.4
155 6   75  6   30  13.9
208 6   8   1   3   13.5
97  6   8   2   8   13.5
325 6   8   3   13  14.3
235 6   8   4   18  14.1
112 6   8   5   23  13.6
188 6   8   6   28  14.1

随机和嵌套随机效应被视为因素。固定效应 F1 的值是 0、1、2、3、4 和 6。固定效应 F2 的值是 0、8、25 和 75。我将固定效应视为连续的，而不是有序的，因为我会喜欢识别因变量 CT 的单调单向变化而不是上下变化。

我之前使用 lme4 包将数据作为混合模型进行分析：

library(lme4)

m1 <- lmer(CT ~ F1*F2 + (1|R/NR) +
offset(LOG_TOT), data = df, verbose=FALSE)

随后在 multcomp 包中使用 glht 进行采用公式方法的事后分析：

library(multcomp)

glht_fixed1 <- glht(m1, linfct = c(
"F1 == 0",
"F1 + 8*F1:F2 == 0",
"F1 + 25*F1:F2 == 0",
"F1 + 75*F1:F2 == 0",
"F1 + (27)*F1:F2 == 0"))

glht_fixed2 <- glht(m1, linfct = c(
"F2 + 1*F1:F2 == 0",
"F2 + 2*F1:F2 == 0",
"F2 + 3*F1:F2 == 0",
"F2 + 4*F1:F2 == 0",
"F2 + 6*F1:F2 == 0",
"F2 + (3.2)*F1:F2 == 0"))

glht_omni <- glht(m1)

这是相应的负二项式 glmmTMB 模型，我现在更喜欢它：

library(glmmTMB)

m2 <- glmmTMB(CT ~ F1*F2 + (1|R/NR) + 
offset(LOG_TOT), data = df, verbose=FALSE, family="nbinom2")

根据 Ben Bolker 的建议 ( https://stat.ethz.ch/pipermail/r-sig-mixed-models/2017q3/025813.html )，使用 glmmTMB 进行事后测试的最佳方法是使用 lsmeans (?或其更新的等价物 emmeans）。

我听从了本的建议，跑了

source(system.file("other_methods","lsmeans_methods.R",package="glmmTMB"))

然后我可以在 glmmTMB 对象上使用 emmeans。例如，

as.glht(emmeans(m2,~(F1 + 27*F1:F2)))

General Linear Hypotheses

Linear Hypotheses:
Estimate
3.11304347826087, 21 == 0 -8.813

但这似乎不正确。我也可以将 F1 和 F2 更改为因数，然后试试这个：

as.glht(emmeans(m2,~(week + 27*week:conc)))

     General Linear Hypotheses

Linear Hypotheses:
           Estimate
0, 0 == 0    -6.721
1, 0 == 0    -6.621
2, 0 == 0    -6.342
3, 0 == 0    -6.740
4, 0 == 0    -6.474
6, 0 == 0    -6.967
0, 8 == 0    -6.694
1, 8 == 0    -6.651
2, 8 == 0    -6.227
3, 8 == 0    -6.812
4, 8 == 0    -6.371
6, 8 == 0    -6.920
0, 25 == 0   -6.653
1, 25 == 0   -6.648
2, 25 == 0   -6.282
3, 25 == 0   -6.766
4, 25 == 0   -6.338
6, 25 == 0   -6.702
0, 75 == 0   -6.470
1, 75 == 0   -6.642
2, 75 == 0   -6.091
3, 75 == 0   -6.531
4, 75 == 0   -5.762
6, 75 == 0   -6.612

但是，再一次，我不确定如何让这个输出符合我的意愿。如果哪位好心人能告诉我如何正确地将 glht 和 linfct 中的公式使用转移到 glmmTMB 的 emmeans 场景中，我将不胜感激。我已经阅读了所有的手册和小插曲，直到我脸色发青（或者至少感觉是这样），但我仍然不知所措。在我的辩护中（罪魁祸首？）我是一个统计专家，如果我在这里提出一个答案非常明显的问题，我深表歉意。

glht 软件和事后测试直接使用 glmmADMB 包，但 glmmADMB 比 glmmTMB 慢 10 倍。我需要多次运行此分析，每次运行 300,000 个负二项式混合模型示例，因此速度至关重要。

非常感谢您的建议和帮助！

Answer 1

的第二个参数 ( specs)与emmeans中的参数不同，因此不能以相同的方式使用它。您必须按照预期的方式调用它。该函数将结果转换为对象，但实际上没有必要这样做，因为摘要会产生类似的结果。linfctglhtemmeans()as.glht()glhtemmeans

我认为您试图获得的结果可以通过

emmeans(m2, ~ F2, at = list(F2 = c(0, 8, 25, 75)))

（使用带有预测变量的原始模型作为定量变量）。这将计算调整后的均值保持F1在其平均值和每个指定的值F2。

请查看emmeans(). 此外，还有许多提供解释和示例的小插曲——从https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html开始。

Answer 2

根据我优秀的统计顾问的建议，我认为下面的解决方案提供了我之前使用 glht 和 linfct 获得的结果。

F1 的斜率是通过使用对比度和 emmeans 在 F2 的各个水平上计算的，以计算 F1 的两个值之间的因变量差异，这些 F1 值相隔一个单位（即 c(0,1)）。（因为回归是线性的，所以 F1 的两个值是任意的，只要它们被一个单位分开，例如 c(3,4)）。F2 的斜坡反之亦然。

因此，F1 在 F2 = 0、8、25、75 和 27 处的斜率（27 是 F2 的平均值）：

contrast(emmeans(m1, specs="F1", at=list(F1=c(0,1), F2=0)),list(c(-1,1)))
        (above equivalent to: summary(m1)$coefficients$cond["F1",])
        contrast(emmeans(m1, specs="F1", at=list(F1=c(0,1), F2=8)),list(c(-1,1)))
        contrast(emmeans(m1, specs="F1", at=list(F1=c(0,1), F2=25)),list(c(-1,1)))
        contrast(emmeans(m1, specs="F1", at=list(F1=c(0,1), F2=75)),list(c(-1,1)))
        contrast(emmeans(m1, specs="F1", at=list(F1=c(0,1), F2=27)),list(c(-1,1)))

和 F2 在 F1 = 1、2、3、4、6 和 3.2 时的斜率（3.2 是 F1 的平均值，不包括零值）：

contrast(emmeans(m1, specs="F2", at=list(F2=c(0,1), F1=0)),list(c(-1,1)))
(above equivalent to: summary(m1)$coefficients$cond["F2",])
contrast(emmeans(m1, specs="F2", at=list(F2=c(0,1), F1=1)),list(c(-1,1)))
contrast(emmeans(m1, specs="F2", at=list(F2=c(0,1), F1=2)),list(c(-1,1)))
contrast(emmeans(m1, specs="F2", at=list(F2=c(0,1), F1=3)),list(c(-1,1)))
contrast(emmeans(m1, specs="F2", at=list(F2=c(0,1), F1=4)),list(c(-1,1)))
contrast(emmeans(m1, specs="F2", at=list(F2=c(0,1), F1=6)),list(c(-1,1)))
contrast(emmeans(m1, specs="F2", at=list(F2=c(0,1), F1=3.2)),list(c(-1,1)))

F1 = 0 和 F2 = 0 时 F1 和 F2 斜率的交互作用

contrast(emmeans(m1, specs=c("F1","F2"), at=list(F1=c(0,1),F2=c(0,1))),list(c(1,-1,-1,1)))
(above equivalent to: summary(m1)$coefficients$cond["F1:F2",])

从 提供的结果emmGrid对象中contrast()，可以根据需要挑选出斜率estimate的估计值 ( )、估计斜率的标准差 ( SE)、估计斜率与零假设斜率z.ratio(emmGrid从estimate除以SE) 和相应的 P 值 (p.value由emmGridas计算2*pnorm(-abs(z.ratio))。

例如：

contrast(emmeans(m1, specs="F1", at=list(F2=c(0,1), F1=0)),list(c(-1,1)))

产量：

NOTE: Results may be misleading due to involvement in interactions
 contrast    estimate          SE df z.ratio p.value
 c(-1, 1) 0.001971714 0.002616634 NA   0.754  0.4511

后记在 1.25 年后添加：

以上给出了正确的解决方案，但正如 Russell Lenth 指出的那样，使用 emtrends 更容易获得答案。但是，我选择这个答案是正确的，因为在展示如何使用 emmeans 计算斜率以在自变量变化 1 时找到预测因变量的结果变化时可能有一些教学价值。