Casual Inference Assignment 4

library(boot)
library(caret)
PEACEKEEPING WORKOUT (based on  King, Gary;Zeng, Langche, 2007, 
                      "Replication data for: When Can History be Our Guide? 
                      The Pitfalls of Counterfactual Inference", 
                      https://hdl.handle.net/1902.1/DXRXCFAWPK, 
                      Harvard Dataverse, V4, 
                      UNF:3:DaYlT6QSX9r0D50ye+tXpA== [fileUNF] )

peace <- read.csv("https://course-resources.minerva.kgi.edu/uploaded_files/mke/00086677-3767/peace.csv")
# extract relevant columns
peace <- peace[, 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)
peace <- peace[c(-19, -47), ]
attach(peace)

#variables based on the table 2
glm1 <- glm(pbs2s3~ wartype + logcost + wardur + factnum +
              factnum2 + trnsfcap + develop + exp + decade +
              treaty + untype4, data=peace, family = "binomial")

#empty data frames for original model
control <- c()
treat <- c()
store1 <- c()

#mean of the values that we are considering for the model
#based on "we compute the marginal
#effects of UN peacekeeping operations as a function of the duration of the civil war,
#holding constant all other variables at their means"

typem <- mean(wartype)
costm <- mean(logcost)
factm <-mean(factnum)
fact2m <- mean(factnum2)
trnsm <- mean(trnsfcap)
develm <- mean(develop)
expm <- mean(exp)
decm <- mean(decade)
treatm <- mean(treaty)

#defining two observations (one control and one treatment) with mean values of the variables to predict to different models


#we loop for each 315 which is the number of month for treatment effect according to figure 8

for(i in 1:315){ 
  control <- data.frame(typem,costm,i,factm,fact2m,trnsm,0,treatm,
                     develm,expm,decm)
  treat <- data.frame(typem,costm,i,factm,fact2m,trnsm,1,treatm,
                      develm,expm,decm)
  names(control) <- c('wartype','logcost','wardur','factnum','factnum2','trnsfcap','untype4','treaty','develop',
                   'exp','decade')
  names(treat) <- c('wartype','logcost','wardur','factnum','factnum2','trnsfcap','untype4','treaty','develop',
                    'exp','decade')
  pre_treat <- mean(predict(glm1, newdata = treat,type = 'response'))
  pre_control <- mean(predict(glm1, newdata = control, type = 'response'))
  store1[i] <- pre_treat - pre_control
}



#model with the interaction term

glm2 <- glm(pbs2s3 ~ wartype + logcost + wardur + factnum +
              factnum2 + trnsfcap + develop + exp + decade +
              treaty + untype4 + wardur*untype4 , data=peace, family = "binomial")

#loop for each month for treatment effect with interaction term
store2 <- c()
for(i in 1:315){
  control <- data.frame(typem,costm,i,factm,fact2m,trnsm,0,treatm,
                     develm,expm,decm)
  treat <- data.frame(typem,costm,i,factm,fact2m,trnsm,1,treatm,
                      develm,expm,decm)
  names(control) <- c('wartype','logcost','wardur','factnum','factnum2','trnsfcap','untype4','treaty','develop',
                   'exp','decade')
  names(treat) <- c('wartype','logcost','wardur','factnum','factnum2','trnsfcap','untype4','treaty','develop',
                    'exp','decade')
  pre_treat <- mean(predict(glm2, newdata = treat,type = 'response'))
  pre_control <- mean(predict(glm2, newdata = control, type = 'response'))
  store2[i] <- pre_treat - pre_control
}

#plottting the store 1
plot(store1, type = "l", lty = "dotted",ylim = c(0,0.8),xlim = c(0,315),axes=FALSE,
     xlab = "",ylab="",
     sub = "FIG. 8. Causal Effect of Multidimensional UN Peacekeeping Operations", font=2)
par(new=TRUE)
#plottting the store 2
plot(store2,type = "l", ylim = c(0,0.8),axes = FALSE,
     ylab = "Marginal effects of UN peacekeeping operations",
     xlab="Duration of wars in months",cex.lab=1)
#defining the axis and labels    
axis(1, seq(5,315,15), cex.axis=0.5, font=2)
axis(2, seq(0.0,0.8,0.1), cex.axis=0.5, font=2)
axis(3, seq(5,315,15), cex.axis=0.5, font=2)
axis(4, seq(0.0,0.8,0.1), cex.axis=0.5, font=2)
box()
text(75,0.2,"Model with Interaction Term",cex=.8)
text(205,0.45,"Dotted:Original Model",cex=.8)


#Question 4:

library(Matching)
peace <- read.csv("https://course-resources.minerva.kgi.edu/uploaded_files/mke/00086677-3767/peace.csv",stringsAsFactors = FALSE)
# remove rows with missing data (there are better ways to handle missing data)
peace <- peace[c(-19,-47,-4,-16,-84,-93,-98), ]


#reassigning the uncint pbs2l, pbs5l values based on the codebook to define Tr accoring to Q3
peace[peace$uncint == "None", "uncint"] <- 0
peace[peace$uncint == "Observer", "uncint"] <-  2 
peace[peace$uncint == "PKO", "uncint"] <-  3 
peace[peace$uncint == "Enforcement", "uncint"] <- 4 
peace[peace$pbs2l == "Failure", "pbs2l"] <-  0
peace[peace$pbs2l == "Success", "pbs2l"] <- 1 
peace[peace$pbs5l == "Failure", "pbs5l"] <-  0
peace[peace$pbs5l == "Success", "pbs5l"] <-  1

