MandlaSibanda196/gist:ca6c2ddcef84e67c65d0462bd5ecb8d9

## gistfile1.txt

## loading necessary packages

library(Zelig)
library("ZeligChoice")
library(xtable)
library(reshape)
library(apsrtable)
library(stargazer)
library(rms)
library(rgenoud)
library(Matching)
library(quantreg)
library(foreign)


########################################################
## Loading the Data
########################################################

## This loads three data sets
## 1. Gender cases
## 2. Cleaned by-judge data
## 3. Non-gender cases, downloaded from Kuersten and Haire's website
## (Commented out; their data are available from their website)

# ## Gender cases only
women.cases <- read.csv("glynn_sen_daughters_by_case_1.csv", stringsAsFactors = FALSE)

# ## Removing male plaintiffs:
women.cases <- subset(women.cases, femplaintiff == 1)

women.cases <- subset(women.cases, area == "employment" | area == "Title IX" | area == "pregnancy" | area == "abortion" | area == "reproductive rights")

women.cases$area <- factor(women.cases$area, levels = c("abortion","employment","pregnancy","reproductive rights","Title IX"))

# ## All cases, including non-gender cases (from Ashlyn Kuersten/Susan Haire's coding)
# cases <- read.csv("cases.for.analysis.csv")

judge.means <- read.csv("glynn_sen_daughters_by_judge.csv", stringsAsFactors = FALSE)

all <- subset(judge.means, girls != "NaN") # subsets judge
# data to those for
# whom we have fertility data

######################################################
## Table 1: Number of children and girls
######################################################
aa <- table(all$republican, all$child)
bb <- table(all$republican, all$girls)

## and now for the table:

xtable(aa)
xtable(bb)

########################################################
## Table 2: Judge Demographics
########################################################

dems <- subset(all, republican == 0)
reps <- subset(all, republican == 1)
women <- subset(all, woman == 1)
men <- subset(all, woman == 0)


mean.kids <-  cbind(mean(na.omit(all$child)),
                    mean(na.omit(dems$child)),
                    mean(na.omit(reps$child)),
                    mean(na.omit(women$child)),
                    mean(na.omit(men$child))
)

mean.girls <- cbind(mean(na.omit(all$girls)),
                    mean(na.omit(dems$girls)),
                    mean(na.omit(reps$girls)),
                    mean(na.omit(women$girls)),
                    mean(na.omit(men$girls))
)

prop.zero <- cbind(prop.table(table(na.omit(all$child)))[1],
                   prop.table(table(na.omit(dems$child)))[1],
                   prop.table(table(na.omit(reps$child)))[1],
                   prop.table(table(na.omit(women$child)))[1],
                   prop.table(table(na.omit(men$child)))[1]
)


prop.one <- cbind(prop.table(table(na.omit(all$child)))[2],
                  prop.table(table(na.omit(dems$child)))[2],
                  prop.table(table(na.omit(reps$child)))[2],
                  prop.table(table(na.omit(women$child)))[2],
                  prop.table(table(na.omit(men$child)))[2]
)

prop.two <- cbind(prop.table(table(na.omit(all$child)))[3],
                  prop.table(table(na.omit(dems$child)))[3],
                  prop.table(table(na.omit(reps$child)))[3],
                  prop.table(table(na.omit(women$child)))[3],
                  prop.table(table(na.omit(men$child)))[3]
)

prop.three <- cbind(prop.table(table(na.omit(all$child)))[4],
                    prop.table(table(na.omit(dems$child)))[4],
                    prop.table(table(na.omit(reps$child)))[4],
                    prop.table(table(na.omit(women$child)))[4],
                    prop.table(table(na.omit(men$child)))[4]
)

prop.four <- cbind(prop.table(table(na.omit(all$child)))[5],
                   prop.table(table(na.omit(dems$child)))[5],
                   prop.table(table(na.omit(reps$child)))[5],
                   prop.table(table(na.omit(women$child)))[5],
                   prop.table(table(na.omit(men$child)))[5]
)

