如何用dplyr对多分组数据应用自定义water.model函数？

问题背景 #

你已经构建了包含多位置、多年份的水文数据集，定义了依赖water.update内部函数的water.model函数——单个位置+年份的子数据运行该函数时，能正常生成WAT、RO、DR三列结果。现在需要通过dplyr的group_by(loc.id, year)批量处理所有分组数据，最终得到包含新增列的完整数据集。

你的数据集生成代码：

big.data <- data.frame(loc.id = rep(1:3, each = 10*3), year = rep(rep(1981:1983, each = 10),times = 3), day = rep(1:10, times = 3*3), CN = rep(c(50,55,58), each = 10*3), top.FC = rep(c(72,76,80),each = 10*3), DC = rep(c(0.02,0.5,0.8), each = 10*3), WAT0 = rep(c(20,22,26), each = 10*3), Precp = sample(1:100,90, replace = T), ETo = sample(1:10,90, replace = T))

你的模型函数定义：

water.model <- function(dat){
  top.FC <- unique(dat$top.FC)
  dat$WAT <- -9.9
  dat$RO <- -9.9
  dat$DR <- -9.9
  dat$WAT[1] <- top.FC/2 # WAT.i is a constant
  dat$RO[1] <- NA
  dat$DR[1] <- NA
  for(d in 1:(nrow(dat)-1)){
    dat[d + 1,10:12] <- water.update(WAT0 = dat$WAT[d], RAIN.i = dat$Precp[d + 1], ETo.i = dat$ETo[d + 1], CN = unique(dat$CN), DC = unique(dat$DC), top.FC = unique(dat$top.FC))
  }
  return(dat)
}
water.update <- function(WAT0, RAIN.i, ETo.i, CN, DC, top.FC){
  S = 25400/CN - 254; IA = 0.2*S
  if (RAIN.i > IA) {
    RO = (RAIN.i - 0.2 * S)^2/(RAIN.i + 0.8 * S)
  } else {
    RO = 0
  }
  if (WAT0 + RAIN.i - RO > top.FC) {
    DR = DC * (WAT0 + RAIN.i - RO - top.FC)
  } else {
    DR = 0
  }
  dWAT = RAIN.i - RO - DR - ETo.i
  WAT1 = WAT0 + dWAT
  WAT1 <- ifelse(WAT1 < 0, 0, WAT1)
  return(list(WAT1,RO,DR))
}

解决方案 #

在dplyr中，我们可以用group_by()结合group_modify()或reframe()实现分组批量处理，以下是两种可行方法：

方法1：使用group_modify（兼容性更广） #

library(dplyr)

# 分组并批量应用模型函数
final_result <- big.data %>%
  group_by(loc.id, year) %>%
  group_modify(~ water.model(.x)) %>% # .x代表每个分组的子数据框
  ungroup() # 取消分组，得到普通数据框格式

方法2：使用reframe（dplyr 1.1.0+版本推荐） #

如果你的dplyr版本在1.1.0及以上，reframe()是更简洁的选择：

final_result <- big.data %>%
  group_by(loc.id, year) %>%
  reframe(water.model(pick(everything()))) %>%
  ungroup()

关键说明 #

• 两种方法都会自动将每个分组的处理结果合并成完整数据集，最终的final_result包含所有原始列，以及新增的WAT、RO、DR列。
• 运行前需确保已加载dplyr包（library(dplyr)）。
• group_modify中的.x参数会自动传入每个分组的子数据框，直接对接你的water.model函数即可。

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