Data source

Startup

Packages, data, recode.

options(digits = 2)
library(pacman)
p_load(kirkegaard, dplyr, haven, rms, ggeffects)
theme_set(theme_bw())

#load data
pew = haven::read_spss("data/Pew-Research-Center-2014-U.S.-Religious-Landscape-Study/Dataset - Pew Research Center 2014 Religious Landscape Study National Telephone Survey - Version 1.1 - December 1 2016.sav")

#table of variables
var_table = data_frame(
  number = seq_along(pew),
  code = names(pew)
)
## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
#function to recode
spss_labels = function(x) {
  purrr::map2_chr(x, names(x), function(x, x_orig) {
  #get label, replace illegal characters
  x_label = attr(x, "label")
  x_label_nice = str_replace_all(x_label, pattern = "[ \\-\\/]", replacement = "_") %>% 
    str_replace_all(pattern = "[\\(\\)]", replacement = "")
  
  #if no label, use orignal
  if (are_equal(x_label_nice, character())) return(x_orig)
  x_label_nice
  })
}


#save meaningful names to var table
var_table$name = spss_labels(pew)

#get rid of labelled type
spss_rm_labelled = function(x) {
  #does object have labels?
  has_labels = has_attr(x, "labels")
  
  #if labels, convert to factor
  if (has_labels) {
    y = as_factor(x)
  } else {
    #convert to ordinary numeric
    y = as.vector(x)
  }
  
  y
}

#convert each variable
pew = df_colFunc(pew, spss_rm_labelled)

#function to note type
spss_type = function(x) {
  
  map_chr(x, function(x) {
      #factor?
  if (is.factor(x)) {
    y = "Nominal: " + str_c(levels(x), collapse = "; ")
    return(y)
  }
  
  #con
  return("Continuous")
  })

}

var_table$type = spss_type(pew)

Big govt support

Policy % in favor by group.

#recode
pew %<>% mutate(
  race = racecmb,
  big_govt = qb20 == "Bigger government, more services",
  #relabel
  education = factor(educ, labels = levels(educ) %>% str_match(pattern = "^[^,\\(]+"))
)

#race
GG_group_means(pew, var = "big_govt", groupvar = "race", subgroupvar = "SEX") +
  scale_y_continuous("% who support bigger governments with more services", labels = scales::percent) +
  xlab("Self-identified race/ethnicity")

ggsave("figures/big_govt_race.png")
## Saving 7 x 5 in image
#income
GG_group_means(pew, var = "big_govt", groupvar = "income", subgroupvar = "SEX") +
  scale_y_continuous("% who support bigger governments with more services", labels = scales::percent) +
  theme(axis.text.x = element_text(size = 7, angle = -30, vjust = .9)) +
  xlab("Income group")

ggsave("figures/big_govt_income.png")
## Saving 7 x 5 in image
#education
GG_group_means(pew, var = "big_govt", groupvar = "education", subgroupvar = "SEX") +
  scale_y_continuous("% who support bigger governments with more services", labels = scales::percent) +
  theme(axis.text.x = element_text(size = 6)) +
  xlab("Highest educational attainment")

ggsave("figures/big_govt_education.png")
## Saving 7 x 5 in image
#fertility
GG_group_means(pew, var = "big_govt", groupvar = "fertrec", subgroupvar = "SEX") +
  scale_y_continuous("% who support bigger governments with more services", labels = scales::percent) +
  xlab("Number of kids")

ggsave("figures/big_govt_fert.png")
## Saving 7 x 5 in image

Regression models

models = list(
  "nothing" = lrm(big_govt ~ SEX, data = pew),
  "race" = lrm(big_govt ~ SEX * race, data = pew),
  "race & age" = lrm(big_govt ~ SEX * race + agerec, data = pew),
  "race & age & education & income" = lrm(big_govt ~ SEX * race + agerec + education + income, data = pew),
  "race & age & education & income & fertility" = lrm(big_govt ~ SEX * race + agerec + education + income + fertrec, data = pew)
)

#summarize effect size of being female
models_sum = map2_df(models, names(models), function(m, n) {
  tibble(
    adjust_for = n,
    Female_effect = coef(m)[2],
    predict = list(ggeffect(m, "SEX"))
  )
})
models_sum
## # A tibble: 5 x 3
##   adjust_for                                  Female_effect predict         
##   <chr>                                               <dbl> <list>          
## 1 nothing                                             0.452 <df[,6] [2 × 6]>
## 2 race                                                0.501 <df[,6] [2 × 6]>
## 3 race & age                                          0.554 <df[,6] [2 × 6]>
## 4 race & age & education & income                     0.547 <df[,6] [2 × 6]>
## 5 race & age & education & income & fertility         0.567 <df[,6] [2 × 6]>
#plot effect sizes
#ggeffect(models[[1]], "SEX") %>% plot() + xlab("Sex") + ylab("Probability of supporting bigger government")

#but ggeffects plot only supports results from 1 model at a time, so we need to manually combine predictions
models_sum %>% 
  unnest(predict) %>% 
  ggplot(aes(x, predicted, color = adjust_for)) +
  geom_point(position = position_dodge(width = 1)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(width = 1)) +
  xlab("Sex") + scale_y_continuous("Probability of supporting bigger government", labels = scales::percent)

GG_save("figures/model_predictions.png")