####################################################################
## author:    Michael L. Wicki
## contact:   wimi@ethz.ch, ETH Zurich
## file name: mun_ref_analysis.R
## Context:   ISTP mobility lab - The Decline of Hometown Glory
## started:   2018-01-26
## Summary:   Structural equation model & linear fixed-effects model
####################################################################

#Copyright (C) 2020 Michael Wicki

#This R-Script is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#any later version.

#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.

#You should have received a copy of the GNU General Public License
#along with this program.  If not, see <http://www.gnu.org/licenses/>.

rm(list = ls())
getwd()
setwd("Z:\\People\\Wicki Michael\\Mobility Lab\\Analyse")
set.seed(42)

################################ load packages ################################ 

library(foreign)
library(Hmisc)
library(psych)
library(ggplot2)
library(plyr)
library(lme4)
library(taRifx)
library(epiDisplay)
library(corrplot)
require(stats)
library(robustlmm)
library(stargazer)
library(readxl)
library(scales)
library(reshape2)
library(MuMIn)
library(cowplot)
library(car)
library(xtable)
require(nnet)
library(reshape2)
library(lmtest)
library(multiwayvcov)


load("replicationData.Rda")

################################ Graphs Paper ################################ 

######### Violin plot
datasettotal$all <- 1
datm <- melt(datasettotal,id.vars='all', measure.vars=c('nattMM1','nattMM2','nattMM4'))
datm$variable <- factor(datm$variable, levels=c('nattMM1','nattMM2','nattMM4'), labels=c("amalgamation", "inter-municipal cooperation", "cantonal delegation"))

dodge <- position_dodge(width = 1)
vpl2 <- ggplot(datm, aes(x=variable, y=value)) +
  geom_violin(position=dodge, fill = "grey") +
  geom_boxplot(width=.1, position = dodge) +
  ggtitle("") +
  ylab("reform support") +
  theme_bw() +
  scale_y_continuous(breaks=c(1:11), labels=c(1:11),limits=c(1,11)) +
#  scale_x_continuous() +
  theme(plot.title = element_text(size=16),
        axis.text.x=element_text(size=16, angle = 35, hjust = 1),
        axis.text.y=element_text(size=16),
        axis.title.y=element_text(size=16),
        axis.title.x = element_blank(),
        panel.grid.minor = element_blank())
vpl2

ggsave(filename = "./figs/figure2.png", vpl2)


############# Social network geography and support plots with boxplots
datasettotal$decileelAkm2 <- with(datasettotal, cut(elAkm2_w_original, 
                                breaks=quantile(elAkm2_w_original, probs=seq(0,1, by=0.1), na.rm=TRUE), 
                                include.lowest=TRUE),
                                labels=c("1","2","3","4","5","6","7","8","9","10"))

datasettotal$decilefSpaOverb <- with(datasettotal, cut(fSpaOverb, 
                                                    breaks=quantile(fSpaOverb, probs=seq(0,1, by=0.1), na.rm=TRUE), 
                                                    include.lowest=TRUE),
                                     labels=c("1","2","3","4","5","6","7","8","9","10"))


pl22 <- ggplot(data=subset(datasettotal, !is.na(decileelAkm2)), aes(x = decileelAkm2, y = nattMM1)) +
  geom_boxplot() +
  ggtitle("") +
  xlab("") +
  ylab("amalgamation") +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text=element_text(size=12),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
pl22

pl23 <- ggplot(data=subset(datasettotal, !is.na(decileelAkm2)), aes(x = decileelAkm2, y = nattMM2)) + 
  geom_boxplot() +
  ggtitle("") +
  xlab("") +
  ylab("inter-municipal cooperation") +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text=element_text(size=12),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
pl23

pl24 <- ggplot(data=subset(datasettotal, !is.na(decileelAkm2)), aes(x = decileelAkm2, y = nattMM4)) + 
  geom_boxplot() +
  scale_x_discrete(labels=c("1","2","3","4","5","6","7","8","9","10")) +
  ggtitle("") +
  xlab("Individual social network geography [deciles]") +
  ylab("cantonal delegation") +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text=element_text(size=12),
        panel.grid.minor = element_blank())
pl24

pl25 <- ggplot(data=subset(datasettotal, !is.na(decilefSpaOverb)), aes(x = decilefSpaOverb, y = nattMM1)) + 
  geom_boxplot() +
  ggtitle("") +
  xlab("") +
  ylab("") +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text=element_text(size=12),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())
pl25

pl26 <- ggplot(data=subset(datasettotal, !is.na(decilefSpaOverb)), aes(x = decilefSpaOverb, y = nattMM2)) + 
  geom_boxplot() +
  ggtitle("") +
  xlab("") +
  ylab("") +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text=element_text(size=12),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())
pl26