prop.five <- cbind(prop.table(table(na.omit(all$child)))[6],
                   prop.table(table(na.omit(dems$child)))[6],
                   prop.table(table(na.omit(reps$child)))[6],
                   prop.table(table(na.omit(women$child)))[6],
                   prop.table(table(na.omit(men$child)))[6]
)

aa <- table(na.omit(all$child))
plus6.all <- sum(aa[7:length(aa)])/sum(aa)

bb <- table(na.omit(dems$child))
plus6.dems <- sum(bb[7:length(bb)])/sum(bb)

cc <- table(na.omit(reps$child))
plus6.reps <- sum(cc[7:length(cc)])/sum(cc)

dd <- table(na.omit(women$child))
plus6.women <- sum(dd[7:length(dd)])/sum(dd)

ee <- table(na.omit(men$child))
plus6.men <- sum(ee[7:length(ee)])/sum(ee)

prop.six.or.greater <- cbind(plus6.all, plus6.dems, plus6.reps, plus6.women, plus6.men)

mean.female <- cbind(mean(na.omit(all$woman)),
                     mean(na.omit(dems$woman)),
                     mean(na.omit(reps$woman)),
                     mean(na.omit(women$woman)),
                     mean(na.omit(men$woman))
)

mean.rep <- cbind(mean(na.omit(all$republican)),
                  mean(na.omit(dems$republican)),
                  mean(na.omit(reps$republican)),
                  mean(na.omit(women$republican)),
                  mean(na.omit(men$republican))
)


mean.white <- cbind(mean(na.omit(all$race == 1)),
                    mean(na.omit(dems$race == 1)),
                    mean(na.omit(reps$race == 1)),
                    mean(na.omit(women$race == 1)),
                    mean(na.omit(men$race == 1))
)

mean.yearb <- cbind(mean(na.omit(all$yearb)),
                    mean(na.omit(dems$yearb)),
                    mean(na.omit(reps$yearb)),
                    mean(na.omit(women$yearb)),
                    mean(na.omit(men$yearb))
)


no_judges <- cbind(nrow(all), nrow(dems), nrow(reps), nrow(women), nrow(men))

demographic_table <- rbind(mean.kids, mean.girls, prop.zero, prop.one,
                           prop.two, prop.three, prop.four, prop.five, prop.six.or.greater, mean.female, mean.rep, mean.white, mean.yearb, no_judges)

colnames(demographic_table) <- c("All", "Democrats", "Republicans","Women", "Men")
rownames(demographic_table) <- c("Mean No. Children", "Mean No. Girls",
                                 "Proportion who have 0 children","1 children",
                                 "2 children", "3 children",
                                 "4 children", "5 Children", "6 Children or More", "Proportion Female", "Proportion Republican",
                                 "Proportion White", "Mean Year Born", "N")
xtable(demographic_table, digits = 2, caption = "Demographics of U.S. Court of Appeal Judges who voted on gender-related cases (1996-2002)", label = "t:statsgender", align = "l|ccccc")

########################################################
## Calculating the Weights (Number of Cases)
########################################################

no_cases <- matrix(data = NA, nrow = nrow(judge.means), ncol = 1)
for(i in 1:length(no_cases)){
  no_cases[i] <- nrow(women.cases[which(women.cases$name == judge.means$name[i]),])
}

judge.means <- cbind(judge.means, no_cases)

## total number of cases we are working with

sum(judge.means$no_cases) 	#  2,674 as reported in the paper

##########################################################
##########################################################
########## Calculating the Outcome Var ###################
##########################################################

no_liberalvote <- matrix(data = NA, nrow = nrow(judge.means), ncol = 1)
for(i in 1:length(no_liberalvote)){
  stuff <- women.cases[which(women.cases$name == judge.means$name[i]),]
  no_liberalvote[i] <- nrow(subset(stuff, vote == 2 | vote == 3))
}

