ATP Analysis.R

library(haven)
library(stargazer)
library(tidyverse)

#Obviously you'll need to download the ATP data (linked in the post) and put them in the right directory relative to this script.
atp47 <- read_sav("W47_Apr19/ATP W47.sav")
atp104 <- read_sav("W104_Mar22/ATP W104.sav")

#Convert the spectrum of agree/disagree to binary variable, and drop the non-responses
atp47 %>%
  mutate(F_AGECAT = as_factor(F_AGECAT),
         across(starts_with("ME_SPPRT"), ~case_when(.x == 99 ~ NA,
                                                    .x %in% 1:2 ~ 1,
                                                    .x %in% 3:4 ~ 0))) -> atp47

#get summary stats for each question
atp47 %>%
  as_survey_design(id = QKEY, weights = WEIGHT_W47) %>%
  group_by(F_AGECAT) %>%
  summarise(a = survey_mean(ME_SPPRTa_W47, na.rm = T),
            b = survey_mean(ME_SPPRTb_W47, na.rm = T),
            c = survey_mean(ME_SPPRTc_W47, na.rm = T),
            d = survey_mean(ME_SPPRTd_W47, na.rm = T)) %>%
  pivot_longer(cols = !F_AGECAT) %>%
  mutate(question = gsub("_se", "", name),
         statistic = ifelse(endsWith(name, "_se"), "se", "p")) %>%
  select(-name) %>%
  pivot_wider(names_from = statistic) %>%
  mutate(Year = 2019) -> atp47_me_responses

#same analysis, but for the 104 wave
atp104 %>%
  mutate(F_AGECAT = as_factor(F_AGECAT),
         across(starts_with("ME_SPPRT"), ~case_when(.x == 99 ~ NA,
                                                    .x %in% 1:2 ~ 1,
                                                    .x %in% 3:4 ~ 0))) -> atp104

atp104 %>%
  as_survey_design(id = QKEY, weights = WEIGHT_W104) %>%
  group_by(F_AGECAT) %>%
  summarise(a = survey_mean(ME_SPPRT_a_W104, na.rm = T),
            b = survey_mean(ME_SPPRT_b_W104, na.rm = T),
            c = survey_mean(ME_SPPRT_c_W104, na.rm = T),
            d = survey_mean(ME_SPPRT_d_W104, na.rm = T)) %>%
  pivot_longer(cols = !F_AGECAT) %>%
  mutate(question = gsub("_se", "", name),
         statistic = ifelse(endsWith(name, "_se"), "se", "p")) %>%
  select(-name) %>%
  pivot_wider(names_from = statistic) %>%
  mutate(Year = 2022) -> atp104_me_responses

#combine the two summary stats data frames to make the chart
bind_rows(atp47_me_responses, atp104_me_responses) %>%
  filter(!F_AGECAT %in% c('DK/REF', "Refused")) %>%
  mutate(question = case_when(question == 'a' ~ "The Israeli Government",
                              question == 'b' ~ "The Palestinian Government",
                              question == 'c' ~ "The Israeli People",
                              question == 'd' ~ "The Palestinian People")) %>%
  ggplot(aes(x=F_AGECAT, y = p, fill = factor(Year))) +
  geom_col(position = position_dodge(width = 0.5), width = 0.5, color = 'black') +
  geom_errorbar(aes(ymin = p - 1.96 * se, ymax = p + 1.96 * se), position = position_dodge(width = 0.5), width = 0.25) +
  facet_wrap(~question) +
  scale_y_continuous(labels = scales::percent, breaks = seq(0, .8, .2))+
  theme(panel.grid.major.x = element_blank(), 
        panel.border = element_rect(fill = NA),
        legend.position = 'bottom',
        legend.title = element_blank()) +
  labs(x = "Age Group", y = "Share Somewhat/Very Favorable", title = "Do you have a favorable or unfavorable opinion of _____")

ggsave("img/plot1.png", width = 8, height = 5, bg = 'white')

