#####################################################################
#                                                                   #
# Electoral systems and party systems in Germany on the local level #
#                                                                   #
#                         Replication Script                        #
#                                                                   #
#####################################################################

getwd()
setwd("")
.libPaths("...")

library(openxlsx)
library(ggpubr)
library(sjPlot)
library(lme4)

#Access data
data <- read.xlsx("...\\ESPSGerLocLvl.xlsx", sheet = "Results")


#Recode Länder 
data$Land_c <- 0
data$Land_c[data$Land==1] <- "BW"
data$Land_c[data$Land==2] <- "BY"
data$Land_c[data$Land==3] <- "BE"
data$Land_c[data$Land==4] <- "BB"
data$Land_c[data$Land==5] <- "HB"
data$Land_c[data$Land==6] <- "HH"
data$Land_c[data$Land==7] <- "HE"
data$Land_c[data$Land==8] <- "MV"
data$Land_c[data$Land==9] <- "NI"
data$Land_c[data$Land==10] <- "NW"
data$Land_c[data$Land==11] <- "RP"
data$Land_c[data$Land==12] <- "SL"
data$Land_c[data$Land==13] <- "SN"
data$Land_c[data$Land==14] <- "ST"
data$Land_c[data$Land==15] <- "SH"
data$Land_c[data$Land==16] <- "TH"
data$East <- 0
data$East[data$Land_c=="MV"] <- 1
data$East[data$Land_c=="ST"] <- 1
data$East[data$Land_c=="TH"] <- 1
data$East[data$Land_c=="SN"] <- 1
data$East[data$District=="Treptow-Köpenick"] <- 1
data$East[data$District=="Pankow"] <- 1
data$East[data$District=="Marzahn-Hellersdorf"] <- 1
data$East[data$District=="Lichtenberg"] <- 1

#rename # of parties
data$N <- data$`(raw).no..of.parties`

#Figure 1
ggplot(data, aes(x=reorder(Land_c, LSI, FUN=median), y= LSI))+ylab("LSI")+geom_boxplot()+xlab("Länder")+theme(axis.title.x=element_text(face="italic"))

#Figure 2
ggplot(data, aes(x=reorder(Land_c, N, FUN=median), y= N))+ylab("N")+geom_boxplot()+xlab("Länder")+theme(axis.title.x=element_text(face="italic"))
ggplot(data, aes(x=reorder(Land_c, ENP, FUN=median), y= ENP))+ylab("ENP")+geom_boxplot()+xlab("Länder")+theme(axis.title.x=element_text(face="italic"))
ggplot(data, aes(x=reorder(Land_c, LRI, FUN=median), y= LRI))+ylab("LRI")+geom_boxplot()+xlab("Länder")+theme(axis.title.x=element_text(face="italic"))

#Table 3
summary(data$LHI)
summary(data$LSI)
summary(data$N)
summary(data$ENP)
summary(data$LRI)

#Table 4
data$threshold.dummy <- 0
data$threshold.dummy[data$legal.threshold!=0] <-1
data_sub <- subset(data, data$Land_c!="HH")
LSI_m0 <- lmer(LSI~1+(1|Land), data=data_sub)
LSI_m1 <- lmer(LSI~as.factor(Electoral.system)+(1|Land), data=data_sub)
LSI_m2 <- lmer(LSI~threshold.dummy+Assembly.size+(1|Land), data=data_sub)
LSI_m3 <- lmer(LSI~Eff_Threshold+(1|Land), data=data_sub)
LSI_m4 <- lmer(LSI~as.factor(Electoral.system)+threshold.dummy+Assembly.size+(1|Land), data=data_sub)
LSI_m5 <- lmer(LSI~as.factor(Electoral.system)+Eff_Threshold+(1|Land), data=data_sub)

tab_model(LSI_m0, LSI_m1,LSI_m2,LSI_m3,LSI_m4,LSI_m5, show.ci = F, collapse.se = T, p.style = "stars")

#Table 5 (N) 
N_m0 <- lmer(N~1+(1|Land), data=data_sub)
N_m1 <- lmer(N~as.factor(Electoral.system)+(1|Land), data=data_sub)
N_m2 <- lmer(N~threshold.dummy+Assembly.size+(1|Land), data=data_sub)
N_m3 <- lmer(N~Eff_Threshold+(1|Land), data=data_sub)
N_m4 <- lmer(N~as.factor(Electoral.system)+threshold.dummy+Assembly.size+(1|Land), data=data_sub)
N_m5 <- lmer(N~as.factor(Electoral.system)+Eff_Threshold+(1|Land), data=data_sub)

tab_model(N_m0, N_m1,N_m2,N_m3,N_m4,N_m5, show.ci = F, collapse.se = T, p.style = "stars")

