按树木与年份分组应用GAM模型预测树木生长阶段

核心思路 #

通过分组嵌套+自定义处理函数的方式，对Year和Tree的每个组合单独拟合GAM模型并生成预测结果，替代全局模型的拟合方式。以下是基于mgcv、dplyr和purrr的完整实现代码：

步骤1：加载依赖包 #

library(mgcv)
library(dplyr)
library(purrr)

步骤2：定义分组处理函数 #

该函数接收单个分组的数据（某一年某一棵树），完成模型拟合、预测结果生成与整合：

fit_gam_by_group <- function(group_data) {
  # 提取当前分组的年份和树ID
  current_year <- unique(group_data$Year)
  current_tree <- unique(group_data$Tree)
  
  # 为当前分组生成DOY序列（可选择全局范围或分组内范围，这里用分组内的DOY区间）
  doy_seq <- seq(min(group_data$DOY), max(group_data$DOY), 1)
  new_data <- data.frame(DOY = doy_seq)
  
  # 拟合三个生长阶段的GAM模型
  enlarging_mod <- gam(Enlarging ~ s(DOY), data = group_data, 
                       family = quasipoisson, gamma = 1, min.sp = 0.01)
  thickening_mod <- gam(Thickening ~ s(DOY), data = group_data, 
                        family = quasipoisson, gamma = 1, min.sp = 0.01)
  mature_mod <- gam(Mature ~ s(DOY), data = group_data, 
                    family = quasipoisson, gamma = 1, min.sp = 0.01)
  
  # 生成预测结果
  pred_enlarging <- predict(enlarging_mod, newdata = new_data, type = "response")
  pred_thickening <- predict(thickening_mod, newdata = new_data, type = "response")
  pred_mature <- predict(mature_mod, newdata = new_data, type = "response")
  
  # 整合结果并添加分组标识
  result_df <- data.frame(
    Year = current_year,
    Tree = current_tree,
    DOY = doy_seq,
    enlarging_gam = pred_enlarging,
    thickening_gam = pred_thickening,
    mature_gam = pred_mature
  )
  
  return(result_df)
}

步骤3：执行分组预测 #

使用dplyr的分组嵌套功能，对每个Year-Tree组合应用上述函数：

# 按Year和Tree分组并嵌套数据
grouped_data <- xylo %>%
  group_by(Year, Tree) %>%
  nest()

# 对每个分组应用函数并整合结果
final_pred <- grouped_data %>%
  mutate(prediction = map(data, fit_gam_by_group)) %>%
  select(-data) %>%
  unnest(prediction)

关键说明 #

• DOY序列范围：上述代码使用单个分组内的DOY区间生成预测序列，若需要统一使用全局DOY范围（即所有样本的最小到最大DOY），只需将函数内的doy_seq替换为全局序列：

  # 全局DOY序列（提前在函数外定义）
  global_doy <- seq(min(xylo$DOY), max(xylo$DOY), 1)
  # 在函数内使用
  new_data <- data.frame(DOY = global_doy)
  

• 模型稳定性：若部分分组的样本量过小导致GAM拟合报错，可在函数内添加判断逻辑（如if(nrow(group_data) < 5) return(NULL)），跳过无效分组。
• 结果验证：可通过filter(final_pred, Year == 2012, Tree == 25)查看单棵树单年份的预测结果，确认与手动拟合的一致性。

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