pl27 <- ggplot(data=subset(datasettotal, !is.na(decilefSpaOverb)), aes(x = decilefSpaOverb, y = nattMM4)) + 
  geom_boxplot() +
  ggtitle("") +
  scale_x_discrete(labels=c("1","2","3","4","5","6","7","8","9","10")) +
  xlab("Under-bounded administrative space [deciles]") +
  ylab("") +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text=element_text(size=12),
        panel.grid.minor = element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())
pl27

p_completebp <- ggdraw() +
  draw_plot(pl22, x = 0, y = .67, width = .5, height = .33) +
  draw_plot(pl25, x = .5, y = .67, width = .5, height = .33) +
  draw_plot(pl23, x = 0, y = 0.34, width = .5, height = 0.33) +
  draw_plot(pl26, x = .5, y = 0.34, width = .5, height = 0.33) +
  draw_plot(pl24, x = 0, y = 0, width = .5, height = .33) +
  draw_plot(pl27, x = .5, y = 0, width = .5, height = .33) +
  draw_plot_label(label = c("A", "B", "C", "D", "E", "F"), size = 15,
                  x = c(0, 0.5, 0, 0.5, 0, 0.5), y = c(1, 1, 0.67, 0.67, 0.34, 0.34))
p_completebp

ggsave(filename = "./figs/figure4.png", p_completebp,  width = 8, height = 10)


###functional vs. administrative space
pl6 <- ggplot(datasettotal, aes(x = fSpaOverb, y = munkm2)) + geom_point() +
  ggtitle("") +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("Under-bounded administrative space [log]") +
  ylab("municipality size [km2]") +
  theme_bw() +
  theme(text = element_text(size=17),
        axis.text=element_text(size=14),
        panel.grid.minor = element_blank()) +
  annotation_logticks(sides = "b")
pl6

pl7 <- ggplot(datasettotal, aes(x = munkm2, y = elAkm2_w_original)) + geom_point() +
  ggtitle("") +
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  geom_smooth(method='lm') +
  xlab("municipality size [km2]") +
  ylab("Individual social network geography [log(km2)]") +
  theme_bw() +
  theme(text = element_text(size=17),
        axis.text=element_text(size=14),
        panel.grid.minor = element_blank()) +
  annotation_logticks(sides = "lb")
pl7

summary(lm(log(datasettotal$elAkm2_w_original + 1) ~ log(datasettotal$munkm2 + 1)))




