在R中结合dplyr group_by与approxfun/approx实现生长阶段插值分析

问题背景 #

现有数据集记录了不同年份（Year）、不同树木（Tree）在部分年积日（DOY）下，Enlarging（膨大期）和Maturing（成熟期）两个细胞生长阶段的细胞数量。采样存在间隔：前一次采样细胞数为0，后一次已超过2。需要完成两个目标：

1. 对每个Year-Tree分组，生成覆盖全年1-365天的细胞数插值数据集
2. 计算每个Year-Tree分组下，两个生长阶段的起始DOY（细胞数>1）和结束DOY（细胞数<1）

数据集示例：

library(dplyr)

df <- data.frame(
  "Year" = c(2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012,
             2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013),
  "Tree" = c(15, 15, 15, 15, 15, 22, 22, 22, 22, 22, 41, 41, 41, 41, 41,
             53, 53, 53, 53, 53),
  "DOY" = c(65, 97, 125, 177, 214, 65, 97, 125, 177, 214, 61, 99, 118, 166, 221,
            61, 99, 118, 166, 221),
  "Enlarging" = c(0, 2, 4, 5, 0, 0, 3, 6, 3, 0, 0, 5, 4, 4, 0, 0, 4, 7, 5, 0),
  "Maturing" = c(0, 0, 3, 7, 0, 0, 0, 3, 4, 0, 0, 3, 6, 8, 0, 0, 0, 4, 7, 0)
) %>%
  mutate(Year = as.factor(Year),
         Tree = as.factor(Tree),
         DOY = as.numeric(DOY),
         Enlarging = as.numeric(Enlarging),
         Maturing = as.numeric(Maturing))

print(df)

解决方案1：生成全年插值数据集 #

使用dplyr的分组功能结合approxfun，对每个Year-Tree分组构建线性插值函数，预测全年1-365天的细胞数：

library(dplyr)

# 生成全年插值数据
interpolated_full <- df %>%
  group_by(Year, Tree) %>%
  summarise(
    # 生成全年1-365的DOY序列
    DOY = seq(1, 365, by = 1),
    # 对Enlarging阶段进行线性插值
    Enlarging_interp = approxfun(DOY, Enlarging, method = "linear")(DOY),
    # 对Maturing阶段进行线性插值
    Maturing_interp = approxfun(DOY, Maturing, method = "linear")(DOY),
    .groups = "drop"
  )

# 查看结果示例
head(interpolated_full)

说明：

• approxfun默认使用线性插值，适合处理采样点间的连续变化
• 若需要更平滑的插值，可将method参数改为"spline"或"constant"

解决方案2：直接计算各生长阶段的起止DOY #

无需生成全年数据，直接通过相邻采样点的线性插值，计算细胞数等于1时的DOY，作为阶段的起始/结束点：

library(dplyr)

# 定义函数：计算单个阶段的起止DOY
get_stage_dates <- function(doy_vec, count_vec) {
  # 找到所有细胞数从0到>1的区间（起始点）
  start_indices <- which(diff(count_vec > 1) == 1)
  start_doys <- sapply(start_indices, function(i) {
    # 线性插值计算count=1时的DOY
    approx(x = count_vec[c(i, i+1)], y = doy_vec[c(i, i+1)], xout = 1)$y
  })
  
  # 找到所有细胞数从>1到0的区间（结束点）
  end_indices <- which(diff(count_vec > 1) == -1)
  end_doys <- sapply(end_indices, function(i) {
    approx(x = count_vec[c(i, i+1)], y = doy_vec[c(i, i+1)], xout = 1)$y
  })
  
  # 返回起止DOY，若没有则返回NA
  data.frame(
    Start_DOY = ifelse(length(start_doys) > 0, start_doys, NA),
    End_DOY = ifelse(length(end_doys) > 0, end_doys, NA)
  )
}

# 按Year-Tree分组，计算两个阶段的起止日期
stage_dates <- df %>%
  group_by(Year, Tree) %>%
  summarise(
    # 处理Enlarging阶段
    bind_rows(
      get_stage_dates(DOY, Enlarging) %>% mutate(Stage = "Enlarging"),
      get_stage_dates(DOY, Maturing) %>% mutate(Stage = "Maturing")
    ),
    .groups = "drop"
  )

# 查看结果
print(stage_dates)

说明：

• 函数get_stage_dates通过diff(count_vec > 1)识别阶段的起始和结束区间
• 利用approx函数在区间内插值得到细胞数等于1的精确DOY
• 若某个阶段无有效起止点（如全年细胞数始终为0），则返回NA

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