如何用facet_wrap为Obs数据集添加两条含/不含异常值的回归线？

我有一组按处理组拆分的图表（每个处理对应一个facet），每个图表包含观测数据（Obs）和模型输出数据（Model）。观测数据中存在疑似无效的异常值，因此需要为每个facet中的Obs数据集添加两条回归线：一条包含所有数据点，一条剔除异常值；Model数据集无需添加回归线。由于不同站点的处理组数量不固定，希望找到通用的适配方案。

示例代码与当前效果 #

df <- data.frame(year=factor(c(2000,2001,2002,2003,2004,2005,
                               2000,2001,2002,2003,2004,2005,
                               2000,2001,2002,2003,2004,2005,
                               2000,2001,2002,2003,2004,2005)),
                 treatment_code=c("T1","T1","T1","T1","T1","T1",
                                  "T1","T1","T1","T1","T1","T1",
                                  "T2","T2","T2","T2","T2","T2",
                                  "T2","T2","T2","T2","T2","T2"),
                 value=c(9,10,8.5,7.5,22,10.5,
                         11,9,12,9,10,11.5,
                         8,11,9.5,12,10,10.5,
                         7,9,11,10,12,11.5),
                 model=c("Obs","Obs","Obs","Obs","Obs","Obs",
                         "A","A","A","A","A","A",
                         "Obs","Obs","Obs","Obs","Obs","Obs",
                         "A","A","A","A","A","A"))

gS_calib <- df %>%
  ggplot(aes(x=year, y=value, color=model, show.legend=TRUE)) +
  geom_point(show.legend=TRUE) +
  xlab("Year") +
  facet_wrap(~treatment_code, ncol=length(unique(df$treatment_code))) 

gS_calib

当前图表效果：

完全可以实现，核心是先预处理生成Obs数据的两种回归拟合结果，再结合facet_wrap自动适配不同处理组数量，具体步骤如下：

1. 编写通用异常值剔除函数 #

这里用常用的IQR法（四分位距）识别并剔除异常值，你也可以根据需求换成标准差法或自定义阈值：

remove_outliers <- function(data, group_col, value_col) {
  data %>%
    group_by({{group_col}}) %>%
    # 标记异常值
    mutate(is_outlier = ifelse({{value_col}} < quantile({{value_col}}, 0.25) - 1.5*IQR({{value_col}}) |
                                {{value_col}} > quantile({{value_col}}, 0.75) + 1.5*IQR({{value_col}}),
                              TRUE, FALSE)) %>%
    # 剔除异常值
    filter(!is_outlier) %>%
    ungroup()
}

2. 生成两种回归拟合数据 #

分别对所有Obs数据和剔除异常值后的Obs数据做分组回归，用broom包生成拟合值：

library(broom)
library(dplyr)
library(ggplot2)

# 提取仅Obs的数据
obs_data <- df %>% filter(model == "Obs")

# 全数据的回归拟合（包含异常值）
fit_full <- obs_data %>%
  group_by(treatment_code) %>%
  do(augment(lm(value ~ as.numeric(year), data = .))) %>%
  mutate(fit_type = "全数据回归线")

# 剔除异常值后的回归拟合
obs_no_outliers <- remove_outliers(obs_data, treatment_code, value)
fit_clean <- obs_no_outliers %>%
  group_by(treatment_code) %>%
  do(augment(lm(value ~ as.numeric(year), data = .))) %>%
  mutate(fit_type = "剔除异常值回归线")

# 合并两种拟合结果
fit_data <- bind_rows(fit_full, fit_clean)

3. 绘制最终图表 #

把原始点和两条回归线整合到分面中，facet_wrap会自动根据处理组数量调整布局：

gS_calib <- df %>%
  ggplot(aes(x = year, y = value, color = model)) +
  # 绘制所有原始数据点
  geom_point(show.legend = TRUE) +
  # 添加全数据回归线（仅针对Obs）
  geom_line(data = fit_data %>% filter(fit_type == "全数据回归线"),
            aes(y = .fitted, color = fit_type), linetype = "solid") +
  # 添加剔除异常值后的回归线（仅针对Obs）
  geom_line(data = fit_data %>% filter(fit_type == "剔除异常值回归线"),
            aes(y = .fitted, color = fit_type), linetype = "dashed") +
  xlab("年份") +
  # 自动适配处理组数量的分面
  facet_wrap(~treatment_code) +
  # 自定义颜色区分不同数据类型
  scale_color_manual(values = c("Obs" = "black", "A" = "blue", 
                                "全数据回归线" = "red", "剔除异常值回归线" = "darkred")) +
  labs(color = "数据/拟合类型")

gS_calib

关键优势 #

• 完全通用：不管处理组数量多少，facet_wrap都会自动调整分面布局，无需手动指定列数；
• 灵活可改：异常值剔除逻辑、回归线样式（颜色、线型）都可以根据需求自定义；
• 分组独立拟合：每个分面的回归线都是基于该处理组内的Obs数据生成，结果准确。

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