############################################################# Percent outside municipality ###########################################################
lprcnt1 <- lmer(nattMM1 ~ prcntMun * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lprcnt2 <- lmer(nattMM2 ~ prcntMun * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lprcnt3 <- lmer(nattMM4 ~ prcntMun * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)

r.squaredGLMM(lprcnt1)
r.squaredGLMM(lprcnt2)
r.squaredGLMM(lprcnt3)

tablelcontacts <- stargazer(lprcnt1, lprcnt2, lprcnt3,
                          type="html", 
                          single.row = TRUE,
                          title="Table A 4: Predictors of municipal reform support (Mixed-Effects Linear Regression)", 
                          omit=c("educ","inc"), 
                          align=TRUE,
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          digits = 3,
                          dep.var.labels = c("Amalgamation", "Inter-municipal cooperation", "Cantonal delegation"),
                          covariate.labels = c("Mun. contacts [share]", "Rel. municipality attachment", "Commuting [log]", "Service satisfaction", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity", "L-R self placement", "Rel. governm. trust", "Tax rate", "Population [log]", "Mun. contacts*rel. municipality attachment", "Commuting*rel. municipality attachment", "Intercept"),
                          object.names = FALSE,
                          notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Conditional R2", "16.6", "10.2", "11.2"),
                                         c("Municipalities", "164", "165", "164")),
                          notes = "Table entries are maximum-likelihood estimates with estimated standard errors in parantheses. Education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                          out='tablelcontacts.html')


# nattMM5 
##SNG
lfemSNG5 <- lmer(nattMM5 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemOver5 <- lmer(nattMM5 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)


r.squaredGLMM(lfemSNG5)
r.squaredGLMM(lfemOver5)


tablelargSNG <- stargazer(lfemSNG5, lfemOver5,
                           type="html", 
                           single.row = TRUE,
                           title="Predictors of municipal reform policy support (Mixed-Effects Linear Regression)", 
                           omit=c("educ","inc"), 
                           align=TRUE,
                           star.char = c("+", "*", "**", "***"),
                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                           digits = 3,
                           dep.var.labels = c("Identity loss"),
                           covariate.labels = c("SNG [log]", "SNG over-bound [log]", "Rel. municipality attachment", "Commuting [log]", "Service satisfaction", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity", "L-R self placement", "Rel. governm. trust", "Political interest (mun.)", "Tax rate", "Population [log]","SNG*rel. municipality attachment", "SNG over-bound*rel. municipality attachment", "Commuting*rel. municipality attachment", "Intercept"),
                           object.names = FALSE,
                           notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Conditional R2", "9.8", "11.4"),
                                         c("Municipalities", "164", "164")),
                           notes = "Table entries are maximum-likelihood estimates with estimated standard errors in parentheses. Education and income were included in the regression, but were not significant and are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                           out='tablelidentityloss.html')




################################# Main Models Papers ######################
########## Amalgamation ############
lfemMrOver2 <- lmer(nattMM1 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log | gemeinde), data=datasettotal, REML = FALSE)
lfemMrSng2 <- lmer(nattMM1 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
rlfemMrOver2 <- rlmer(nattMM1 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 | gemeinde), data=datasettotal, REML = FALSE)
rlfemMrSng2 <- rlmer(nattMM1 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 | gemeinde), data=datasettotal, REML = FALSE)

summary(lfemMrOver2)
summary(lfemMrSng2)

lfemMrOver2vif <- vif(lfemMrOver2)
lfemMrSng2vif <- vif(lfemMrSng2)
stargazer(lfemMrOver2vif, out='lfemMrOver2vif.html')
stargazer(lfemMrSng2vif, out='lfemMrSng2vif.html')

########## inter-municipal cooperation #############
lfemMmOver2 <- lmer(nattMM2 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemMmSng2 <- lmer(nattMM2 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
rlfemMmOver2 <- rlmer(nattMM2 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 | gemeinde), data=datasettotal, REML = FALSE)
rlfemMmSng2 <- rlmer(nattMM2 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 | gemeinde), data=datasettotal, REML = FALSE)

summary(rlfemMmOver2)
summary(rlfemMmSng2)

lfemMmOver2vif <- vif(lfemMmOver2)
lfemMmSng2vif <- vif(lfemMmSng2)
stargazer(lfemMmOver2vif, out='lfemMmOver2vif.html')
stargazer(lfemMmSng2vif, out='lfemMmSng2vif.html')


########## cantonal delegation #############
lfemCdOver2 <- lmer(nattMM4 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log | gemeinde), data=datasettotal, REML = FALSE)
lfemCdSng2 <- lmer(nattMM4 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log | gemeinde), data=datasettotal, REML = FALSE)
rlfemCdOver2 <- rlmer(nattMM4 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 | gemeinde), data=datasettotal, REML = FALSE)
rlfemCdSng2 <- rlmer(nattMM4 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 | gemeinde), data=datasettotal, REML = FALSE)

summary(rlfemCdOver2)
summary(rlfemCdSng2)

lfemCdOver2vif <- vif(lfemCdOver2)
lfemCdSng2vif <- vif(lfemCdSng2)
stargazer(lfemCdOver2vif, out='lfemCdOver2vif.html')
stargazer(lfemCdSng2vif, out='lfemCdSng2vif.html')


#Model w/out robust estimators
r.squaredGLMM(lfemMrOver2)
r.squaredGLMM(lfemMrSng2)
r.squaredGLMM(lfemMmOver2)
r.squaredGLMM(lfemMmSng2)
r.squaredGLMM(lfemCdOver2)
r.squaredGLMM(lfemCdSng2)