#Analyzing the four possible forms of support
atp104 %>%
  mutate(lean_gov = case_when(ME_SPPRT_a_W104 == 1 & ME_SPPRT_b_W104 == 0 ~ "Lean Israel",
                              ME_SPPRT_a_W104 == 0 & ME_SPPRT_b_W104 == 1 ~ "Lean Palestine",
                              ME_SPPRT_a_W104 == 1 & ME_SPPRT_b_W104 == 1 ~ "Equal favor",
                              ME_SPPRT_a_W104 == 0 & ME_SPPRT_b_W104 == 0 ~ "Equal disfavor"),
         lean_people = case_when(ME_SPPRT_c_W104 == 1 & ME_SPPRT_d_W104 == 0 ~ "Lean Israel",
                                 ME_SPPRT_c_W104 == 0 & ME_SPPRT_d_W104 == 1 ~ "Lean Palestine",
                                 ME_SPPRT_c_W104 == 1 & ME_SPPRT_d_W104 == 1 ~ "Equal favor",
                                 ME_SPPRT_c_W104 == 0 & ME_SPPRT_d_W104 == 0 ~ "Equal disfavor"),) %>%
  filter(!is.na(lean_gov),
         !is.na(lean_people),
         F_AGECAT != "Refused") %>%
  as_survey_design(id = QKEY, weights = WEIGHT_W104) -> atp104_lean_sd

#Is this tidy? No! But it works.
bind_rows(
  mutate(summarise(group_by(atp104_lean_sd, F_AGECAT, lean_gov), p = survey_prop()), Question = "Government"),
  mutate(summarise(group_by(atp104_lean_sd, F_AGECAT, lean_people), p = survey_prop()), Question = "People")  
  ) %>%
  mutate(lean = ifelse(is.na(lean_gov), lean_people, lean_gov)) %>%
  ggplot(aes(x=F_AGECAT, y=p, fill = lean)) + 
  geom_col(position = position_dodge(width = 1), color = 'black') + 
  geom_errorbar(aes(ymin = p - 1.96 * p_se, ymax = p + 1.96 * p_se), position = position_dodge(width = 1), width = 0.5) +
  facet_grid(rows = vars(Question), cols = vars(lean)) +
  scale_y_continuous(labels = scales::percent) +
  theme(panel.grid.major.x = element_blank(),
        panel.border = element_rect(fill = NA),
        legend.position = 'none',
        legend.title = element_blank()) + 
  labs(x = "Age Group", y = "", title = "Favorability Balance in 2022 ATP Responses")

ggsave("img/plot2.png", width = 8, height = 5, bg = 'white')

#Create the variable measuring support for Israel but not Palestine
atp47 <- filter(atp47, !F_AGECAT %in% c("DK/REF", "Refused")) %>% 
  mutate(Support_Only_Isr_Gov47 = ifelse(ME_SPPRTa_W47 - ME_SPPRTb_W47 == 1, 1, 0),
         Support_Only_Isr_People47 = ifelse(ME_SPPRTc_W47 - ME_SPPRTd_W47 == 1, 1, 0)) %>%
  #Also a rename to make the regression table easier to read
  rename(F_PARTY_FINAL = PARTY_W47) 

atp104 <- filter(atp104, !F_AGECAT == "Refused") %>%
  mutate(Support_Only_Isr_Gov104 = ifelse(ME_SPPRT_a_W104 - ME_SPPRT_b_W104 == 1, 1, 0),
         Support_Only_Isr_People104 = ifelse(ME_SPPRT_c_W104 - ME_SPPRT_d_W104 == 1, 1, 0))

#Generate the government regression table
stargazer(lm(Support_Only_Isr_Gov47 ~ F_AGECAT, atp47, weights = WEIGHT_W47),
          lm(Support_Only_Isr_Gov47 ~ F_AGECAT + factor(F_RACECMB) + factor(F_HISP), atp47, weights = WEIGHT_W47), 
          lm(Support_Only_Isr_Gov47 ~ F_AGECAT + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), atp47, weights = WEIGHT_W47),
          lm(Support_Only_Isr_Gov104 ~ F_AGECAT, atp104, weights = WEIGHT_W104),
          lm(Support_Only_Isr_Gov104 ~ F_AGECAT + factor(F_RACECMB) + factor(F_HISP), atp104, weights = WEIGHT_W104),
          lm(Support_Only_Isr_Gov104 ~ F_AGECAT  + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), atp104, weights = WEIGHT_W104), 
          covariate.labels = c("30-49", "50-64", "65+", "Black", "Asian", "Mixed", "Other Race", "Refused Race", "Hispanic", "Refused Ethnicity", "Conservative", "Moderate", "Liberal", "Very Liberal", "Refused Ideo", "Democrat", "Independent", 'Other Party', "Refused Party"), 
          dep.var.labels = c("Lean Israeli Gov't (2019)", "Lean Israeli Gov't (2022)"), 
          model.numbers = F, column.labels = paste("(", as.character(rep(seq(1, 3, 1), 2)), ")", sep = ""),
          out = 'table1.html', single.row = T, keep.stat = 'n', report = 'vc*'
)

