时间戳

首页 > 解决方案 > 如何按组获取 lm() 系数和置信区间？

r - 如何按组获取 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, 2, 0.1)-df$v2
df$a2 <- df$v1*rnorm(1000, 1, 0.2)+df$v2*rnorm(1000, 2, 0.2)

数据按 分组country，我想计算lm()许多变量（即a1、a2、v1、 ）的系数和置信区间 (CI v2)。为此，我想使用group_byfromdplyr来按国家/地区分组，并使用 afor 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？

标签： rfunctionfor-loopdplyr

解决方案

我建议使用新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, 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>

推荐阅读