NOTE: DON’T WORRY ABOUT FIGURES LOOKING “SQUISHED.” THAT IS JUST AN ARTIFACT OF THE WEBSITE.
Preparation
Show the code
# librarieslibrary(tidyverse)library(haven)library(survey)library(patchwork)library(sjPlot)library(MASS)library(margins)# load datasetwd("~/Downloads")data <-read_dta('barkerdata2.dta')data21 <-read_dta('ces2021.dta')# uniform aesthesticsdavid_func <-function() {list(theme_bw(),theme(panel.grid =element_blank(),axis.text =element_text(size =10),axis.title =element_text(size =12),legend.text =element_text(size =10),legend.title =element_text(size =12)) )}# Creating a few variablesdata <- data %>%mutate(pidfac_d =case_when(pid3lean01 ==0~1, T ~0),pidfac_i =case_when(pid3lean01 ==0.5~1, T ~0),pidfac_r =case_when(pid3lean01 ==1~1, T ~0))# Creating survey design objects to include (1) sampling weights, and (2) clustered SEs at the state levelmydesign <-svydesign(ids =~inputstate, weights =~weight, data = data)data21a <- data21 %>%drop_na(teamweight)mydesign21 <-svydesign(ids =~inputstate, weights =~teamweight, data = data21a)
Chapter #2
Figure 2.1
NOTE: THE ORIGINAL NUMBERS IN-TEXT DID NOT APPEAR TO USE SURVEY-WEIGHTS?? WHEN I RUN THE CORS W/ WEIGHTS, THE NUMBERS VARY SLIGHTLY
Show the code
library(weights)# The code is very long because because we need to bootstrap the confidence intervals with weights :)# Core ----------------------------------------------------# for Ds/Rscor1 <- data %>%drop_na(intel901march23newest, antiintelcore601newest) %>%filter(pid3lean01 !=0.5) %>%group_by(pid3lean01) %>% dplyr::summarize(cor = weights::wtd.cor(x = intel901march23newest, y = antiintelcore601newest,weight = weight,bootse = T,bootn =1000)) %>%mutate(type ='Core anti-intellectualism')# Full samplecor2 <- data %>%drop_na(intel901march23newest, antiintelcore601newest) %>% dplyr::summarize(cor = weights::wtd.cor(x = intel901march23newest, y = antiintelcore601newest,weight = weight,bootse = T,bootn =1000)) %>%mutate(pid3lean01 =0.5,type ='Core anti-intellectualism')# indignant ----------------------------------------------------# for Ds/Rscor3 <- data %>%drop_na(intel901march23newest, antiintpop2pluspolitic3_501) %>%filter(pid3lean01 !=0.5) %>%group_by(pid3lean01) %>% dplyr::summarize(cor = weights::wtd.cor(x = intel901march23newest, y = antiintpop2pluspolitic3_501,weight = weight,bootse = T,bootn =1000)) %>%mutate(type ='Indignant anti-intellectualism')# Full samplecor4 <- data %>%drop_na(intel901march23newest, antiintpop2pluspolitic3_501) %>% dplyr::summarize(cor = weights::wtd.cor(x = intel901march23newest, y = antiintpop2pluspolitic3_501,weight = weight,bootse = T,bootn =1000)) %>%mutate(pid3lean01 =0.5,type ='Indignant anti-intellectualism')# combine into one df + calculate CIs ----------------------------------cors <-bind_rows(cor1, cor2, cor3, cor4)cors$cor1 <- cors$cor[, 1] # corecors$cor2 <- cors$cor[, 3] # boostrapped secors2 <-data.frame(pid3lean01 = cors$pid3lean01, cor = cors$cor1, se = cors$cor2, type = cors$type)cors2 <- cors2 %>%mutate(min = cor -1.96*se, # CIsmax = cor +1.96*se) %>%mutate(group =factor(pid3lean01,levels =c(0.5, 0, 1),labels =c('Republicans', 'Democrats', 'Full sample')),cor =round(cor,2))# time to plot ---------------------------------------------------------cors2 %>%ggplot(aes(y = cor, x =factor(type, levels =c('Indignant anti-intellectualism','Core anti-intellectualism')), fill = group,ymin = min, ymax = max)) +geom_col(color ='black', position =position_dodge()) +geom_errorbar(position =position_dodge(width =0.9), width =0.2) +geom_text(aes(label =paste(cor), vjust =ifelse(cor >=0, 0, 0.5), hjust =ifelse(cor >=0, 0, 1.9)), position =position_dodge(width =0.9)) +scale_fill_grey(start =0.8, end =0.3) +labs(fill ='Sample',x ='', y ='Pearson\'s Correlation with Intellectual Identity') +david_func() +coord_flip() +ylim(-1, .15) +guides(fill =guide_legend(reverse =TRUE))
Chapter #3
Figure 3.1
Show the code
# Models (using separate binomial probits because the I can't find a way to run multinomial probits while using complex design elements, e.g., clusters/weights)R <-svyglm(pidfac_r ~ justbachelorsdegree + graddegree + black + hispanic + male + age01 + faminc_missmedsub01 + churchatt601, design = mydesign, family =quasibinomial(link ="probit"))I <-svyglm(pidfac_i ~ justbachelorsdegree + graddegree + black + hispanic + male + age01 + faminc_missmedsub01 + churchatt601, design = mydesign, family =quasibinomial(link ="probit"))D <-svyglm(pidfac_d ~ justbachelorsdegree + graddegree + black + hispanic + male + age01 + faminc_missmedsub01 + churchatt601, design = mydesign, family =quasibinomial(link ="probit"))# Function to extract predicted probabilitiespred_func <-function(reg, bach, grad, pid, educ) { predict_data <-data.frame(# desired levels of educationjustbachelorsdegree = bach,graddegree = grad,# covariates fixed at their menasblack =weighted.mean(data$black, na.rm =TRUE, w = data$weight),hispanic =weighted.mean(data$hispanic, na.rm =TRUE, w = data$weight),male =weighted.mean(data$male, na.rm =TRUE), w = data$weight,age01 =weighted.mean(data$age01, na.rm =TRUE, w = data$weight),faminc_missmedsub01 =weighted.mean(data$faminc_missmedsub01, na.rm =TRUE, w = data$weight),churchatt601 =weighted.mean(data$churchatt601, na.rm =TRUE, w = data$weight) )# tidying preds <-predict(reg, newdata = predict_data, type ="response") preds <-as.data.frame(preds) preds$pid <- pid preds$educ <- educreturn(preds)}# Creating dataframe of predictionsd1 <-pred_func(D, 0, 0, 'Democrat', '< Bachelor\'s degree') d2 <-pred_func(D, 1, 0, 'Democrat', 'Bachelor\'s degree') d3 <-pred_func(D, 1, 1, 'Democrat', 'Graduate degree') i1 <-pred_func(I, 0, 0, 'Independent', '< Bachelor\'s degree') i2 <-pred_func(I, 1, 0, 'Independent', 'Bachelor\'s degree') i3 <-pred_func(I, 1, 1, 'Independent', 'Graduate degree') r1 <-pred_func(R, 0, 0, 'Republican', '< Bachelor\'s degree') r2 <-pred_func(R, 1, 0, 'Republican', 'Bachelor\'s degree') r3 <-pred_func(R, 1, 1, 'Republican', 'Graduate degree') f3.1<-bind_rows(d1, d2, d3, i1, i2, i3, r1, r2, r3)f3.1<- f3.1%>%mutate(ymax = response + SE*1.96,ymin = response - SE*1.96)# Plotf3.1%>%ggplot(aes(x = educ, y = response, group = pid)) +geom_line() +geom_point() +geom_errorbar(aes(ymin = ymin, ymax = ymax), width =0.1) +david_func() +facet_wrap(~pid) +labs(x ='Educational Attainment', y ='Probability of PID') +theme(axis.text.x =element_text(angle =45, hjust =1)) +scale_y_continuous(labels = scales::percent,limits =c(0,1))
Figure 6.4 -WAITING FOR NICK (estimates lightly off)
Show the code
# library(lavaan)# # # Define the path model with observed variables# model <- '# # Regression paths# antisemitismindex201 ~ biblebornchurchpca01 + pid3lean01 + bachelorsdegree + postgrad + white + male + age01 + famincmeansub01# antiintelpolpluspop401 ~ antisemitismindex201 + biblebornchurchpca01 + pid3lean01 + bachelorsdegree + postgrad + white + male + age01 + famincmeansub01# '# # # Fit the model# fit <- sem(model, data = data21a, cluster = "inputstate", sampling.weights = 'teamweight')# # # Get the summary of the model fit# summary(fit)#