#Table 6 (ENP) 
ENP_m0 <- lmer(ENP~1+(1|Land), data=data_sub)
ENP_m1 <- lmer(ENP~as.factor(Electoral.system)+(1|Land), data=data_sub)
ENP_m2 <- lmer(ENP~threshold.dummy+Assembly.size+(1|Land), data=data_sub)
ENP_m3 <- lmer(ENP~Eff_Threshold+(1|Land), data=data_sub)
ENP_m4 <- lmer(ENP~as.factor(Electoral.system)+threshold.dummy+Assembly.size+(1|Land), data=data_sub)
ENP_m5 <- lmer(ENP~as.factor(Electoral.system)+Eff_Threshold+(1|Land), data=data_sub)

tab_model(ENP_m0, ENP_m1,ENP_m2,ENP_m3,ENP_m4,ENP_m5, show.ci = F, collapse.se = T, p.style = "stars")

#Table 7 (LRI) 
LRI_m0 <- lmer(LRI~1+(1|Land), data=data_sub)
LRI_m1 <- lmer(LRI~as.factor(Electoral.system)+(1|Land), data=data_sub)
LRI_m2 <- lmer(LRI~threshold.dummy+Assembly.size+(1|Land), data=data_sub)
LRI_m3 <- lmer(LRI~Eff_Threshold+(1|Land), data=data_sub)
LRI_m4 <- lmer(LRI~as.factor(Electoral.system)+threshold.dummy+Assembly.size+(1|Land), data=data_sub)
LRI_m5 <- lmer(LRI~as.factor(Electoral.system)+Eff_Threshold+(1|Land), data=data_sub)

tab_model(LRI_m0, LRI_m1,LRI_m2,LRI_m3,LRI_m4,LRI_m5, show.ci = F, collapse.se = T, p.style = "stars")


#Appendix 2
tab_corr(as.data.frame(cbind(data$LHI, data$LSI, data$N, data$ENP,data$LRI)),
         var.labels = c("LHI", "LSI", "N","ENP", "LRI"),
         show.p = F, 
         p.numeric = F, 
         fade.ns = F, 
         digits = 3,
         triangle = "upper")

#Appendix 3
LHI_m0 <- lmer(LHI~1+(1|Land), data=data_sub)
LHI_m1 <- lmer(LHI~as.factor(Electoral.system)+(1|Land), data=data_sub)
LHI_m2 <- lmer(LHI~threshold.dummy+Assembly.size+(1|Land), data=data_sub)
LHI_m3 <- lmer(LHI~Eff_Threshold+(1|Land), data=data_sub)
LHI_m4 <- lmer(LHI~as.factor(Electoral.system)+threshold.dummy+Assembly.size+(1|Land), data=data_sub)
LHI_m5 <- lmer(LHI~as.factor(Electoral.system)+Eff_Threshold+(1|Land), data=data_sub)

tab_model(LHI_m0, LHI_m1,LHI_m2,LHI_m3,LHI_m4,LHI_m5, show.ci = F, collapse.se = T, p.style = "stars")

#Cooks Distance
library(HLMdiag)
CD <- as.data.frame(cbind(cooks.distance(LSI_m0, level="Land"),
                          cooks.distance(LSI_m1, level="Land"),
                          cooks.distance(LSI_m2, level="Land"),
                          cooks.distance(LSI_m3, level="Land"),
                          cooks.distance(LSI_m4, level="Land"),
                          cooks.distance(LSI_m5, level="Land"),
                          cooks.distance(N_m0, level="Land"),
                          cooks.distance(N_m1, level="Land"),
                          cooks.distance(N_m2, level="Land"),
                          cooks.distance(N_m3, level="Land"),
                          cooks.distance(N_m4, level="Land"),
                          cooks.distance(N_m5, level="Land"),
                          cooks.distance(ENP_m0, level="Land"),
                          cooks.distance(ENP_m1, level="Land"),
                          cooks.distance(ENP_m2, level="Land"),
                          cooks.distance(ENP_m3, level="Land"),
                          cooks.distance(ENP_m4, level="Land"),
                          cooks.distance(ENP_m5, level="Land"),
                          cooks.distance(LRI_m0, level="Land"),
                          cooks.distance(LRI_m1, level="Land"),
                          cooks.distance(LRI_m2, level="Land"),
                          cooks.distance(LRI_m3, level="Land"),
                          cooks.distance(LRI_m4, level="Land"),
                          cooks.distance(LRI_m5, level="Land")))

CD <- t(CD)
CD <- as.data.frame(CD)
tab_df(CD, )

#
t.test(data$ENP[data$Land_c=="BW"|data$Land_c=="SN"|data$Land_c=="ST"], data$ENP[data$Land_c=="HE"|data$Land_c=="NW"])
t.test(data$LSI[data$Land_c=="BW"|data$Land_c=="SN"|data$Land_c=="ST"], data$LSI[data$Land_c=="HE"|data$Land_c=="NW"])
t.test(data$N[data$Land_c=="BW"|data$Land_c=="SN"|data$Land_c=="ST"], data$N[data$Land_c=="HE"|data$Land_c=="NW"])
t.test(data$LRI[data$Land_c=="BW"|data$Land_c=="SN"|data$Land_c=="ST"], data$LRI[data$Land_c=="HE"|data$Land_c=="NW"])