Created
November 23, 2019 13:50
-
-
Save t-flora/18aa2eaade5270848a757ad0c4a21382 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
################################## | |
################################## | |
######## CAUSAL INFERENCE ######## | |
################################## | |
################################## | |
#Setting up libraries | |
library(rbounds) | |
library(Matching) | |
library(arm) | |
set.seed(8008) | |
################# | |
### Problem 1 ### | |
################# | |
# Part A. | |
# The code in this part tries to find match balance before doing the matching itself. | |
# It takes the output of “genout”, which calls GenMatch and is thus a weight matrix, to evaluate match balance, even though no matching is being performed by GenMatch. | |
# Part B | |
# The arguments passed to the GenMatch function are not the same as the ones passed to Match() later on. | |
# Specifically, there is no “caliper” argument in Match, which defaults it to NULL, and the “estimand” argument would also be a problem were it not defaulted to “ATT” or if we had another estimand in GenMatch. | |
# Part C | |
# Y is defined only as “re78”, without specifying the dataset the variable is in. | |
# This means that will be a traceback error. | |
################# | |
### Problem 2 ### | |
################# | |
foo <- read.csv("https://course-resources.minerva.kgi.edu/uploaded_files/mke/00086677-3767/peace.csv") | |
foo2 <- foo[, c(6:8, 11:16, 99, 50, 114, 49, 63, 136, 109, 126, 48, 160, 142, 10)] | |
# remove 2 rows with missing data (there are better ways to handle missing data) | |
foo2 <- na.omit(foo2) | |
# take a peek at the data set (identify the columns) | |
head(foo2) | |
summary(foo2) | |
glm.orig <- glm(pbs2s3 ~ wartype + logcost + wardur + factnum + factnum2 + | |
trnsfcap + develop + exp + decade + treaty + untype4, | |
data = foo2, family = 'binomial') | |
glm.interac <- glm(pbs2s3 ~ wartype + logcost + wardur + factnum + factnum2 + | |
trnsfcap + develop + exp + decade + treaty + untype4 + | |
I(foo2$exp*foo2$untype4) + I(foo2$wardur*foo2$logcost), | |
data = foo2, family = 'binomial') | |
# Geographic variables were not taken into account in Gary King's model of interest | |
mean.wartype <- mean(foo2$wartype) | |
mean.logdead <- mean(foo2$logcost) | |
# In this case, all variables but wardur (which we plot) are held at their means | |
# mean.wardur <- mean(foo2$wardur) | |
mean.factnum <- mean(foo2$factnum) | |
mean.factnum2 <- mean(foo2$factnum2) | |
mean.trnsfcap <- mean(foo2$trnsfcap) | |
mean.develop <- mean(foo2$develop) | |
mean.exp <- mean(foo2$exp) | |
mean.decade <- mean(foo2$decade) | |
mean.treaty <- mean(foo2$treaty) | |
# Hold interaction terms at their means as well | |
mean.exptreat <- mean(foo2$exp*foo2$untype4) | |
mean.durlogcost <- mean(foo2$wardur*foo2$logcost) | |
# How many wars in the data set had peacebuilding missions? | |
# mean.untype4 <- mean(foo2$untype4) | |
# There is only one observation with 'wardur' above 315 months, at 600 months | |
# We discard the observation as it is an outlier too far from the range of the remaining data | |
length(which(foo2$wardur>315)) | |
logits <- function(Xs, coef) { | |
logit <- coef[1] + sum(coef[2:length(coef)]*Xs) | |
return(exp(logit) / (1 + exp(logit))) | |
} | |
# Create empty vectors to store war durations | |
# Each vector has a length of 315 as that is | |
# the maximum duration plotted in King (2007) | |
wardur.origtreat <- rep(NA, 315) | |
wardur.origctrl <- rep(NA, 315) | |
for (k in 1:315) { | |
Xs.orig.treat <- c(mean.wartype, mean.logdead, k, mean.factnum, mean.factnum2, | |
mean.trnsfcap, mean.develop, mean.exp, mean.decade, mean.treaty, 1) | |
Xs.orig.ctrl <- c(mean.wartype, mean.logdead, k, mean.factnum, mean.factnum2, | |
mean.trnsfcap, mean.develop, mean.exp, mean.decade, mean.treaty, 0) | |
wardur.origtreat[k] <- logits(Xs.orig.treat, coef(glm.orig)) | |
wardur.origctrl[k] <- logits(Xs.orig.ctrl, coef(glm.orig)) | |
} | |
original.y <- wardur.origtreat - wardur.origctrl | |
wardur.interac.treat <- rep(NA, 315) | |
wardur.interac.ctrl <- rep(NA, 315) | |
for (i in 1:315) { | |
Xs.interac.treat <- c(mean.wartype, mean.logdead, i, mean.factnum, mean.factnum2, | |
mean.trnsfcap, mean.develop, mean.exp, mean.decade, mean.treaty, 1, 1*mean.exp, | |
i*mean.logdead) | |
Xs.interac.ctrl <- c(mean.wartype, mean.logdead, i, mean.factnum, mean.factnum2, | |
mean.trnsfcap, mean.develop, mean.exp, mean.decade, mean.treaty, 0, 0*mean.exp, i*mean.logdead) | |
wardur.interac.treat[i] <- logits(Xs.interac.treat, coef(glm.interac)) | |
wardur.interac.ctrl[i] <- logits(Xs.interac.ctrl, coef(glm.interac)) | |
} | |
interac_y <- wardur.interac.treat - wardur.interac.ctrl | |
plot(1:315, original.y, xlab = "Duration of wars in months", ylab="Marginal effects of UN peacekeeping operations", | |
type = "l", ylim = c(0, 0.8), main = 'Causal effect of multidimensional UN peacekeeping operations') | |
lines(1:315, interac_y, col = "blue", ylim = c(0, 0.8)) | |
legend(1, 95, legend=c("Original model", "Model with interaction terms"), | |
col=c("blue", "red"), lty=1:2, cex=0.8) | |
################# | |
### Problem 3 ### | |
################# | |
## If we are to only select the variables to be used in problems 3 & 4, we can use the | |
# the following code to only pick those columns of the dataset | |
## Find 'uncint', 'pbs2l', and 'pbs5l | |
# which(foo=='Enforcement', arr.ind = T) | |
# which(colnames(foo)=='pbs2l') | |
# which(colnames(foo)=='pbs5l') | |
## Keep 'uncint', 'pbs2l', and 'pbs5l | |
foo3 <- foo[, c(6:8, 11:16, 99, 50, 34, 35, 52, 114, 49, 63, 136, 109, 126, 48, 160, 142, 10)] | |
## Remove rows with missing data | |
foo3 <- na.omit(foo3) | |
# Create a vector full of zeroes of the same length | |
# as there are data points for 'uncint' | |
Tr <- rep(0, length(foo3$uncint)) | |
Tr[which(foo3$uncint != "None")] <- 1 | |
foo3$Tr <- Tr | |
################# | |
### Problem 4 ### | |
################# | |
## Logistic regressions | |
# Check balance for logistic regressions | |
mbinit <- MatchBalance(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), | |
data = foo3, nboots=500) | |
# From mbinit, we can see that there is quite a lot of difference in logcost | |
# between the treatment and control groups, with treatment group wars having more deaths | |
# There are even more noticeable differences to keep in mind in regard to ineraction terms | |
# between the human cost and the number of factions in the conflict | |
glm2l <- glm(pbs2l ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost) + Tr, | |
data = foo3, family = 'binomial') | |
summary(glm2l) | |
glm5l <- glm(pbs5l ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost) + Tr, | |
data = foo3, family = 'binomial') | |
summary(glm5l) | |
# Below, there are calculations of treatment effect estimates by considering | |
# the unit-level difference in outcomes obtained by simulating counterfactuals | |
# This was not included in the table in the responses as it is part of | |
# last year's answer key and is not my own work, but it is worth it | |
# to compare the estimates and the coefficients | |
# Two years | |
# foo.counter_factual <- foo3 | |
# foo.counter_factual$Tr <- 1 - foo3$Tr | |
# counter.factuals <- predict(glm2l, newdata=foo.counter_factual, type="response") | |
# unit_treat_effects <- rep(NA, nrow(foo3)) | |
# mask <- foo3$Tr == 1 | |
# unit_treat_effects[mask] <- glm2l$fitted.values[mask] - counter.factuals[mask] | |
# unit_treat_effects[!mask] <- counter.factuals[!mask] - glm2l$fitted.values[!mask] | |
# mean(unit_treat_effects) | |
# Out[1] = 0.09374585 | |
# Five years | |
# foo5.counter_factual <- foo3 | |
# foo5.counter_factual$Tr <- 1 - foo3$Tr | |
# counter.factuals5 <- predict(glm5l, newdata=foo5.counter_factual, type="response") | |
# unit_treat_effects5 <- rep(NA, nrow(foo3)) | |
# mask <- foo3$Tr == 1 | |
# unit_treat_effects5[mask] <- glm5l$fitted.values[mask] - counter.factuals5[mask] | |
# unit_treat_effects5[!mask] <- counter.factuals5[!mask] - glm5l$fitted.values[!mask] | |
# mean(unit_treat_effects5) | |
# Out[2] = 0.126225 | |
### Propensity Score Matching | |
glmps <- glm(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), data = foo3, family = 'binomial') | |
summary(glmps) | |
pscores <- glmps$fitted.values | |
Y2l <- foo3$pbs2l | |
Y5l <- foo3$pbs5l | |
m2lps <- Match(Y=Y2l, Tr=Tr, X=pscores, M=1, BiasAdjust = T, replace = T) | |
summary(m2lps) | |
m2lps$est.noadj | |
mb2lps <- MatchBalance(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), data = foo3, | |
match.out = m2lps, nboots=500) | |
m5lps <- Match(Y=Y5l, Tr=Tr, X=pscores, M=1, BiasAdjust = T, replace = T) | |
summary(m5lps) | |
m5lps$est.noadj | |
mb5lps <- MatchBalance(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), data = foo3, | |
match.out = m5lps, nboots=500) | |
### Genetic Matching w/o Propensity Scores | |
attach(foo3) | |
X4 <- cbind(wartype, logcost, foo3$wardur, factnum, | |
factnum2, trnsfcap, develop, exp, decade, treaty, I(foo3$logcost*foo3$wardur), I(foo3$factnum*foo3$logcost)) | |
Xps <- cbind(glmps$fitted.values, wartype, logcost, foo3$wardur, factnum, | |
factnum2, trnsfcap, develop, exp, decade, treaty, I(foo3$logcost*foo3$wardur), I(foo3$factnum*foo3$logcost)) | |
detach(foo3) | |
genout <- GenMatch(Tr=Tr, X=X4, M=1, replace = T, pop.size=200, max.generations=10, wait.generations=25) | |
genout2 <- GenMatch(Tr=Tr, X=X4, M=1, replace = T, pop.size=200, max.generations=10, wait.generations=25) | |
match.nops2l <- Match(Y=Y2l, Tr=Tr, X=X4, M=1, BiasAdjust = T, Weight.matrix = genout, replace = T) | |
summary(match.nops2l) | |
match.nops2l$est.noadj | |
match.nops5l <- Match(Y=Y5l, Tr=Tr, X=X4, M=1, BiasAdjust = T, Weight.matrix = genout2, replace = T) | |
summary(match.nops5l) | |
match.nops5l$est.noadj | |
mbnops2l <- MatchBalance(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), data = foo3, | |
match.out = match.nops2l, nboots=500) | |
mbnops5l <- MatchBalance(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), data = foo3, | |
match.out = match.nops5l, nboots=500) | |
### Genetic Matching w/ Propensity Scores | |
genoutps <- GenMatch(Tr=Tr, X=Xps, M=1, | |
pop.size=200, max.generations=10, wait.generations=25, replace = T) | |
genoutps2 <- GenMatch(Tr=Tr, X=Xps, M=1, | |
pop.size=200, max.generations=10, wait.generations=25, replace = T) | |
gm2lps <- Match(Y=Y2l, Tr=Tr, X=Xps, M=1, BiasAdjust = T, Weight.matrix = genoutps, replace = T) | |
summary(gm2lps) | |
gm2lps$est.noadj | |
gmb2lps <- MatchBalance(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), data = foo3, | |
match.out = gm2lps, nboots=500) | |
gm5lps <- Match(Y=Y5l, Tr=Tr, X=Xps, M=1, BiasAdjust = T, Weight.matrix = genoutps2, replace = T) | |
summary(gm5lps) | |
gm5lps$est.noadj | |
gmb5lps <- MatchBalance(Tr ~ wartype + logcost + wardur + factnum + | |
trnsfcap + develop + exp + decade + treaty + | |
I(foo3$logcost*foo3$wardur) + I(foo3$factnum*foo3$logcost), data = foo3, | |
match.out = gm5lps, nboots=500) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment