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

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

我有一个这样的数据集：

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,0.1)-df$v2
df$a2 <- df$v1*rnorm(1000,1,0.2)+df$v2*rnorm(1000,0.2)

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

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

for loop

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

解决方法

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

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,0.1) - df$v2
df$a2 <- df$v1 * rnorm(1000,1,0.2) + df$v2 * rnorm(1000,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>

,

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

df <- data.frame(
  country = rep(c("W",1))

df$a1 <- rnorm(1000,0.1) + df$v1*rnorm(1000,0.1)-df$v2
df$a2 <- df$v1*rnorm(1000,0.2)+df$v2*rnorm(1000,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]],".l")
          name.u <- paste0(atts[[i]],".u")
          mod <- lm(paste0(atts[[i]],"~",socs[[j]]),.)
          assign(name.b,summary(mod)$coefficients[socs[[j]],1])
          cis <- confint(mod,level=0.95)
          assign(name.l,cis[socs[j],1])
          assign(name.u,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)