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

Question

我正在使用一个分层模型（lmerMod），它有一个结果和两个分类变量，

cat1

和

cat2

。我想使用

ggplot2::

制作交互图，可视化模型系数并考虑两个特定的复杂性：

成对比较：我需要在

cat2

（1 和 2）的每个级别内相互比较

cat1

（控制、T1、T2、T3）每个级别的系数。 统计显着性：我想使用星号或一些类似的符号来强调这些比较的统计显着性。例如

ggsignif

或

ggpval

包，或

ggpubr::stat_compare_means()

函数。

要求：

• 几何形状应为

  pointrange

  。

• cat1

  应位于 x 轴上。

• cat2

  应定义色彩美感。
• 我想使用

  lmerTest

  包进行模型拟合（以获得正确的 p 值）。

挑战：

• 对于大型数据集（N=46,000，G1=4,000，G2=40），使用

  ggeffects::hypothesis_test()

  需要很长时间。
• 我很难找到带有括号的

  ggplot

  解决方案（我猜重要的星星可以与

  geom_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)

Answer 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()

生成以下图：