Emmeans 没有从模型中给我正确的调整方法

Posted 2023-02-18

【中文标题】Emmeans 没有从模型中给我正确的调整方法

【英文标题】：Emmeans does not give me the correct adjusted means from the model

【发布时间】：2021-11-07 11:19:27

【问题描述】：

我使用 emmeans 从我的线性混合效应回归模型中得出调整后的均值，但结果似乎不正确。我想绘制模型拟合和各个数据点的调整值，但结果看起来很奇怪：

估计调整后的平均值对于课程 A 似乎太高，而在课程 C 上太低。在我的线性混合效应回归中，我预测后测与前测作为协变量以及组和课程的主要影响和交互.因为我对 Course 和不同的测试条件进行了重复测量，所以我为 Course 和 School 包括了一个随机截距。使用 emmeans 我得到以下估计：

# model fit
CI_post <- lmer(
  post.diff ~ 
    pre.diff +
    group * course 
  + (1|bib) 
  + (1|school), 
  data = dat, 
  REML = FALSE)

#estimated adjusted means
emmeans(CI_post, specs = c("course", "group"),lmer.df = "satterthwaite")

# Results
 course group       emmean    SE   df lower.CL upper.CL
 A      blocked      0.311 0.191 6.65  -0.1452    0.768
 B      blocked      0.649 0.180 5.38   0.1954    1.102
 C      blocked      1.141 0.195 7.28   0.6847    1.598
 A      interleaved  0.189 0.194 7.15  -0.2666    0.645
 B      interleaved  0.497 0.179 5.31   0.0451    0.949
 C      interleaved  1.046 0.191 6.72   0.5907    1.502

我绘制的正是这些值，我认为这些值是不正确的。有人可以帮助我获得正确的估计调整均值吗？

看了this后，我怀疑是因为pre.diff是一个固定值的错误？

ref_grid(CI_post)

#result
'emmGrid' object with variables:
    pre.diff = 1.5065
    group = blocked, interleaved
    course = A, B, C

编辑 按照 Lenth 的建议，我尝试了： post.diff.adj = post.diff + b * (1.506 - pre.diff)，这给了我下图：

它看起来更好，更正确。我使用了模型中的模型回归系数：

Fixed effects:
                          Estimate Std. Error        df t value             Pr(>|t|)    
(Intercept)               -0.66087    0.18158   5.58701  -3.639             0.012280 *  
pre.diff                   0.64544    0.06178 130.60667  10.448 < 0.0000000000000002 ***
groupinterleaved          -0.12209    0.15189  65.38709  -0.804             0.424431    
courseB                    0.33714    0.09703 131.63603   3.475             0.000693 ***
courseC                    0.82993    0.16318 151.09201   5.086           0.00000107 ***
groupinterleaved:courseB  -0.02922    0.11777 101.47596  -0.248             0.804563    
groupinterleaved:courseC   0.02692    0.11763 100.29319   0.229             0.819435 

然后我在我的 tibble 中使用了计算它：

dat <- dat %>%
  mutate(adjustedMean = (post.diff) + (0.6454358 * (1.506 - pre.diff)))

然后我用 ggplot 绘制它：

CI_post_plot <- ggplot(dat, aes(x = interaction(group, course), y = adjustedMean)) +
  geom_point(aes(color=group), size=1.5, position=position_jitter(width=0.1), alpha=0.7)+
  scale_y_continuous(name = "Time substracted from straight gliding time (sec.)", breaks = seq(-2, 6, 1)) +
  theme_pubr()+
  theme(legend.position="none",
        axis.title.x=element_blank()) +
  geom_hline(aes(yintercept=0), linetype = "dashed", size=0.2) + 
  scale_x_discrete(labels = c("Blocked\nCourse A", "Interleaved\nCourse A", "Blocked\nCourse B", "Interleaved\nCourse B", "Blocked\nCourse C",  "Interleaved\nCourse C")) 

CI_post_plot <- CI_post_plot + 
  geom_point(data = estmarg_mean, aes(x=interaction(group, course), y=emmean, group=group), size=2.5) +
  geom_errorbar(data = estmarg_mean, aes(x= interaction(group, course), y = emmean, ymin = lower.CL,ymax = upper.CL), width=0.1)

https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html

【问题讨论】：

仅仅因为它们与观察到的均值不匹配并不意味着它们不正确。您在模型中有一个协变量 pre.diff，并且 EMM 使用该协变量平均值（大约 1.5）的预测。如果事实上pre.diff 与这些因素有关，那可能会对 EMM 产生很大影响。

干杯。这是最好的方法还是有更好的方法？

我想说，如果获得调整均值的统计目标是正确的，那么这是获得它们的好方法。但是该图显示了未经调整的数据和调整后的平均值。我想也许你可以通过post.diff.adj = b * (1.506 - pre.diff) 创建调整后的响应值，其中b 是拟合模型中pre.diff 的回归系数。这会减去估计的协变量效应，并在平均 pre.diff 值处添加协变量效应。

PS 如果你试试这个，我很想看看修改后的情节；也许您可以将其添加到您的帖子中。

太棒了。我建议在第二个的 y 轴标签上添加“调整”

【参考方案1】：

回顾一些 cmets，OP 中的第二个图显示了调整后的响应值和调整后的均值 (AKA EMM)。这个粗略的草图显示了这一点： 大多数实验设计文本将显示如何调整均值的类似图片：该模型适合每个处理的平行线；这些线穿过各自数据云的中心。调整后的均值是协变量平均值的估计值。

为了获得调整后的数据，我们对每个数据点做同样的事情；显示了几个具有代表性的点。调整后的数据是每个数据点在平均线上的投影，投影的路径与回归线平行。

几年前The American Statistician 中有一篇关于“对齐数据”的文章，与此类似。我不记得作者，也没有在快速搜索中找到它。这也与一些回归文本中讨论的“分量加残差”图有关。其基本思想是从数据中去除干扰变量的估计影响，或者等效地获得模型残差并添加非干扰影响。