library(PSAgraphics)
data(lindner)
dim(lindner)
attach(lindner)
head(lindner)
table(abcix)

###############
##basic check##  are they related?
###############
chisq.test(lifepres, abcix)


###########
##MatchIt##   Week 2 lecture
###########
library(MatchIt)  ## try full matching
m2full = matchit(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc,data = lindner, method = "full")

m2full.dat = match.data(m2full)
attach(m2full.dat)

boxplot(distance ~ abcix) # propensity score overlap
# can also do overlapping histogram as in week 1
  
summary(m2full)   # check balance improvement
summary(m2full, standardize = T)  # look at standardized mean diffs < .1
plot(summary(m2full, standardize = T))  # see plot, balance improvement

head(m2full.dat)

#####################
##do some inference##
#####################
library(lme4) #compare outcomes over subclasses, log(cardbill) outcome
mfullL.lmer = lmer(log(cardbill) ~ abcix + (1 + abcix|subclass), data = m2full.dat)
summary(mfullL.lmer)
confint(mfullL.lmer)
exp(confint(mfullL.lmer))


############
##optmatch##
############
library(optmatch)
d <- lindner
treat_var <- "abcix"
covars <- c("stent","height","female","diabetic","acutemi","ejecfrac","ves1proc")

opt_1 <- fullmatch(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc,data = d)
summary(opt_1)

###pulling out matches###
sets <- unique(opt_1)


#Use this block to get set sizes
set_sizes <- matrix(NA,nrow=length(sets),ncol=1) 
for(i in 1:length(sets)){ 
  set_sizes[i,] <- table(opt_1)[which(names(table(opt_1))==sets[i])]
}
rownames(set_sizes) <- sets
set_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  set_sizes_column <- ifelse(opt_1==sets[i], set_sizes[i], set_sizes_column)
}


#Use this block to get number of treated units in each set
treat_sizes <- matrix(NA,nrow=length(sets),ncol=1)  
for(i in 1:length(sets)){ 
  treat_sizes[i,] <- sum((opt_1==sets[i])*(d[,treat_var]==1))
}
rownames(treat_sizes) <- sets
treat_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  treat_sizes_column <- ifelse(opt_1==sets[i], treat_sizes[i], treat_sizes_column)
}
  

#Create a dataframe that stores averaged values within set, within treatment type
d_sets <- matrix(NA,nrow=2*length(sets),ncol=length(covars)+4)
colnames(d_sets) <- c(treat_var,covars,"set_size","num_treat","weights")
for(i in 1:length(sets)){
  treated_temp <- d[as.logical((d[,treat_var]==1)*(opt_1==sets[i])),covars]
  control_temp <- d[as.logical((d[,treat_var]==0)*(opt_1==sets[i])),covars]
  d_sets[i             ,covars] <- apply(control_temp,2,mean)
  d_sets[i+length(sets),covars] <- apply(treated_temp,2,mean)
    
}
d_sets[,treat_var] <- c(rep(1,length(sets)),rep(0,length(sets)))
d_sets[,"set_size"]  <- rep(set_sizes  ,2)
d_sets[,"num_treat"] <- rep(treat_sizes,2)


###Choosing whether to weight to ATE, ATT, or ATC###

#for ATE
#d_sets[,"weights"] <- d_sets[,"set_size"]/sum(set_sizes)  

#for ATT
d_sets[,"weights"] <- d_sets[,"num_treat"]/sum(treat_sizes)

#for ATC
#d_sets[,"weights"] <- (d_sets[,"set_size"]-d_sets[,"num_treat"])/(sum(set_sizes)-sum(treat_sizes))


