Search code examples
rtidyversenested-loopsglm

Combing tidyverse + survey [R]: How to use svyglm in Nest-Map-Unnest-Chain?


I am currently struggling to run weighted regression models on multiple variables in R.

When using (non-weighted) glm, I was successful by running the following:

mtcars_1 <- mtcars %>%
  nest(-gear)%>%
  mutate(model_0      = map(data, ~ glm(vs ~ drat, family = "binomial", data = .)))%>%
  mutate(model_0_tidy = map(model_0, tidy))%>%
  select(gear, model_0_tidy)%>%
  ungroup()%>%
  unnest(model_0_tidy)

That is I receive the following:

# A tibble: 6 x 6
   gear term        estimate std.error statistic p.value
  <dbl> <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1     4 (Intercept)  -15.3       22.6     -0.677   0.499
2     4 drat           4.26       5.76     0.740   0.459
3     3 (Intercept)   -3.91       7.39    -0.529   0.597
4     3 drat           0.801      2.32     0.345   0.730
5     5 (Intercept)    5.20      14.4      0.362   0.718
6     5 drat          -1.71       3.77    -0.453   0.651

However, when I would like to weight my observations and thus use svyglm from the survey-package, nesting does not work.

This was my approach:

design_0 <- svydesign(ids=~0, data = mtcars, weights = mtaars$wt)

mtcars_2 <- mtcars%>%
  nest(-gear)%>%
  mutate(model_1 = map(data, ~ svyglm(vs ~ drat, family = quasibinomial(logit), design = design_0, data = .)))%>%
  mutate(model_1_tidy = map(model_1, tidy))%>%
  select(gear, model_1_tidy)%>%
  ungroup()%>%
  unnest(model_1_tidy)

# If suggested that wt serves as frequency weight

# Outcome

   gear term        estimate std.error statistic p.value
  <dbl> <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1     4 (Intercept)    -8.12      3.88     -2.09  0.0451
2     4 drat            2.12      1.07      1.99  0.0554
3     3 (Intercept)    -8.12      3.88     -2.09  0.0451
4     3 drat            2.12      1.07      1.99  0.0554
5     5 (Intercept)    -8.12      3.88     -2.09  0.0451
6     5 drat            2.12      1.07      1.99  0.0554

Estimates for each type of gear (that is 3,4,5) turns out to be the same.

It appears as if nesting was essentially ignored here.

Are there any solutions for combining svyglm with nest-map-unnest? Or will I have to look for other, less comfortable ways?

Thank you!


Solution

  • try to do it this way

      mtcars%>%
      nest(-gear) %>% 
      mutate(design = map(data, ~ svydesign(ids=~0, data = .x, weights = ~ wt)),
             model = map(.x = design,
                         .f = ~ svyglm(vs ~ drat,
                                       family = quasibinomial(logit),
                                       design = .x))) %>% 
      mutate(model_tidy = map(model, tidy)) %>% 
      select(gear, model_tidy)%>%
      ungroup()%>%
      unnest(model_tidy)