glht() 和 lsmeans() 在 lmer() 模型中找不到对比

Question

我有以下情况：

我的固定效应模型在称为“作曲家”的参与者组中找到了 Relation_PenultimateLast 的主效应。因此，我想找出 Relation_PenultimateLast 的水平在统计上与其他水平有何不同。

f.e.model.composers = lmer(Score ~ Relation_PenultimateLast + (1|TrajectoryType) + (1|StimulusType) + (1|Relation_FirstLast) + (1|LastPosition), data=datasheet.complete.composers)

摘要（f.e.model.composers）

Random effects:
 Groups             Name        Variance Std.Dev.
 TrajectoryType     (Intercept) 0.005457 0.07387 
 LastPosition       (Intercept) 0.036705 0.19159 
 Relation_FirstLast (Intercept) 0.004298 0.06556 
 StimulusType       (Intercept) 0.019197 0.13855 
 Residual                       1.318116 1.14809 
Number of obs: 2200, groups:  
TrajectoryType, 25; LastPosition, 8; Relation_FirstLast, 4; StimulusType, 4

Fixed effects:
                         Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)               2.90933    0.12476 14.84800  23.320 4.15e-13 ***
Relation_PenultimateLast  0.09987    0.02493 22.43100   4.006 0.000577 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我必须对我的 lmer() 模型进行 Tukey 比较。 现在，我找到了两种比较 Relation_PenultimateLast 级别的方法（我在这里找到了它们：https://stats.stackexchange.com/questions/237512/how-to-perform-post-hoc-test-on-lmer-model）：

summary(glht(f.e.model.composers, linfct = mcp(Relation_PenultimateLast = "Tukey")), test = adjusted("holm"))

和

lsmeans(f.e.model.composers, list(pairwise ~ Relation_PenultimateLast), adjust = "holm")

这些不起作用。 前者报道：

Variable(s) ‘Relation_PenultimateLast’ of class ‘integer’ is/are not contained as a factor in ‘model’

后者：

Relation_PenultimateLast   lsmean        SE  df lower.CL upper.CL
                      2.6 3.168989 0.1063552 8.5 2.926218  3.41176

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

    $` of contrast`
     contrast  estimate SE df z.ratio p.value
     (nothing)   nonEst NA NA      NA      NA

有人可以帮我理解为什么我会得到这个结果吗？

Answer 1

首先，重要的是要意识到您安装的模型不合适。它使用Relation_PenultimateLast 作为数字预测器；因此，它适合其值 1、2、3 和 4 的线性趋势，而不是作为一个因素对每个水平的单独估计。我还想知道，鉴于您展示的情节，为什么Test 在模型中不是；它看起来应该是（再次作为一个因素，而不是一个数字预测器）。我建议您获得一些统计咨询帮助，以检查您在研究中是否使用了适当的模型。或许你可以给统计学专业的研究生一些实际应用的基础——一个双赢的主张。

要将Relation_PenultimateLast 建模为因子，一种方法是在模型公式中将其替换为factor(Relation_PenultimateLast)。这适用于lsmeans()，但不适用于glht()。更好的方法可能是在数据集中进行更改：

datasheet.complete.composers = transform(datasheet.complete.composers,
    Relation_PenultimateLast = factor(Relation_PenultimateLast))
f.e.model.composers = lmer(...) ### (as before, assuming Test isn't needed)

（顺便说一句，你一定是一个比我更好的打字员；我会使用更短的名字，虽然我很赞成使用信息丰富的名字。）

（注意：f.e.model.composers 是否应该建议一个固定效应模型？它不是一个；它是一个混合模型。再次，顾问......）

lsmeans 包注定要被弃用，所以我建议你使用它的延续，emmeans 包：

library(emmeans)
emmeans(f.e.model.composers, pairwise ~ Relation_PenultimateLast)

我建议对这个应用程序使用默认的"tukey" 调整而不是 Holm。

如果确实Test 应该在模型中，那么看起来您需要包含交互；所以它会是这样的：

model.composers = lmer(Score ~ Relation_PenultimateLast * factor(Test) + ...)

### A plot like the one shown, but based on the model predictions:
emmip(model.composers, Relation_PenultimateLast ~ Test) 

### Estimates and comparisons of Relation_PenultimateLast for each Test:
emmeans(model.composers, pairwise ~ Relation_PenultimateLast | Test)