Skip to content

Instantly share code, notes, and snippets.

@jongbinjung
Created March 2, 2018 17:35
Show Gist options
  • Save jongbinjung/d7dbc0df2a23d1581f70a3af974615dd to your computer and use it in GitHub Desktop.
Save jongbinjung/d7dbc0df2a23d1581f70a3af974615dd to your computer and use it in GitHub Desktop.
results <- foreach(i=1:1000, .combine = bind_rows) %:%
foreach(method=c("binomial", "quasibinomial"), .combine = bind_rows) %dopar% {
beta <- 0.1
df <- data.frame(id = rep(1:1000, 2)) %>%
mutate(m = rnorm(2000, 0, 1),
x = rnorm(2000, m, 2),
e = rnorm(2000, 0, beta / 2),
y = inv_logit(1 + beta*x+e)) %>%
group_by(id) %>%
mutate(p = rbeta(1, 2, 2),
w = c(p[1], 1-p[1])) %>%
ungroup()
df_ <- df %>%
group_by(id) %>%
summarize_at(vars(x,y), funs(weighted.mean(., w)))
bind_rows(
broom::tidy(glm(y ~ x, family = method, data = df, weights = w)) %>%
mutate(type = "weighted", method = method),
broom::tidy(glm(y ~ x, family = method, data = df_)) %>%
mutate(type = "aggregated", method = method)
) %>%
filter(term == "x")
}
results %>%
mutate(lower = estimate - 2*std.error,
upper = estimate + 2*std.error,
good = beta < upper & beta > lower) %>%
group_by(type, method) %>%
summarize(coverage = mean(good),
bias = mean(estimate) - beta,
estimated_se = mean(std.error),
empirical_se = sd(estimate))
#> type method coverage bias estimated_se empirical_se
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 aggregated binomial 1.000 -4.589e-04 0.0415065 0.0007334
#> 2 aggregated quasibinomial 0.901 -4.251e-04 0.0007210 0.0007770
#> 3 weighted binomial 1.000 -6.743e-05 0.0322919 0.0005430
#> 4 weighted quasibinomial 0.929 -3.488e-05 0.0005052 0.0005649
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment