#---------------------------------------------------------------------------------#
### information: 
#---------------------------------------------------------------------------------#
### authors: björn bremer (bremer@mpifg.de), reto bürgisser (buergisser@ipz.uzh.ch)
### last updated: 12 December 2021

#---------------------------------------------------------------------------------#
### description
#---------------------------------------------------------------------------------#

### this file replicates the analysis of the following paper:
### "Do Citizens Care About Government Debt? Evidence From Survey Experiments on Budgetary Priorities"
### this file focuses on the results from the SPLIT SAMPLE EXPERIMENT shown in the MAIN TEXT

#---------------------------------------------------------------------------------#

set.seed(1234567)

#--------------------------------------------------------
### load libraries ###
#--------------------------------------------------------

library(grid)
library(gridExtra)
library(ggplot2)
library(ggpubr)
library(rstudioapi)
library(margins)
library(plyr)
library(dplyr)

#--------------------------------------------------------
### set working directory and load data ###
#--------------------------------------------------------

setwd(dirname(getActiveDocumentContext()$path)) # set working directory to the location of this file
load("../data/df_split_ejpr_replication.Rdata") # load the file

setwd("../figurestables/")

##########################################################################################
### figure 1a: predicted pverage support for fiscal consolidation by treatment, pooled ###
##########################################################################################

#--------------------------------------------------------
### calculate OLS regressions ###
#--------------------------------------------------------

## pension
debt_1 <- lm(debt ~ as.factor(split), data = df_split)
summary(debt_1)

debt_2 <- lm(debt ~ as.factor(split) + age + female + married + children + education_c + occup_gh + income + unemployed + retired + union + party_fam_rec, data = df_split)
summary(debt_2)

debt_3 <- lm(debt ~ as.factor(split) + age + female + married + children + education_c + occup_gh + income + unemployed + retired + union + party_fam_rec + country, data = df_split)
summary(debt_3)

#--------------------------------------------------------
### predict values ###
#--------------------------------------------------------

## debt
pre_debt_1 <- prediction(debt_1, df_split)
pre_debt_2 <- prediction(debt_2, df_split)
pre_debt_3 <- prediction(debt_3, df_split)

# calculate average predicted values by split
sum_debt_1 <- ddply(pre_debt_1, c("split"), summarise,
                    mean = mean(fitted, na.rm=TRUE),
                    mean_se = mean(se.fitted, na.rm=TRUE))
sum_debt_1$model <- "M1"

sum_debt_2 <- ddply(pre_debt_2, c("split"), summarise,
                    mean = mean(fitted, na.rm=TRUE),
                    mean_se = mean(se.fitted, na.rm=TRUE))
sum_debt_2$model <- "M2"

sum_debt_3 <- ddply(pre_debt_3, c("split"), summarise,
                    mean = mean(fitted, na.rm=TRUE),
                    mean_se = mean(se.fitted, na.rm=TRUE))
sum_debt_3$model <- "M3"

# combine all estimates in a single df
sum_debt_all <- rbind(sum_debt_1, sum_debt_2, sum_debt_3)

#--------------------------------------------------------
### plot the graph ###
#--------------------------------------------------------

plot_pre_debt <- ggplot(sum_debt_all,aes(as.factor(split), mean, colour = model)) +
   geom_errorbar(aes(ymin=mean-(1.96*mean_se), ymax=mean+(1.96*mean_se)), width=0, position=position_dodge(width=0.5)) + 
   geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
   scale_y_continuous(limits=c(3.5, 8.5), breaks = round(seq(3.5, 8.5, by = 1), 1)) +
   ylab("Mean support for consolidation (0-10 scale)") +
   theme_bw() +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  
   scale_x_discrete(name = "Treatment", breaks = c(1,2,3), labels = c("Control", "Lower spending", "Higher taxes")) +      
   scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE")) +
   scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE"), start = 0.2, end = 0.8, aesthetics = "colour")
plot_pre_debt
ggsave(plot_pre_debt, width = 14.5, height = 12, units = c("cm"), file ="figure1a.eps")


##########################################################################################
### figure 1b: predicted average support for fiscal consolidation by treatment, pooled ###
##########################################################################################

# create new debt binary support variable 
df_split_bin <- df_split[which (!is.na(df_split$split)),]

df_split_bin$debt_bin <- ifelse(df_split$debt <= 5, 0,
                                ifelse(df_split$debt > 5, 1, NA))

df_split_bin$debt_bin <- df_split_bin$debt_bin * 100 

# calculate means and se overall
mean_debt_bin <- ddply(df_split_bin, c("split"), summarise,
                       N = sum(!is.na(debt_bin)),
                       mean = mean(debt_bin,  na.rm=TRUE),
                       sd = sd(debt_bin,  na.rm=TRUE),
                       se = sd / sqrt(N)
)

mean_debt_bin$split <- as.factor(mean_debt_bin$split)

# bar plot with CIs overall
barplot_mean_debt_bin <- ggplot(mean_debt_bin, aes(x = split, y = mean)) +
   geom_hline(yintercept = 50, colour = "gray50") + 
   geom_bar(stat="identity", width=0.6) +
   geom_errorbar(aes(ymin=mean-(1.96*se), ymax=mean+(1.96*se)), width=0.3, color="grey20") + 
   scale_y_continuous(limits=c(0,80), breaks = round(seq(0, 80, by = 10), 1)) +
   ylab("Support for fiscal consolidation (in %)") +
   scale_x_discrete(name = "Treatment", labels = c("Control", "Lower spending", "Higher taxes")) +  
   theme_bw() +
   theme(axis.text.x = element_text(angle = 45, hjust = 1))  
