如何为数据框每行高效添加线性模型的截距、斜率、R²、p值列？

问题背景 #

你有一个包含1000个雨量站30年xy坐标与年降水量的数据框（每行对应一个站点），需要为每个站点拟合整个记录期的线性模型，并将回归的截距、斜率、R平方、p值作为新列添加到原数据框中。

解决方案 #

以下是高效实现的步骤，基于tidyverse和broom包，适合处理大规模数据：

• 步骤1：安装并加载必要的包
  tidyverse用于数据转换和分组计算，broom用于提取回归模型的结构化结果：

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

• 步骤2：将宽格式数据转换为长格式
  原始数据是宽格式（每列对应年份），线性回归需要长格式（每行对应单个站点的单年份数据）：

  df_long <- df %>%
    pivot_longer(cols = starts_with("Year."), 
                 names_to = "Year", 
                 values_to = "Precipitation") %>%
    mutate(Year = as.numeric(str_remove(Year, "Year."))) # 提取纯年份数值
  

• 步骤3：分组拟合回归并提取统计量
  按每个站点的x和y坐标分组，拟合Precipitation ~ Year的线性模型，然后提取所需统计量：

  reg_stats <- df_long %>%
    group_by(x, y) %>%
    nest() %>% # 按站点嵌套数据
    mutate(
      # 为每个站点拟合线性模型
      model = map(data, ~ lm(Precipitation ~ Year, data = .x)),
      # 提取截距、斜率及其p值
      coefs = map(model, tidy),
      # 提取R平方等模型整体统计量
      glance_stats = map(model, glance)
    ) %>%
    # 展开系数数据并重塑为宽格式
    unnest(coefs) %>%
    pivot_wider(names_from = term, values_from = c(estimate, p.value)) %>%
    # 展开模型整体统计量
    unnest(glance_stats) %>%
    # 选择并重命名需要的列
    select(
      x, y,
      intercept = estimate_(Intercept),
      slope = estimate_Year,
      r_squared = r.squared,
      slope_p_value = p.value_Year
    )
  

• 步骤4：合并结果到原始数据框
  将提取的回归统计量与原始数据框按x和y合并，得到最终数据：

  final_df <- df %>% left_join(reg_stats, by = c("x", "y"))
  

关键说明 #

• 这种嵌套+映射的方式相比循环更高效，处理1000个站点的30年数据毫无压力；
• broom包的tidy()和glance()函数能将回归结果转换为整洁的数据框，方便后续合并；
• 如果不需要额外包，也可以用apply系列函数，但代码可读性和维护性会大幅下降。

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