如何从R语言ggplot的中位曲线获取代数方程？

你用stat_summary(fun.y = "median")生成的中位线，本质是每个x值对应的y中位数点连接成的折线，它本身没有预设的代数方程（比如线性、多项式），而是离散点的拟合结果。要得到可复用的表达式，需要对这些离散的中位线数据点进行模型拟合，以下是具体方案：

1. 正确提取中位线的原始数据 #

你之前的提取逻辑可以优化，通过特征筛选确保拿到的是stat_summary生成的中位线数据：

# 先把你的ggplot对象保存为gg
gg <- df %>% ggplot(aes(
  x = flowDayShifted,
  y = flowScaled,
  group = interaction(clusterFourier, as.factor(peak_number)),
  colour = as.factor(year)
)) + geom_line() +
  facet_wrap(~ clusterFourier,  ncol = 3) +
  stat_summary(
    fun.y = "median",
    aes(group = clusterFourier),
    color = "black",
    geom = "line",
    size = 1.2
  )

# 提取黑色中位线的图层数据，按聚类分组
plot_data <- ggplot_build(gg)$data
median_data <- plot_data[[which(sapply(plot_data, function(x) "group" %in% colnames(x) && any(x$color == "black")))]]
# 按cluster拆分数据，方便后续分别处理
cluster_medians <- split(median_data, median_data$group)

2. 拟合模型并导出表达式 #

洪水过程线通常是非线性的，根据需求可以选择两种拟合方式：

方式一：多项式回归（得到明确的代数方程） #

根据曲线的复杂度选择合适的多项式阶数（比如3阶、5阶），拟合后能直接写出多项式方程：

# 以第一个聚类为例，拟合5阶多项式
cluster1_poly <- lm(y ~ poly(x, 5, raw = TRUE), data = cluster_medians[[1]])
# 查看方程系数
summary(cluster1_poly)

输出结果中的系数可以直接整理为多项式方程：y = a₀ + a₁x + a₂x² + ... + a₅x⁵

方式二：loess局部回归（保留曲线特性，生成预测函数） #

如果不需要简洁的显式方程，只需要能根据x计算对应y值，loess更适合拟合非线性的洪水过程线：

# 拟合loess模型
cluster1_loess <- loess(y ~ x, data = cluster_medians[[1]], span = 0.75)
# 生成可直接调用的预测函数
predict_median <- function(x_val) {
  predict(cluster1_loess, newdata = data.frame(x = x_val))
}
# 示例：预测x=5时的中位y值
predict_median(5)

3. 验证拟合效果 #

可以把拟合曲线和原中位线画在一起，确认拟合精度：

# 生成拟合用的x序列
x_seq <- seq(min(cluster_medians[[1]]$x), max(cluster_medians[[1]]$x), length.out = 100)
# 用多项式模型预测y值
y_fit <- predict(cluster1_poly, newdata = data.frame(x = x_seq))

# 绘图对比原中位线（黑色）和拟合曲线（红色虚线）
ggplot() +
  geom_line(data = cluster_medians[[1]], aes(x, y), color = "black", size = 1.2) +
  geom_line(aes(x_seq, y_fit), color = "red", linetype = "dashed")

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