######################################################################################################################
#                                                                                                                    #
###########################################   Replication script   ###################################################                                                  #
#                                                                                                                    #
######################################################################################################################
########################################## written by Kay Grunenberg #################################################


# set working directory
getwd()
setwd("F:/")

# loading the data into the R environment
data <- read.table("MMP_dataset.csv", header = TRUE, sep = ";", dec = ",")

# subset data  
data1 <- subset(data,national==1)
data2 <- subset(data,post1990==1)
data3 <- subset(data,c(national==1 & post1990==1))

# Testing through t.test for each variable in the full data set
x1 <- t.test(data$lsi[data$bobw==1],data$lsi[data$bobw==0])
x1
x2 <- t.test(data$enps[data$bobw==1],data$enps[data$bobw==0])
x2
x3 <- t.test(data$thresholdPR[data$bobw==1],data$thresholdPR[data$bobw==0])
x3
x4 <- t.test(data$prmedianmagnitude[data$bobw==1],data$prmedianmagnitude[data$bobw==0])
x4
x5 <- t.test(data$smd_share[data$bobw==1],data$smd_share[data$bobw==0])
x5
x6 <- t.test(data$minority_districts[data$bobw==1],data$minority_districts[data$bobw==0])
x6
x7 <- t.test(data$established[data$bobw==1],data$established[data$bobw==0])
x7
x8 <- t.test(data$ethnic[data$bobw==1],data$ethnic[data$bobw==0])
x8

# Testing only national data
x11 <- t.test(data1$lsi[data1$bobw==1],data1$lsi[data1$bobw==0])
x11
x12 <- t.test(data1$enps[data1$bobw==1],data1$enps[data1$bobw==0])
x12
x13 <- t.test(data1$thresholdPR[data1$bobw==1],data1$thresholdPR[data1$bobw==0])
x13
x14 <- t.test(data1$prmedianmagnitude[data1$bobw==1],data1$prmedianmagnitude[data1$bobw==0])
x14
x15 <- t.test(data1$smd_share[data1$bobw==1],data1$smd_share[data1$bobw==0])
x15
x16 <- t.test(data1$minority_districts[data1$bobw==1],data1$minority_districts[data1$bobw==0])
x16
x17 <- t.test(data1$established[data1$bobw==1],data1$established[data1$bobw==0])
x17
x18 <- t.test(data1$ethnic[data1$bobw==1],data1$ethnic[data1$bobw==0])
x18

# Testing only data after 1990
x21 <- t.test(data2$lsi[data2$bobw==1],data2$lsi[data2$bobw==0])
x21
x22 <- t.test(data2$enps[data2$bobw==1],data2$enps[data2$bobw==0])
x22
x23 <- t.test(data2$thresholdPR[data2$bobw==1],data2$thresholdPR[data2$bobw==0])
x23
x24 <- t.test(data2$prmedianmagnitude[data2$bobw==1],data2$prmedianmagnitude[data2$bobw==0])
x24
x25 <- t.test(data2$smd_share[data2$bobw==1],data2$smd_share[data2$bobw==0])
x25
x26 <- t.test(data2$minority_districts[data2$bobw==1],data2$minority_districts[data2$bobw==0])
x26
x27 <- t.test(data2$established[data2$bobw==1],data2$established[data2$bobw==0])
x27
x28 <- t.test(data2$ethnic[data2$bobw==1],data2$ethnic[data2$bobw==0])
x28

# Testing only national data after 1990
x31 <- t.test(data3$lsi[data3$bobw==1],data3$lsi[data3$bobw==0])
x31
x32 <- t.test(data3$enps[data3$bobw==1],data3$enps[data3$bobw==0])
x32
x33 <- t.test(data3$thresholdPR[data3$bobw==1],data3$thresholdPR[data3$bobw==0])
x33
x34 <- t.test(data3$prmedianmagnitude[data3$bobw==1],data3$prmedianmagnitude[data3$bobw==0])
x34
x35 <- t.test(data3$smd_share[data3$bobw==1],data3$smd_share[data3$bobw==0])
x35
x36 <- t.test(data3$minority_districts[data3$bobw==1],data3$minority_districts[data3$bobw==0])
x36
x37 <- t.test(data3$established[data3$bobw==1],data3$established[data3$bobw==0])
x37
x38 <- t.test(data3$ethnic[data3$bobw==1],data3$ethnic[data3$bobw==0])
x38

