如何按单株树木分组应用GAM函数分析树木生长阶段

数据背景 #

拥有树木生长数据库，包含10年共60株树木（每年6株）的日序（DOY）数据，记录三个生长阶段（Enlarging、Thickening、Mature）的细胞数量，简化数据结构如下：

Year Tree DOY Enlarging Thickening Mature
  2012  25  80    0         0          0
  2012  25  87    1         0          0
  2012  25  94    4         0          0
  2012  25  103   5         1          0
  2012  25  111   3         3          0
  2012  25  119   1         4          1
  2012  25  127   1         5          3
  2012  30  80    0         0          0
  2012  30  87    2         0          0
  ...

需求：按单株树木-年度分组拟合GAM模型，获取各生长阶段起止时间（Enlarging起始：细胞数>1；结束：细胞数<0），并绘制生长曲线。

问题重现 #

使用dplyr::do()分组建模后，提取拟合值和调用predict()时出现报错：

• 提取拟合值时因各组数据行数不一致，无法合并为数据框
• predict()无法直接作用于存储模型的数据框

解决方案 #

基于tidyverse（dplyr + purrr + ggplot2）实现分组建模、结果提取与可视化：

1. 加载依赖包 #

library(tidyverse)
library(mgcv)

2. 分组嵌套并拟合GAM模型 #

用nest()将每组数据打包，再用map()为每组拟合三个阶段的GAM模型：

# 按Tree和Year分组嵌套数据
nested_df <- df %>%
  group_by(Tree, Year) %>%
  nest() %>%
  ungroup()

# 为每组拟合三个生长阶段的GAM模型
model_df <- nested_df %>%
  mutate(
    gam_enlarging = map(data, ~ gam(Enlarging ~ s(DOY), data = ., 
                                    family = quasipoisson, gamma = 1, min.sp = 0.01)),
    gam_thickening = map(data, ~ gam(Thickening ~ s(DOY), data = ., 
                                     family = quasipoisson, gamma = 1, min.sp = 0.01)),
    gam_mature = map(data, ~ gam(Mature ~ s(DOY), data = ., 
                                 family = quasipoisson, gamma = 1, min.sp = 0.01))
  )

3. 合并拟合值到原始数据 #

用map2()将每组模型的拟合值与对应原始数据绑定，避免行数不匹配问题：

# 为每个阶段添加拟合值，并合并所有数据
df_with_fitted <- model_df %>%
  # 为Enlarging阶段添加拟合值
  mutate(data_enlarging = map2(data, gam_enlarging, ~ mutate(.x, fitted_enlarging = fitted(.y)))) %>%
  # 为Thickening阶段添加拟合值
  mutate(data_thickening = map2(data, gam_thickening, ~ mutate(.x, fitted_thickening = fitted(.y)))) %>%
  # 为Mature阶段添加拟合值
  mutate(data_mature = map2(data, gam_mature, ~ mutate(.x, fitted_mature = fitted(.y)))) %>%
  # 合并所有阶段数据
  select(Tree, Year, data_enlarging) %>%
  unnest(data_enlarging) %>%
  left_join(unnest(select(model_df, Tree, Year, data_thickening)), 
            by = c("Tree", "Year", "DOY")) %>%
  left_join(unnest(select(model_df, Tree, Year, data_mature)), 
            by = c("Tree", "Year", "DOY"))

4. 计算生长阶段起止时间 #

按每组筛选符合条件的DOY，确定各阶段起止（注：实际细胞数不会为负，可根据需求调整结束阈值，如<0.5）：

stage_dates <- df_with_fitted %>%
  group_by(Tree, Year) %>%
  summarise(
    # Enlarging阶段起止
    enlarging_start = min(DOY[fitted_enlarging > 1], na.rm = TRUE),
    enlarging_end = max(DOY[fitted_enlarging < 0], na.rm = TRUE),
    # Thickening阶段起止
    thickening_start = min(DOY[fitted_thickening > 0], na.rm = TRUE),
    thickening_end = max(DOY[fitted_thickening > 0], na.rm = TRUE),
    # Mature阶段起止
    mature_start = min(DOY[fitted_mature > 0], na.rm = TRUE),
    mature_end = max(DOY[fitted_mature > 0], na.rm = TRUE)
  ) %>%
  ungroup() %>%
  # 替换无匹配值的NA
  mutate(across(c(enlarging_start:mature_end), ~ replace_na(., NA)))

5. 绘制生长曲线 #

用ggplot2绘制每组的原始数据点与拟合曲线：

# Enlarging阶段曲线
ggplot(df_with_fitted, aes(x = DOY)) +
  geom_point(aes(y = Enlarging), alpha = 0.5) +
  geom_line(aes(y = fitted_enlarging), color = "#2c7fb8") +
  facet_wrap(~ Tree + Year, scales = "free_x") +
  labs(title = "Enlarging Stage Growth Curves",
       x = "Day of Year (DOY)",
       y = "Cell Count") +
  theme_bw()

# Thickening阶段曲线
ggplot(df_with_fitted, aes(x = DOY)) +
  geom_point(aes(y = Thickening), alpha = 0.5) +
  geom_line(aes(y = fitted_thickening), color = "#7fbf7b") +
  facet_wrap(~ Tree + Year, scales = "free_x") +
  labs(title = "Thickening Stage Growth Curves",
       x = "Day of Year (DOY)",
       y = "Cell Count") +
  theme_bw()

# Mature阶段曲线
ggplot(df_with_fitted, aes(x = DOY)) +
  geom_point(aes(y = Mature), alpha = 0.5) +
  geom_line(aes(y = fitted_mature), color = "#d73027") +
  facet_wrap(~ Tree + Year, scales = "free_x") +
  labs(title = "Mature Stage Growth Curves",
       x = "Day of Year (DOY)",
       y = "Cell Count") +
  theme_bw()

报错原因说明 #

1. 拟合值提取报错：各组数据的DOY观测数量不一致，直接合并所有拟合值会导致行数不匹配。使用nest()+map2()绑定每组数据与对应拟合值，可避免此问题。
2. predict()报错：predict()仅能作用于单个GAM模型对象，无法直接处理存储模型的数据框。通过map2()为每组模型传入对应数据，可实现批量预测。

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