求助：基于tidyverse按jar_camp计算甲烷-时间二阶回归的R²值

嘿，作为tidyverse新手，这个需求其实很好实现！咱们可以结合dplyr的分组功能和broom包（专门用来把统计模型结果转成整洁数据框）来完成，步骤清晰又符合tidy风格。

第一步：加载所需包 #

首先确保你安装了tidyverse和broom，如果没装先运行安装命令：

install.packages(c("tidyverse", "broom"))

然后加载包：

library(tidyverse)
library(broom)

第二步：核心代码实现 #

咱们的思路是：按jar_camp分组，对每个组拟合二阶多项式回归模型，提取模型的R²值，再把这个值合并回原数据集的每一行（同一个jar_camp组的所有行都对应同一个R²）。

完整代码如下（假设你的数据集名为df，替换成你实际的数据集名称即可）：

df_with_r2 <- df %>%
  # 按jar_camp分组并嵌套数据，每个组的数据打包成列表元素
  group_nest(jar_camp) %>%
  # 对每个组的嵌套数据拟合二阶多项式模型
  mutate(model = map(data, ~ lm(ch4_umol ~ poly(stamp, 2, raw = TRUE), data = .x))) %>%
  # 从模型中提取R²值
  mutate(r_squared = map_dbl(model, ~ glance(.x)$r.squared)) %>%
  # 展开嵌套数据，把R²合并回原数据集
  unnest(data) %>%
  # 取消分组（可选，根据后续需求调整）
  ungroup()

代码细节解释 #

• group_nest(jar_camp)：把每个jar_camp组的数据打包成列表列，避免写循环，完美适配tidyverse的向量化操作逻辑；
• map(data, ~ lm(...))：用purrr::map对每个组的数据拟合线性模型，poly(stamp, 2, raw=TRUE)表示用原始二次项（而非正交多项式），后续如果要计算甲烷生成速率（浓度随时间的变化率）会更直观；
• map_dbl(model, ~ glance(.x)$r.squared)：用broom::glance提取模型的关键统计量（包括R²），map_dbl把结果转成数值列，方便合并；
• unnest(data)：把嵌套的数据集展开，让原始数据的每一行都带上对应组的R²值。

扩展：计算甲烷生成速率 #

如果你还需要计算每个时间点的甲烷生成速率，二次函数形式为ch4_umol = a + b*stamp + c*stamp²，速率就是导数rate = b + 2*c*stamp，可以在上面的代码基础上继续拓展：

df_with_r2_and_rate <- df_with_r2 %>%
  # 从模型中提取系数a、b、c
  mutate(coefs = map(model, ~ tidy(.x)$estimate)) %>%
  # 将系数拆分到单独列
  mutate(
    a = map_dbl(coefs, ~ .x[1]),
    b = map_dbl(coefs, ~ .x[2]),
    c = map_dbl(coefs, ~ .x[3])
  ) %>%
  # 计算每个时间点的甲烷生成速率
  mutate(ch4_rate = b + 2*c*as.numeric(stamp)) %>%
  # 移除不需要的中间列（可选）
  select(-model, -coefs)

这样你就得到了带R²和甲烷生成速率的完整数据集啦！

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