barplot_mean_debt_bin
ggsave(barplot_mean_debt_bin, width = 12, height = 12, units = c("cm"), file ="figure1b.eps")


###########################################################################
### figure 2a: support for fiscal consolidation by trade-off and income ###
###########################################################################

# add labels
df_split_bin$income_brackets_three <- factor(df_split_bin$income_brackets_three,
                                         levels = c("1" , "2", "3"),
                                         labels = c("Low", "Middle", "High"))

# calculate means and se by income #
mean_debt_income_bin <- ddply(df_split_bin, c("split", "income_brackets_three"), summarise,
                              N = sum(!is.na(debt_bin)),
                              mean = mean(debt_bin,  na.rm=TRUE),
                              sd = sd(debt_bin,  na.rm=TRUE),
                              se = sd / sqrt(N)
)

mean_debt_income_bin$split <- as.factor(mean_debt_income_bin$split)

# bar plot with CIs by income
barplot_mean_debt_income_bin <- ggplot(mean_debt_income_bin[which (!is.na(mean_debt_income_bin$income_brackets_three)),],
                                       aes(split, mean, fill = income_brackets_three)) +
   geom_hline(yintercept = 50, colour = "gray50") + 
   geom_bar(aes(fill = income_brackets_three),  stat="identity", position=position_dodge(), colour="black", size=.3) +
   geom_errorbar(aes(ymin=mean-(1.41*se), ymax=mean+(1.41*se)), width=.4, position=position_dodge(.9)) +
   scale_y_continuous(limits=c(0,90), breaks = round(seq(0, 90, by = 10), 1)) +
   ylab("Support for fiscal consolidation (in %)") +
   scale_x_discrete(name = "Treatment", labels = c("Control", "Lower spending", "Higher taxes")) +  
   theme_bw() +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
   scale_fill_manual(values=c("gray20","gray50","gray80")) +
   labs(fill = "Income") + 
   theme(legend.position="bottom")
barplot_mean_debt_income_bin
ggsave(barplot_mean_debt_income_bin, width = 12, height = 12, units = c("cm"), file ="figure2a.eps")


##################################################################################
### figure 2b: support for fiscal consolidation by trade-off and partisanship  ###
##################################################################################

# calculate means and se by party_lr #
mean_debt_party_bin <- ddply(df_split_bin, c("split", "party_lr"), summarise,
                             N = sum(!is.na(debt_bin)),
                             mean = mean(debt_bin,  na.rm=TRUE),
                             sd = sd(debt_bin,  na.rm=TRUE),
                             se = sd / sqrt(N)
)

mean_debt_party_bin$split <- as.factor(mean_debt_party_bin$split)

# bar plot with CIs by partisanship (with abstention)
barplot_mean_debt_party_bin2 <- ggplot(mean_debt_party_bin[which (mean_debt_party_bin$party_lr == "Abstention"  |  mean_debt_party_bin$party_lr == "Left" | mean_debt_party_bin$party_lr == "Right"),],
                                       aes(split, mean, fill = party_lr)) +
   geom_hline(yintercept = 50, colour = "gray50") + 
   geom_bar(aes(fill = party_lr), stat="identity", position=position_dodge(), colour="black", size=.3) +
   geom_errorbar(aes(ymin=mean-(1.41*se), ymax=mean+(1.41*se)), width=.4, position=position_dodge(.9)) +
   scale_y_continuous(limits=c(0,90), breaks = round(seq(0, 90, by = 10), 1)) +
   ylab("Support for fiscal consolidation (in %)") +
   scale_x_discrete(name = "Treatment", labels = c("Control", "Lower spending", "Higher taxes")) +  
   theme_bw() +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
   scale_fill_manual(values=c("gray20","gray50","gray80")) +
   labs(fill = "Partisanship") + 
   theme(legend.position="bottom")
barplot_mean_debt_party_bin2
ggsave(barplot_mean_debt_party_bin2, width = 12, height = 12, units = c("cm"), file ="figure2b.eps")


############################################################################
### Figure 3: support for fiscal consolidation by trade-off and country  ###
############################################################################

# calculate means and se by country 
mean_debt_country_bin <- ddply(df_split_bin, c("split", "country"), summarise,
                               N = sum(!is.na(debt_bin)),
                               mean = mean(debt_bin,  na.rm=TRUE),
                               sd = sd(debt_bin,  na.rm=TRUE),
                               se = sd / sqrt(N)
)

mean_debt_country_bin$split <- as.factor(mean_debt_country_bin$split)

# bar plot with CIs by country
barplot_mean_debt_country_bin <- ggplot(mean_debt_country_bin, aes(split, mean, fill = country)) +
   geom_hline(yintercept = 50, colour = "gray50") + 
   geom_bar(aes(fill = country), stat="identity", position=position_dodge(), colour="black", size=.3) +
   geom_errorbar(aes(ymin=mean-(1.41*se), ymax=mean+(1.41*se)), width=.4, position=position_dodge(.9)) +
   scale_y_continuous(limits=c(0,90), breaks = round(seq(0,90, by = 10), 1)) +
   ylab("Support for fiscal consolidation (in %)") +
   scale_x_discrete(name = "Treatment", labels = c("Control", "Lower spending", "Higher taxes")) +  
   theme_bw() +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
   scale_fill_manual(values=c("gray20","gray40","gray60","gray80")) +
   labs(fill = "Country") + 
   theme(legend.position="bottom")
barplot_mean_debt_country_bin
ggsave(barplot_mean_debt_country_bin, width = 20, height = 12, units = c("cm"), file ="figure3.eps")
