R语言如何使用dplyr和broom对mpg与其余变量逐一拟合线性模型

给定R内置mtcars数据集，前6行内容如下：

> head(mtcars)
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

需求 #

以mpg为因变量，数据集中其余所有变量分别作为单一自变量，逐一拟合一元线性模型，覆盖从mpg ~ cyl到mpg ~ carb的全部单变量回归，要求基于dplyr生态实现。

现有代码的问题 #

最初编写的代码如下：

library(broom)
 library(tidyr)
 library(dplyr)
 lm(mpg ~ ., data = mtcars) %>% 
 # glance() # 运行失败
 tidy()

这段代码实际拟合的是包含所有自变量的多元线性模型，返回结果如下：

# A tibble: 11 × 5
   term        estimate std.error statistic p.value
   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
 1 (Intercept)  12.3      18.7        0.657  0.518 
 2 cyl          -0.111     1.05      -0.107  0.916 
 3 disp          0.0133    0.0179     0.747  0.463 
 4 hp           -0.0215    0.0218    -0.987  0.335 
 5 drat          0.787     1.64       0.481  0.635 
 6 wt           -3.72      1.89      -1.96   0.0633
 7 qsec          0.821     0.731      1.12   0.274 
 8 vs            0.318     2.10       0.151  0.881 
 9 am            2.52      2.06       1.23   0.234 
10 gear          0.655     1.49       0.439  0.665 
11 carb         -0.199     0.829     -0.241  0.812 

单独运行cyl对mpg的一元回归代码，得到的结果为：

> tidy(lm(mpg ~ cyl, data = mtcars))
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    37.9      2.07      18.3  8.37e-18
2 cyl            -2.88     0.322     -8.92 6.11e-10

两种方式得到的cyl变量系数从-2.88变为-0.11，差异显著，原代码不符合需求。

原代码的核心问题是公式里的.代表自动代入除因变量外的所有字段，所以拟合的是多元回归模型，系数是控制其他变量后的条件效应，和一元回归的边际效应天然存在差异，不是代码报错。
要批量拟合所有一元回归，可以通过构造嵌套数据框逐行建模实现，代码如下：

library(dplyr)
library(tidyr)
library(broom)

# 提取除mpg外的所有自变量名
names(mtcars)[-1] %>%
  tibble(predictor = .) %>%
  mutate(
    # 逐个生成单变量回归公式
    formu = paste("mpg ~", predictor),
    # 逐公式拟合线性模型
    model = lapply(formu, function(f) lm(as.formula(f), data = mtcars)),
    # 提取模型系数结果
    coef_res = lapply(model, tidy)
    # 若需要R²、AIC等模型整体指标，可添加下一行
    # perf_res = lapply(model, glance)
  ) %>%
  # 展开结果表，得到所有单变量回归的汇总结果
  unnest(coef_res)

运行后得到的结果中，每个自变量对应的系数、标准误、p值和单独运行对应一元回归的结果完全一致，例如cyl的系数为-2.88，和单独运行lm(mpg ~ cyl)的输出匹配。

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