如何在R中用循环按单株逐年拟合生长阶段逻辑回归并生成日预测

问题说明 #

现有数据集包含10年共60棵树（每年6棵）的细胞生长阶段观测数据：

• 核心字段：Year（年份）、Tree（树木编号）、DOY（年积日）、DOY2（DOY平方项）
• 响应变量：Enlarging、Thickening、Maturing三个生长阶段的存在(1)/缺失(0)标记

需求：

• 对Enlarging和Thickening阶段，拟合包含DOY + DOY2的二项式GLM逻辑回归
• 对Maturing阶段，拟合仅包含DOY的二项式GLM逻辑回归
• 按每棵树每年为一组（共60组）批量完成模型拟合，并生成全年1-365天的预测概率

实现方案 #

核心步骤 #

1. 构建分组标识：将Year和Tree组合为唯一分组键，确保每组对应单棵树单一年份的观测
2. 定义模型公式：针对不同生长阶段预设对应的回归公式
3. 批量拟合与预测：通过循环遍历所有分组，完成模型训练并生成全年预测
4. 整合结果：将所有分组的预测结果合并为统一数据集，便于后续分析

完整代码实现 #

# 示例数据集（替换为你的实际数据）
df <- data.frame(
  Year = c(2012, 2012, 2012, 2012, 2012, 2012, 2012,
           2012, 2012, 2012, 2013, 2013, 2013,
           2013, 2013),
  Tree = c(15, 15, 15, 15, 15, 22, 22, 22, 22, 22, 41, 41,
           41, 41, 41),
  DOY = c(65, 97, 125, 177, 214, 65, 97, 125, 177, 214,
          61, 99, 118, 166, 221),
  DOY2 = c(4225, 9409, 15625, 31329, 45796, 4225, 9409, 15625,
           31329, 45796, 3721, 9801, 13924, 27556, 48841),
  Enlarging = c(0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0),
  Thickening = c(0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0),
  Maturing = c(0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0)
)

# 1. 创建分组标识：Year + Tree的组合
df$group <- paste(df$Year, df$Tree, sep = "_")
unique_groups <- unique(df$group)

# 2. 定义各阶段的模型公式
model_formulas <- list(
  Enlarging = formula(Enlarging ~ DOY + DOY2),
  Thickening = formula(Thickening ~ DOY + DOY2),
  Maturing = formula(Maturing ~ DOY)
)

# 3. 生成全年预测用的DOY数据集
new_doy <- data.frame(DOY = 1:365, DOY2 = (1:365)^2)

# 4. 批量拟合模型并生成预测
all_predictions <- list()

for (grp in unique_groups) {
  # 提取当前分组的观测数据
  group_data <- subset(df, group == grp)
  
  # 拆分年份和树木编号
  grp_split <- strsplit(grp, "_")[[1]]
  grp_year <- as.integer(grp_split[1])
  grp_tree <- as.integer(grp_split[2])
  
  # 遍历每个生长阶段拟合模型并预测
  stage_preds <- list()
  for (stage in names(model_formulas)) {
    # 拟合GLM模型
    glm_model <- glm(model_formulas[[stage]], 
                     family = binomial(link = "logit"), 
                     data = group_data)
    
    # 生成预测概率
    preds <- predict(glm_model, newdata = new_doy, type = "response")
    
    # 整理为数据框
    stage_preds[[stage]] <- data.frame(
      Year = grp_year,
      Tree = grp_tree,
      DOY = new_doy$DOY,
      Stage = stage,
      Predicted_Prob = preds
    )
  }
  
  # 合并当前分组的所有阶段预测
  group_preds <- do.call(rbind, stage_preds)
  all_predictions[[grp]] <- group_preds
}

# 5. 整合所有分组的预测结果
final_predictions <- do.call(rbind, all_predictions)

# 查看结果前几行
head(final_predictions)

代码说明 #

• 分组逻辑：通过Year_Tree的字符串组合确保每组是独立的树木年度观测，避免交叉干扰
• 模型区分：用列表存储不同阶段的公式，统一遍历逻辑，减少重复代码
• 预测输出：将预测结果整理为包含年份、树木编号、DOY、阶段和预测概率的结构化数据，方便后续可视化或统计分析
• 效率优化：基础循环结合列表存储结果，在60组数据下运行效率足够；若数据量更大，可替换为purrr包的向量化操作进一步提速

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