lib_vote_share <- no_liberalvote/no_cases
lib_vote_share

judge.means <- cbind(judge.means, no_liberalvote, lib_vote_share)
judge.means <- subset(judge.means, girls != "NaN")

## treated group:
treated.means <- subset(judge.means, girls > 0)
## control group:
control.means <- subset(judge.means, girls = 0)

control.means$child

##### histograms ####
hist(treated.means$child)
hist(control.means$child)
hist(treated.means$girls)
hist(control.means$girls)


A <- summary(treated.means$girls)
#ANSWER MEAN = 1.71
B <- summary(control.means$girls)
#ANSWER MEAN = 1.237

treatment_effect <-  A - B
treatment_effect
#ANSWER
#treatment effect = 0.4733

lmA1 <- lm(judge.means$lib_vote_share ~ judge.means$girls)
summary(lmA1)
confint(lmA1, level=0.95)

#ANSWER
#(Intercept)        0.35648869 0.45768779
#judge.means$girls -0.01027287 0.05224176
# p-value: 0.1872 THEREFORE NOT STATISTICALLY SIGNIFICANT

#######################################################
## Figure 1 DEMOCRAT AND REPUBLICANS
########################################################

## just women:
women.means <- subset(judge.means, woman == 1)

## just men:
men.means <- subset(judge.means, woman == 0)

## just republicans:
rep.means <- subset(judge.means, republican == 1)

## just democrats
dem.means <- subset(judge.means, republican == 0)

#pdf(file= "lib_votes.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(density(judge.means$lib_vote_share),
     xlim = c(-.4,1.4), ylim = c(0,2.5),
     ylab = "", xlab = "Proportion of Cases Decided in a Feminist Direction",
     yaxt = "n",
     bty = "n",
     main = "",
     col = "black", lwd = 2)
lines(density(rep.means$lib_vote_share),
      col = "firebrick", lwd = 2, lty  = 2)
lines(density(dem.means$lib_vote_share),
      col = "dodgerblue", lwd = 2, lty  = 3)
abline(v = .5, col = "grey50", lty = 2)
text(x = .5, y = 2.4, "Less Feminist", col = "grey50", pos = 2, cex = 0.9)
text(x = .5, y = 2.4, "More Feminist", col = "grey50", pos = 4, cex = 0.9)
text(x = .25, y = 1.7, "Republicans", pos = 2, cex = 0.9)
text(x = .7, y = 1, "Democrats", pos = 4, cex = 0.9)
text(x = .075, y = .6, "All", pos = 4, cex = 0.9)
#dev.off()

#######################################################
## Figure 2 CONTROL AND TREATMENT COMPARISON
########################################################

#pdf(file= "lib_votes.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(density(judge.means$lib_vote_share),
     xlim = c(-.4,1.4), ylim = c(0,2.5),
     ylab = "", xlab = "Proportion of Cases Decided in a Feminist Direction",
     yaxt = "n",
     bty = "n",
     main = "",
     col = "black", lwd = 2)
lines(density(treated.means$lib_vote_share),
      col = "firebrick", lwd = 2, lty  = 2)
lines(density(control.means$lib_vote_share),
      col = "dodgerblue", lwd = 2, lty  = 3)
abline(v = .5, col = "grey50", lty = 2)
text(x = .5, y = 2.4, "Less Feminist", col = "grey50", pos = 2, cex = 0.9)
text(x = .5, y = 2.4, "More Feminist", col = "grey50", pos = 4, cex = 0.9)
text(x = .25, y = 1.7, "With NO girls", pos = 2, cex = 0.9)
text(x = .7, y = 1, "With girls", pos = 4, cex = 0.9)
text(x = .075, y = .6, "All", pos = 4, cex = 0.9)
#dev.off()

##########################################################
###################### MATCHING ########################
##########################################################

judge.means$girls_A <- ifelse(judge.means$girls>0, 1, 0)

