quasibinomial_investigation.R

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