如何修改GetExtrema函数，筛选满足特定导数条件的曲线极值？

原函数GetExtrema通过二阶差分符号变化识别局部极值，但无法筛选出满足"极值点前后0.25年导数符号符合趋势"的有效极值。以下是修改方案，通过计算一阶导数并验证区间内的导数符号，精准筛选目标极值点。

实现思路 #

1. 用相邻预测点的斜率近似计算曲线的一阶导数；
2. 对每个候选极值点，锁定其前0.25年和后0.25年的年龄范围；
3. 验证区间内导数符号是否符合要求：
   • 有效最大值：前区间导数全部为正，后区间导数全部为负；
   • 有效最小值：前区间导数全部为负，后区间导数全部为正；
4. 仅保留满足条件的极值点（若数据存在少量噪声，可将"全部符合"改为"多数符合"的宽松判断）。

修改后的完整代码 #

# 加载包
library(tidyverse)
library(ggeffects)
library(mgcv)

# 读取数据（请确保bmi_long.csv文件在工作目录中）
dat <- read.csv("bmi_long.csv") %>% 
  select(age, bmi)

# 拟合GAM模型
model <- gam(
  bmi ~ s(sqrt(age), bs = 'ps'), 
  method = "REML", data = dat
)

# 生成预测数据
pred <- predict_response(model, terms = "age [all]") %>% 
  as.data.frame() %>% 
  select(x, predicted) %>% 
  rename(y = predicted)

# 绘制预测曲线
with(pred, plot(x, y, type = "l"))

# 定义修改后的GetExtrema函数
GetExtrema <- function(x, y, window = 0.25) {
  # 计算一阶导数（相邻点斜率近似）
  dy_dx <- diff(y) / diff(x)
  # 导数对应的年龄为相邻年龄的中点
  dx_x <- (x[-1] + x[-length(x)]) / 2
  
  # 识别初始局部极值点索引
  turn <- diff(sign(diff(y)))
  max_candidates <- which(turn == -2) + 1L
  min_candidates <- which(turn == 2) + 1L
  
  # 筛选有效最大值
  valid_max <- sapply(max_candidates, function(ind) {
    x0 <- x[ind]
    # 截取极值点前后0.25年范围内的导数
    pre_dx <- dy_dx[dx_x > (x0 - window) & dx_x < x0]
    post_dx <- dy_dx[dx_x > x0 & dx_x < (x0 + window)]
    
    # 验证导数符号趋势（若需宽松判断，可改用mean(sign(pre_dx)) > 0）
    all(sign(pre_dx) == 1) && all(sign(post_dx) == -1)
  })
  
  # 筛选有效最小值
  valid_min <- sapply(min_candidates, function(ind) {
    x0 <- x[ind]
    pre_dx <- dy_dx[dx_x > (x0 - window) & dx_x < x0]
    post_dx <- dy_dx[dx_x > x0 & dx_x < (x0 + window)]
    
    all(sign(pre_dx) == -1) && all(sign(post_dx) == 1)
  })
  
  # 提取最终有效极值点
  list(
    min.x = x[min_candidates[valid_min]],
    min.y = y[min_candidates[valid_min]],
    max.x = x[max_candidates[valid_max]],
    max.y = y[max_candidates[valid_max]]
  )
}

# 调用函数
GetExtrema(pred$x, pred$y)

代码说明 #

• dy_dx：用相邻预测点的斜率近似一阶导数，dx_x对应导数的中间年龄；
• 极值筛选逻辑：通过截取目标区间内的导数，严格验证符号趋势，避免局部小波动干扰；
• window参数：默认值为0.25年，可根据需求灵活调整时间窗口大小。

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