我有一个这样的数据集:
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()系数和置信区间(CIs)。为此，我想使用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 ?
我建议使用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>
我在另一个问题中看到，这些类型的转换来创建新的变量可以使用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)