如何按组获取 lm() 系数和置信区间？

Question

我有一个这样的数据集：

df <- data.frame(country = rep(c("W", "Q"), 500),
                 v1 = rnorm(1000, 10, 2),
                 v2 = rnorm(1000, 8, 1))

df$a1 <- rnorm(1000, 2, 0.1) + df$v1*rnorm(1000, 2, 0.1)-df$v2
df$a2 <- df$v1*rnorm(1000, 1, 0.2)+df$v2*rnorm(1000, 2, 0.2)

数据按country分组，我想计算许多变量（即a1、a2、v1、v2）的lm()系数和置信区间（CI）。为此，我想使用dplyr 中的group_by 按国家/地区分组，并使用for loop 指定我要为其计算beta 系数和CI 的变量。

到目前为止，我已经设法创建了一个几乎可以做到这一点的函数，但是当应用于分组数据集时，它只返回第一个系数和计算的 CI。

betas <- function(dat, atts, socs){
    for (i in 1:length(atts)) {
      for (j in 1:length(socs)) {
        mod <- lm(paste0(atts[[i]], "~", socs[[j]]), dat)
        mod.s <- summary(mod)$coefficients
        cis <- confint(mod, socs[[j]], level=0.95)
        dat <- dat %>% mutate(!!paste0(atts[[i]], ".", socs[[j]], ".b") := mod.s[socs[[j]], 1],
                 !!paste0(atts[[i]], ".", socs[[j]], ".l") := cis[socs[j], 1],
                 !!paste0(atts[[i]], ".", socs[[j]], ".u") := cis[socs[j], 2])
      }
    }
  return(dat)
}

df1 <- df %>% group_by(country) %>% betas(., atts = c("a1", "a2"), socs = c("v1", "v2"))

#To see how the output is all the same.
hist(df1$a2.v2.b)

如何计算按国家分组的数据集的多对变量的系数和 CI？

Answer 1

我建议使用dplyr 中的新grouped_ 函数来执行此操作。这是为每个分组变量级别运行所有感兴趣的回归模型的简洁方法：

set.seed(123)
library(dplyr)

# data
df <- data.frame(
  country = rep(c("W", "Q"), 500),
  v1 = rnorm(1000, 10, 2),
  v2 = rnorm(1000, 8, 1)
)

df$a1 <- rnorm(1000, 2, 0.1) + df$v1 * rnorm(1000, 2, 0.1) - df$v2
df$a2 <- df$v1 * rnorm(1000, 1, 0.2) + df$v2 * rnorm(1000, 2, 0.2)


dplyr::bind_rows(
  # a1 as IV
  df %>%
    group_by(country) %>%
    group_modify(~ parameters::model_parameters(stats::lm(cbind(v1, v2) ~ a1, data = .x))),
  # a2 as IV
  df %>%
    group_by(country) %>%
    group_modify(~ parameters::model_parameters(stats::lm(cbind(v1, v2) ~ a2, data = .x)))
)
#> # A tibble: 16 x 11
#> # Groups:   country [2]
#>    country Parameter   Coefficient      SE    CI  CI_low  CI_high     t df_error
#>    <chr>   <chr>             <dbl>   <dbl> <dbl>   <dbl>    <dbl> <dbl>    <int>
#>  1 Q       (Intercept)      3.83   0.103    0.95  3.63    4.03    37.2       498
#>  2 Q       a1               0.442  0.00714  0.95  0.428   0.456   61.9       498
#>  3 Q       (Intercept)      8.73   0.154    0.95  8.43    9.03    56.8       498
#>  4 Q       a1              -0.0470 0.0107   0.95 -0.0680 -0.0261  -4.41      498
#>  5 W       (Intercept)      3.59   0.118    0.95  3.36    3.82    30.5       498
#>  6 W       a1               0.460  0.00797  0.95  0.445   0.476   57.8       498
#>  7 W       (Intercept)      8.43   0.163    0.95  8.11    8.75    51.7       498
#>  8 W       a1              -0.0302 0.0110   0.95 -0.0519 -0.00858 -2.74      498
#>  9 Q       (Intercept)      3.17   0.487    0.95  2.22    4.13     6.51      498
#> 10 Q       a2               0.259  0.0185   0.95  0.223   0.295   14.0       498
#> 11 Q       (Intercept)      4.38   0.245    0.95  3.90    4.86    17.9       498
#> 12 Q       a2               0.142  0.00929  0.95  0.124   0.161   15.3       498
#> 13 W       (Intercept)      2.83   0.494    0.95  1.86    3.80     5.73      498
#> 14 W       a2               0.279  0.0186   0.95  0.243   0.316   15.0       498
#> 15 W       (Intercept)      4.17   0.243    0.95  3.69    4.65    17.1       498
#> 16 W       a2               0.146  0.00919  0.95  0.128   0.164   15.9       498
#> # ... with 2 more variables: p <dbl>, Response <chr>

Answer 2

我在另一个问题中看到，这些类型的创建新变量的转换可以使用来自dplyr 的do({}) 来完成。以下似乎可以解决问题：

df <- data.frame(
  country = rep(c("W", "Q"), 500),
  v1 = rnorm(1000, 10, 2),
  v2 = rnorm(1000, 8, 1))

df$a1 <- rnorm(1000, 2, 0.1) + df$v1*rnorm(1000, 2, 0.1)-df$v2
df$a2 <- df$v1*rnorm(1000, 1, 0.2)+df$v2*rnorm(1000, 2, 0.2)


betas <- function(dw, atts, socs){
  for (i in 1:length(atts)) {
    for (j in 1:length(socs)) {
      dw1 <- dw %>% group_by(country) %>% do(
        {
          name.b <- paste0(atts[[i]], ".", socs[[j]], ".b")
          name.l <- paste0(atts[[i]], ".", socs[[j]], ".l")
          name.u <- paste0(atts[[i]], ".", socs[[j]], ".u")
          mod <- lm(paste0(atts[[i]], "~", socs[[j]]), .)
          assign(name.b, summary(mod)$coefficients[socs[[j]], 1])
          cis <- confint(mod, socs[[j]], level=0.95)
          assign(name.l, cis[socs[j], 1])
          assign(name.u, cis[socs[j], 2])
          dat <- data.frame(., get(name.b), get(name.l), get(name.u))
          colnames(dat)[(length(names(dat))-2):length(names(dat))] <- c(name.b, name.l, name.u)
          dat
        })
    }
  }
  dw1
}


df1 <- betas(df, c("a1", "a2"), c("v1", "v2"))

hist(df1$a2.v2.b)