#treatment
Tr <- rep(0, length(peace$uncint))
Tr[which(as.numeric(peace$uncint) != 0 & as.numeric(peace$uncint) != 1)] <- 1

#logistic regression

#year 2 
#Treatment effects without matching
glm1 <- glm(as.numeric(peace$pbs2l) ~ Tr + wartype + logcost + wardur + factnum + factnum2 + 
              trnsfcap + treaty + develop + exp + decade, family = binomial, data = peace)
#matchbalance
mb_y2_glm <- MatchBalance(as.numeric(peace$pbs2l) ~ Tr + wartype + logcost + wardur + factnum + factnum2 + 
                            trnsfcap + treaty + develop + exp + decade, data=peace, nboots=500)
#Tmt effect (Tr estimate)
summary(glm1)

#year 5 
glm2 <- glm(as.numeric(peace$pbs5l) ~ Tr + wartype + logcost + wardur + factnum + factnum2 + 
              trnsfcap + treaty + develop + exp + decade, family = binomial, data = peace)
#matchbalance
mb_y5_glm <- MatchBalance(as.numeric(peace$pbs5l) ~ Tr + wartype + logcost + wardur + factnum + factnum2 + 
                            trnsfcap + treaty + develop + exp + decade, data=peace, nboots=500)
#Tmt effect (Tr estimate)
summary(glm2)

#Propensity score matching 

#model for propensity score
PS.model<- glm(Tr ~ wartype + logcost + wardur + 
                        factnum + factnum2 + trnsfcap + treaty + develop +
                        exp + decade , family = binomial, data = peace)
mout_prop<- Match(Tr=Tr, X=PS.model$fitted)
mb1  <- MatchBalance(Tr ~  wartype + logcost + wardur + factnum + 
                      factnum2 + trnsfcap + 
                      treaty + develop + exp + decade, 
                    data=peace, match.out = mout_prop, nboots=500)


#year 2
mout_propy_2 <- Match(Y=as.numeric(peace$pbs2l), Tr=Tr, X=PS.model$fitted, BiasAdjust = TRUE)
summary(mout_propy_2)
#Tmt effect
signif(mout_propy_2$est)
signif(mout_propy_2$est.noadj)
#matchbalance
mb2  <- MatchBalance(Tr ~  wartype + logcost + wardur + factnum + 
                      factnum2 + trnsfcap + treaty + develop + exp + decade, data=peace, 
                    match.out = mout_propy_2, nboots=500)

#year 5
mout_propy_5 <- Match(Y=as.numeric(peace$pbs5l), Tr=Tr, X=PS.model$fitted, BiasAdjust = TRUE)
summary(mout_propy_5)
#Tmt effect
signif(mout_propy_5$est)
signif(mout_propy_5$est.noadj)
#Matchbalance
mb3  <- MatchBalance(Tr ~  wartype + logcost + wardur + factnum + 
                       factnum2 + trnsfcap + treaty + develop + exp + decade, data=peace, 
                     match.out = mout_propy_5, nboots=500)

#GenMatch

X = cbind(peace$wartype, peace$logcost, peace$wardur, 
          peace$factnum2, peace$trnsfcap, peace$treaty, peace$develop, peace$exp, peace$decade)

genout = GenMatch(Tr = Tr, X = X, estimand = "ATT", M = 1, pop.size = 200, max.generations = 25,
                  wait.generations = 7, nboots = 200)

#2 year
gen_mouty2 <- Match(Y=as.numeric(peace$pbs2l), Tr=Tr, X=X, Weight.matrix = genout, BiasAdjust = TRUE)
summary(gen_mouty2)
#tmt effect
signif(gen_mouty2$est)
signif(gen_mouty2$est.noadj)
#Matchbalance
mb_y2  <- MatchBalance(Tr ~  wartype + logcost + wardur + factnum + 
                      factnum2 + trnsfcap + treaty + develop + exp + decade, data=peace, 
                    match.out = gen_mouty2, nboots=500)

#5 years
gen_mouty5 <- Match(Y=as.numeric(peace$pbs5l), Tr=Tr, X=X, Weight.matrix = genout, BiasAdjust = TRUE)
summary(gen_mouty5)
#tmt effect
signif(gen_mouty5$est)
signif(gen_mouty5$est.noadj)
#Matchbalance
mb_y5 <- MatchBalance(Tr ~  wartype + logcost + wardur + factnum + 
                      factnum2 + trnsfcap + treaty + develop + exp + decade, data=peace, 
                    match.out = gen_mouty5, nboots=500)

#table: 

tab<-matrix(c("NA*",0.5996839,
  "NA**(0.0020)", "NA*", 0.8233143, "NA**(0.0027)",0.363588
  ,0.363636
  ,"NA**(0.006)",
  0.393891,
  0.393939,
  "NA**(0.01)",0.116015, "NA**(0.0606061)",
  0.182, 
  0.146318,
  "NA**(0.0909091)",
  0.204
),ncol=3,byrow=TRUE)
 rownames(tab)<-c(
   "Logistic regression/len success 2 years",
   "Logistic regression/len success 5 years",
   
   
   "P- score matching/len success 2 years",
   "P- score matching/len success 5 years",
   
   
   "GenMatch/len success 2 years",
   "GenMatch/len success 5 years"
 )
colnames(tab)<-c("Tmt effect (bias adj)","Tmt effect (no bias adj)",
  "P-value (from MatchBalance)"
)
 tab <- as.table(tab)
 tab