R：创建几个回归输出的可发布表

Question

我一直在为一篇论文分析数据，现在已经获得了多元线性回归的结果。然而，R 提供的摘要并不是真正适合在最终论文中发表。我还以几种不同的方式指定了一个变量，以展示结果的稳健性。

我可以在 R 中创建一个漂亮的、可导出的表，其中包含变量名称（理想情况下还允许以更具信息性的方式命名变量）、估计值、标准误差、稳健的标准误差 p 值以及理想情况下的显着性指标?举例说明：

我有这样的总结输出：

Residuals:
    Min      1Q  Median      3Q     Max 
-50.868  -4.644   1.583   7.054  20.490 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)  
(Intercept)               4.710e+01  1.848e+01   2.549   0.0136 *
Var1                     -8.588e-01  2.201e+00  -0.390   0.6979  
Var2                      2.486e+00  1.055e+00   2.357   0.0220 *
log(specification1)       3.376e+00  2.152e+00   1.569   0.1223  
Var4                     -3.651e-04  2.797e-04  -1.305   0.1971  
Var5                      4.809e+00  2.654e+00   1.812   0.0753 .
Var6                     -8.706e+00  6.972e+00  -1.249   0.2170  
Var7                     -8.172e+00  5.755e+00  -1.420   0.1612  
Var8                     -3.276e+00  7.067e+00  -0.463   0.6448  
Var9                     -1.477e+01  7.849e+00  -1.882   0.0650 .

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

和

Residuals:
    Min      1Q  Median      3Q     Max 
-48.881  -5.699   0.956   8.947  17.888 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)  
(Intercept)               4.258e+01  1.750e+01   2.405   0.0195 *
Var1                      4.298e-01  2.120e+00   0.200   0.8421  
Var2                      5.179e+00  1.027e+00   2.122   0.0271 *
log(specification 2)      2.050e+00  9.435e-01   2.173   0.0338 *
Var4                     -1.420e-04  2.261e-04  -1.513   0.1356  
Var5                      4.584e+00  2.511e+00   1.826   0.0730 .

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我想找一张看起来像这样的桌子：

             Model1                                             Model2
Intercept    Estimate  Std.Error  p robust_Std.Error robust_p  Estimate  Std.Error  p  robust ...
Var1
Var2
Var3
Var4
Var5
Var6
Var7
Var8
Var9

列中当然包含估计值。有没有可以很好地做到这一点的函数/包？

提前致谢

Answer 1

我建议你使用broom 包，像这样：

fit1 <- lm(mpg ~ ., mtcars)
broom::tidy(fit1)

# # A tibble: 11 x 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 

它将从lm函数的输出中提取tibble。

---

如果您有多个模型，并且想将所有小标题与常用术语一起设置，则可以这样处理：

为您的模型创建一个x 列表。

fit1 <- lm(mpg ~ cyl + disp + gear, mtcars)
fit2 <- lm(mpg ~ cyl + hp   + drat, mtcars)

x <- list(fit1, fit2)

你可以使用这个解决方案：

library(purrr)
library(dplyr)
library(stringr)

# set names for the list
names(x) <- paste("Model", seq_along(x), sep = "_")

# tidy them up
x <- map(x, broom::tidy)

# set the list names at the beginning of each column
x <- imap(x, ~set_names(.x, paste(.y, names(.x), sep = "_")))

# rename each term column as "term"
x <- map(x, ~rename_with(.x, str_replace, pattern = ".*term", replacement = "term"))

# join them all together
reduce(x, full_join, by = "term")

它返回您要求的输出：

# A tibble: 6 x 9
  term        Model1_estimate Model1_std.error Model1_statistic Model1_p.value Model2_estimate Model2_std.error Model2_statistic Model2_p.value
  <chr>                 <dbl>            <dbl>            <dbl>          <dbl>           <dbl>            <dbl>            <dbl>          <dbl>
1 (Intercept)         34.0              4.76              7.13    0.0000000925         22.5              7.99               2.82        0.00880
2 cyl                 -1.59             0.724            -2.20    0.0366               -1.36             0.735             -1.85        0.0747 
3 disp                -0.0200           0.0109           -1.83    0.0774               NA               NA                 NA          NA      
4 gear                 0.158            0.910             0.174   0.863                NA               NA                 NA          NA      
5 hp                  NA               NA                NA      NA                    -0.0288           0.0153            -1.88        0.0704 
6 drat                NA               NA                NA      NA                     2.84             1.52               1.87        0.0725 

如果您的列表有两个以上的模型，则代码将是稳定的。