##Create Table 1
weighted.means.apply <- function(x){weighted.mean(x,d_sets[1:length(sets),"weights"])}
out <- matrix(NA,nrow=(length(covars)+1),ncol=4)
out[,1] <- c(treat_var,covars)
out[,2] <- apply(d_sets[1:length(sets)                   ,c(treat_var,covars)],2,weighted.means.apply)
out[,3] <- apply(d_sets[(length(sets)+1):(2*length(sets)),c(treat_var,covars)],2,weighted.means.apply)
sd.all  <- apply(d[,c(treat_var,covars)],2,sd)
out[,4] <- (as.numeric(out[,2])-as.numeric(out[,3]))/sd.all
colnames(out) <- c("variables","mean treated","mean control","smd")

write.csv(out,file="table1_att.csv",row.names=FALSE)









############
##optmatch##    adding a caliper
############
ppty <- glm(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc
            , family = binomial(),
            data = d)
ppty.dist <- match_on(ppty)
mhd.pptyc <- caliper(ppty.dist, width = 1) +
  match_on(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc,data = d)


opt_2 <- fullmatch(mhd.pptyc, data = d)
summary(opt_2)

###pulling out matches###
sets <- unique(opt_2)

#remove units matched to sinks
d_matched <- d[!is.na(opt_2),]
sets <- unique(opt_2[!is.na(opt_2)])
opt_2 <- opt_2[!is.na(opt_2)]

#Use this block to get set sizes
set_sizes <- matrix(NA,nrow=length(sets),ncol=1) 
for(i in 1:length(sets)){ 
  set_sizes[i,] <- table(opt_2)[which(names(table(opt_2))==sets[i])]
}
rownames(set_sizes) <- sets
set_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  set_sizes_column <- ifelse(opt_2==sets[i], set_sizes[i], set_sizes_column)
}


#Use this block to get number of treated units in each set
treat_sizes <- matrix(NA,nrow=length(sets),ncol=1)  
for(i in 1:length(sets)){ 
  treat_sizes[i,] <- sum((opt_2==sets[i])*(d_matched[,treat_var]==1))
}
rownames(treat_sizes) <- sets
treat_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  treat_sizes_column <- ifelse(opt_2==sets[i], treat_sizes[i], treat_sizes_column)
}


#Create a dataframe that stores averaged values within set, within treatment type
d_sets <- matrix(NA,nrow=2*length(sets),ncol=length(covars)+4)
colnames(d_sets) <- c(treat_var,covars,"set_size","num_treat","weights")
for(i in 1:length(sets)){
  treated_temp <- d_matched[as.logical((d_matched[,treat_var]==1)*(opt_2==sets[i])),covars]
  control_temp <- d_matched[as.logical((d_matched[,treat_var]==0)*(opt_2==sets[i])),covars]
  d_sets[i             ,covars] <- apply(control_temp,2,mean)
  d_sets[i+length(sets),covars] <- apply(treated_temp,2,mean)
  
}
d_sets[,treat_var] <- c(rep(1,length(sets)),rep(0,length(sets)))
d_sets[,"set_size"]  <- rep(set_sizes  ,2)
d_sets[,"num_treat"] <- rep(treat_sizes,2)


###Choosing whether to weight to ATE, ATT, or ATC###

#for ATE
d_sets[,"weights"] <- d_sets[,"set_size"]/sum(set_sizes)  

#for ATT
#d_sets[,"weights"] <- d_sets[,"num_treat"]/sum(treat_sizes)

#for ATC
#d_sets[,"weights"] <- (d_sets[,"set_size"]-d_sets[,"num_treat"])/(sum(set_sizes)-sum(treat_sizes))


##Create Table 1
weighted.means.apply <- function(x){weighted.mean(x,d_sets[1:length(sets),"weights"])}
out <- matrix(NA,nrow=(length(covars)+1),ncol=4)
out[,1] <- c(treat_var,covars)
out[,2] <- apply(d_sets[1:length(sets)                   ,c(treat_var,covars)],2,weighted.means.apply)
out[,3] <- apply(d_sets[(length(sets)+1):(2*length(sets)),c(treat_var,covars)],2,weighted.means.apply)
sd.all  <- apply(d[,c(treat_var,covars)],2,sd)
out[,4] <- (as.numeric(out[,2])-as.numeric(out[,3]))/sd.all
colnames(out) <- c("variables","mean treated","mean control","smd")

