reproduce_data_analysis.Rmd

---
title: "Reproduce Data Analysis"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
```

```{r load packages, message=FALSE, warning=FALSE}
# # Ensure that pacman is installed for package management and loading.
# if (!require("pacman")) install.packages("pacman")
# # for data reading wrangling and visualization
pacman::p_load(tidyverse) 
# for working directories
pacman::p_load(here) 
# # for cross tabulation and data cleaning
# pacman::p_load(janitor) 
# for working with strings
pacman::p_load(glue) 
# For randomized inference, also loads randomizr and estimatr 
pacman::p_load(ri2) 
# for marginal effects from lineal regressions
pacman::p_load(margins)
# Tests for linear regression models
pacman::p_load(lmtest)
pacman::p_load(car)
# Tables
pacman::p_load(kableExtra)
# for updated ggplot2 theme
pacman::p_load(hrbrthemes)
# for updated ggplot2 colorblind-friendly scheme
pacman::p_load(ggthemes)
# theme_set(hrbrthemes::theme_ipsum())
pacman::p_load(reshape2)
# for plotting of covariate balance
pacman::p_load(cobalt)
# for matching only
pacman::p_load(MatchIt)
library(stargazer)
library(usmap) 
library(ggplot2)  
library(maps)  
#install.packages("mapproj") 
library(mapproj)
#install.packages("viridis")
library(viridis)
library(dplyr)
library(ggthemes) 
library(sjPlot)
library(sjmisc)
library(here)
```

```{r download data}
df <- read.csv(here('', 'rla_clean_5_12.csv'))
dim(df)

```

```{r graph of income}
# make bar graph of income
p2 <- ggplot(data = subset(df, income_clean != "NA"), aes(income_clean)) + 
  geom_bar(binwidth=3, fill="#1c3ca8cf", color="#e9ecef", alpha=0.9) +
  ggtitle("Income of respondents") +
theme(
  plot.background = element_rect(fill = "white"),
  panel.background = element_rect(fill = "white"),
  axis.line.x = element_line(color = "grey")
  ) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), plot.title = element_text(size=10))

#ggsave("income.png", plot = p2)
p2

p3 <- ggplot(data = df, aes(treatment_group)) +
  geom_bar(binwidth=3, fill="#1c3ca8cf", color="#e9ecef", alpha=0.9) +
  ggtitle("Treatment groups of respondents") +
theme(
  plot.background = element_rect(fill = "white"),
  panel.background = element_rect(fill = "white"),
  axis.line.x = element_line(color = "grey")
  ) +
  xlab("Treatment Group") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), plot.title = element_text(size=10))

p3

#ggsave("treatment_group.png", plot = p3)
```

```{r map of respondents}
# load US state map data
us_states <- map_data("state")  
head(us_states)  

# aggregate number of respondents per state from your survey data
bystate <- df %>% 
  group_by(state) %>% 
  dplyr::summarise(
    numbperstate = n()
  )

# merge RLA data w/maps data

bystate$state <- tolower(bystate$state)

# fix errors
bystate[bystate$state=="new hampshir",]$state <- "new hampshire"
bystate[bystate$state=="minnesotta",]$state <- "minnesota"
bystate[bystate$state=="west virgina",]$state <- "west virginia"

us_states_df <- us_states %>% 
  left_join(bystate, by = c("region" = "state"))

map <- ggplot(data = us_states_df,
             mapping = aes(x = long, y = lat,
                           group = group, fill = numbperstate))+
  geom_polygon(color = "gray90", size = 0.1) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  scale_fill_continuous(type = "viridis", name = "Number of respondents \n per state") +
  theme_map() +
  theme(legend.position="bottom",
        axis.line=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank(),
        panel.background=element_blank(),
        panel.border=element_blank(),
        panel.grid=element_blank(),
        text = element_text(size=14))

map

#ggsave("map.png", plot = map)

```


```{r balance_plots}
df$age <- as.numeric(df$age)
# exclude data with missing covariates
no_na_df <- df %>% drop_na('age', 'state', 'income_clean', 'party', 'gender_female', 'parent_yes', 'race_hispanic', 'edu_4college', 'edu_hs')

love.plot(treatment_group_num ~ age + gender_female + parent_yes + state + income_clean + race_hispanic + edu_4college + edu_hs + dv_pre_state_conf + dv_pre_national_conf + party,
                 data = no_na_df, 
          binary = "std", 
          limits = c(-.5, .5),
          thresholds = c(m = .1)) # threshold determines the dotted lines on the graph

m.out <- matchit(dummy_treat ~ age + gender_female + parent_yes + state + income_clean + race_hispanic + edu_4college + edu_hs + dv_pre_state_conf + dv_pre_national_conf + party,
                 data = no_na_df, 
                method = "nearest", replace = TRUE) # perform matching
# plot covariate balance with and without matching
love.plot(m.out, binary = "std", thresholds = c(m = .1))

```

```{r compare tests with attention check failures and 2SD of duration}
# people within 2 SD of duration
sd2 = sd(df$duration_sec)*2

duration_sd2_df <- no_na_df %>% 
  filter(abs(duration_sec - mean(duration_sec)) < sd2)

# only removes 15 people so would not sig affect results
rm(duration_sd2_df)

# people who passed attention check
atten_check_df <- no_na_df %>% 
  filter(attention_pass == 1)

love.plot(treatment_group_num ~ age + gender_female + parent_yes + state + income_clean + race_hispanic + edu_4college + edu_hs + dv_pre_state_conf + dv_pre_national_conf + party,
                 data = atten_check_df, 
          binary = "std", 
          limits = c(-.5, .5),
          thresholds = c(m = .1)) # threshold determines the dotted lines on the graph

m.out <- matchit(dummy_treat ~ age + gender_female + parent_yes + state + income_clean + race_hispanic + edu_4college + edu_hs + dv_pre_state_conf + dv_pre_national_conf + party,
                 data = atten_check_df, 
                method = "nearest", replace = TRUE) # perform matching
# plot covariate balance with and without matching
love.plot(m.out, binary = "std", thresholds = c(m = .1))

# results look pretty similar for those who failed and did not fail the attention check
```

```{r lm_ftest_balance}
# regression of covariates on treatment assignment variable
balance_lm <- lm(treatment_group_num ~ age + gender_female + parent_yes + state + income_clean + race_hispanic + edu_4college + edu_hs + party + dv_pre_state_conf + dv_pre_national_conf, data = no_na_df) # factor(region) encodes the string variable as a factor for analysis


# Test whether all coefficients from the balancce_lm regression are equal to zero
# using heteroskedasticity-robust standard errors, denoted by hc2
car::linearHypothesis(balance_lm, c("gender_female = 0", "age = 0", "parent_yes = 0", "race_hispanic = 0", "edu_hs = 0", "dv_pre_state_conf = 0", "dv_pre_national_conf = 0"),
  test = "F", white.adjust = "hc2", singular.ok = TRUE)

summary(balance_lm)
```

```{r balance_table}
# Define a function to generate a balance check table using two-sample t-test 
# Parameters: "cov_list" includes covariates, "treat" indicates the treatment, "alpha" is the significance level
t_table <- function(data, cov_list, treat,
                    alpha = 0.05) {
  # for each covariate in cov_list, apply a function to conduct the t-test
  out <- lapply(cov_list, function(x) {
      tres <- t.test(data[[x]] ~ data[[treat]]) # t-test result between control and treatment groups' covariate x
    # capture the group means for both control and treatment, and the p-value
    c(mean_control = as.numeric(tres$estimate[1]), mean_treat = as.numeric(tres$estimate[2]), p_value = tres$p.value)    
  })
  
  # save results by binding the results for all above covariates in cov_list by row (rbind)
  out <- as.data.frame(do.call(rbind, out))
  # combine covaraite names and results by column (cbind)
  out <- cbind(covariate = cov_list, out)
  # get ride of non-important strings in the names
  names(out) <- gsub("[^0-9A-Za-z_]", "", names(out)) # gsub(pattern, replacement, x) replaces pattern in x with replacement
  
  # code presentation of p-value according to whether they are below a threshold, for e.g., 0.001
  out$p_value <- ifelse(out$p_value < 0.001,
    "<0.001",
    round(out$p_value, 3) # rounding of p values to 3 decimal places
  )
  return(out)
}

# Apply above t-t_table function to our dataset and selected covariates
balance_table <- t_table(
  data = df,
  c("gender_female", "parent_yes", "race_hispanic", "race_white", "race_asian", "race_black", "edu_hs", "edu_4college", "trust_federal", "trust_state"),
  "dummy_treat"
)

# output the balance table
balance_table %>% 
  knitr::kable() %>%
  kableExtra::kable_styling(bootstrap_options = "striped")
```

```{r region_balance}
# balance table for Democrats
balance_dems <- t_table(
  data = subset(df, party=="Democrat"), # subsetting data to those with "party" equal to "Democrat"
  c("gender_female", "parent_yes", "race_hispanic", "race_white", "race_asian", "race_black", "edu_hs", "edu_4college", "trust_federal", "trust_state"),
  "dummy_treat"
)

# balance table for Republicans
balance_repubs <- t_table(
  data = subset(df, party=="Republican"), # subsetting data to those with "party" equal to "Republican"
  c("gender_female", "parent_yes", "race_hispanic", "race_white", "race_asian", "race_black", "edu_hs", "edu_4college", "trust_federal", "trust_state"),
  "dummy_treat"
)

balance_indpt_other <- t_table(
  data = subset(df, party=="Independent" | party == "Other Party"), # subsetting data to those with "party" equal to "Republican"
  c("gender_female", "parent_yes", "race_hispanic", "race_white", "race_asian", "race_black", "edu_hs", "edu_4college", "trust_federal", "trust_state"),
  "dummy_treat"
)

# combine the 4 t-tables, excluding first column of covariate names for the latter 3, so that only one covariate name column is there
balance_party = cbind(balance_dems, balance_repubs %>% select(2:4), balance_indpt_other %>% select(2:4))

# Format: adding in headers
x <- knitr::kable(balance_party, digits = 2) %>% kable_styling()
# add in a header to label what we're cross-tabulating with
add_header_above(x, c('', 'Democrats'=3, 'Republicans'=3, 'Other and Independent'=3))


```

```{r basic regression with dummy treat}
df[c(6:10,12,19:21,23,27,29,30,32,35:42,44,45)] <- lapply(df[c(6:10,12,19:21,23,27,29,30,32,35:42,44,45)], factor)

atten_check_df$treatment_group <- relevel(factor(atten_check_df$treatment_group), ref = "Control")

# run regressions to estimate treatment effect with robust standard errors

# run the regression for state confidence level 
reg1 <- lm(dv_post_state_conf ~ dummy_treat + dv_pre_state_conf, data = df)

# run the regression for national confidence level
reg2 <- lm(dv_post_national_conf ~ dummy_treat + dv_pre_national_conf, data = df)

# run the regression with controls for state level
reg3 <- lm(dv_post_state_conf ~ dummy_treat + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)

# run the regression with controls for national level
reg4 <- lm(dv_post_national_conf ~ dummy_treat + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

# compare with those who did not fail attention check
reg1_att <- lm(dv_post_state_conf ~ dummy_treat + dv_pre_state_conf, data = atten_check_df)

reg2_att <- lm(dv_post_national_conf ~ dummy_treat + dv_pre_national_conf, data = atten_check_df)

reg3_att <- lm(dv_post_state_conf ~ dummy_treat + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = atten_check_df)

reg4_att <- lm(dv_post_national_conf ~ dummy_treat + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = atten_check_df)

# results do not vary by including those who did or did not fail the attention check

# compare by controlling for attention check
reg1_att_con <- lm(dv_post_state_conf ~ dummy_treat + dv_pre_state_conf + attention_pass, data = df)

# results look the same

# run the regression for national confidence level
reg2_att_con <- lm(dv_post_national_conf ~ dummy_treat + dv_pre_national_conf + attention_pass, data = df)

# results look different, treatment is statistic sig with control for attention pass

# run the regression with controls for state level
reg3_att_con <- lm(dv_post_state_conf ~ dummy_treat + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + attention_pass + trust_state, data = df)

# results look the same for state level w/ and w/out attention pass

# run the regression with controls for national level
reg4_att_con <- lm(dv_post_national_conf ~ dummy_treat + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + attention_pass + trust_federal, data = df)

# results look the same here too
```

```{r regression with treatment groups}
# run the regression with treatment groups
df$treatment_group <- relevel(factor(df$treatment_group), ref = "Control")

#state level with treatment groups
reg5 <- lm(dv_post_state_conf ~ treatment_group + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)

# national with treatment groups
reg6 <- lm(dv_post_national_conf ~ treatment_group + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

# compared with those who did not fail the attention check

reg5_att <- lm(dv_post_state_conf ~ treatment_group + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = atten_check_df)

# results are not different for those did and did not fail the attention check at state level 

reg6_att <- lm(dv_post_national_conf ~ treatment_group + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

# results do not vary much at the national level for those who fail the attention check

# compare controlling for the attention check

#state level with treatment groups
reg5_att_con <- lm(dv_post_state_conf ~ treatment_group + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + attention_pass + trust_state, data = df)

# looks the same at the state level

# national with treatment groups
reg6_att_con <- lm(dv_post_national_conf ~ treatment_group + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + attention_pass + trust_federal, data = df)

# looks the same at the national level

# we can get rid of the extra regressions now
rm(reg1_att, reg1_att_con,reg2_att, reg2_att_con,reg3_att, reg3_att_con, reg4_att, reg4_att_con,reg5_att, reg5_att_con,reg6_att, reg6_att_con)
rm(atten_check_df)

```

```{r secondary effect}
reg_state_seceffect <- lm(seceffect ~ treatment_group + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)

```

```{r treatment plot dummy_treat}
# # for the state level no controls
# lm_state_dummy1 <- bind_rows(list(tidy(reg1))) %>% 
#   filter((term %in% c("dummy_treat"))) 
# 
# lm_natl_dummy1 <- bind_rows(list(tidy(reg2))) %>% 
#   filter((term %in% c("dummy_treat"))) 
# 
# lm_state_dummy2 <- bind_rows(list(tidy(reg3))) %>% 
#   filter((term %in% c("dummy_treat"))) 
# 
# lm_natl_dummy2 <- bind_rows(list(tidy(reg4))) %>% 
#   filter((term %in% c("dummy_treat"))) 
# 
# combined_te <- bind_rows(lm_state_dummy1, lm_natl_dummy1, lm_state_dummy2, lm_natl_dummy2)
# 
# combined_te$intervention <- c("State Effect", "National Effect", "State Effect w/controls", "National effects w/controls")
# 
# ggplot(combined_te, aes(x = term,y = estimate)) +
#   facet_grid(cols = vars(intervention), scales = 'free_x', space = 'free_x') + # transform into facet grid by intervention
#   geom_point(position=position_dodge(width=0.5)) +
#   geom_errorbar(aes(x=term,ymin = estimate - 1.96 * std.error,ymax = estimate + 1.96 * std.error),width = .1,position=position_dodge(width=0.5)) + 
#   xlab('Treatment Effect') +
#   ylab('Estimate') +
#   theme(plot.title = element_text(hjust = 0.5),
#         axis.text.x = element_text(angle = 45, vjust = 1, hjust =1)) + 
#   geom_hline(yintercept = 0, colour = 'grey60', linetype = 2) +
#   coord_cartesian(ylim=c(-0.5, 0.5)) + 
#   ggtitle('Post-Test Treatment Effect (conditioned on pre-test)') +
#   theme(legend.position = "none",
#         panel.grid = element_blank(),
#         axis.text.x = element_blank())

```

```{r treatment plot treatment_group}
labs <- c("Bipartisan", "Hand Count", "Local", "Loser", "RLA Percentage", "Soup")

# for the state level with controls
lm_state <- bind_rows(list(tidy(reg5))) %>% 
  filter((term %in% c("treatment_groupControl", "treatment_groupBipartisan", "treatment_groupHandcount", "treatment_groupLocal", "treatment_groupLoser", "treatment_groupRL_percentage", "treatment_groupSoup"))) 

p4 <- ggplot(lm_state, aes(x = term,y =estimate)) +
  geom_point(position=position_dodge(width=0.5)) +
  geom_errorbar(aes(x=term,ymin = estimate - 1.96 * std.error,ymax = estimate + 1.96 * std.error),width = .1,position=position_dodge(width=0.5), color = "steelblue2") + 
  xlab('Treatment') +
  ylab('Estimate') +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust =1)) +
  geom_hline(yintercept = 0, colour = 'grey60', linetype = 2) + 
    coord_cartesian(ylim=c(-.75, .75)) + 
  ggtitle('Post-Test State Outcome w/ controls (conditioned on pre-test)') +
    scale_x_discrete(labels = labs)


# for the national level

lm_national <- bind_rows(list(tidy(reg6))) %>% 
  filter((term %in% c("treatment_groupControl", "treatment_groupBipartisan", "treatment_groupHandcount", "treatment_groupLocal", "treatment_groupLoser", "treatment_groupRL_percentage", "treatment_groupSoup"))) 

# plot treatment effect
p5 <- ggplot(lm_national, aes(x = term,y =estimate)) +
  geom_point(position=position_dodge(width=0.5)) +
  geom_errorbar(aes(x=term,ymin = estimate - 1.96 * std.error,ymax = estimate + 1.96 * std.error),width = .1,position=position_dodge(width=0.5), color = "firebrick") + 
  xlab('Treatment') +
  ylab('Estimate') +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust =1)) + 
  geom_hline(yintercept = 0, colour = 'grey60', linetype = 2) +
  coord_cartesian(ylim=c(-.75, .75)) + 
  ggtitle('Post-Test National Outcome w/ controls (conditioned on pre-test)') +
    scale_x_discrete(labels = labs)

# for the state level secondary effect with controls
lm_state_sec <- bind_rows(list(tidy(reg_state_seceffect))) %>% 
  filter((term %in% c("treatment_groupControl", "treatment_groupBipartisan", "treatment_groupHandcount", "treatment_groupLocal", "treatment_groupLoser", "treatment_groupRL_percentage", "treatment_groupSoup"))) 

sec_effect <- ggplot(lm_state_sec, aes(x = term,y =estimate)) +
  geom_point(position=position_dodge(width=0.5)) +
  geom_errorbar(aes(x=term,ymin = estimate - 1.96 * std.error,ymax = estimate + 1.96 * std.error),width = .1,position=position_dodge(width=0.5), color = "steelblue2") + 
  xlab('Treatment') +
  ylab('Estimate') +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust =1)) +
  geom_hline(yintercept = 0, colour = 'grey60', linetype = 2) + 
    coord_cartesian(ylim=c(-.5, 1)) + 
  ggtitle('Secondary Effect State Outcome w/ controls (conditioned on pre-test)') +
  scale_x_discrete(labels = labs)

p4
p5
sec_effect
# ggsave("state_treat.png", plot = p4)
# ggsave("national_treat.png", plot = p5)
# ggsave("sec_effect.png", plot = sec_effect)
```

```{r difference_in_means_simple}
# simple difference in means for those treated with a message
with(df, mean(dv_post_state_conf[dummy_treat == 1]) - mean(dv_post_state_conf[dummy_treat == 0])) #with(data, expr, …) evaluates the expr on the data

with(df, mean(dv_post_national_conf[dummy_treat == 1]) - mean(dv_post_national_conf[dummy_treat == 0])) #with(data, expr, …) evaluates the expr on the data
```

```{r difference_in_means}
# Use the estimatr package for difference in means estimates
estimatr::difference_in_means(dv_state_treatment_diff ~ dummy_treat, data = df) 
# difference_in_means(y ~ x, data) computes mean(y when x==1) - mean(y when x==0) in data, along with standard errors and p-values from two-sided t-tests

estimatr::difference_in_means(dv_national_treatment_diff ~ dummy_treat, data = df) 

```

```{r difference_in_means_subsample}
# difference in means for women at state level

df$gender_female <- to_factor(df$gender_female)

estimatr::difference_in_means(dv_state_treatment_diff ~ dummy_treat, data = df, subset = gender_female == "1")

# difference in means for women at national level
estimatr::difference_in_means(dv_national_treatment_diff ~ dummy_treat, data = df, subset = gender_female == "1")

# difference in means for age at state level
df$age <- to_numeric(df$age)
# age < 40
estimatr::difference_in_means(dv_state_treatment_diff ~ dummy_treat, data = df, subset = age < 40)
# age > 40
estimatr::difference_in_means(dv_state_treatment_diff ~ dummy_treat, data = df, subset = age > 40)

# difference in means for age at national level
# age > 40
estimatr::difference_in_means(dv_national_treatment_diff ~ dummy_treat, data = df, subset = age > 40)

# difference in means for HS education at state level
estimatr::difference_in_means(dv_state_treatment_diff ~ dummy_treat, data = df, subset = edu_hs == "1")

# difference in means for HS education at national level
estimatr::difference_in_means(dv_national_treatment_diff ~ dummy_treat, data = df, subset = edu_hs== "1")

# difference in means for college education at state level
estimatr::difference_in_means(dv_state_treatment_diff ~ dummy_treat, data = df, subset = edu_4college == "1")

# difference in means for college education at national level
estimatr::difference_in_means(dv_national_treatment_diff ~ dummy_treat, data = df, subset = edu_4college == "1")

```

```{r gender interaction}
# Gender interaction

