使用 ggplot2 进行分层模型 （lmerMod） 的交互作用图，并突出显示统计显着性

Interaction plot using ggplot2 for hierarchical model (lmerMod) with statistical significance highlighting

提问人：CompSocialSciR 提问时间：9/21/2023 更新时间：10/5/2023 访问量：56

问：

我正在使用一个分层模型 （lmerMod），该模型具有一个结果和两个分类变量，以及 .我想制作一个交互作用图，用它来可视化模型系数并考虑两个特定的复杂性：cat1cat2ggplot2::

成对比较：我需要比较每个级别（对照、T1、T2、T3）在（1 和 2）的每个级别内的系数。 统计意义：我想使用星号或类似符号来强调这些比较的统计意义。例如，或包，或函数。cat2cat1ggsignifggpvalggpubr::stat_compare_means()

要求：

  • 几何图形应为 .pointrange
  • cat1应位于 x 轴上。
  • cat2应该定义颜色美学。
  • 我想使用该包进行模型拟合（以获得正确的 p 值）。lmerTest

挑战：

  • 对于大型数据集（N=46,000，G1=4,000，G2=40），使用需要很长时间。ggeffects::hypothesis_test()
  • 我为拥有括号的解决方案而苦苦挣扎（意义星星可以是带有的地方，我猜）ggplotgeom_text()ifelse()

可重现的例子：

下面是使用玩具数据集的简单 R 示例：

# Load required packages
library(lmerTest)
library(ggplot2)

# Create toy data
set.seed(123)
N <- 1000  # Increased sample size
G1 <- 20   # Increased number of levels for G1
G2 <- 10   # Increased number of levels for G2

df <- data.frame(
  outcome = rnorm(N),
  cat1 = factor(rep(1:2, each = N / 2)),
  cat2 = factor(rep(c("Control", "T1", "T2", "T3"), each = N / 4)),
  G1 = factor(rep(1:G1, each = N / G1)),
  G2 = factor(rep(1:G2, each = N / G2))
) 

# Fit the hierarchical model
fit <- lmerTest::lmer(outcome ~ cat1 * cat2 + (1|G1) + (1|G2), data = df)
summary(fit)

# Attempt to make the plot
interactions::cat_plot(model = fit, pred = cat1, modx = cat2, 
                       interval.geom = "linerange", colors = "Set1")
# ... (I want significance stars and brackets connecting the coefficients: this is where I am stuck)

Categorical by categorical interaction without significant differences

R GGPLOT2 LME4 边际效应

评论

答：

0赞 Robert Long 10/5/2023 #1

首先，请注意，在拟合模型时，数据生成存在两个问题

固定效应模型矩阵有秩缺陷，因此删除了 4 列/系数

边界（奇异）拟合：参见 help（'isSingular'）

第一个问题的出现是因为预测变量中存在完美的多重共线性。这意味着可以使用模型中的其他预测变量来完美预测一个或多个预测变量。发生这种情况时，预测变量矩阵（通常称为设计矩阵）无法反转，并且无法唯一估计模型系数。

第二种可能是由于其中一种随机效应的变化太小。为了这个答案，我将省略 的随机效应。所以，希望我们可以通过以下代码实现你想要的：G2

library(lmerTest)
library(ggplot2)
library(effects)

# Create toy data
set.seed(123)
N <- 1000
G1 <- 20
G2 <- 10

df <- data.frame(
  outcome = rnorm(N),
  cat1 = factor(rep(1:2, each = N / 2)),
  cat2 = factor(sample(c("Control", "T1", "T2", "T3"), N, replace = TRUE)), # Shuffle the cat2 levels randomly
  G1 = factor(rep(1:G1, each = N / G1)),
  G2 = factor(rep(1:G2, each = N / G2))
)

# Fit the hierarchical model with interaction
fit <- lmerTest::lmer(outcome ~ cat1 * cat2 + (1|G1) + (1|G2), data = df)

# Extract effects of interaction
eff <- effect("cat1:cat2", fit)

# Convert effect object to dataframe for plotting
eff_df <- as.data.frame(eff)

# Rename columns for ease of use
names(eff_df) <- c("cat1", "cat2", "fit", "lower", "upper")

# Visualization using ggplot2
ggplot(eff_df, aes(x = cat1, y = fit, color = cat2, group = cat2)) +
  geom_point(aes(shape = cat2), size = 3) +
  geom_line(aes(linetype = cat2)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  labs(y = "Predicted Outcome", title = "Interaction Effect of cat1 and cat2") +
  theme_minimal()

生成以下图：

enter image description here

评论

0赞 CompSocialSciR 10/5/2023
谢谢。数据生成只是为了产生一个可重复的例子。您的可视化效果确实显示了这些更改，但未显示其重要性。在上面的示例中，试图添加带有星号的括号，以突出显示，例如，当 cat1==1 时，从 Cat2==“Control” 到 Cat2==“T1” 是否有统计学显着性增加。例如，ggsignif intro：cran.r-project.org/web/packages/ggsignif/vignettes/intro.html
0赞 Robert Long 10/5/2023
尽量不要担心统计学意义。p 值不是很有用。恕我直言，效果大小更为重要。假设您的效应为 -5，p 值为 0.051，这在 5% 水平上并不显著。将其与效应大小为 -0.1 且 p 值为 0.0499 进行比较，这将是显著的。假设这是一项关于降低 BMI> 30 的人的 BMI 的干预措施的研究。减少 -5 非常有意义，而 -0.1 则完全无关紧要。
0赞 CompSocialSciR 10/5/2023
我从事科学研究工作，我们关心实质性和统计学意义。你不应该排除大的影响只是由于偶然。这在功效不足的研究中很常见。
0赞 Robert Long 10/5/2023
事实上，这是确保样本量计算具有足够功效的充分理由。无论如何，我稍后会看看我是否可以纳入我的答案。ggsignif

上一个：如何绘制交互作用项标准差变化的边际效应？

下一个：具有估计边际均值的障碍模型的解释