如何使用sjPlot可视化协调afex混合效果模型输出

问题描述

我在R中使用了afex::mixed的混合效果模型（通常我使用lme4::lmer，但我读到||表示法不适用于该软件包see here中的分类变量，如下所示：

>str(DF)
'data.frame':   1521 obs. of  3 variables:
 $ p: Factor w/ 100 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
 $ a: Factor w/ 2 levels "Down","Up": 2 2 2 2 2 2 2 2 2 2 ...
 $ y: num  12 0 13 0 0 10 5 0 0 5 ...
>
>
> # fit mixed effects model with afex::mixed
> m1 <- mixed(y ~ a + (a||p),# random slopes and intercepts by participant,no correlation parameters
+             data = DF,+             expand_re = TRUE,+             method = "S",+             return = "merMod")
Contrasts set to contr.sum for the following variables: a,p

我主要感兴趣的是a对y的固定作用，所以我这样检查：

> fixef(m1)
(Intercept)          a1 
   6.837455    4.608073 

我（可能是错误的？）将其解释为该模型预测，当a == "down"时，y将比a == "up"时大4.61。

然后我用sjPlot::plot_model可视化模型，如下所示：

plot_model(m1,type = "pred",terms = "a")

为什么该图似乎显示出比模型的统计输出更大的固定效果？两点之差不应该是4.61吗？如果没有，该图显示什么？

---

我可以使用sjPlot::get_model_data检索要显示的确切数据：

> get_model_data(m1,terms = "a",type = "pred")

# Predicted values of y
# x = a

x | Predicted |   SE | group_col |         95% CI
-------------------------------------------------
1 |     11.45 | 0.40 |         1 | [10.65,12.24]
2 |      2.23 | 0.55 |         1 | [ 1.15,3.31]

Adjusted for:
* re1.a1 = -0.05
*      p = 0 (population-level)

但是我仍然对它与模型输出的关系，特别是a的固定效果感到困惑。

---

我还可以使用sjPlot重新创建afex::afex_plot可视化文件：

> afex_plot(m1,x = "a",mapping = c("color"))
Aggregating data over: p 

这里最有趣的是消息Aggregating data over: p。这是否表示未考虑p的随机效应？如果是这样，这与根本没有引用模型的绘图有何不同？例如，尽管只是一条y ~ x行，但以下内容似乎绘制了相同的值...

ggplot(DF,aes(x = a,y = y,color = a))+
  geom_smooth(aes(group = 1),method = "lm",se = T,color = "black")

解决方法

这是由于变量从 Afex 获取 sum 对比引起的，这是对 lme4 的标准虚拟编码的重要改进。

本质上你说的是对的，你可能错误地解释了模型系数，因为和对比编码可以将你的变量级别“向下”设置为 -1，而将“向上”设置为 1。简单的算术表明这可能是右：6.83 - 4.61 ~ 2.23； 6.83 + 4.61 ~ 11.45。

您可以尝试查看模型中的因子水平如何编码以解决此问题。我自己无法做到这一点，因为您没有提供可复制的示例。

关于和对比的解释：https://rpubs.com/monajhzhu/608609

,

出于好奇，我稍后将尝试对此进行更多研究，但这会产生任何用处吗？

library(dplyr)
library(broom.mixed)

augment(m1,DF) %>%
  group_by(a) %>%
  summarise(with_re = mean(.fitted,na.rm = TRUE),with_re_se = sd(.fitted,na.rm = TRUE) / sqrt(n()),fixed_only = mean(.fixed,fixed_only_se = sd(.fixed,na.rm = TRUE) / sqrt(n()))