#ANSWER - HAD TO CHANGE THE GIRLS VARIABLE TO A LOGICAL VARIABLE OF ONES AND ZEROS AND THEN APPEND IT AS girls_A TO THE MAIN DATA

# Estimate the propensity model
glm1  <- glm(girls_A ~ republican + child + woman
             , family=binomial, data=judge.means)

#save data objects

X  <- glm1$fitted
Y  <- judge.means$lib_vote_share
Tr  <- judge.means$girls_A

# one-to-one matching with replacement (the "M=1" option).
# Estimating the treatment effect on the treated (the "estimand" option defaults to ATT).

rr  <- Match(Y=Y, Tr=Tr, X=X, M=1);
summary(rr)
rr_data <- rr$mdata

#Finding balance
mb  <- MatchBalance(girls_A ~ republican + child + woman
                    ,data=judge.means, match.out=rr, nboots=500)

lmC <- lm(rr_data$Y ~ rr_data$Tr, data = rr_data)
summary(lmC)
coef(lmC)[2]
confint(lmC, level = 0.95)

#ANSWER = PEstimate = 0.356784
#ANSWER = (Intercept) 0.33921008 0.3743573
#rr_data$Tr  0.06294106 0.1126467
#p-value: 6.596e-12

###multivariate matching procedure that some of
###the covariates

Y2 <- judge.means$lib_vote_share
Tr2 <- judge.means$girls_A
X2 <- cbind(judge.means$child,
            judge.means$republican == 1,
            judge.means$republican == 0,
            judge.means$woman == 1,
            judge.means$woman == 0,
            X)

mout.mv <- Match(Y=Y2, Tr=Tr2, X=X2, M=1)
summary(mout.mv)

# Finding Match Balance of multivariate matching
mb.mv <- MatchBalance(girls_A ~ child + republican + woman
                      ,data=judge.means, match.out=mout.mv, nboots=500)

# Estimating estimates with linear regression
mv_data <- mout.mv$mdata
lm.mv <- lm(mv_data$Y ~ mv_data$Tr, data = mv_data)
summary(lm.mv)

coef(lm.mv)[2]

confint(lm.mv, level = 0.95)

#ANSWER = PEstimate = 0.36091
#ANSWER = (Intercept) 0.34317928 0.3786474
#mv_data$Tr  0.05975124 0.1099108
#p-value: 4.734e-11


# Genetic Matching with multivariate###
Y4 <- judge.means$lib_vote_share
Tr4 <- judge.means$girls_A
X_all <- cbind(judge.means$child,
               judge.means$republican == 1,
               judge.means$republican == 0,
               judge.means$woman == 1,
               judge.means$woman == 0)


gout.mv <- GenMatch(Tr=Tr4, X=X_all, M=1, estimand = "ATT", pop.size = 100, max.generations = 20, wait.generations = 7, caliper=0.25)

mout.gen.mv <- Match(Y=Y4, Tr=Tr4, X=X_all, estimand="ATT", M=1, Weight.matrix = gout.mv)

# Finding Match Balance of multivariate matching
mb.gen.mv <- MatchBalance(girls_A ~ child + republican + woman
                          ,data=judge.means, match.out=mout.gen.mv, nboots=500)

# Estimating the effect with linear regression
gen_mv_data <- mout.gen.mv$mdata
lm.gen.mv <- lm(gen_mv_data$Y ~ gen_mv_data$Tr, data=gen_mv_data)
summary(lm.gen.mv)

coef(lm.gen.mv)[2]
#ANSWER TREATMENT EFFECT = 0.360839

