I am trying to run multiple different formulas on nested data. I was able to code most of what I need, but I am not very familiar with functions or dealing with nested data, so I don’t know how to get the information out of the lists now that I have it.

#packages and data
library(tidyverse)
library(broom)
install.packages("AER")
library("AER")
data(Affairs, package="AER")

#creating response variable
Affairs$ynaffair[Affairs$affairs >  0] <- 1
Affairs$ynaffair[Affairs$affairs == 0] <- 0

#factoring some of the variables
Affairs <- Affairs %>% 
  mutate_at(c("affairs", "religiousness", "occupation", "rating", "ynaffair"), as.factor)

#nesting by religious group
by_relig_group <- Affairs %>%
  group_by(religiousness) %>%
  nest()

#Create formulas
gender_formula <- as.formula("ynaffair ~ gender")
age_formula <- as.formula("ynaffair ~ age")
occupation_formula <- as.formula("ynaffair ~ occupation")
formula_list <- list(gender_formula,age_formula, occupation_formula)

# Designate the models and function
model <- function(df, formula_list, index) {
    list(glm(formula = formula_list[[index]] , family=binomial, data = Affairs))
}

#Apply to each data group and formula list by formula number 
df1 <- by_relig_group %>% 
  mutate(model1 = map(data, ~model(., formula_list, 1)), 
         model2 = map(data, ~model(., formula_list, 2))) 

And it looks like the code worked, but now I am not quite sure how to get it out of its current form.

At a smaller level, I was able to get all of this to run

#create nested data frame
by_relig_group <- Affairs %>%
  group_by(religiousness) %>%
  nest()

#create a function to run the linear models
model <- function(df){
  glm(ynaffair~gender + age + yearsmarried + children + religiousness + education + occupation + rating,
               family = binomial,
             data = Affairs)
}

#run the models
by_relig_group <- by_relig_group %>%
  mutate(model = map(data,model))


#unnest the models using the broom package 
by_relig_group_unnest <- by_relig_group %>%
  mutate(tidy = map(model, broom::tidy)) %>% ## I can't figure out how to apply these two lines of code above
  unnest(tidy)

The end goal is to get the frequency of each level(next problem to solve), estimate (which I’ll exponentiate to get the OR), CI, and p-value and export to Excel. If you have any suggestions on how to get CI while I am running this code, I’d be grateful for the advice.