如何为facet_wrap后的ggplot2分面图匹配对应月份的二次拟合函数

我有两个数据框：

• 一个记录1880-2023年各年份每月的全球平均气温；
• 另一个是tidy格式的二次模型参数，包含月份列、截距列，以及beta（一次项）、beta平方（二次项）两个模型估计列。

我需要在按月份使用facet_wrap分面的ggplot2图表中，为每个月份的气温数据叠加对应月份的二次函数曲线。目前已完成数据加载、清洗和二次模型拟合，但geom_function无法为每个分面匹配对应月份的模型参数，求解决办法。

相关代码如下：

加载数据 #

library(tidytuesdayR)
library(tidyverse)

tuesdata <- tidytuesdayR::tt_load(2023, week = 28)
global_temps <- tuesdata$global_temps

数据清洗 #

# 清洗数据
global_temps <- global_temps |> 
  janitor::clean_names()

# 重塑数据
global_temps <- global_temps |> 
    pivot_longer(
        c(jan:dec),
        names_to = "month",
        values_to = "temperature"
    ) |> 
    select(year, month, temperature) |> 
  mutate(month = str_to_title(month))

# 排序月份
month_order <- c(
  "Jan", "Feb", "Mar", "Apr", "May", "Jun", 
  "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"
)

global_temps <- global_temps |> 
  mutate(month = factor(month, levels = month_order))

二次模型拟合 #

# 为温度增长拟合二次项模型
global_temp_quad <- global_temps |> 
  group_by(month) |> 
  nest() |> 
  mutate(model = map(data, ~ lm(temperature ~ poly(year, 2, raw = TRUE), data = .x)))

global_temp_quad <- global_temp_quad |> 
  mutate(coef = map(model, broom::tidy)) |> 
  unnest(coef) |> 
  select(term, estimate)

# 宽表转换
global_temp_quad <- global_temp_quad |> 
  pivot_wider(
    names_from = term,
    values_from = estimate
  ) |> 
  rename(
    intercept = "(Intercept)",
    year = "poly(year, 2, raw = TRUE)1",
    year_2 = "poly(year, 2, raw = TRUE)2"
  ) |> 
  ungroup()

原绘图代码（存在匹配问题） #

ggplot() +
  geom_point(
    data = global_temps,
    aes(x = year, y = temperature), 
    show.legend = FALSE,
    size = 1.5,
    alpha = 1,
    shape = 1
  ) +
  theme_light() +
  labs(
    title = "Global surface temperatures (1890-2023)",
    subtitle = "Temperature is recorded as deviation from the mean temperature for 1890-1950",
    source = "Data from NASA/GISS",
    x = NULL,
    y = "Deviation from mean temperature (Celsius)"
  ) +
  facet_wrap(~ month) +
  geom_function(
    data = global_temp_quad,
    aes(group = interaction(month)),
    color = "red",
    alpha = 0.7,
    show.legend = FALSE,
    fun = function(x) {
      global_temp_quad$intercept +
      global_temp_quad$year * x +
      global_temp_quad$year_2 * x^2
    }
  ) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.ticks = element_blank()
  ) +
  scale_x_continuous(breaks = seq(1890, 2020, 40)) +
  scale_y_continuous(limits = c(-0.5, 1.5))

方法1：生成预测数据并用geom_line绘制（推荐） #

geom_function的fun参数无法直接关联分面的分组参数，最可靠的方式是先基于模型参数生成每个月份的预测数据，再用geom_line绘制曲线，这样能确保每个分面匹配对应月份的函数。

步骤1：生成预测数据 #

# 生成每个月份的预测数据（覆盖数据中的年份范围）
pred_data <- global_temp_quad |>
  rowwise() |>
  mutate(
    # 生成密集的年份序列，让曲线更平滑
    year_pred = list(seq(min(global_temps$year), max(global_temps$year), length.out = 100)),
    # 根据二次模型计算预测温度
    temp_pred = list(intercept + year * year_pred + year_2 * year_pred^2)
  ) |>
  unnest(c(year_pred, temp_pred))

步骤2：修改绘图代码 #

替换原代码中的geom_function部分为geom_line，加载预测数据：

ggplot() +
  geom_point(
    data = global_temps,
    aes(x = year, y = temperature), 
    show.legend = FALSE,
    size = 1.5,
    alpha = 1,
    shape = 1
  ) +
  # 添加对应月份的二次拟合曲线
  geom_line(
    data = pred_data,
    aes(x = year_pred, y = temp_pred),
    color = "red",
    alpha = 0.7
  ) +
  theme_light() +
  labs(
    title = "Global surface temperatures (1890-2023)",
    subtitle = "Temperature is recorded as deviation from the mean temperature for 1890-1950",
    source = "Data from NASA/GISS",
    x = NULL,
    y = "Deviation from mean temperature (Celsius)"
  ) +
  facet_wrap(~ month) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.ticks = element_blank()
  ) +
  scale_x_continuous(breaks = seq(1890, 2020, 40)) +
  scale_y_continuous(limits = c(-0.5, 1.5))

方法2：使用stat_function传递分组参数 #

如果坚持使用函数类图层，可以用stat_function结合args传递每个月份的模型参数，通过group关联分面：

ggplot() +
  geom_point(
    data = global_temps,
    aes(x = year, y = temperature), 
    show.legend = FALSE,
    size = 1.5,
    alpha = 1,
    shape = 1
  ) +
  theme_light() +
  labs(
    title = "Global surface temperatures (1890-2023)",
    subtitle = "Temperature is recorded as deviation from the mean temperature for 1890-1950",
    source = "Data from NASA/GISS",
    x = NULL,
    y = "Deviation from mean temperature (Celsius)"
  ) +
  facet_wrap(~ month) +
  stat_function(
    data = global_temp_quad,
    aes(group = month),
    color = "red",
    alpha = 0.7,
    # 定义带参数的二次函数
    fun = function(x, a, b, c) {
      a + b*x + c*x^2
    },
    # 传递每个月份的模型参数
    args = list(
      a = global_temp_quad$intercept,
      b = global_temp_quad$year,
      c = global_temp_quad$year_2
    )
  ) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.ticks = element_blank()
  ) +
  scale_x_continuous(breaks = seq(1890, 2020, 40)) +
  scale_y_continuous(limits = c(-0.5, 1.5))

内容的提问来源于stack exchange，提问作者Simen Løkken