confint(lm.gen.mv, level = 0.95)
#ANSWER = (Intercept)    0.34315195 0.3785267
#gen_mv_data$Tr 0.05872448 0.1087519

	## loading necessary packages

	library(Zelig)
	library("ZeligChoice")
	library(xtable)
	library(reshape)
	library(apsrtable)
	library(stargazer)
	library(rms)
	library(rgenoud)
	library(Matching)
	library(quantreg)
	library(foreign)


	########################################################
	## Loading the Data
	########################################################

	## This loads three data sets
	## 1. Gender cases
	## 2. Cleaned by-judge data
	## 3. Non-gender cases, downloaded from Kuersten and Haire's website
	## (Commented out; their data are available from their website)

	# ## Gender cases only
	women.cases <- read.csv("glynn_sen_daughters_by_case_1.csv", stringsAsFactors = FALSE)

	# ## Removing male plaintiffs:
	women.cases <- subset(women.cases, femplaintiff == 1)

	women.cases <- subset(women.cases, area == "employment" \| area == "Title IX" \| area == "pregnancy" \| area == "abortion" \| area == "reproductive rights")

	women.cases$area <- factor(women.cases$area, levels = c("abortion","employment","pregnancy","reproductive rights","Title IX"))

	# ## All cases, including non-gender cases (from Ashlyn Kuersten/Susan Haire's coding)
	# cases <- read.csv("cases.for.analysis.csv")

	judge.means <- read.csv("glynn_sen_daughters_by_judge.csv", stringsAsFactors = FALSE)

	all <- subset(judge.means, girls != "NaN") # subsets judge
	# data to those for
	# whom we have fertility data

	######################################################
	## Table 1: Number of children and girls
	######################################################
	aa <- table(all$republican, all$child)
	bb <- table(all$republican, all$girls)

	## and now for the table:

	xtable(aa)
	xtable(bb)

	########################################################
	## Table 2: Judge Demographics
	########################################################

	dems <- subset(all, republican == 0)
	reps <- subset(all, republican == 1)
	women <- subset(all, woman == 1)
	men <- subset(all, woman == 0)


	mean.kids <- cbind(mean(na.omit(all$child)),
	mean(na.omit(dems$child)),
	mean(na.omit(reps$child)),
	mean(na.omit(women$child)),
	mean(na.omit(men$child))
	)

	mean.girls <- cbind(mean(na.omit(all$girls)),
	mean(na.omit(dems$girls)),
	mean(na.omit(reps$girls)),
	mean(na.omit(women$girls)),
	mean(na.omit(men$girls))
	)

	prop.zero <- cbind(prop.table(table(na.omit(all$child)))[1],
	prop.table(table(na.omit(dems$child)))[1],
	prop.table(table(na.omit(reps$child)))[1],
	prop.table(table(na.omit(women$child)))[1],
	prop.table(table(na.omit(men$child)))[1]
	)


	prop.one <- cbind(prop.table(table(na.omit(all$child)))[2],
	prop.table(table(na.omit(dems$child)))[2],
	prop.table(table(na.omit(reps$child)))[2],
	prop.table(table(na.omit(women$child)))[2],
	prop.table(table(na.omit(men$child)))[2]
	)

	prop.two <- cbind(prop.table(table(na.omit(all$child)))[3],
	prop.table(table(na.omit(dems$child)))[3],
	prop.table(table(na.omit(reps$child)))[3],
	prop.table(table(na.omit(women$child)))[3],
	prop.table(table(na.omit(men$child)))[3]
	)

	prop.three <- cbind(prop.table(table(na.omit(all$child)))[4],
	prop.table(table(na.omit(dems$child)))[4],
	prop.table(table(na.omit(reps$child)))[4],
	prop.table(table(na.omit(women$child)))[4],
	prop.table(table(na.omit(men$child)))[4]
	)

	prop.four <- cbind(prop.table(table(na.omit(all$child)))[5],
	prop.table(table(na.omit(dems$child)))[5],
	prop.table(table(na.omit(reps$child)))[5],
	prop.table(table(na.omit(women$child)))[5],
	prop.table(table(na.omit(men$child)))[5]
	)

	prop.five <- cbind(prop.table(table(na.omit(all$child)))[6],
	prop.table(table(na.omit(dems$child)))[6],
	prop.table(table(na.omit(reps$child)))[6],
	prop.table(table(na.omit(women$child)))[6],
	prop.table(table(na.omit(men$child)))[6]
	)

	aa <- table(na.omit(all$child))
	plus6.all <- sum(aa[7:length(aa)])/sum(aa)

	bb <- table(na.omit(dems$child))
	plus6.dems <- sum(bb[7:length(bb)])/sum(bb)

	cc <- table(na.omit(reps$child))
	plus6.reps <- sum(cc[7:length(cc)])/sum(cc)

	dd <- table(na.omit(women$child))
	plus6.women <- sum(dd[7:length(dd)])/sum(dd)

	ee <- table(na.omit(men$child))
	plus6.men <- sum(ee[7:length(ee)])/sum(ee)

	prop.six.or.greater <- cbind(plus6.all, plus6.dems, plus6.reps, plus6.women, plus6.men)

	mean.female <- cbind(mean(na.omit(all$woman)),
	mean(na.omit(dems$woman)),
	mean(na.omit(reps$woman)),
	mean(na.omit(women$woman)),
	mean(na.omit(men$woman))
	)

	mean.rep <- cbind(mean(na.omit(all$republican)),
	mean(na.omit(dems$republican)),
	mean(na.omit(reps$republican)),
	mean(na.omit(women$republican)),
	mean(na.omit(men$republican))
	)


	mean.white <- cbind(mean(na.omit(all$race == 1)),
	mean(na.omit(dems$race == 1)),
	mean(na.omit(reps$race == 1)),
	mean(na.omit(women$race == 1)),
	mean(na.omit(men$race == 1))
	)

	mean.yearb <- cbind(mean(na.omit(all$yearb)),
	mean(na.omit(dems$yearb)),
	mean(na.omit(reps$yearb)),
	mean(na.omit(women$yearb)),
	mean(na.omit(men$yearb))
	)


	no_judges <- cbind(nrow(all), nrow(dems), nrow(reps), nrow(women), nrow(men))

	demographic_table <- rbind(mean.kids, mean.girls, prop.zero, prop.one,
	prop.two, prop.three, prop.four, prop.five, prop.six.or.greater, mean.female, mean.rep, mean.white, mean.yearb, no_judges)

	colnames(demographic_table) <- c("All", "Democrats", "Republicans","Women", "Men")
	rownames(demographic_table) <- c("Mean No. Children", "Mean No. Girls",
	"Proportion who have 0 children","1 children",
	"2 children", "3 children",
	"4 children", "5 Children", "6 Children or More", "Proportion Female", "Proportion Republican",
	"Proportion White", "Mean Year Born", "N")
	xtable(demographic_table, digits = 2, caption = "Demographics of U.S. Court of Appeal Judges who voted on gender-related cases (1996-2002)", label = "t:statsgender", align = "l\|ccccc")

	########################################################
	## Calculating the Weights (Number of Cases)
	########################################################

	no_cases <- matrix(data = NA, nrow = nrow(judge.means), ncol = 1)
	for(i in 1:length(no_cases)){
	no_cases[i] <- nrow(women.cases[which(women.cases$name == judge.means$name[i]),])
	}

	judge.means <- cbind(judge.means, no_cases)

	## total number of cases we are working with

	sum(judge.means$no_cases) # 2,674 as reported in the paper

	##########################################################
	##########################################################
	########## Calculating the Outcome Var ###################
	##########################################################

	no_liberalvote <- matrix(data = NA, nrow = nrow(judge.means), ncol = 1)
	for(i in 1:length(no_liberalvote)){
	stuff <- women.cases[which(women.cases$name == judge.means$name[i]),]
	no_liberalvote[i] <- nrow(subset(stuff, vote == 2 \| vote == 3))
	}

	lib_vote_share <- no_liberalvote/no_cases
	lib_vote_share

	judge.means <- cbind(judge.means, no_liberalvote, lib_vote_share)
	judge.means <- subset(judge.means, girls != "NaN")

	## treated group:
	treated.means <- subset(judge.means, girls > 0)
	## control group:
	control.means <- subset(judge.means, girls = 0)

	control.means$child

	##### histograms ####
	hist(treated.means$child)
	hist(control.means$child)
	hist(treated.means$girls)
	hist(control.means$girls)


	A <- summary(treated.means$girls)
	#ANSWER MEAN = 1.71
	B <- summary(control.means$girls)
	#ANSWER MEAN = 1.237

	treatment_effect <- A - B
	treatment_effect
	#ANSWER
	#treatment effect = 0.4733

	lmA1 <- lm(judge.means$lib_vote_share ~ judge.means$girls)
	summary(lmA1)
	confint(lmA1, level=0.95)

	#ANSWER
	#(Intercept) 0.35648869 0.45768779
	#judge.means$girls -0.01027287 0.05224176
	# p-value: 0.1872 THEREFORE NOT STATISTICALLY SIGNIFICANT

	#######################################################
	## Figure 1 DEMOCRAT AND REPUBLICANS
	########################################################

	## just women:
	women.means <- subset(judge.means, woman == 1)

	## just men:
	men.means <- subset(judge.means, woman == 0)

	## just republicans:
	rep.means <- subset(judge.means, republican == 1)

	## just democrats
	dem.means <- subset(judge.means, republican == 0)

	#pdf(file= "lib_votes.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
	plot(density(judge.means$lib_vote_share),
	xlim = c(-.4,1.4), ylim = c(0,2.5),
	ylab = "", xlab = "Proportion of Cases Decided in a Feminist Direction",
	yaxt = "n",
	bty = "n",
	main = "",
	col = "black", lwd = 2)
	lines(density(rep.means$lib_vote_share),
	col = "firebrick", lwd = 2, lty = 2)
	lines(density(dem.means$lib_vote_share),
	col = "dodgerblue", lwd = 2, lty = 3)
	abline(v = .5, col = "grey50", lty = 2)
	text(x = .5, y = 2.4, "Less Feminist", col = "grey50", pos = 2, cex = 0.9)
	text(x = .5, y = 2.4, "More Feminist", col = "grey50", pos = 4, cex = 0.9)
	text(x = .25, y = 1.7, "Republicans", pos = 2, cex = 0.9)
	text(x = .7, y = 1, "Democrats", pos = 4, cex = 0.9)
	text(x = .075, y = .6, "All", pos = 4, cex = 0.9)
	#dev.off()

	#######################################################
	## Figure 2 CONTROL AND TREATMENT COMPARISON
	########################################################

	#pdf(file= "lib_votes.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
	plot(density(judge.means$lib_vote_share),
	xlim = c(-.4,1.4), ylim = c(0,2.5),
	ylab = "", xlab = "Proportion of Cases Decided in a Feminist Direction",
	yaxt = "n",
	bty = "n",
	main = "",
	col = "black", lwd = 2)
	lines(density(treated.means$lib_vote_share),
	col = "firebrick", lwd = 2, lty = 2)
	lines(density(control.means$lib_vote_share),
	col = "dodgerblue", lwd = 2, lty = 3)
	abline(v = .5, col = "grey50", lty = 2)
	text(x = .5, y = 2.4, "Less Feminist", col = "grey50", pos = 2, cex = 0.9)
	text(x = .5, y = 2.4, "More Feminist", col = "grey50", pos = 4, cex = 0.9)
	text(x = .25, y = 1.7, "With NO girls", pos = 2, cex = 0.9)
	text(x = .7, y = 1, "With girls", pos = 4, cex = 0.9)
	text(x = .075, y = .6, "All", pos = 4, cex = 0.9)
	#dev.off()

	##########################################################
	###################### MATCHING ########################
	##########################################################

	judge.means$girls_A <- ifelse(judge.means$girls>0, 1, 0)

	#ANSWER - HAD TO CHANGE THE GIRLS VARIABLE TO A LOGICAL VARIABLE OF ONES AND ZEROS AND THEN APPEND IT AS girls_A TO THE MAIN DATA

	# Estimate the propensity model
	glm1 <- glm(girls_A ~ republican + child + woman
	, family=binomial, data=judge.means)

	#save data objects

	X <- glm1$fitted
	Y <- judge.means$lib_vote_share
	Tr <- judge.means$girls_A

	# one-to-one matching with replacement (the "M=1" option).
	# Estimating the treatment effect on the treated (the "estimand" option defaults to ATT).

	rr <- Match(Y=Y, Tr=Tr, X=X, M=1);
	summary(rr)
	rr_data <- rr$mdata

	#Finding balance
	mb <- MatchBalance(girls_A ~ republican + child + woman
	,data=judge.means, match.out=rr, nboots=500)

	lmC <- lm(rr_data$Y ~ rr_data$Tr, data = rr_data)
	summary(lmC)
	coef(lmC)[2]
	confint(lmC, level = 0.95)

	#ANSWER = PEstimate = 0.356784
	#ANSWER = (Intercept) 0.33921008 0.3743573
	#rr_data$Tr 0.06294106 0.1126467
	#p-value: 6.596e-12

	###multivariate matching procedure that some of
	###the covariates

	Y2 <- judge.means$lib_vote_share
	Tr2 <- judge.means$girls_A
	X2 <- cbind(judge.means$child,
	judge.means$republican == 1,
	judge.means$republican == 0,
	judge.means$woman == 1,
	judge.means$woman == 0,
	X)

	mout.mv <- Match(Y=Y2, Tr=Tr2, X=X2, M=1)
	summary(mout.mv)

	# Finding Match Balance of multivariate matching
	mb.mv <- MatchBalance(girls_A ~ child + republican + woman
	,data=judge.means, match.out=mout.mv, nboots=500)

	# Estimating estimates with linear regression
	mv_data <- mout.mv$mdata
	lm.mv <- lm(mv_data$Y ~ mv_data$Tr, data = mv_data)
	summary(lm.mv)

	coef(lm.mv)[2]

	confint(lm.mv, level = 0.95)

	#ANSWER = PEstimate = 0.36091
	#ANSWER = (Intercept) 0.34317928 0.3786474
	#mv_data$Tr 0.05975124 0.1099108
	#p-value: 4.734e-11


	# Genetic Matching with multivariate###
	Y4 <- judge.means$lib_vote_share
	Tr4 <- judge.means$girls_A
	X_all <- cbind(judge.means$child,
	judge.means$republican == 1,
	judge.means$republican == 0,
	judge.means$woman == 1,
	judge.means$woman == 0)


	gout.mv <- GenMatch(Tr=Tr4, X=X_all, M=1, estimand = "ATT", pop.size = 100, max.generations = 20, wait.generations = 7, caliper=0.25)

	mout.gen.mv <- Match(Y=Y4, Tr=Tr4, X=X_all, estimand="ATT", M=1, Weight.matrix = gout.mv)

	# Finding Match Balance of multivariate matching
	mb.gen.mv <- MatchBalance(girls_A ~ child + republican + woman
	,data=judge.means, match.out=mout.gen.mv, nboots=500)

	# Estimating the effect with linear regression
	gen_mv_data <- mout.gen.mv$mdata
	lm.gen.mv <- lm(gen_mv_data$Y ~ gen_mv_data$Tr, data=gen_mv_data)
	summary(lm.gen.mv)

	coef(lm.gen.mv)[2]
	#ANSWER TREATMENT EFFECT = 0.360839

	confint(lm.gen.mv, level = 0.95)
	#ANSWER = (Intercept) 0.34315195 0.3785267
	#gen_mv_data$Tr 0.05872448 0.1087519