如何用R语言ggplot绘制含公式标签的多元线性回归与交互效应模型图

核心思路 #

1. 连续变量Z无法直接作为分组维度生成多条回归线，需要先将Z转换为若干有代表性的固定取值（通常取均值、均值±1标准差三个水平），基于这些取值拟合全局多元模型后生成预测值绘制线条
2. 多元回归的公式和R²可以直接从已拟合的lm对象中提取，手动构造标签添加到图表中

前置准备 #

先加载所需依赖包，导入你的数据集：

library(ggplot2)
library(ggeffects)
library(dplyr)

# 导入数据集
hw_data <- structure(list(X = c(18, 19, 20, 17, 8, 15, 14, 16, 18, 14, 16, 
13, 16, 17, 10, 18, 19, 25, 18, 13, 18, 16, 11, 17, 15, 18, 19, 
16, 20, 17, 8, 18, 15, 14, 18, 14, 16, 13, 16), Y = c(15, 13, 
14, 22, 2, 11, 15, 11, 20, 17, 20, 17, 20, 14, 21, 10, 13, 16, 
12, 11, 13, 10, 4, 16, 18, 15, 10, 13, 14, 17, 2, 11, 15, 11, 
20, 17, 14, 7, 16), Z = c(32, 42, 37, 34, 32, 39, 44, 49, 36, 
31, 36, 37, 37, 45, 46, 48, 36, 42, 36, 25, 36, 39, 26, 32, 33, 
38, 33, 44, 46, 34, 32, 39, 44, 49, 36, 31, 36, 37, 37)), class = "data.frame", row.names = c(NA, 
-39L))

# 你已有的模型拟合代码
hw_regression_simple <- lm(Y ~ X, data = hw_data)
hw_regression_two_factors <- lm(Y ~ X + Z, hw_data)
hw_regression_interaction <- lm(Y ~ X * Z, hw_data)

1. 交互效应模型（Y ~ X * Z）可视化 #

该模型下X对Y的影响斜率会随Z的取值变化，绘制三条对应Z不同水平的回归线：

# 生成Z的代表性取值：均值-1SD、均值、均值+1SD
z_mean <- mean(hw_data$Z)
z_sd <- sd(hw_data$Z)
z_levels <- c(z_mean - z_sd, z_mean, z_mean + z_sd)

# 生成模型预测值
pred_interaction <- ggpredict(hw_regression_interaction, terms = c("X [all]", paste0("Z [", paste(z_levels, collapse = ","), "]")))

# 绘图
ggplot() +
  geom_point(data = hw_data, aes(x = X, y = Y, color = Z)) +
  geom_line(data = pred_interaction, aes(x = x, y = predicted, color = group), linewidth = 1) +
  geom_ribbon(data = pred_interaction, aes(x = x, ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.1) +
  # 添加模型公式和R²
  annotate("text", x = min(hw_data$X), y = max(hw_data$Y), hjust = 0,
           label = paste0(
             "Y = ", round(coef(hw_regression_interaction)[1], 2),
             " + ", round(coef(hw_regression_interaction)[2], 2), "X",
             ifelse(coef(hw_regression_interaction)[3] >= 0, " + ", " - "), abs(round(coef(hw_regression_interaction)[3], 2)), "Z",
             ifelse(coef(hw_regression_interaction)[4] >= 0, " + ", " - "), abs(round(coef(hw_regression_interaction)[4], 2)), "X×Z\n",
             "R² = ", round(summary(hw_regression_interaction)$r.squared, 3)
           )) +
  labs(title = "X、Z与Y的交互效应", x = "X", y = "Y") +
  scale_color_viridis_c() +
  scale_fill_viridis_c(guide = "none") +
  theme_minimal()

2. 主效应模型（Y ~ X + Z）可视化 #

该模型下X对Y的影响斜率固定，不同Z水平的回归线为平行关系：

# 生成模型预测值
pred_two_factor <- ggpredict(hw_regression_two_factors, terms = c("X [all]", paste0("Z [", paste(z_levels, collapse = ","), "]")))

# 绘图
ggplot() +
  geom_point(data = hw_data, aes(x = X, y = Y, color = Z)) +
  geom_line(data = pred_two_factor, aes(x = x, y = predicted, color = group), linewidth = 1) +
  geom_ribbon(data = pred_two_factor, aes(x = x, ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.1) +
  # 添加模型公式和R²
  annotate("text", x = min(hw_data$X), y = max(hw_data$Y), hjust = 0,
           label = paste0(
             "Y = ", round(coef(hw_regression_two_factors)[1], 2),
             " + ", round(coef(hw_regression_two_factors)[2], 2), "X",
             ifelse(coef(hw_regression_two_factors)[3] >= 0, " + ", " - "), abs(round(coef(hw_regression_two_factors)[3], 2)), "Z\n",
             "R² = ", round(summary(hw_regression_two_factors)$r.squared, 3)
           )) +
  labs(title = "X、Z与Y的主效应", x = "X", y = "Y") +
  scale_color_viridis_c() +
  scale_fill_viridis_c(guide = "none") +
  theme_minimal()

可选方案：按Z分组单独拟合回归 #

如果不需要全局多元回归的结果，只需要按Z的分组分别拟合简单线性回归，可以先将Z离散化为分组变量，代码如下：

# 将Z按三分位离散化为三组
hw_data <- hw_data %>%
  mutate(Z分组 = cut(Z, breaks = quantile(Z, c(0, 0.33, 0.66, 1)),
                     labels = c("低Z", "中Z", "高Z"), include.lowest = TRUE))

ggplot(hw_data, aes(x = X, y = Y, color = Z分组)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~~~")),
               formula = y ~ x) +
  labs(title = "按Z分组的简单线性回归", x = "X", y = "Y") +
  theme_minimal()

内容的提问来源于stack exchange，提问作者Shawn Hemelstrand