tablelattach <- stargazer(lfemMrOver2, lfemMrSng2, lfemMmOver2, lfemMmSng2, lfemCdOver2, lfemCdSng2,
                           type="html", 
                           single.row = TRUE,
                           title="Table 2: Predictors of municipal reform support (mixed-effects linear regression)", 
                           omit=c("educ","inc"), 
                           align=TRUE,
                           star.char = c("+", "*", "**", "***"),
                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                           digits = 3,
                           dep.var.labels = c("Amalgamation", "Inter-municipal cooperation", "Cantonal delegation"),
                           covariate.labels = c("SNG over-bound [log]", "SNG [log]", "Rel. municipality attachment [RMA]", "Commuting [log]", "Service satisfaction", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity", "L-R self placement", "Rel. governm. trust", "Political interest (mun.)", "Tax rate", "Population [log]", "SNG over-bound*RMA", "SNG*RMA", "Commuting*RMA", "Intercept"),
                           object.names = FALSE,
                           notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Conditional R2", "16.6", "17.4", "10.6", "10.6", "11.5", "11.2"),
                                         c("Municipalities", "164", "164", "165", "165", "164", "164")),
                           notes = "Table entries are maximum-likelihood estimates with estimated standard errors in parentheses. Education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                           out='tablelmain.html')

##Model w/ robust estimators
coeffs <- coef(summary(rlfemMrOver2)) 
p <- pnorm(abs(coeffs[, "t value"]), lower.tail = FALSE) * 2 

coeffsp1 <- cbind(coeffs, "p value" = round(p,3)) 
write.csv(coeffsp1, "rlfemMrOver2.csv")

coeffs <- coef(summary(rlfemMrSng2)) 
p <- pnorm(abs(coeffs[, "t value"]), lower.tail = FALSE) * 2 

coeffsp2 <- cbind(coeffs, "p value" = round(p,3)) 
write.csv(coeffsp2, "rlfemMrSng2.csv")


coeffs <- coef(summary(rlfemMmOver2)) 
p <- pnorm(abs(coeffs[, "t value"]), lower.tail = FALSE) * 2 

coeffsp3 <- cbind(coeffs, "p value" = round(p,3)) 
write.csv(coeffsp3, "rlfemMmOver2.csv")


coeffs <- coef(summary(rlfemMmSng2)) 
p <- pnorm(abs(coeffs[, "t value"]), lower.tail = FALSE) * 2 

coeffsp4 <- cbind(coeffs, "p value" = round(p,3)) 
write.csv(coeffsp4, "rlfemMmSng2.csv")


coeffs <- coef(summary(rlfemCdOver2)) 
p <- pnorm(abs(coeffs[, "t value"]), lower.tail = FALSE) * 2 

coeffsp5 <- cbind(coeffs, "p value" = round(p,3)) 
write.csv(coeffsp5, "rlfemCdOver2.csv")


coeffs <- coef(summary(rlfemCdSng2)) 
p <- pnorm(abs(coeffs[, "t value"]), lower.tail = FALSE) * 2 

coeffsp6 <- cbind(coeffs, "p value" = round(p,3)) 
write.csv(coeffsp6, "rlfemCdSng2.csv")

write.csv(c(coeffsp1, coeffsp2, coeffsp3, coeffsp4, coeffsp5, coeffsp6), "rlfemall.csv")


#check for normality of normally distributed residuals
dev.copy(png,'residualplot1.png')
plot(lfemMrSng2, main="(1)")
dev.off()
dev.copy(png,'residualplot2.png')
plot(lfemMrOver2, main="(2)")
dev.off()
dev.copy(png,'residualplot3.png')
plot(lfemMmSng2, main="(3)")
dev.off()
dev.copy(png,'residualplot4.png')
plot(lfemMmOver2, main="(4)")
dev.off()
dev.copy(png,'residualplot5.png')
plot(lfemCdSng2, main="(1)")
dev.off()
dev.copy(png,'residualplot6.png')
plot(lfemCdOver2, main="(2)")
dev.off()

dev.copy(png,'qqplot.png')
par(mfrow=c(3,2))
qqnorm(resid(lfemMrSng2), main="Normal Q-Q Plot (1)") 
lfemMrSng2qqnorm <- qqline(resid(lfemMrSng2), col = "red")
qqnorm(resid(lfemMrOver2), main="Normal Q-Q Plot (2)")
lfemMrOver2qqnorm <- qqline(resid(lfemMrOver2), col = "red")
qqnorm(resid(lfemMmSng2), main="Normal Q-Q Plot (3)") 
lfemMmSng2qqnorm <- qqline(resid(lfemMmSng2), col = "red")
qqnorm(resid(lfemMmOver2), main="Normal Q-Q Plot (4)")
lfemMmOver2qqnorm <- qqline(resid(lfemMmOver2), col = "red")
qqnorm(resid(lfemCdSng2), main="Normal Q-Q Plot (3)") 
lfemMmSng2qqnorm <- qqline(resid(lfemMmSng2), col = "red")
qqnorm(resid(lfemCdOver2), main="Normal Q-Q Plot (4)")
lfemMmOver2qqnorm <- qqline(resid(lfemMmOver2), col = "red")
dev.off()



qqnorm(resid(lfemMmOver2))
qqnorm(resid(lfemMmSng2))




#check for sensitivity of data
ggplot(data.frame(lev=hatvalues(lfemMrSng2),pearson=residuals(lfemMrSng2,type="pearson")),
       aes(x=lev,y=pearson)) +
  geom_point() +
  theme_bw()

ggplot(data.frame(lev=hatvalues(lfemMrOver2),pearson=residuals(lfemMrOver2,type="pearson")),
       aes(x=lev,y=pearson)) +
  geom_point() +
  theme_bw()


########## relative attachement #############
lfemattSng <- lmer(cmAttach ~ (fSpaOverblog + comm_dist_log) * (satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust + tax_rate + pop_log) + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemattOver <- lmer(cmAttach ~ (elAkm2_w_log + comm_dist_log) * (satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust + tax_rate + pop_log) + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemsatSng <- lmer(satSvc ~ (fSpaOverblog + comm_dist_log) * (cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust + tax_rate + pop_log) + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemsatOver <- lmer(satSvc ~ (elAkm2_w_log + comm_dist_log) * (cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust + tax_rate + pop_log) + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemtrstSng <- lmer(rel_trust ~ (fSpaOverblog + comm_dist_log) * (yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + tax_rate + pop_log) + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemtrstOver <- lmer(rel_trust ~ (elAkm2_w_log + comm_dist_log) * (yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + tax_rate + pop_log) + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)

summary(lfemattSng1)

r.squaredGLMM(lfemattSng)
r.squaredGLMM(lfemattOver)
r.squaredGLMM(lfemsatSng)
r.squaredGLMM(lfemsatOver)
lfemattOvervif <- vif(lfemattOver)
lfemattSngvif <- vif(lfemattSng)
stargazer(lfemattOvervif, out='lfemattOvervif.html')
stargazer(lfemattSngvif, out='lfemattSngvif.html')
lfemsatSngvif <- vif(lfemsatSng)
lfemsatOvervif <- vif(lfemsatOver)
stargazer(lfemsatOvervif, out='lfemsatOvervif.html')
stargazer(lfemsatSngvif, out='lfemsatSngvif.html')


#Model w/out robust estimators
tablelattach <- stargazer(lfemattSng, lfemattOver, lfemsatSng, lfemsatOver,
                          type="html", 
                          single.row = TRUE,
                          title="Table 1: Predictors of attachment and service satisfaction (mixed-effects linear regression)", 
                          omit=c("educ","inc",
                                 "fSpaOverblog:satSvc" , "fSpaOverblog:cmAttach" , "fSpaOverblog:yearsresidence" , "fSpaOverblog:homeowner" ,
                                 "fSpaOverblog:relocD" , "fSpaOverblog:age" , "fSpaOverblog:sexFSmpl" , "fSpaOverblog:swiss" ,
                                 "fSpaOverblog:locciveng" , "fSpaOverblog:relD" , "fSpaOverblog:tax_rate" , "fSpaOverblog:pop_log" ,
                                 "elAkm2_w_log:satSvc" , "elAkm2_w_log:cmAttach" , "elAkm2_w_log:yearsresidence" , "elAkm2_w_log:homeowner" , "elAkm2_w_log:relocD" , "elAkm2_w_log:age" ,
                                 "elAkm2_w_log:sexFSmpl" , "elAkm2_w_log:swiss" , "elAkm2_w_log:locciveng" , "elAkm2_w_log:relD" , "elAkm2_w_log:tax_rate" , "elAkm2_w_log:pop_log" , "comm_dist_log:satSvc" , "comm_dist_log:cmAttach" , "comm_dist_log:yearsresidence" , "comm_dist_log:homeowner" ,
                                 "comm_dist_log:relocD" , "comm_dist_log:age" , "comm_dist_log:sexFSmpl" , "comm_dist_log:swiss" ,
                                 "comm_dist_log:locciveng" , "comm_dist_log:relD" , "comm_dist_log:tax_rate" , "comm_dist_log:pop_log",
                                 "comm_dist_log:rel_trust", "elAkm2_w_log:rel_trust", "fSpaOverblog:rel_trust"), 
                          align=TRUE,
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          digits = 3,
                          dep.var.labels = c("rel. municipal attachment", "Service satisfaction"),
                          covariate.labels = c("SNG over-bound [log]", "SNG [log]", "Commuting [log]", "Service satisfaction", "rel. municipal attachment", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity",  "Trust in local Government", "Tax rate", "Population [log]", "Intercept"),
                          object.names = FALSE,
                          notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Interactions included?", "Yes", "Yes", "Yes", "Yes"), 
                                         c("Conditional R2", "14.8", "14.6", "14.9", "14.8"),
                                         c("Municipalities", "167", "167", "167", "167")),
                          notes = "Table entries are maximum-likelihood estimates with estimated standard errors in parantheses. Interaction terms of the mobility measures with all other variables, education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                          out='tablelattach2.html')

dev.copy(png,'residualplotlfemattSng.png')
plot(lfemattSng, main="(2)")
dev.off()
dev.copy(png,'residualplotlfemattOver.png')
plot(lfemattOver, main="(1)")
dev.off()
dev.copy(png,'residualplotlfemsatSng.png')
plot(lfemsatSng, main="(4)")
dev.off()
dev.copy(png,'residualplotlfemsatOver.png')
plot(lfemsatOver, main="(3)")
dev.off()

################################# Service Satisfaction Models Papers ######################
########## Amalgamation ############
lfemMrOver2 <- lmer(nattMM1 ~ fSpaOverblog * satSvc + comm_dist_log * satSvc + cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemMrSng2 <- lmer(nattMM1 ~ elAkm2_w_log * satSvc + comm_dist_log * satSvc + cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)

summary(lfemMrOver2)
summary(lfemMrSng2)
r.squaredGLMM(lfemMrOver2)
r.squaredGLMM(lfemMrSng2)
lfemMrOver2vif <- vif(lfemMrOver2)
lfemMrSng2vif <- vif(lfemMrSng2)
stargazer(lfemMrOver2vif, out='lfemMmOver2vif.html')
stargazer(lfemMrSng2vif, out='lfemMmSng2vif.html')

########## inter-municipal cooperation #############
lfemMmOver2 <- lmer(nattMM2 ~ fSpaOverblog * satSvc + comm_dist_log * satSvc + cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemMmSng2 <- lmer(nattMM2 ~ elAkm2_w_log * satSvc + comm_dist_log * satSvc + cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)

summary(rlfemMmOver2)
summary(rlfemMmSng2)
r.squaredGLMM(lfemMmOver2)
r.squaredGLMM(lfemMmSng2)
lfemMmOver2vif <- vif(lfemMmOver2)
lfemMmSng2vif <- vif(lfemMmSng2)
stargazer(lfemMmOver2vif, out='lfemMmOver2vif.html')
stargazer(lfemMmSng2vif, out='lfemMmSng2vif.html')


########## cantonal delegation #############
lfemCdOver2 <- lmer(nattMM4 ~ fSpaOverblog * satSvc + comm_dist_log * satSvc + cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)
lfemCdSng2 <- lmer(nattMM4 ~ elAkm2_w_log * satSvc + comm_dist_log * satSvc + cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + tax_rate + pop_log + (1 + tax_rate + pop_log  | gemeinde), data=datasettotal, REML = FALSE)

summary(rlfemCdOver2)
summary(rlfemCdSng2)
r.squaredGLMM(lfemCdOver2)
r.squaredGLMM(lfemCdSng2)
lfemMmOver2vif <- vif(lfemCdOver2)
lfemMmSng2vif <- vif(lfemCdSng2)
stargazer(lfemCdOver2vif, out='lfemCdOver2vif.html')
stargazer(lfemCdSng2vif, out='lfemCdSng2vif.html')


#Model w/out robust estimators
tablelattach <- stargazer(lfemMrOver2, lfemMrSng2, lfemMmOver2, lfemMmSng2, lfemCdOver2, lfemCdSng2,
                          type="html", 
                          single.row = TRUE,
                          title="Table A: Predictors of municipal reform support (mixed-effects linear regression)", 
                          omit=c("educ","inc"), 
                          align=TRUE,
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          digits = 3,
                          dep.var.labels = c("Amalgamation", "Inter-municipal cooperation", "Cantonal delegation"),
                          covariate.labels = c("SNG over-bound [log]", "SNG [log]", "Service satisfaction", "Commuting [log]", "Rel. municipality attachment", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity", "L-R self placement", "Rel. governm. trust", "Political interest (mun.)", "Tax rate", "Population [log]", "SNG over-bound*Service satisfaction", "SNG*Service satisfaction", "Commuting*Service satisfaction", "Intercept"),
                          object.names = FALSE,
                          notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Conditional R2", "17.7", "19.3", "9.4", "10.1", "11.6", "10.1"),
                                         c("Municipalities", "164", "164", "165", "165", "164", "164")),
                          notes = "Table entries are maximum-likelihood estimates with estimated standard errors in parentheses. Education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                          out='tablelservice.html')

############################ Fixed effects Table 1 ##################
fe1attSng <- lm(cmAttach ~ (fSpaOverblog + comm_dist_log) * (satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)
fe1attOver <- lm(cmAttach ~ (elAkm2_w_log + comm_dist_log) * (satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)
fe1satSng <- lm(satSvc ~ (fSpaOverblog + comm_dist_log) * (cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)
fe1satOver <- lm(satSvc ~ (elAkm2_w_log + comm_dist_log) * (cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)

fe1table1 <- stargazer(fe1attSng, fe1attOver, fe1satSng, fe1satOver,
                          type="html", 
                          single.row = TRUE,
                          title="Table A7: Predictors of attachment and service satisfaction (fixed-effects linear regression)", 
                          omit=c("educ","inc", "gemeinde",
                                 "fSpaOverblog:satSvc" , "fSpaOverblog:cmAttach" , "fSpaOverblog:yearsresidence" , "fSpaOverblog:homeowner" ,
                                 "fSpaOverblog:relocD" , "fSpaOverblog:age" , "fSpaOverblog:sexFSmpl" , "fSpaOverblog:swiss" ,
                                 "fSpaOverblog:locciveng" , "fSpaOverblog:relD" , "fSpaOverblog:tax_rate" , "fSpaOverblog:pop_log" ,
                                 "elAkm2_w_log:satSvc" , "elAkm2_w_log:cmAttach" , "elAkm2_w_log:yearsresidence" , "elAkm2_w_log:homeowner" , "elAkm2_w_log:relocD" , "elAkm2_w_log:age" ,
                                 "elAkm2_w_log:sexFSmpl" , "elAkm2_w_log:swiss" , "elAkm2_w_log:locciveng" , "elAkm2_w_log:relD", "comm_dist_log:satSvc" , "comm_dist_log:cmAttach" , "comm_dist_log:yearsresidence" , "comm_dist_log:homeowner" ,
                                 "comm_dist_log:relocD" , "comm_dist_log:age" , "comm_dist_log:sexFSmpl" , "comm_dist_log:swiss" ,
                                 "comm_dist_log:locciveng" , "comm_dist_log:relD" , 
                                 "comm_dist_log:rel_trust", "elAkm2_w_log:rel_trust", "fSpaOverblog:rel_trust"), 
                          align=TRUE,
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          digits = 3,
                          dep.var.labels = c("rel. municipal attachment", "Service satisfaction"),
                          covariate.labels = c("SNG over-bound [log]", "SNG [log]", "Commuting [log]", "Service satisfaction", "rel. municipal attachment", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity",  "Trust in local Government", "Intercept"),
                          object.names = FALSE,
                          notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Interactions included?", "Yes", "Yes", "Yes", "Yes"), 
                                         c("Municipality fixed effects", "Yes", "Yes", "Yes", "Yes"),
                                         c("Municipalities", "167", "167", "167", "167")),
                          notes = "Table entries are maximum-likelihood estimates with estimated standard errors in parantheses. Interaction terms of the mobility measures with all other variables, education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                          out='fe1table1.html')

#w/ clustered standard errors
feattSng <- lm(cmAttach ~ (fSpaOverblog + comm_dist_log) * (satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)
feattOver <- lm(cmAttach ~ (elAkm2_w_log + comm_dist_log) * (satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)
fesatSng <- lm(satSvc ~ (fSpaOverblog + comm_dist_log) * (cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)
fesatOver <- lm(satSvc ~ (elAkm2_w_log + comm_dist_log) * (cmAttach + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + rel_trust) + as.factor(gemeinde), data=datasettotal)

# standard errors clustered by group
fe2attSng <- coeftest(feattSng, cluster.vcov(feattSng, datasettotal$gemeinde))
fe2attOver <- coeftest(feattOver, cluster.vcov(feattOver, datasettotal$gemeinde))
fe2satSng <- coeftest(fesatSng, cluster.vcov(fesatSng, datasettotal$gemeinde))
fe2satOver <- coeftest(fesatOver, cluster.vcov(fesatOver, datasettotal$gemeinde))

#table
fe2table1 <- stargazer(fe2attSng, fe2attOver, fe2satSng, fe2satOver,
                       type="html", 
                       single.row = TRUE,
                       title="Table A8: Predictors of attachment and service satisfaction (fixed-effects linear regression with clustered standard errors)", 
                       omit=c("educ","inc", "gemeinde",
                              "fSpaOverblog:satSvc" , "fSpaOverblog:cmAttach" , "fSpaOverblog:yearsresidence" , "fSpaOverblog:homeowner" ,
                              "fSpaOverblog:relocD" , "fSpaOverblog:age" , "fSpaOverblog:sexFSmpl" , "fSpaOverblog:swiss" ,
                              "fSpaOverblog:locciveng" , "fSpaOverblog:relD" , "fSpaOverblog:tax_rate" , "fSpaOverblog:pop_log" ,
                              "elAkm2_w_log:satSvc" , "elAkm2_w_log:cmAttach" , "elAkm2_w_log:yearsresidence" , "elAkm2_w_log:homeowner" , "elAkm2_w_log:relocD" , "elAkm2_w_log:age" ,
                              "elAkm2_w_log:sexFSmpl" , "elAkm2_w_log:swiss" , "elAkm2_w_log:locciveng" , "elAkm2_w_log:relD", "comm_dist_log:satSvc" , "comm_dist_log:cmAttach" , "comm_dist_log:yearsresidence" , "comm_dist_log:homeowner" ,
                              "comm_dist_log:relocD" , "comm_dist_log:age" , "comm_dist_log:sexFSmpl" , "comm_dist_log:swiss" ,
                              "comm_dist_log:locciveng" , "comm_dist_log:relD" , 
                              "comm_dist_log:rel_trust", "elAkm2_w_log:rel_trust", "fSpaOverblog:rel_trust"), 
                       align=TRUE,
                       star.char = c("+", "*", "**", "***"),
                       star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                       digits = 3,
                       dep.var.labels = c("rel. municipal attachment", "Service satisfaction"),
                       covariate.labels = c("SNG over-bound [log]", "SNG [log]", "Commuting [log]", "Service satisfaction", "rel. municipal attachment", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity",  "Trust in local Government", "Intercept"),
                       object.names = FALSE,
                       notes.append = FALSE, notes.align = "l",
                       add.lines=list(c("Interactions included?", "Yes", "Yes", "Yes", "Yes"), 
                                      c("Municipality fixed effects", "Yes", "Yes", "Yes", "Yes"),
                                      c("Clustered standard errors", "Yes", "Yes", "Yes", "Yes"),
                                      c("Municipalities", "167", "167", "167", "167")),
                       notes = "Table entries are maximum-likelihood estimates with estimated standard errors clustered at the municipality level in parantheses. Interaction terms of the mobility measures with all other variables, education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                       out='fe2table1.html')


############################ Fixed effects Table 2 ##################
fe1mm1over <- lm(nattMM1 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
fe1mm1sng <- lm(nattMM1 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
fe1mm2over <- lm(nattMM2 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
fe1mm2sng <- lm(nattMM2 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
fe1mm4over <- lm(nattMM4 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
fe1mm4sng <- lm(nattMM4 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)

#table
tablefe1 <- stargazer(fe1mm1over, fe1mm1sng, fe1mm2over, fe1mm2sng, fe1mm4over, fe1mm4sng,
                          type="html", 
                          single.row = TRUE,
                          title="Table A9: Predictors of municipal reform support (fixed-effects linear regression)", 
                          omit=c("educ","inc", "gemeinde"), 
                          align=TRUE,
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          digits = 3,
                          dep.var.labels = c("Amalgamation", "Inter-municipal cooperation", "Cantonal delegation"),
                          covariate.labels = c("SNG over-bound [log]", "SNG [log]", "Rel. municipality attachment [RMA]", "Commuting [log]", "Service satisfaction", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity", "L-R self placement", "Rel. governm. trust", "Political interest (mun.)", "SNG over-bound*RMA", "SNG*RMA", "Commuting*RMA", "Intercept"),
                          object.names = FALSE,
                          notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Municipality fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                            c("Municipalities", "164", "164", "165", "165", "164", "164")),
                          notes = "Table entries are maximum-likelihood estimates with estimated standard errors in parentheses. Education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                          out='fixedeffects.html')

#w/ clustered standard errors
femm1over <- lm(nattMM1 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
femm1sng <- lm(nattMM1 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
femm2over <- lm(nattMM2 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
femm2sng <- lm(nattMM2 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
femm4over <- lm(nattMM4 ~ fSpaOverblog * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)
femm4sng <- lm(nattMM4 ~ elAkm2_w_log * cmAttach + comm_dist_log * cmAttach + satSvc + yearsresidence + homeowner + relocD + inc + educ + age + sexFSmpl + swiss + locciveng + relD + npolLR + rel_trust + polIntL + as.factor(gemeinde), data=datasettotal)

# standard errors clustered by group
fe2mm1over <- coeftest(femm1over, cluster.vcov(femm1over, datasettotal$gemeinde))
fe2mm1sng <- coeftest(femm1sng, cluster.vcov(femm1sng, datasettotal$gemeinde))
fe2mm2over <- coeftest(femm2over, cluster.vcov(femm2over, datasettotal$gemeinde))
fe2mm2sng <- coeftest(femm2sng, cluster.vcov(femm2sng, datasettotal$gemeinde))
fe2mm4over <- coeftest(femm4over, cluster.vcov(femm4over, datasettotal$gemeinde))
fe2mm4sng <- coeftest(femm4sng, cluster.vcov(femm4sng, datasettotal$gemeinde))

tablefe2 <- stargazer(fe2mm1over, fe2mm1sng, fe2mm2over, fe2mm2sng, fe2mm4over, fe2mm4sng,
                          type="html", 
                          single.row = TRUE,
                          title="Table A10: Predictors of municipal reform support (fixed-effects linear regression with clustered standard errors)", 
                          omit=c("educ","inc", "gemeinde"), 
                          align=TRUE,
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          digits = 3,
                          dep.var.labels = c("Amalgamation", "Inter-municipal cooperation", "Cantonal delegation"),
                          covariate.labels = c("SNG over-bound [log]", "SNG [log]", "Rel. municipality attachment [RMA]", "Commuting [log]", "Service satisfaction", "Tenure", "Homeowner [dummy]", "Relocation [dummy]", "Age", "Female [dummy]", "Swiss [dummy]", "Local civic engagement", "Religiosity", "L-R self placement", "Rel. governm. trust", "Political interest (mun.)", "SNG over-bound*RMA", "SNG*RMA", "Commuting*RMA", "Intercept"),
                          object.names = FALSE,
                          notes.append = FALSE, notes.align = "l",
                          add.lines=list(c("Municipality fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                                         c("Clustered standard errors", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                                         c("Municipalities", "164", "164", "165", "165", "164", "164")),
                          notes = "Table entries are maximum-likelihood estimates with estimated standard errors clustered at the municipality level in parentheses. Education and income were included in the regression, but are not presented in this table. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                          out='fixedeffectsclustered.html')