## title: "Swiss_Analysis"
## date: "2022-12-16"

library(texreg)
library(readxl)
library(ggplot2)
library(dplyr)
library(devtools)
library(psy)
require(psych)
library(data.table)
library(nlme)
library(visreg)
library(apaTables)
library(ncf)
library(pgirmess)
library(ape)
library(plm)
library(visreg)
library(readstata13)
library(estimatr)
library(texreg)
library(robustbase)
library(readstata13)
library(reghelper)
library(extrafont)
rm(list=ls())


## Set Directory
setwd("...")


## Open file and generate new variables
uschdata1 <- read.dta13("Swiss.dta")

uschdata1$VotingResult2 <- uschdata1$VotingResult2*100 ## to move from percentages to absolute numbers
uschdata1$RFDummyDiff_2 <- 2-uschdata1$RFDummyDiff ## recoding so that 1 = RFD; 0 = no RFD


## Figure 2
uschdata1$RFDummyDiff_2 <- factor(uschdata1$RFDummyDiff_2, levels = c(0,1), labels = c("No Regulatory Focus Difference", "Regulatory Focus Difference"))

# Remove the referendum with the miscoded dew point and relative humidity
uschdata1 <- uschdata1[uschdata1$Election != 50,]

random <- plm(VotingResult2 ~ windspeed*RFDummyDiff_2 + Turnout + Pressure_c + Temp + DewPoint + Precipitation + CloudCover_c + Unemployment + Education + Rural + Age + Income + Race, data=uschdata1, index=c("State", "Election"), model="random")
summary(random)
loadfonts(quiet = TRUE)
pdf(file="Wind_Voting_Figure2.pdf", width = 6, height = 4, family = "Times")
visreg(random, "windspeed", "RFDummyDiff_2", jitter = T, xlab = "Wind Speed (in km/h)", ylab = "Prevention-Oriented Campaign Vote (in %)")
dev.off()

# Remove each remaining referendum in turn and see what that does to the interaction coefficient.
# Removing "Election 50" (referendum 24) at the end doesn't do anything as we already removed it,
#  but it allows us to show the coefficient for that case (N=23) in the table.
cat("Removed\tB\tSE\tp\n")
for (ref in 27:50) {
  df <- uschdata1[uschdata1$Election != ref,]
  m22 <- plm(VotingResult2 ~ windspeed*RFDummyDiff_2 + Turnout + Pressure_c + Temp + DewPoint + Precipitation + CloudCover_c + Unemployment + Education + Rural + Age + Income + Race, data=df, index=c("State", "Election"), model="random")
  coef <- (summary(m22)$coefficients)
  b <- coef[16, 1]
  se <- coef[16, 2]
  p <- coef[16, 4]
  cat(ref, "\t", sprintf("%.3f", b), "\t", sprintf("%.3f", se), "\t", sprintf("%.3f", p), "\n", sep="")
}