如何在R中使用GAM方法对多维时间序列进行三步插值

核心思路 #

针对每个geom_id的时间序列，用GAM的时间平滑项拟合变量趋势，再预测填补缺失值。无需局限于“三步时间步长”的局部插值，GAM的全局平滑能更好捕捉时间序列的整体趋势；若需局部约束，可调整平滑项参数实现。

步骤1：数据预处理 #

先确保时间列格式正确并按分组排序：

library(mgcv)
library(dplyr)
library(lubridate)

# 转换时间格式（根据你的实际日期格式调整，比如ymd_hms）
df <- df %>%
  mutate(acqdate = ymd(acqdate)) %>%
  arrange(geom_id, acqdate) # 按分组和时间排序

步骤2：单变量插值示例（以b02为例） #

按geom_id分组拟合GAM，预测填补NA：

df_interp <- df %>%
  group_by(geom_id) %>%
  do({
    # 拟合GAM：用三次回归样条平滑时间变量
    gam_model <- gam(b02 ~ s(acqdate, bs = "cr"), data = .)
    # 预测所有行（包括含NA的行）
    .$b02_filled <- predict(gam_model, newdata = .)
    .
  }) %>%
  ungroup() %>%
  # 用预测值替换原NA
  mutate(b02 = ifelse(is.na(b02), b02_filled, b02)) %>%
  select(-b02_filled)

步骤3：批量处理多变量（b02至b06） #

用循环批量完成所有目标变量的插值：

# 定义需要插值的变量列表
target_vars <- c("b02", "b03", "b04", "b05", "b06")

df_final <- df
for (var in target_vars) {
  df_final <- df_final %>%
    group_by(geom_id) %>%
    do({
      dat <- .
      # 动态生成公式，引用当前变量
      gam_formula <- formula(paste(var, "~ s(acqdate, bs = 'cr')"))
      gam_model <- gam(gam_formula, data = dat)
      # 生成预测列名
      pred_col <- paste0(var, "_filled")
      dat[[pred_col]] <- predict(gam_model, newdata = dat)
      dat
    }) %>%
    ungroup() %>%
    # 替换NA并删除临时列
    mutate(!!sym(var) := ifelse(is.na(!!sym(var)), !!sym(pred_col), !!sym(var))) %>%
    select(-all_of(pred_col))
}

关键参数说明 #

• s(acqdate, bs = "cr")：使用三次回归样条拟合时间趋势，适合常规时间序列；若需更灵活的拟合，可换用bs = "tp"（薄板样条）或bs = "ps"（P样条）。
• 若需加入空间效应（经纬度），可扩展公式为：b02 ~ s(acqdate, bs = "cr") + s(lon, lat, bs = "tp")，同时拟合时间和空间趋势。
• 调整平滑项自由度k：比如s(acqdate, bs = "cr", k = 8)，k越小拟合越平滑，越大越贴近数据，可通过gam.check(gam_model)验证拟合效果。

关于gam_interp()的使用 #

gam_interp()主要用于对已拟合的GAM模型在新的时间点（更细的时间粒度）插值，而非填补现有时间步的NA。如果需要生成三步时间步长的新序列，示例：

# 针对单个geom_id的模型
sample_dat <- df %>% filter(geom_id == "xxx")
gam_model <- gam(b02 ~ s(acqdate, bs = "cr"), data = sample_dat)

# 生成间隔为3天的新时间点（根据你的"三步"定义调整）
new_dates <- seq(min(sample_dat$acqdate), max(sample_dat$acqdate), by = "3 days")
# 插值预测
interp_results <- gam_interp(gam_model, newdata = data.frame(acqdate = new_dates))

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