R语言按分组做线性回归时如何为结果添加各组R²值

问题描述 #

需要在R中按照shop、art两个分组变量的类别组合分别拟合线性回归模型，使用的示例数据集如下：

df=structure(list(shop = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L), .Label = c("a", "c"), class = "factor"), art = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("b", "d"), class = "factor"), 
    Y = c(177L, 122L, 175L, 140L, 201L, 202L, 279L, 253L, 236L, 
    137L, 166L, 241L, 195L, 221L, 238L, 203L, 254L, 219L, 101L, 
    157L, 188L, 219L, 267L, 126L, 291L, 239L, 230L), x1 = c(1L, 
    0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 
    0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L), x2 = c(0L, 1L, 
    1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 
    1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L), x3 = c(0L, 0L, 0L, 
    1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 
    1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L), x4 = c(0L, 0L, 1L, 1L, 
    0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 
    0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L), x5 = c(0L, 0L, 1L, 1L, 0L, 
    0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 
    1L, 0L, 0L, 1L, 1L, 1L, 0L), x6 = c(0L, 1L, 0L, 0L, 1L, 1L, 
    0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 
    1L, 1L, 1L, 1L, 0L, 1L), x7 = c(1L, 1L, 0L, 0L, 1L, 0L, 0L, 
    0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 
    0L, 1L, 1L, 1L, 0L), x8 = c(0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 
    1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 
    0L, 1L, 0L, 1L), x9 = c(1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 
    0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 
    1L, 1L, 0L)), .Names = c("shop", "art", "Y", "x1", "x2", 
"x3", "x4", "x5", "x6", "x7", "x8", "x9"), class = "data.frame", row.names = c(NA, 
-27L))

数据集字段说明：

• Y：数值型因变量
• x1~x9：数值型自变量
• shop、art：分组因子变量

现有实现代码如下，运行后可输出分组回归系数，但缺少模型R²值：

library(tidyverse)
library(broom)
my_lm <- function(df) {
  lm(Y ~ ., data = df)
}
df %>% 
  group_by(art, shop) %>% 
  nest() %>%
  mutate(fit = map(data, my_lm),
         tidy = map(fit, tidy)) %>%
  select(-fit, - data) %>%
  unnest()

需求为在输出结果中新增R-square字段，让同一art+shop分组下的所有回归系数行都附带该组模型的R²值。

修改方案 #

借助broom::glance()函数提取模型整体拟合指标，将R²值和对应组的系数表按行绑定即可，修改后的完整代码：

library(tidyverse)
library(broom)

my_lm <- function(df) {
  lm(Y ~ ., data = df)
}

df %>% 
  group_by(art, shop) %>% 
  nest() %>%
  mutate(
    fit = map(data, my_lm),
    # 提取回归系数表
    coef_df = map(fit, tidy),
    # 提取模型拟合指标，仅保留R²字段
    r2_df = map(fit, ~glance(.x) %>% select(r.squared))
  ) %>%
  # 同组R²绑定到所有系数行
  mutate(final_res = map2(coef_df, r2_df, cbind)) %>%
  select(-fit, -data, -coef_df, -r2_df) %>%
  unnest(final_res)

代码说明 #

• broom::glance()用于输出线性模型的整体优度指标，除R²外还可提取调整R²、F统计量、AIC、BIC、残差标准差等指标，按需在select步骤保留对应字段即可。
• map2配合cbind可以将单值的R²按行补全，匹配同组所有系数行，完全符合输出格式要求。

示例数据中每个分组样本量仅11~16个，自变量共9个，拟合得到的R²可能存在过拟合偏高的情况，属于样本量与自变量数量匹配的正常现象，与代码逻辑无关。

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