#### preparing and finishing the presentation of the results

# identifying t.test results for each subset
co1 <- cbind(which(x1$estimate == max(x1$estimate)), which(x2$estimate == max(x2$estimate)),
             which(x3$estimate == max(x3$estimate)), which(x4$estimate == max(x4$estimate)),
             which(x5$estimate == max(x5$estimate)), which(x6$estimate == max(x6$estimate)),
             which(x7$estimate == max(x7$estimate)), which(x8$estimate == max(x8$estimate)))
co2 <- cbind(which(x11$estimate == max(x11$estimate)), which(x12$estimate == max(x12$estimate)),
             which(x13$estimate == max(x13$estimate)), which(x14$estimate == max(x14$estimate)),
             which(x15$estimate == max(x15$estimate)), which(x16$estimate == max(x16$estimate)),
             which(x17$estimate == max(x17$estimate)), which(x18$estimate == max(x18$estimate)))
co3 <- cbind(which(x21$estimate == max(x21$estimate)), which(x22$estimate == max(x22$estimate)),
             which(x23$estimate == max(x23$estimate)), which(x24$estimate == max(x24$estimate)),
             which(x25$estimate == max(x25$estimate)), which(x26$estimate == max(x26$estimate)),
             which(x27$estimate == max(x27$estimate)), which(x28$estimate == max(x28$estimate)))
co4 <- cbind(which(x31$estimate == max(x31$estimate)), which(x32$estimate == max(x32$estimate)),
             which(x33$estimate == max(x33$estimate)), which(x34$estimate == max(x34$estimate)), 
             which(x35$estimate == max(x35$estimate)),which(x36$estimate == max(x36$estimate)),
             which(x37$estimate == max(x37$estimate)), which(x38$estimate == max(x38$estimate)))

# identifying p-value for each subset
p1 <- c(x1$p.value,x2$p.value,x3$p.value,x4$p.value,x5$p.value,x6$p.value,x7$p.value,x8$p.value)
p2 <- c(x11$p.value,x12$p.value,x13$p.value,x14$p.value,x15$p.value,x16$p.value,x17$p.value,x18$p.value)
p3 <- c(x21$p.value,x22$p.value,x23$p.value,x24$p.value,x25$p.value,x26$p.value,x27$p.value,x28$p.value)
p4 <- c(x31$p.value,x32$p.value,x33$p.value,x34$p.value,x35$p.value,x36$p.value,x37$p.value,x38$p.value)

p1
p2
p3
p4

# preparing and generating results for table
c1 <- 0
c2 <- 0
c3 <- 0
c4 <- 0

