#############################################################################
##                                                                         ## 
##    Chapter 2: Covariate Selection and Propensity Score Estimation       ##
##                                                                         ##
#############################################################################


############################################################################
#                 Logit Model to Estimate Propensity Scores                #
############################################################################

# Import the psdemo.csv dataset to R
psdemo <- read.csv("D:/Sage Little Green Book/website materials/psdemo.csv") 

# Run the logitic regression model to predic propensity scores
logitMod <- glm(s_treatment ~ student_gender+ s_grade+ S_CLIMATE_COMMUNITY+ S_CLIMATE_SCHOOLSAFETY+S_CONFLICTRES_AGGRESSIVE+ S_LEARNING_RECESSEFFECT+ S_LEARNING_ENGAGEMENT+ S_RECESS_ORGANIZED+ S_YOUTHDEV_INTERACTIONS+ S_YOUTHDEV_PEERNONCONFLICT + S_PHYSICAL_SELFCONCEPT, data=psdemo, family=binomial(link="logit"))

# predicted scores; 
predicted <- plogis(predict(logitMod, psdemo))

# or for predicted scores
predicted <- predict(logitMod, psdemo, type="response")  

####################################################################################
#                 Classification Tree for Propensity Scores                        #
####################################################################################

#Coming soon!

####################################################################################
#                 Classification Tree for Propensity Scores                        #
####################################################################################

#Coming soon!




#############################################################################
##                                                                         ## 
##              Chapter 3: Propensity Score Adjustment Methods             ##
##                                                                         ##
#############################################################################

############################################################
#                 Data Preparation                         #
############################################################

# Install MatchIt and related packages (do it once)
install.packages("MatchIt")
install.packages("optmatch")
library(MASS)
library(MatchIt)
library(optmatch)

# Import the psdemo.csv dataset into R

# For PC
psdemo <- read.csv("C:/Users/.../psdemo.csv")

#For Mac
psdemo <- read.csv("/Users/.../psdemo.csv")

# View psdemo dataset
psdemo

#############################################################
#               Analyses before Matching                    #
#############################################################

# T-test on the matched data
t.test(S_CLIMATE_RECESSSAFETY ~ s_treatment, data = psdemo)

# Regression on the matched s_treatmentdata
mrl <- lm(S_CLIMATE_RECESSSAFETY ~ s_treatment+ S_CONFLICTRES_RELATIONSHIPS+ S_CONFLICTRES_AGGBELIEF+ S_LEARNING_SPORTSEFFECT+ S_RECESS_ENJOYMENT+ S_YOUTHDEV_PEERCONFLICT, data = psdemo)
summary(mrl)


############################################################
#                    PSM Matching                          #
############################################################

# Run optimal matching
om.out <- matchit(s_treatment ~ S_CONFLICTRES_RELATIONSHIPS+ S_CONFLICTRES_AGGBELIEF+ S_LEARNING_SPORTSEFFECT+ S_RECESS_ENJOYMENT+ S_YOUTHDEV_PEERCONFLICT, data = psdemo, method = "optimal", distance = "logit", ratio = 1)




#############################################################################
##                                                                         ## 
##              Chapter 4: Propensity Score Adjustment Methods             ##
##                                                                         ##
#############################################################################


##########################################################################################################
#            #4.5a Example for checking the balance of covariates before and after PS matching           # 
##########################################################################################################

# Checking Balance 
summary(om.out)
plot(om.out, type="hist")
plot(summary(om.out, standardize = T), interactive = T)


############################################################
#               Analyses after Matching                    #
############################################################

# Write the matched data in a dataframe for after matching analysis
mpsdemo = match.data(om.out) 

# T-test on the matched data
t.test(S_CLIMATE_RECESSSAFETY ~ s_treatment, data = mpsdemo)

# Regression on the matched data
mr2 <- lm(S_CLIMATE_RECESSSAFETY ~ s_treatment+ S_CONFLICTRES_RELATIONSHIPS+ S_CONFLICTRES_AGGBELIEF+ S_LEARNING_SPORTSEFFECT+ S_RECESS_ENJOYMENT+ S_YOUTHDEV_PEERCONFLICT, data = mpsdemo)
summary(mr2)

############################################################
#                 Sensitivity analysis                     #
############################################################

# Installation of packages to prepare for sensitivity analysis

