如何使用扫帚在整洁的输出中包含来自多个模型的置信区间？答案

【问题标题】：How to include confidence intervals from multiple models in tidy output using broom?如何使用扫帚在整洁的输出中包含来自多个模型的置信区间？
【发布时间】：2019-12-27 14:00:22
【问题描述】：

我正在尝试使用 broom::tidy 从许多线性模型中输出一些结果，包括置信区间，但输出似乎只包括第一个模型的置信区间。

线性模型具有相同的预测变量，但响应不同。

考虑以下示例：

library(tidyverse)
library(broom)

# Create toy dataframe.

df <- tibble(
  x = sample(100, replace = TRUE),
  y1 = runif(100),
  y2 = rnorm(100)
)


# Fit linear models, each with x as predictor and y1 and y2 respectively as responses.

my_models <- lm(
  cbind(y1, y2) ~ x,
  data = df
)


# Output results as a tidy tibble.

tidy(my_models, conf.int = TRUE)


# Check confidence intervals with other function.

confint(my_models)

函数tidy(my_models, conf.int = TRUE)返回以下内容：

> tidy(my_models, conf.int = TRUE)
# A tibble: 4 x 8
  response term          estimate std.error statistic       p.value  conf.low conf.high
  <chr>    <chr>            <dbl>     <dbl>     <dbl>         <dbl>     <dbl>     <dbl>
1 y1       (Intercept)  0.370      0.0572      6.47   0.00000000392  0.256      0.483  
2 y1       x            0.00176    0.000949    1.86   0.0663        -0.000121   0.00365
3 y2       (Intercept) -0.0252     0.215      -0.117  0.907          0.256      0.483  
4 y2       x            0.0000574  0.00357     0.0161 0.987         -0.000121   0.00365

请注意，截距和x 的置信区间的边界是两个模型（或响应）。我希望它们会有所不同。

比较函数confint(my_models)的输出：

> confint(my_models)
                       2.5 %      97.5 %
y1:(Intercept)  0.2562157921 0.483051716
y1:x           -0.0001209424 0.003646348
y2:(Intercept) -0.4520961653 0.401713738
y2:x           -0.0070326154 0.007147456

正如预期的那样，这里的边界不同。这也是我对tidy(my_models, conf.int = TRUE) 的预期结果。由于包括y1 作为响应的模型的边界在两个函数中是相同的，我假设tidy 仅输出第一个模型的置信区间。所以我想知道我在这里做错了什么？

【问题讨论】：

标签： r lm broom

【解决方案1】：

这是旧版本扫帚中多响应线性模型的reported issue：

library(broom)
packageVersion("broom")
[1] ‘0.5.4’

mod <- lm(cbind(mpg, disp) ~ wt, mtcars)
tidy(mod, conf.int = TRUE)
# A tibble: 4 x 8
  response term        estimate std.error statistic  p.value conf.low conf.high
  <chr>    <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 mpg      (Intercept)    37.3      1.88      19.9  8.24e-19    33.5      41.1 
2 mpg      wt             -5.34     0.559     -9.56 1.29e-10    -6.49     -4.20
3 disp     (Intercept)  -131.      35.7       -3.67 9.33e- 4    33.5      41.1 
4 disp     wt            112.      10.6       10.6  1.22e-11    -6.49     -4.20

升级到最新版本就可以了：

library(broom)
packageVersion("broom")
[1] ‘0.7.0’

mod <- lm(cbind(mpg, disp) ~ wt, mtcars)
tidy(mod, conf.int = TRUE)
# A tibble: 4 x 8
  response term        estimate std.error statistic  p.value conf.low conf.high
  <chr>    <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 mpg      (Intercept)    37.3      1.88      19.9  8.24e-19    33.5      41.1
2 mpg      wt             -5.34     0.559     -9.56 1.29e-10    -6.49     -4.20
3 disp     (Intercept)  -131.      35.7       -3.67 9.33e- 4  -204.      -58.2
4 disp     wt            112.      10.6       10.6  1.22e-11    90.8     134.

【讨论】：