for (i in 1:8) {
  
  c1[i] <- ifelse(co1[i]==1,ifelse(p1[i]<0.001,"b***",ifelse(p1[i]<0.01,"b**",ifelse(p1[i]<0.05,"b*",ifelse(p1[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p1[i]<0.01,"o**",ifelse(p1[i]<0.05,"o*",ifelse(p1[i]<0.1,"o+","o")))))
  c2[i] <- ifelse(co2[i]==1,ifelse(p2[i]<0.001,"b***",ifelse(p2[i]<0.01,"b**",ifelse(p2[i]<0.05,"b*",ifelse(p2[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p2[i]<0.01,"o**",ifelse(p2[i]<0.05,"o*",ifelse(p2[i]<0.1,"o+","o")))))
  c3[i] <- ifelse(co3[i]==1,ifelse(p3[i]<0.001,"b***",ifelse(p3[i]<0.01,"b**",ifelse(p3[i]<0.05,"b*",ifelse(p3[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p3[i]<0.01,"o**",ifelse(p3[i]<0.05,"o*",ifelse(p3[i]<0.1,"o+","o")))))
  c4[i] <- ifelse(co4[i]==1,ifelse(p4[i]<0.001,"b***",ifelse(p4[i]<0.01,"b**",ifelse(p4[i]<0.05,"b*",ifelse(p4[i]<0.1,"b+","b")))),
                  ifelse(p4[i]<0.001,"o***",ifelse(p4[i]<0.01,"o**",ifelse(p4[i]<0.05,"o*",ifelse(p4[i]<0.1,"o+","o")))))
  
}

c1
c2
c3
c4

# generating first raw table
matte1 <- rbind(c1,c2,c3,c4)
matte1

# adding labels
rownames(matte1) <- c("Full Set","National cases","Post 1990","National Post 1990")
colnames(matte1) <- c("mean LSI | ","mean ENP | ", "mean PR threshold | ", "mean PR district magnitude | " , 
                      "mean share SMD seats | " , "share minority districts | " , "share election est. democracies|",
                      "ethnic fractionalization")
matte1

# export first table with stargazer
install.packages("stargazer")
library(stargazer)

stargazer(matte1,type="html",out="models1.htm")


######################################################################################################################

# redifining data1 and data3 to include Scotland and Wales
data1 <- subset(data,national!=0)
data3 <- subset(data,c(national!=0 & post1990==1))

# Executing tests again for data including only national data but also Scotland and Wales 
x11 <- t.test(data1$lsi[data1$bobw==1],data1$lsi[data1$bobw==0])
x11
x12 <- t.test(data1$enps[data1$bobw==1],data1$enps[data1$bobw==0])
x12
x13 <- t.test(data1$thresholdPR[data1$bobw==1],data1$thresholdPR[data1$bobw==0])
x13
x14 <- t.test(data1$prmedianmagnitude[data1$bobw==1],data1$prmedianmagnitude[data1$bobw==0])
x14
x15 <- t.test(data1$smd_share[data1$bobw==1],data1$smd_share[data1$bobw==0])
x15
x16 <- t.test(data1$minority_districts[data1$bobw==1],data1$minority_districts[data1$bobw==0])
x16
x17 <- t.test(data1$established[data1$bobw==1],data1$established[data1$bobw==0])
x17
x18 <- t.test(data1$ethnic[data1$bobw==1],data1$ethnic[data1$bobw==0])
x18

# Executing tests again for data including only national data and Scotland and Wales AFTER 1990
x31 <- t.test(data3$lsi[data3$bobw==1],data3$lsi[data3$bobw==0])
x31
x32 <- t.test(data3$enps[data3$bobw==1],data3$enps[data3$bobw==0])
x32
x33 <- t.test(data3$thresholdPR[data3$bobw==1],data3$thresholdPR[data3$bobw==0])
x33
x34 <- t.test(data3$prmedianmagnitude[data3$bobw==1],data3$prmedianmagnitude[data3$bobw==0])
x34
x35 <- t.test(data3$smd_share[data3$bobw==1],data3$smd_share[data3$bobw==0])
x35
x36 <- t.test(data3$minority_districts[data3$bobw==1],data3$minority_districts[data3$bobw==0])
x36
x37 <- t.test(data3$established[data3$bobw==1],data3$established[data3$bobw==0])
x37
x38 <- t.test(data3$ethnic[data3$bobw==1],data3$ethnic[data3$bobw==0])
x38

#### preparing and finishing the presentation of the results

# identifying t.test results for each subset
co5 <- cbind(which(x11$estimate == max(x11$estimate)), which(x12$estimate == max(x12$estimate)),
             which(x13$estimate == max(x13$estimate)), which(x14$estimate == max(x14$estimate)),
             which(x15$estimate == max(x15$estimate)), which(x16$estimate == max(x16$estimate))[1],
             which(x17$estimate == max(x17$estimate)), which(x18$estimate == max(x18$estimate)))
co6 <- cbind(which(x31$estimate == max(x31$estimate)), which(x32$estimate == max(x32$estimate)),
             which(x33$estimate == max(x33$estimate)), which(x34$estimate == max(x34$estimate)),
             which(x35$estimate == max(x35$estimate)), which(x36$estimate == max(x36$estimate)),
             which(x37$estimate == max(x37$estimate)), which(x38$estimate == max(x38$estimate)))

# identifying the p-value of each test
p5 <- c(x11$p.value,x12$p.value,x13$p.value,x14$p.value,x15$p.value,x16$p.value,x17$p.value,x18$p.value)
p6 <- c(x31$p.value,x32$p.value,x33$p.value,x34$p.value,x35$p.value,x36$p.value,x37$p.value,x38$p.value)

# preparing and generating results for table
c5 <- 0
c6 <- 0

for (i in 1:8) {
  
  c5[i] <- ifelse(co5[i]==1,ifelse(p1[i]<0.001,"b***",ifelse(p1[i]<0.01,"b**",ifelse(p1[i]<0.05,"b*",ifelse(p1[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p1[i]<0.01,"o**",ifelse(p1[i]<0.05,"o*",ifelse(p1[i]<0.1,"o+","o")))))
  c6[i] <- ifelse(co6[i]==1,ifelse(p2[i]<0.001,"b***",ifelse(p2[i]<0.01,"b**",ifelse(p2[i]<0.05,"b*",ifelse(p2[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p2[i]<0.01,"o**",ifelse(p2[i]<0.05,"o*",ifelse(p2[i]<0.1,"o+","o")))))

}

c5
c6

# generating the second table
matte2 <- rbind(c5,c6)
matte2

rownames(matte2) <- c("National cases (with S+W)","National Post 1990 (with S+W)")
colnames(matte2) <- c("mean LSI | ","mean ENP | ", "mean PR threshold | ", "mean PR district magnitude | " , 
                      "mean share SMD seats | " , "share minority districts | " , "share election est. democracies|",
                      "ethnic fractionalization")
matte2

# export second table
stargazer(matte2,type="html",out="models2.htm")


######################################################################################################################
######################################################################################################################

# redefining full set and subset to exclude cases of Manipulation (manipulation == 0)
datan <- data_epsr_nref
data <- subset(datan,manipulation==0)
data1 <- subset(data,national==1)
data2 <- subset(data,post1990==1)
data3 <- subset(data,c(national==1 & post1990==1))

# executing t-tests with the new full data set
x1 <- t.test(data$lsi[data$bobw==1],data$lsi[data$bobw==0])
x1
x2 <- t.test(data$enps[data$bobw==1],data$enps[data$bobw==0])
x2
x3 <- t.test(data$thresholdPR[data$bobw==1],data$thresholdPR[data$bobw==0])
x3
x4 <- t.test(data$prmedianmagnitude[data$bobw==1],data$prmedianmagnitude[data$bobw==0])
x4
x5 <- t.test(data$smd_share[data$bobw==1],data$smd_share[data$bobw==0])
x5
x6 <- t.test(data$minority_districts[data$bobw==1],data$minority_districts[data$bobw==0])
x6
x7 <- t.test(data$established[data$bobw==1],data$established[data$bobw==0])
x7
x8 <- t.test(data$ethnic[data$bobw==1],data$ethnic[data$bobw==0])
x8

# executing t-tests for the first subset only including national cases WITHOUT Scotland and Wales
x11 <- t.test(data1$lsi[data1$bobw==1],data1$lsi[data1$bobw==0])
x11
x12 <- t.test(data1$enps[data1$bobw==1],data1$enps[data1$bobw==0])
x12
x13 <- t.test(data1$thresholdPR[data1$bobw==1],data1$thresholdPR[data1$bobw==0])
x13
x14 <- t.test(data1$prmedianmagnitude[data1$bobw==1],data1$prmedianmagnitude[data1$bobw==0])
x14
x15 <- t.test(data1$smd_share[data1$bobw==1],data1$smd_share[data1$bobw==0])
x15
x16 <- t.test(data1$minority_districts[data1$bobw==1],data1$minority_districts[data1$bobw==0])
x16
x17 <- t.test(data1$established[data1$bobw==1],data1$established[data1$bobw==0])
x17
x18 <- t.test(data1$ethnic[data1$bobw==1],data1$ethnic[data1$bobw==0])
x18

# executing t-tests for the third subset only including cases AFTER 1990
x21 <- t.test(data2$lsi[data2$bobw==1],data2$lsi[data2$bobw==0])
x21
x22 <- t.test(data2$enps[data2$bobw==1],data2$enps[data2$bobw==0])
x22
x23 <- t.test(data2$thresholdPR[data2$bobw==1],data2$thresholdPR[data2$bobw==0])
x23
x24 <- t.test(data2$prmedianmagnitude[data2$bobw==1],data2$prmedianmagnitude[data2$bobw==0])
x24
x25 <- t.test(data2$smd_share[data2$bobw==1],data2$smd_share[data2$bobw==0])
x25
x26 <- t.test(data2$minority_districts[data2$bobw==1],data2$minority_districts[data2$bobw==0])
x26
x27 <- t.test(data2$established[data2$bobw==1],data2$established[data2$bobw==0])
x27
x28 <- t.test(data2$ethnic[data2$bobw==1],data2$ethnic[data2$bobw==0])
x28

# executing t-tests for the fourth subset only including all national cases after 1990 except Scotland and Wales
x31 <- t.test(data3$lsi[data3$bobw==1],data3$lsi[data3$bobw==0])
x31
x32 <- t.test(data3$enps[data3$bobw==1],data3$enps[data3$bobw==0])
x32
x33 <- t.test(data3$thresholdPR[data3$bobw==1],data3$thresholdPR[data3$bobw==0])
x33
x34 <- t.test(data3$prmedianmagnitude[data3$bobw==1],data3$prmedianmagnitude[data3$bobw==0])
x34
x35 <- t.test(data3$smd_share[data3$bobw==1],data3$smd_share[data3$bobw==0])
x35
x36 <- t.test(data3$minority_districts[data3$bobw==1],data3$minority_districts[data3$bobw==0])
x36
x37 <- t.test(data3$established[data3$bobw==1],data3$established[data3$bobw==0])
x37
x38 <- t.test(data3$ethnic[data3$bobw==1],data3$ethnic[data3$bobw==0])
x38

#### preparing and finishing the presentation of the results

# identifying t.test results for each subset
co7 <- cbind(which(x1$estimate == max(x1$estimate)), which(x2$estimate == max(x2$estimate)),
             which(x3$estimate == max(x3$estimate)), which(x4$estimate == max(x4$estimate)),
             which(x5$estimate == max(x5$estimate)), which(x6$estimate == max(x6$estimate)),
             which(x7$estimate == max(x7$estimate)), which(x8$estimate == max(x8$estimate)))
co8 <- cbind(which(x11$estimate == max(x11$estimate)), which(x12$estimate == max(x12$estimate)),
             which(x13$estimate == max(x13$estimate)), which(x14$estimate == max(x14$estimate)),
             which(x15$estimate == max(x15$estimate)), which(x16$estimate == max(x16$estimate)),
             which(x17$estimate == max(x17$estimate)), which(x18$estimate == max(x18$estimate)))
co9 <- cbind(which(x21$estimate == max(x21$estimate)), which(x22$estimate == max(x22$estimate)),
             which(x23$estimate == max(x23$estimate)), which(x24$estimate == max(x24$estimate)), 
             which(x25$estimate == max(x25$estimate)), which(x26$estimate == max(x26$estimate)),
             which(x27$estimate == max(x27$estimate)), which(x28$estimate == max(x28$estimate)))
co10 <- cbind(which(x31$estimate == max(x31$estimate)), which(x32$estimate == max(x32$estimate)),
              which(x33$estimate == max(x33$estimate)), which(x34$estimate == max(x34$estimate)), 
              which(x35$estimate == max(x35$estimate)),which(x36$estimate == max(x36$estimate)),
              which(x37$estimate == max(x37$estimate)), which(x38$estimate == max(x38$estimate)))

# identifying the p-value of each test
p7 <- c(x1$p.value,x2$p.value,x3$p.value,x4$p.value,x5$p.value,x6$p.value,x7$p.value,x8$p.value)
p8 <- c(x11$p.value,x12$p.value,x13$p.value,x14$p.value,x15$p.value,x16$p.value,x17$p.value,x18$p.value)
p9 <- c(x21$p.value,x22$p.value,x23$p.value,x24$p.value,x25$p.value,x26$p.value,x27$p.value,x28$p.value)
p10 <- c(x31$p.value,x32$p.value,x33$p.value,x34$p.value,x35$p.value,x36$p.value,x37$p.value,x38$p.value)

# preparing the results for presentation
c7 <- 0
c8 <- 0
c9 <- 0
c10 <- 0

for (i in 1:8) {
  
  c7[i] <- ifelse(co7[i]==1,ifelse(p1[i]<0.001,"b***",ifelse(p1[i]<0.01,"b**",ifelse(p1[i]<0.05,"b*",ifelse(p1[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p1[i]<0.01,"o**",ifelse(p1[i]<0.05,"o*",ifelse(p1[i]<0.1,"o+","o")))))
  c8[i] <- ifelse(co8[i]==1,ifelse(p2[i]<0.001,"b***",ifelse(p2[i]<0.01,"b**",ifelse(p2[i]<0.05,"b*",ifelse(p2[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p2[i]<0.01,"o**",ifelse(p2[i]<0.05,"o*",ifelse(p2[i]<0.1,"o+","o")))))
  c9[i] <- ifelse(co9[i]==1,ifelse(p3[i]<0.001,"b***",ifelse(p3[i]<0.01,"b**",ifelse(p3[i]<0.05,"b*",ifelse(p3[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p3[i]<0.01,"o**",ifelse(p3[i]<0.05,"o*",ifelse(p3[i]<0.1,"o+","o")))))
  c10[i] <- ifelse(co10[i]==1,ifelse(p4[i]<0.001,"b***",ifelse(p4[i]<0.01,"b**",ifelse(p4[i]<0.05,"b*",ifelse(p4[i]<0.1,"b+","b")))),
                  ifelse(p4[i]<0.001,"o***",ifelse(p4[i]<0.01,"o**",ifelse(p4[i]<0.05,"o*",ifelse(p4[i]<0.1,"o+","o")))))
  
}

c7
c8
c9
c10

# generating third result table
matte3 <- rbind(c7,c8,c9,c10)
matte3

rownames(matte3) <- c("Full Set","National cases","Post 1990","National Post 1990")
colnames(matte3) <- c("mean LSI | ","mean ENP | ", "mean PR threshold | ", "mean PR district magnitude | " , 
                      "mean share SMD seats | " , "share minority districts | " , "share election est. democracies|",
                      "ethnic fractionalization")
matte3

# export third table
stargazer(matte3,type="html",out="models3.htm")


######################################################################################################################

# redefining subsets 1 and 3 to include Scotland and Wales
data1 <- subset(data,national!=0)
data3 <- subset(data,c(national!=0 & post1990==1))

# executing t-tests for the redefined first subset which includes all national cases as well as Scotland and Wales
x11 <- t.test(data1$lsi[data1$bobw==1],data1$lsi[data1$bobw==0])
x11
x12 <- t.test(data1$enps[data1$bobw==1],data1$enps[data1$bobw==0])
x12
x13 <- t.test(data1$thresholdPR[data1$bobw==1],data1$thresholdPR[data1$bobw==0])
x13
x14 <- t.test(data1$prmedianmagnitude[data1$bobw==1],data1$prmedianmagnitude[data1$bobw==0])
x14
x15 <- t.test(data1$smd_share[data1$bobw==1],data1$smd_share[data1$bobw==0])
x15
x16 <- t.test(data1$minority_districts[data1$bobw==1],data1$minority_districts[data1$bobw==0])
x16
x17 <- t.test(data1$established[data1$bobw==1],data1$established[data1$bobw==0])
x17
x18 <- t.test(data1$ethnic[data1$bobw==1],data1$ethnic[data1$bobw==0])
x18

# executing t-tests for the redefined third subset including all national cases after 1990 as well as Scotland and
#   Wales
x31 <- t.test(data3$lsi[data3$bobw==1],data3$lsi[data3$bobw==0])
x31
x32 <- t.test(data3$enps[data3$bobw==1],data3$enps[data3$bobw==0])
x32
x33 <- t.test(data3$thresholdPR[data3$bobw==1],data3$thresholdPR[data3$bobw==0])
x33
x34 <- t.test(data3$prmedianmagnitude[data3$bobw==1],data3$prmedianmagnitude[data3$bobw==0])
x34
x35 <- t.test(data3$smd_share[data3$bobw==1],data3$smd_share[data3$bobw==0])
x35
x36 <- t.test(data3$minority_districts[data3$bobw==1],data3$minority_districts[data3$bobw==0])
x36
x37 <- t.test(data3$established[data3$bobw==1],data3$established[data3$bobw==0])
x37
x38 <- t.test(data3$ethnic[data3$bobw==1],data3$ethnic[data3$bobw==0])
x38

#### preparing and finishing the presentation of the results

# identifying t.test results for each subset 
co11 <- cbind(which(x11$estimate == max(x11$estimate)), which(x12$estimate == max(x12$estimate)),
              which(x13$estimate == max(x13$estimate)), which(x14$estimate == max(x14$estimate)),
              which(x15$estimate == max(x15$estimate)), which(x16$estimate == max(x16$estimate)),
              which(x17$estimate == max(x17$estimate)), which(x18$estimate == max(x18$estimate)))
co12 <- cbind(which(x31$estimate == max(x31$estimate)), which(x32$estimate == max(x32$estimate)),
              which(x33$estimate == max(x33$estimate)), which(x34$estimate == max(x34$estimate)),
              which(x35$estimate == max(x35$estimate)),which(x36$estimate == max(x36$estimate)),
              which(x37$estimate == max(x37$estimate)),which(x38$estimate == max(x38$estimate)))

# identifying the p-value of each test
p11 <- c(x11$p.value,x12$p.value,x13$p.value,x14$p.value,x15$p.value,x16$p.value,x17$p.value,x18$p.value)
p12 <- c(x31$p.value,x32$p.value,x33$p.value,x34$p.value,x35$p.value,x36$p.value,x37$p.value,x38$p.value)

# preparing the results for presentation 
c11 <- 0
c12 <- 0

for (i in 1:8) {
  
  c11[i] <- ifelse(co11[i]==1,ifelse(p1[i]<0.001,"b***",ifelse(p1[i]<0.01,"b**",ifelse(p1[i]<0.05,"b*",ifelse(p1[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p1[i]<0.01,"o**",ifelse(p1[i]<0.05,"o*",ifelse(p1[i]<0.1,"o+","o")))))
  c12[i] <- ifelse(co12[i]==1,ifelse(p2[i]<0.001,"b***",ifelse(p2[i]<0.01,"b**",ifelse(p2[i]<0.05,"b*",ifelse(p2[i]<0.1,"b+","b")))),
                  ifelse(p1[i]<0.001,"o***",ifelse(p2[i]<0.01,"o**",ifelse(p2[i]<0.05,"o*",ifelse(p2[i]<0.1,"o+","o")))))

}

c11
c12

# generating the fourth result table
matte4 <- rbind(c11,c12)
matte4

rownames(matte4) <- c("National cases (with S+W)","National Post 1990 (with S+W)")
colnames(matte4) <- c("mean LSI | ","mean ENP | ", "mean PR threshold | ", "mean PR district magnitude | " , 
                      "mean share SMD seats | " , "share minority districts | " , "share election est. democracies|",
                      "ethnic fractionalization")
matte4

# export fourth table
stargazer(matte4,type="html",out="models4.htm")


######################################################################################################################


# generating the final table, which brings together all previos results
matte <- rbind(c("","","","","","","",""),matte1,matte2,
               c("","","","","","","",""),matte3,matte4)

rownames(matte) <- c("With Manipulation","Full Set","National cases","Post 1990","National Post 1990",
                     "National cases (with S+W)","National Post 1990 (with S+W)","Without Manipulation",
                     "Full Set","National cases","Post 1990","National Post 1990",
                     "National cases (with S+W)","National Post 1990 (with S+W)")

# export the final table
stargazer(matte,type="html",out="all_models.htm")