如何在ggplot分面图中添加对应lmer模型的方程与R²值

核心思路是按分面变量（species）分组计算对应模型的回归方程和R²，再将结果关联到分面中，用geom_text/geom_label添加到对应子图。以下分两种常见场景给出实现方案：

场景1：基于全局混合模型提取分组参数 #

如果你已经构建了包含species与自变量交互项的全局混合模型，直接从模型中提取每个species的系数和R²即可：

步骤1：加载依赖包 #

library(lme4)
library(ggplot2)
library(dplyr)
library(MuMIn)  # 用于计算混合模型的R²
library(ggeffects)  # 用于生成模型预测值（可选）

步骤2：定义分组计算函数 #

该函数从全局模型中提取指定species的截距、斜率，以及边际/条件R²：

get_group_eqn <- function(data, global_model) {
  current_sp <- unique(data$species)
  coefs <- fixef(global_model)
  # 依据因子编码调整系数计算（这里假设treatment编码，参考水平为第一个物种）
  if (current_sp == levels(data$species)[1]) {
    intercept <- coefs["(Intercept)"]
    slope <- coefs["x"]  # 替换为你的自变量名
  } else {
    intercept <- coefs["(Intercept)"] + coefs[paste0("species", current_sp)]
    slope <- coefs["x"] + coefs[paste0("x:species", current_sp)]
  }
  # 计算混合模型的R²（边际=固定效应解释方差；条件=固定+随机效应解释方差）
  r2_vals <- r.squaredGLMM(global_model)
  # 格式化方程和R²文本
  eqn_text <- sprintf("y = %.2f + %.2fx", intercept, slope)
  r2_text <- sprintf("R²(边际) = %.3f\nR²(条件) = %.3f", r2_vals[1], r2_vals[2])
  # 返回分组结果
  tibble(species = current_sp, eqn = eqn_text, r2 = r2_text)
}

步骤3：生成分组统计数据 #

# 替换为你的全局模型和数据集
global_model <- lmer(y ~ x * species + (1|random_effect), data = your_data)

# 按species分组计算方程与R²
group_stats <- your_data %>%
  group_by(species) %>%
  summarize(get_group_eqn(cur_data(), global_model)) %>%
  ungroup()

步骤4：绘制分面图并添加文本 #

ggplot(your_data, aes(x = x, y = y)) +
  geom_point(alpha = 0.6) +
  # 添加全局模型的预测线与置信区间
  geom_line(data = ggpredict(global_model, terms = c("x", "species")), 
            aes(y = predicted), color = "#E64B35", linewidth = 1) +
  geom_ribbon(data = ggpredict(global_model, terms = c("x", "species")), 
              aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, fill = "#E64B35") +
  # 分面
  facet_wrap(~species) +
  # 添加分组的方程与R²（放在子图右上角）
  geom_label(data = group_stats, aes(x = Inf, y = Inf, label = paste(eqn, r2, sep = "\n")),
             hjust = 1, vjust = 1, inherit.aes = FALSE, size = 3.5) +
  theme_bw()

场景2：每个分面单独拟合混合模型 #

如果需要为每个species单独构建混合模型，调整分组函数即可：

步骤1：定义单组模型计算函数 #

get_single_sp_eqn <- function(data) {
  current_sp <- unique(data$species)
  # 为当前分组拟合混合模型
  single_model <- lmer(y ~ x + (1|random_effect), data = data)
  coefs <- fixef(single_model)
  r2_vals <- r.squaredGLMM(single_model)
  # 格式化文本
  eqn_text <- sprintf("y = %.2f + %.2fx", coefs["(Intercept)"], coefs["x"])
  r2_text <- sprintf("R²(边际) = %.3f\nR²(条件) = %.3f", r2_vals[1], r2_vals[2])
  tibble(species = current_sp, eqn = eqn_text, r2 = r2_text)
}

步骤2：生成分组统计数据 #

group_stats_single <- your_data %>%
  group_by(species) %>%
  summarize(get_single_sp_eqn(cur_data())) %>%
  ungroup()

步骤3：绘图（预测线需对应单组模型） #

如果要绘制单组模型的预测线，需要先生成每个分组的预测数据，再传入ggplot：

# 生成每个分组的预测数据
pred_data <- your_data %>%
  group_by(species) %>%
  group_modify(~{
    model <- lmer(y ~ x + (1|random_effect), data = .x)
    new_x <- seq(min(.x$x), max(.x$x), length.out = 100)
    pred <- predict(model, newdata = data.frame(x = new_x, random_effect = unique(.x$random_effect)[1]), interval = "confidence")
    tibble(x = new_x, predicted = pred[,1], conf.low = pred[,2], conf.high = pred[,3])
  }) %>%
  ungroup()

# 绘图
ggplot(your_data, aes(x = x, y = y)) +
  geom_point(alpha = 0.6) +
  geom_line(data = pred_data, aes(y = predicted), color = "#4DBBD5", linewidth = 1) +
  geom_ribbon(data = pred_data, aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, fill = "#4DBBD5") +
  facet_wrap(~species) +
  geom_label(data = group_stats_single, aes(x = Inf, y = Inf, label = paste(eqn, r2, sep = "\n")),
             hjust = 1, vjust = 1, inherit.aes = FALSE, size = 3.5) +
  theme_bw()

注意事项 #

• 因子编码会影响系数计算，若使用sum编码（如contr.sum），需调整get_group_eqn中的系数提取逻辑。
• 若只需展示R²的一种类型（如边际R²），直接删除对应文本即可。

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