#state level interaction for gender
reg_state_fem <- lm(dv_post_state_conf ~ treatment_group*gender_female + dv_pre_state_conf + age + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)


# national level interaction for gender
reg_natl_fem <- lm(dv_post_national_conf ~ treatment_group*gender_female + dv_pre_national_conf + age + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

```

```{r age interaction}
# Age interaction

#state level 
reg_state_age <- lm(dv_post_state_conf ~ treatment_group*age + gender_female + dv_pre_state_conf + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)


# national level 
reg_natl_age <- lm(dv_post_national_conf ~ treatment_group*age + gender_female + dv_pre_national_conf + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

```

```{r educ interaction}
# Education interaction

# state only HS 
reg_state_edu_hs <- lm(dv_post_state_conf ~ treatment_group*edu_hs + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + party + trust_state, data = df)

# national only HS
reg_natl_edu_hs <- lm(dv_post_national_conf ~ treatment_group*edu_hs + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + party + trust_federal, data = df)

# state college or above
reg_state_edu <- lm(dv_post_state_conf ~ treatment_group*edu_4college + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_hs + party + trust_state, data = df)

# national college or above
reg_natl_edu <- lm(dv_post_national_conf ~ treatment_group*edu_4college + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_hs + party + trust_federal, data = df)

```

```{r income interaction}
# income interaction

# state 
reg_state_inc <- lm(dv_post_state_conf ~ treatment_group*income_num +  dv_pre_state_conf + age + gender_female + parent_yes + state + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)


# national
reg_natl_inc <- lm(dv_post_national_conf ~ treatment_group*income_num + dv_pre_national_conf + age + gender_female + parent_yes + state + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

```

```{r always_never believer interaction}
# always vs never believer interaction
df$always_believer <- to_factor(df$always_believer)
df$never_believer <- to_factor(df$never_believer)

# state never believer
reg_state_never <- lm(dv_post_state_conf ~ treatment_group*never_believer + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)

# national never believer
reg_natl_never <- lm(dv_post_national_conf ~ treatment_group*never_believer + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

# state always believer
reg_state_always <- lm(dv_post_state_conf ~ treatment_group*always_believer + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + trust_state, data = df)

# national always believer
reg_natl_always <- lm(dv_post_national_conf ~ treatment_group*always_believer + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)
```

```{r political party interaction}
# political party interaction

# state level
df$party <- to_factor(df$party)

reg_state_party <- lm(dv_post_state_conf ~ treatment_group*party + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + trust_state, data = df)

# national level
reg_natl_party <- lm(dv_post_national_conf ~ treatment_group*party +  dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + trust_federal, data = df)

```

```{r race interaction terms}
# race interaction
# state level race_hispanic Y/N

reg_state_hisp <- lm(dv_post_state_conf ~ treatment_group*race_hispanic + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)

# national level race_hispanic Y/N
reg_natl_hisp <- lm(dv_post_national_conf ~ treatment_group*race_hispanic + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_white + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

# state level race_white
reg_state_white <- lm(dv_post_state_conf ~ treatment_group*race_white + race_hispanic + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_black + race_asian + edu_4college + edu_hs + party + trust_state, data = df)

# national level race_white 
reg_natl_white <- lm(dv_post_national_conf ~ treatment_group*race_white + race_hispanic + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_black + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

# state level race_black
reg_state_black <- lm(dv_post_state_conf ~ treatment_group*race_black + race_hispanic + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_white + race_asian + edu_4college + edu_hs + party + trust_state, data = df)

# national level race_black
reg_natl_black <- lm(dv_post_national_conf ~ treatment_group*race_black + race_hispanic + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_white + race_asian + edu_4college + edu_hs + party + trust_federal, data = df)

```

```{r filter out extreme views pre test}
# state level middle views ( > 1 and < 10)
reg_middle_views_state <- lm(dv_post_state_conf ~ treatment_group + dv_pre_state_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party, data = df, subset = (dv_pre_state_conf < 10 & dv_pre_state_conf > 1))

# national middle views ( > 1 and < 10)
reg_middle_views_natl <- lm(dv_post_national_conf ~ treatment_group + dv_pre_national_conf + age + gender_female + parent_yes + state + income_clean + race_hispanic + race_white + race_black + race_asian + edu_4college + edu_hs + party, data = df, subset = (dv_pre_state_conf < 10 & dv_pre_state_conf > 1))

```

```{r visuals on party politics}
# Visual on trust in state government broken down by political party

state_trust = df %>% group_by(trust_state, party) %>% summarise(number = n())
state_trust = subset(state_trust, state_trust$party == "Democrat" | state_trust$party == "Republican")
state_trust = state_trust[-c(11,12),]

state_trust$percent = ifelse(state_trust$party == "Democrat", (state_trust$number/nrow(df[df$party == "Democrat", ])*100), state_trust$number/nrow(df[df$party == "Republican",])*100)

state_trust$percent = round(state_trust$percent, 0)

p6 <- ggplot(state_trust, aes(x=trust_state, y = percent, fill = party))+
  geom_bar(position = "Dodge", stat = "identity", width = 0.5) +
  ggtitle("Trust in State Government by Political Party")+
  theme_bw(base_size =16)+
  geom_text(aes(label=percent), family="serif", position= position_dodge(.5),vjust=-1.25, size=6)+
  theme(axis.text=element_text(size=16, family="serif"),
        axis.title=element_text(size=16, family="serif"),
        legend.title = element_blank(),
        legend.text= element_text(size=16, family="serif"), legend.position = "bottom")+
  ylim(0,100)+
  ylab("Percentage")+
  xlab("Trust")+
  scale_fill_manual(breaks = c("Democrat","Republican"), values = c("blue", "red"))

# Visual for trust in federal government based on political party

federal_trust = df %>% group_by(trust_federal, party) %>% summarise(number = n())
federal_trust = subset(federal_trust, federal_trust$party == "Democrat" | federal_trust$party == "Republican")
federal_trust = federal_trust[-c(11,12),]

federal_trust$percent = ifelse(federal_trust$party == "Democrat", federal_trust$number/nrow(df[df$party == "Democrat", ])*100, federal_trust$number/nrow(df[df$party == "Republican",])*100)

federal_trust$percent = round(federal_trust$percent, 0)

p7 <- ggplot(federal_trust, aes(x=trust_federal, y = percent, fill = party))+
  geom_bar(position = "Dodge", stat = "identity", width = 0.5) +
  ggtitle("Trust in Federal Government by Political Party")+
  theme_bw(base_size =16)+
  geom_text(aes(label=percent), family="serif", position= position_dodge(.5),vjust=-1.25, size=6)+
  theme(axis.text=element_text(size=16, family="serif"),
        axis.title=element_text(size=16, family="serif"),
        legend.title = element_blank(),
        legend.text= element_text(size=16, family="serif"), legend.position = "bottom")+
  ylim(0,100)+
  ylab("Percentage")+
  xlab("Trust")+
  scale_fill_manual(breaks = c("Democrat","Republican"), values = c("blue", "red"))

```

```{r visuals on always/never believer confidence}

# Build visual for always believers confidence

always_state_pre_df = df %>% group_by(always_believer, dv_pre_state_conf) %>% summarise(number = n())
always_state_post_df = df %>% group_by(always_believer, dv_post_state_conf) %>% summarise(number = n())

always_state_pre_df = always_state_pre_df[-c(1:10),]
always_state_post_df = always_state_post_df[-c(1:10),]
colnames(always_state_pre_df)[2] = c("confidence")
colnames(always_state_post_df)[2] = c("confidence")


always_state_pre_df$percent = round(always_state_pre_df$number/nrow(df[df$always_believer == 1,])*100,2)
always_state_post_df$percent = round(always_state_post_df$number/nrow(df[df$always_believer == 1,])*100,2)

p8 <- ggplot(mapping = aes(x = confidence))+
  geom_bar(data = always_state_pre_df, aes(x = confidence-0.2, y = percent, fill = "Confidence Pre-Treatment"), stat = "identity", width = 0.2) +
  geom_bar(data = always_state_post_df, aes(x = confidence, y = percent, fill = "Confidence Post-Treatment"), stat = "identity", width = 0.2) +
  ggtitle("Distribution of Confidence in Always Believers")+
  theme_bw(base_size =16)+
  #geom_text(aes(label=percent), family="serif", position= position_dodge(.5),vjust=-1.25, size=6)+
  theme(axis.text=element_text(size=16, family="serif"),
        axis.title=element_text(size=16, family="serif"),
        legend.title = element_blank(),
        legend.text= element_text(size=16, family="serif"), legend.position = "bottom")+
  ylim(0,50)+
  ylab("Percentage")+
  xlab("Confidence in Vote Count")+
  scale_fill_manual(values = c("black", "red"))

# Build visual for never believers confidence

never_state_pre_df = df %>% group_by(never_believer, dv_pre_state_conf) %>% summarise(number = n())
never_state_post_df = df %>% group_by(never_believer, dv_post_state_conf) %>% summarise(number = n())

never_state_pre_df = never_state_pre_df[-c(1:10),]
never_state_post_df = never_state_post_df[-c(1:10),]
colnames(never_state_pre_df)[2] = c("confidence")
colnames(never_state_post_df)[2] = c("confidence")


never_state_pre_df$percent = round(never_state_pre_df$number/nrow(df[df$never_believer == 1,])*100,2)
never_state_post_df$percent = round(never_state_post_df$number/nrow(df[df$never_believer == 1,])*100,2)

p9 <- ggplot(mapping = aes(x = confidence))+
  geom_bar(data = never_state_pre_df, aes(x = confidence-0.2, y = percent, fill = "Confidence Pre-Treatment"), stat = "identity", width = 0.2) +
  geom_bar(data = never_state_post_df, aes(x = confidence, y = percent, fill = "Confidence Post-Treatment"), stat = "identity", width = 0.2) +
  ggtitle("Distribution of Confidence in Never Believers")+
  theme_bw(base_size =16)+
  #geom_text(aes(label=percent), family="serif", position= position_dodge(.5),vjust=-1.25, size=6)+
  theme(axis.text=element_text(size=16, family="serif"),
        axis.title=element_text(size=16, family="serif"),
        legend.title = element_blank(),
        legend.text= element_text(size=16, family="serif"), legend.position = "bottom")+
  ylim(0,50)+
  ylab("Percentage")+
  xlab("Confidence in Vote Count")+
  scale_fill_manual(values = c("black", "red"))
```

```{r visuals on always/never beleiver and party}

belief_by_party_df = df %>% group_by(accuracy_2016, party) %>% summarise(number2016 = n())
belief_party_2 = df %>% group_by(accuracy_2020, party) %>% summarise(number2020 = n())
belief_party_2 = belief_party_2[-c(1:3), ]

belief_by_party_df = merge(belief_by_party_df, belief_party_2, by = "party")
belief_by_party_df = subset(belief_by_party_df, belief_by_party_df$party == "Democrat" | belief_by_party_df$party == "Republican")
  
state_trust = df %>% group_by(trust_state, party) %>% summarise(number = n())
state_trust = subset(state_trust, state_trust$party == "Democrat" | state_trust$party == "Republican")
state_trust = state_trust[-c(11,12),]

state_trust$percent = ifelse(state_trust$party == "Democrat", (state_trust$number/nrow(df[df$party == "Democrat", ])*100), state_trust$number/nrow(df[df$party == "Republican",])*100)

state_trust$percent = round(state_trust$percent, 0)

```

```{r histograms of confidence for overall population}
# state confidence levels

p10 <- ggplot()+
  geom_histogram(data = df, aes(x = dv_pre_state_conf-0.2, color = "Confidence Pre-Treatment"), alpha = 0.2, binwidth = 0.2) +
  geom_histogram(data = df, aes(x = dv_post_state_conf, color = "Confidence Post-Treatment"),  alpha = 0.2, binwidth = 0.2) +
  ggtitle("Confidence in Vote Count (State Level)")+
  theme_bw(base_size =16)+
  #geom_text(aes(label=percent), family="serif", position= position_dodge(.5),vjust=-1.25, size=6)+
  theme(axis.text=element_text(size=16, family="serif"),
        axis.title=element_text(size=16, family="serif"),
        legend.title = element_blank(),
        legend.text= element_text(size=16, family="serif"), legend.position = "bottom")+
  ylab("Frequency")+
  xlab("Confidence")+
  scale_color_manual(values = c("black", "orange"))

# national confidence levels

p11 <- ggplot()+
  geom_histogram(data = df, aes(x = dv_pre_national_conf-0.2, color = "Confidence Pre-Treatment"), alpha = 0.2, binwidth = 0.2) +
  geom_histogram(data = df, aes(x = dv_post_national_conf, color = "Confidence Post-Treatment"),  alpha = 0.2, binwidth = 0.2) +
  ggtitle("Confidence in Vote Count (National Level)")+
  theme_bw(base_size =16)+
  #geom_text(aes(label=percent), family="serif", position= position_dodge(.5),vjust=-1.25, size=6)+
  theme(axis.text=element_text(size=16, family="serif"),
        axis.title=element_text(size=16, family="serif"),
        legend.title = element_blank(),
        legend.text= element_text(size=16, family="serif"), legend.position = "bottom")+
  ylab("Frequency")+
  xlab("Confidence")+
  scale_color_manual(values = c("black", "orange"))

```

```{r summary table of treatment effects}

# create subsets for each treatment

control_df = subset(df, df$treatment_group == "Control")
bipartisan_df = subset(df, df$treatment_group == "Bipartisan")
handcount_df = subset(df, df$treatment_group == "Handcount")
local_df = subset(df, df$treatment_group == "Local")
loser_df = subset(df, df$treatment_group == "Loser")
rla_df = subset(df, df$treatment_group == "RL_percentage")
soup_df = subset(df, df$treatment_group == "Soup")

# get average treatment effect for each treatment group for state elections
avg_state_control = round(mean(control_df$dv_state_treatment_diff), 2)
avg_state_bipartisan = round(mean(bipartisan_df$dv_state_treatment_diff), 2)
avg_state_handcount = round(mean(handcount_df$dv_state_treatment_diff), 2)
avg_state_local = round(mean(local_df$dv_state_treatment_diff), 2)
avg_state_loser = round(mean(loser_df$dv_state_treatment_diff), 2)
avg_state_rla = round(mean(rla_df$dv_state_treatment_diff), 2)
avg_state_soup = round(mean(soup_df$dv_state_treatment_diff), 2)

# get average treatment effect for each treatment group for national elections
avg_national_control = round(mean(control_df$dv_national_treatment_diff),2)
avg_national_bipartisan = round(mean(bipartisan_df$dv_national_treatment_diff) ,2)
avg_national_handcount = round(mean(handcount_df$dv_national_treatment_diff),2)
avg_national_local = round(mean(local_df$dv_national_treatment_diff),2)
avg_national_loser = round(mean(loser_df$dv_national_treatment_diff),2)
avg_national_rla = round(mean(rla_df$dv_national_treatment_diff),2)
avg_national_soup = round(mean(soup_df$dv_national_treatment_diff),2)

# get percent increase for each treatment group for state elections

inc_state_control = round(nrow(control_df[control_df$dv_state_treatment_diff > 0, ])/nrow(control_df)*100, 2)
inc_state_bipartisan = round(nrow(bipartisan_df[bipartisan_df$dv_state_treatment_diff > 0, ])/nrow(bipartisan_df)*100, 2)
inc_state_handcount = round(nrow(handcount_df[handcount_df$dv_state_treatment_diff > 0, ])/nrow(handcount_df)*100, 2)
inc_state_local = round(nrow(local_df[local_df$dv_state_treatment_diff > 0, ])/nrow(local_df)*100, 2)
inc_state_loser = round(nrow(loser_df[loser_df$dv_state_treatment_diff > 0, ])/nrow(loser_df)*100, 2)
inc_state_rla = round(nrow(rla_df[rla_df$dv_state_treatment_diff > 0, ])/nrow(rla_df)*100, 2)
inc_state_soup = round(nrow(soup_df[soup_df$dv_state_treatment_diff > 0, ])/nrow(soup_df)*100, 2)

# get percent increase for each treatment group for national elections

inc_national_control = round(nrow(control_df[control_df$dv_national_treatment_diff > 0, ])/nrow(control_df)*100, 2)
inc_national_bipartisan = round(nrow(bipartisan_df[bipartisan_df$dv_national_treatment_diff > 0, ])/nrow(bipartisan_df)*100, 2)
inc_national_handcount = round(nrow(handcount_df[handcount_df$dv_national_treatment_diff > 0, ])/nrow(handcount_df)*100, 2)
inc_national_local = round(nrow(local_df[local_df$dv_national_treatment_diff > 0, ])/nrow(local_df)*100, 2)
inc_national_loser = round(nrow(loser_df[loser_df$dv_national_treatment_diff > 0, ])/nrow(loser_df)*100, 2)
inc_national_rla = round(nrow(rla_df[rla_df$dv_national_treatment_diff > 0, ])/nrow(rla_df)*100, 2)
inc_national_soup = round(nrow(soup_df[soup_df$dv_national_treatment_diff > 0, ])/nrow(soup_df)*100, 2)

# get percent decrease for each treatment group for state elections

dec_state_control = round(nrow(control_df[control_df$dv_state_treatment_diff < 0, ])/nrow(control_df)*100, 2)
dec_state_bipartisan = round(nrow(bipartisan_df[bipartisan_df$dv_state_treatment_diff < 0, ])/nrow(bipartisan_df)*100, 2)
dec_state_handcount = round(nrow(handcount_df[handcount_df$dv_state_treatment_diff < 0, ])/nrow(handcount_df)*100, 2)
dec_state_local = round(nrow(local_df[local_df$dv_state_treatment_diff < 0, ])/nrow(local_df)*100, 2)
dec_state_loser = round(nrow(loser_df[loser_df$dv_state_treatment_diff < 0, ])/nrow(loser_df)*100, 2)
dec_state_rla = round(nrow(rla_df[rla_df$dv_state_treatment_diff < 0, ])/nrow(rla_df)*100, 2)
dec_state_soup = round(nrow(soup_df[soup_df$dv_state_treatment_diff < 0, ])/nrow(soup_df)*100, 2)

# get percent decrease for each treatment group for national elections

dec_national_control = round(nrow(control_df[control_df$dv_national_treatment_diff < 0, ])/nrow(control_df)*100, 2)
dec_national_bipartisan = round(nrow(bipartisan_df[bipartisan_df$dv_national_treatment_diff < 0, ])/nrow(bipartisan_df)*100, 2)
dec_national_handcount = round(nrow(handcount_df[handcount_df$dv_national_treatment_diff < 0, ])/nrow(handcount_df)*100, 2)
dec_national_local = round(nrow(local_df[local_df$dv_national_treatment_diff < 0, ])/nrow(local_df)*100, 2)
dec_national_loser = round(nrow(loser_df[loser_df$dv_national_treatment_diff < 0, ])/nrow(loser_df)*100, 2)
dec_national_rla = round(nrow(rla_df[rla_df$dv_national_treatment_diff < 0, ])/nrow(rla_df)*100, 2)
dec_national_soup = round(nrow(soup_df[soup_df$dv_national_treatment_diff < 0, ])/nrow(soup_df)*100, 2)

Treatment = c("Number Treated", "Average State Effect", "Percent Increase for State", "Percent Decrease for State", "Average National Effect", "Percent Increase for National", "Percent Decrease for National")

Control = c(nrow(control_df), avg_state_control, inc_state_control, dec_state_control, avg_national_control, inc_national_control, dec_national_control)

Bipartisan = c(nrow(bipartisan_df), avg_state_bipartisan, inc_state_bipartisan, dec_state_bipartisan, avg_national_bipartisan, inc_national_bipartisan, dec_national_bipartisan)

Handcount = c(nrow(handcount_df), avg_state_handcount, inc_state_handcount, dec_state_handcount, avg_national_handcount, inc_national_handcount, dec_national_handcount)

Local = c(nrow(local_df), avg_state_local, inc_state_local, dec_state_local, avg_national_local, inc_national_local, dec_national_local)

Loser = c(nrow(loser_df), avg_state_loser, inc_state_loser, dec_state_loser, avg_national_loser, inc_national_loser, dec_national_loser)

Risk_Limit = c(nrow(rla_df), avg_state_rla, inc_state_rla, dec_state_rla, avg_national_rla, inc_national_rla, dec_national_rla)

Soup = c(nrow(soup_df), avg_state_soup, inc_state_soup, dec_state_soup, avg_national_soup, inc_national_soup, dec_national_soup)

summary_table_df = rbind(Treatment, Control, Bipartisan, Handcount, Local, Loser, Risk_Limit, Soup)
stargazer(summary_table_df, summary = FALSE )
```

```{r age visualization}
# visualize age distribution

p10 <- ggplot()+
  geom_histogram(data = df, aes(x = age), binwidth = 1, fill = "#1c3ca8cf") +
  scale_fill_viridis() +
  #ggtitle("Age of Respondents") +
theme(
  plot.background = element_rect(fill = "white"),
  panel.background = element_rect(fill = "white"),
  axis.line.x = element_line(color = "grey")
  ) +
  theme(axis.text.x = element_text(vjust = 0.5, hjust=1), plot.title = element_text(size=10)) +
  xlab("Age") +
  geom_vline(color="red",linetype=2, xintercept = mean(df$age,na.rm=T)) + 
  geom_text(aes(x = 50, y = 140, label = "Mean = 44", color = ))

p10 

#ggsave("age.png", plot = p10, width = 3, height = 2, units = "in")

```