write.csv(out,file="table1_ate_caliper.csv",row.names=FALSE)




############
##optmatch##    constraining the match sizes
############
ppty <- glm(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc
            , family = binomial(),
            data = d)
ppty.dist <- match_on(ppty)
mhd.pptyc <- caliper(ppty.dist, width = 1) +
  match_on(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc,data = d)


opt_3 <- fullmatch(mhd.pptyc, min.controls = 0, max.controls = 3, data = d)
opt_3 <- pairmatch(mhd.pptyc, data = d)
summary(opt_3)

###pulling out matches###
sets <- unique(opt_3)

#remove units matched to sinks
d_matched <- d[!is.na(opt_3),]
sets <- unique(opt_3[!is.na(opt_3)])
opt_3 <- opt_3[!is.na(opt_3)]

#Use this block to get set sizes
set_sizes <- matrix(NA,nrow=length(sets),ncol=1) 
for(i in 1:length(sets)){ 
  set_sizes[i,] <- table(opt_3)[which(names(table(opt_3))==sets[i])]
}
rownames(set_sizes) <- sets
set_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  set_sizes_column <- ifelse(opt_3==sets[i], set_sizes[i], set_sizes_column)
}


#Use this block to get number of treated units in each set
treat_sizes <- matrix(NA,nrow=length(sets),ncol=1)  
for(i in 1:length(sets)){ 
  treat_sizes[i,] <- sum((opt_3==sets[i])*(d_matched[,treat_var]==1))
}
rownames(treat_sizes) <- sets
treat_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  treat_sizes_column <- ifelse(opt_3==sets[i], treat_sizes[i], treat_sizes_column)
}


#Create a dataframe that stores averaged values within set, within treatment type
d_sets <- matrix(NA,nrow=2*length(sets),ncol=length(covars)+4)
colnames(d_sets) <- c(treat_var,covars,"set_size","num_treat","weights")
for(i in 1:length(sets)){
  treated_temp <- d_matched[as.logical((d_matched[,treat_var]==1)*(opt_3==sets[i])),covars]
  control_temp <- d_matched[as.logical((d_matched[,treat_var]==0)*(opt_3==sets[i])),covars]
  d_sets[i             ,covars] <- apply(control_temp,2,mean)
  d_sets[i+length(sets),covars] <- apply(treated_temp,2,mean)
  
}
d_sets[,treat_var] <- c(rep(1,length(sets)),rep(0,length(sets)))
d_sets[,"set_size"]  <- rep(set_sizes  ,2)
d_sets[,"num_treat"] <- rep(treat_sizes,2)


###Choosing whether to weight to ATE, ATT, or ATC###

#for ATE
d_sets[,"weights"] <- d_sets[,"set_size"]/sum(set_sizes)  

#for ATT
#d_sets[,"weights"] <- d_sets[,"num_treat"]/sum(treat_sizes)

#for ATC
#d_sets[,"weights"] <- (d_sets[,"set_size"]-d_sets[,"num_treat"])/(sum(set_sizes)-sum(treat_sizes))


##Create Table 1
weighted.means.apply <- function(x){weighted.mean(x,d_sets[1:length(sets),"weights"])}
out <- matrix(NA,nrow=(length(covars)+1),ncol=4)
out[,1] <- c(treat_var,covars)
out[,2] <- apply(d_sets[1:length(sets)                   ,c(treat_var,covars)],2,weighted.means.apply)
out[,3] <- apply(d_sets[(length(sets)+1):(2*length(sets)),c(treat_var,covars)],2,weighted.means.apply)
sd.all  <- apply(d[,c(treat_var,covars)],2,sd)
out[,4] <- (as.numeric(out[,2])-as.numeric(out[,3]))/sd.all
colnames(out) <- c("variables","mean treated","mean control","smd")

write.csv(out,file="table1_ate_caliper.csv",row.names=FALSE)