#Generate the people regression table
stargazer(lm(Support_Only_Isr_People47 ~ F_AGECAT, atp47, weights = WEIGHT_W47),
          lm(Support_Only_Isr_People47 ~ F_AGECAT + factor(F_RACECMB) + factor(F_HISP), atp47, weights = WEIGHT_W47), 
          lm(Support_Only_Isr_People47 ~ F_AGECAT + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), atp47, weights = WEIGHT_W47),
          lm(Support_Only_Isr_People104 ~ F_AGECAT, atp104, weights = WEIGHT_W104),
          lm(Support_Only_Isr_People104 ~ F_AGECAT + factor(F_RACECMB) + factor(F_HISP), atp104, weights = WEIGHT_W104),
          lm(Support_Only_Isr_People104 ~ F_AGECAT  + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), atp104, weights = WEIGHT_W104), 
          covariate.labels = c("30-49", "50-64", "65+", "Black", "Asian", "Mixed", "Other Race", "Refused Race", "Hispanic", "Refused Ethnicity", "Conservative", "Moderate", "Liberal", "Very Liberal", "Refused Ideo", "Democrat", "Independent", 'Other Party', "Refused Party"), 
          dep.var.labels = c("Lean Israeli People (2019)", "Lean Israeli People (2022)"), 
          model.numbers = F, column.labels = paste("(", as.character(rep(seq(1, 3, 1), 2)), ")", sep = ""), 
          out = 'table2.html', single.row = T, keep.stat = 'n', report = 'vc*'
)

#Means of two time online measures, for diagnostic purposes
atp47 %>%
  as_survey_design(id = QKEY, weights = WEIGHT_W47) %>%
  summarise(TIMEONLINE1_W47 = survey_mean(TIMEONLINE1_W47, na.rm = T), 
            TIMEONLINE2_W47 = survey_mean(TIMEONLINE2_W47, na.rm = T))

#Collapse the two time online measures together
atp47 %>%
  mutate(TO = ifelse(is.na(TIMEONLINE1_W47), TIMEONLINE2_W47, TIMEONLINE1_W47),
         T_FLAG = ifelse(is.na(TIMEONLINE1_W47), 1, 0)) -> atp47_time

#Online regression
stargazer(lm(Support_Only_Isr_Gov47 ~ F_AGECAT + TIMEONLINE1_W47 + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47_time, weights =  WEIGHT_W47),
          lm(Support_Only_Isr_Gov47 ~ F_AGECAT + TIMEONLINE2_W47 + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47_time, weights =  WEIGHT_W47),
          lm(Support_Only_Isr_Gov47 ~ F_AGECAT + TO + T_FLAG + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47_time, weights =  WEIGHT_W47),
          lm(Support_Only_Isr_Gov47 ~ F_AGECAT + SITES1_W47 + SITES2_W47 + SITES3_W47 + SITES4_W47 + SITES5_W47 + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47, weights =  WEIGHT_W47),
          lm(Support_Only_Isr_People47 ~ F_AGECAT + TIMEONLINE1_W47 + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47_time, weights =  WEIGHT_W47),
          lm(Support_Only_Isr_People47 ~ F_AGECAT + TIMEONLINE2_W47 + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47_time, weights =  WEIGHT_W47),
          lm(Support_Only_Isr_People47 ~ F_AGECAT + TO + T_FLAG + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47_time, weights =  WEIGHT_W47),
          lm(Support_Only_Isr_People47 ~ F_AGECAT + SITES1_W47 + SITES2_W47 + SITES3_W47 + SITES4_W47 + SITES5_W47 + factor(F_RACECMB) + factor(F_HISP) + factor(F_IDEO) + factor(F_PARTY_FINAL), data = atp47, weights =  WEIGHT_W47),
          omit = c("RACECMB", "HISP", "IDEO", "PARTY"),
          covariate.labels = c("30-49", "50-64", "65+", "Time Online (No Email)", "Time Online", "Time Online (Combined)", "Time Question Choice Flag", "Used Facebook", "Used Instagram", "Used YouTube", "Used 'FizzyPress'", "Used 'Doromojo'"),
          dep.var.labels = c("Lean Israeli Government (2019)", "Lean Israeli People (2019)"),
          model.numbers = F, column.labels = as.character(rep(seq(1, 4, 1), 2)), 
          out = 'table3.html', single.row = T, keep.stat = 'n', report = 'vc*')