install.packages("Matching", dependencies = T)
install.packages("rgenoud")
library("Matching")
library(rbounds)

# Estimate propensity score
ps <- glm(s_treatment ~ S_CONFLICTRES_RELATIONSHIPS+ S_CONFLICTRES_AGGBELIEF+ S_LEARNING_SPORTSEFFECT+ S_RECESS_ENJOYMENT+ S_YOUTHDEV_PEERCONFLICT,
          family = binomial, data = psdemo)

## ATT using independent sample t-test on matched data##
t.test(S_CLIMATE_RECESSSAFETY ~ s_treatment, data = mpsdemo)

## ATT using independent sample t-test on matched original data for comparison with the results from matched data##
t.test(S_CLIMATE_RECESSSAFETY ~ s_treatment, data = psdemo)

##Regression adjustment##
mr2 <- lm(S_CLIMATE_RECESSSAFETY ~ s_treatment+ S_CONFLICTRES_RELATIONSHIPS+ S_CONFLICTRES_AGGBELIEF+ S_LEARNING_SPORTSEFFECT+ S_RECESS_ENJOYMENT+ S_YOUTHDEV_PEERCONFLICT, data = mpsdemo)
summary(mr2)

##Stratification using MatchIt##
stra.out <- matchit(s_treatment ~ student_gender+ s_grade+ S_CLIMATE_COMMUNITY+ S_CLIMATE_SCHOOLSAFETY+S_CONFLICTRES_AGGRESSIVE+ S_LEARNING_RECESSEFFECT+ S_LEARNING_ENGAGEMENT+ S_RECESS_ORGANIZED+ S_YOUTHDEV_INTERACTIONS+ S_YOUTHDEV_PEERNONCONFLICT + S_PHYSICAL_SELFCONCEPT, data = psdemo, method = "subclass")

strpsdemos = match.data(stra.out)

##ANOVA test on treatment X strata on outcome after sub-classification ##
res.aovstr <- aov(S_CLIMATE_RECESSSAFETY ~ s_treatment + subclass, data = strpsdemos)
summary(res.aovstr)

###############Using PS Weight for ATE in R######################

pswt_demo <- read.csv("D:/Sage Little Green Book/website materials/psdemo.csv") 
pswt_demo
attach(pswt_demo)
pswt_demo$PS_ATE <- (s_treatment/PS) + ((1- s_treatment)/(1- PS))
detach(pswt_demo)
summary(lm(S_CLIMATE_RECESSSAFETY~s_treatment, data = pswt_demo, weights = PS_ATE))

#################Using PS Weight for ATT in R###################

pswt_demo <- read.csv("D:/Sage Little Green Book/website materials/psdemo.csv") 
pswt_demo
attach(pswt_demo)
pswt_demo$PS_ATT <- s_treatment + ((1- s_treatment)*PS)/(1- PS)
detach(pswt_demo)
summary(lm(S_CLIMATE_RECESSSAFETY ~s_treatment, data = pswt_demo, weights = PS_ATT))


############################################################
#                 Sensitivity analysis                     #
############################################################

# Installation of packages to prepare for sensitivity analysis

install.packages("Matching", dependencies = T)
install.packages("rgenoud")
library("Matching")
library(rbounds)

T <- psdemo$s_treatment
ps <- glm(s_treatment ~ student_gender+ s_grade+ S_CLIMATE_COMMUNITY+ S_CLIMATE_SCHOOLSAFETY+S_CONFLICTRES_AGGRESSIVE+ S_LEARNING_RECESSEFFECT+ S_LEARNING_ENGAGEMENT+ S_RECESS_ORGANIZED+ S_YOUTHDEV_INTERACTIONS+ S_YOUTHDEV_PEERNONCONFLICT + S_PHYSICAL_SELFCONCEPT, family = binomial, data=psdemo)
O <- psdemo$S_CLIMATE_RECESSSAFETY
m <- Match(Y = O, Tr = T, X = ps$fitted, replace = F)

##Rosenbaum Sensitivity Test for Wilcoxon Signed Rank P-Value## 
psens(m, Gamma = 2, GammaInc = 0.1) 

##Rosenbaum Sensitivity Test for Hodges-Lehmann Point Estimate to ATT## 
hlsens(m, Gamma = 2, GammaInc = 0.1, 0.1)