############
##optmatch##    what if there are more "treated" than "controls"?
############
#d <- lindner

d[,treat_var] <- 1 - d[,treat_var]

ppty <- glm(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc
            , family = binomial(),
            data = d)
ppty.dist <- match_on(ppty)
mhd.pptyc <- caliper(ppty.dist, width = 1) +
  match_on(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc,data = d)


opt_4 <- fullmatch(mhd.pptyc, min.controls = 0, max.controls = 3, data = d)
opt_4 <- pairmatch(mhd.pptyc, data = d)
summary(opt_4)

###pulling out matches###
sets <- unique(opt_4)

#remove units matched to sinks
d_matched <- d[!is.na(opt_4),]
sets <- unique(opt_4[!is.na(opt_4)])
opt_4 <- opt_4[!is.na(opt_4)]

#Use this block to get set sizes
set_sizes <- matrix(NA,nrow=length(sets),ncol=1) 
for(i in 1:length(sets)){ 
  set_sizes[i,] <- table(opt_4)[which(names(table(opt_4))==sets[i])]
}
rownames(set_sizes) <- sets
set_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  set_sizes_column <- ifelse(opt_4==sets[i], set_sizes[i], set_sizes_column)
}


#Use this block to get number of treated units in each set
treat_sizes <- matrix(NA,nrow=length(sets),ncol=1)  
for(i in 1:length(sets)){ 
  treat_sizes[i,] <- sum((opt_4==sets[i])*(d_matched[,treat_var]==1))
}
rownames(treat_sizes) <- sets
treat_sizes_column <- matrix(NA,nrow=dim(d)[1],ncol=1)
for(i in 1:length(sets)){
  treat_sizes_column <- ifelse(opt_4==sets[i], treat_sizes[i], treat_sizes_column)
}


#Create a dataframe that stores averaged values within set, within treatment type
d_sets <- matrix(NA,nrow=2*length(sets),ncol=length(covars)+4)
colnames(d_sets) <- c(treat_var,covars,"set_size","num_treat","weights")
for(i in 1:length(sets)){
  treated_temp <- d_matched[as.logical((d_matched[,treat_var]==1)*(opt_4==sets[i])),covars]
  control_temp <- d_matched[as.logical((d_matched[,treat_var]==0)*(opt_4==sets[i])),covars]
  d_sets[i             ,covars] <- apply(control_temp,2,mean)
  d_sets[i+length(sets),covars] <- apply(treated_temp,2,mean)
  
}
d_sets[,treat_var] <- c(rep(1,length(sets)),rep(0,length(sets)))
d_sets[,"set_size"]  <- rep(set_sizes  ,2)
d_sets[,"num_treat"] <- rep(treat_sizes,2)


###Choosing whether to weight to ATE, ATT, or ATC###

#for ATE
d_sets[,"weights"] <- d_sets[,"set_size"]/sum(set_sizes)  

#for ATT
#d_sets[,"weights"] <- d_sets[,"num_treat"]/sum(treat_sizes)

#for ATC
#d_sets[,"weights"] <- (d_sets[,"set_size"]-d_sets[,"num_treat"])/(sum(set_sizes)-sum(treat_sizes))


##Create Table 1
weighted.means.apply <- function(x){weighted.mean(x,d_sets[1:length(sets),"weights"])}
out <- matrix(NA,nrow=(length(covars)+1),ncol=4)
out[,1] <- c(treat_var,covars)
out[,2] <- apply(d_sets[1:length(sets)                   ,c(treat_var,covars)],2,weighted.means.apply)
out[,3] <- apply(d_sets[(length(sets)+1):(2*length(sets)),c(treat_var,covars)],2,weighted.means.apply)
sd.all  <- apply(d[,c(treat_var,covars)],2,sd)
out[,4] <- (as.numeric(out[,2])-as.numeric(out[,3]))/sd.all
colnames(out) <- c("variables","mean treated","mean control","smd")

write.csv(out,file="table1_ate_caliper.csv",row.names=FALSE)