gmatchv1.R

# Load imports
library(rbounds)
library(rgenoud)

# Load data
data(lalonde)
attach(lalonde)

# Create data objects
Y <- lalonde$re78
Tr <- lalonde$treat
X <- cbind(age, educ, black, hisp, married, nodegr, re74, re75, u74, u75)
BalanceMat <- cbind(age, I(age^2), educ, I(educ^2), black,
                    hisp, married, nodegr, re74 , I(re74^2), re75, I(re75^2),
                    u74, u75, I(re74*re75), I(age*nodegr), I(educ*re74), I(educ*75))

# Genetic matching 
gen1 <- GenMatch(Tr=Tr, X=X, BalanceMat=BalanceMat, pop.size=50,
                 data.type.int=FALSE, print=1, replace=FALSE)

# Matching
mgen1 <- Match(Y=Y, Tr=Tr, X=X, Weight.matrix=gen1, replace=FALSE)
summary(mgen1)

# Sensitivity Analysis
p <- psens(mgen1, Gamma=1.5, GammaInc=.1)
p
1 - p$bounds$`Upper bound`
hlsens(mgen1, Gamma=1.5, GammaInc=.1, .1)

# Sensitivity optim function
my.fitfunc <- function(matches, BM) {
  
  # my.fitfunc requires that the last column of the BM matrix contain the outcomes
  # hence, it is VERY important not to run standard pval balance genetic matching
  # with the same BM!
  
  # NOTE: psens does NOT acount for 1:many matching!
  index.treated <- matches[,1]
  index.control <- matches[,2]
  
  # Get treated and control unit outcomes
  trt <-  BM[, ncol(BM)][index.treated]
  ctrl <- BM[, ncol(BM)][index.control]
  
  # SECOND OPTION: Maximize gamma associated with
  # losing statistical significance
  
  # Compute sensitivity
  psens.out <- psens(trt, ctrl, Gamma=1.5, GammaInc=.1)
  
  # We want to minimize the largest p-value. Thus,
  pvals <- psens.out$bounds$`Upper bound`
  
  return(sort(pvals, decreasing=TRUE))
}

# Demonstration
BM <- cbind(BalanceMat, re78)
df <- data.frame(trt = mgen1$index.treated, ctrl = mgen1$index.control)
my.fitfunc(df, BM)           # With matches from Match()
my.fitfunc(gen1$matches, BM) # Withc matches from GenMatch()

# Genetic Optimization
genout <- GenMatch(Tr=Tr, X=X, BalanceMat=BM, pop.size=505,
                   print=1, ties=FALSE, wait.generations = 30, fit.func = my.fitfunc)

mout <- Match(Y=Y, Tr=treat, X=X, ties=FALSE,
              estimand="ATT", Weight.matrix=genout)
summary(mout)

# INVESTIGATE WHETHER ties=false influence whether genout and mout influence pvals

psens(mout, Gamma=1.5, GammaInc=0.1)
hlsens(mout, Gamma=1.5, GammaInc=0.1)

#Let's determine if balance has actually been obtained on the variables of interest
mb <- MatchBalance(treat~age +educ+black+ hisp+ married+
                     nodegr + u74 + u75 + re75 + re74 + I(re74 * re75),
                   match.out = mout, nboots = 500, ks = TRUE)