1: # -----------------------------------------------------------------------------
   2: # AAAS - Designing Evaluation with Comparison Groups - Webinar - Part 2
   3: # Propensity Score Matching and Weighting
   4: # -----------------------------------------------------------------------------
   5: 
   6: 
   7: # -----------------------------------------------------------------------------
   8: # Download and Load required packages
   9: # -----------------------------------------------------------------------------
  10: # To download remove # in next line
  11: #install.packages(c("MatchIt","cobalt","tidyverse","janitor","haven", "here","WeightIt","marginaleffects","writexl","readxl"))
  12: 
  13: # Load required packages
  14: library(MatchIt) # For matching methods
  15: library(cobalt) # For covariate balance diagnostics
  16: library(tidyverse) # For data manipulation
  17: library(janitor) # For cleaning data
  18: library(haven) # For reading Stata data files
  19: library(here) # For reproducible paths
  20: library(WeightIt) # For propensity score weighting
  21: library(marginaleffects) # For calculating marginal effects and robust standard errors
  22: library(writexl) # To Write excel files 
  23: library(readxl) # To Read in excel files
  24: 
  25: 
  26: # -----------------------------------------------------------------------------
  27: # Set up working directory and load the data
  28: # -----------------------------------------------------------------------------
  29: 
  30: # Set the working directory using setwd() to where you have your file saved.
  31: # This ensures that the current working directory is set correctly and is where you will read in files. 
  32: setwd("C:/Users/mrodgers/AIR/S-STEM REC - Documents/General/Webinars/QED Webinars/R-code/")
  33: 
  34: # Read in the data using a relative path for reproducibility
  35: example_dat <- read_excel("AAAS_pt2_ex_rdata.xlsx")
  36: 
  37: # View the initial rows of the data to check it's loaded correctly
  38: View(example_dat)
  39: names(example_dat)
  40: # [1] "s_ID"        "treat"       "nonwhite"    "age"         "male"        "sesp_pell"   "engin_major" "actcomp"     "collgpa_yr1" 
  41: 
  42: # Display the frequency distribution of the treatment variable
  43: tabyl(example_dat$treat)
  44: 
  45: 
  46: 
  47: 
  48: ####################
  49: # Matching Example #
  50: ####################
  51: 
  52: # -----------------------------------------------------------------------------
  53: # Initial Balance Assessment (No Matching or Weighting)
  54: # -----------------------------------------------------------------------------
  55: 
  56: # Here we assess the baseline balance between treatment and control groups
  57: # without any matching or weighting applied.
  58: 
  59: # Fit a basic model to assess baseline balance
  60: m.out <- matchit(
  61:   treat ~ nonwhite + age + male + sesp_pell + engin_major + actcomp,
  62:   data = example_dat,
  63:   method = NULL # Using NULL method means no matching
  64: )
  65: 
  66: # Generate summary statistics to assess initial balance between groups
  67: summary(m.out)
  68: 
  69: # Visualize the balance across all covariates before matching
  70: plot(summary(m.out))
  71: 
  72: # -----------------------------------------------------------------------------
  73: # Propensity Score (PS) Matching
  74: # -----------------------------------------------------------------------------
  75: 
  76: # Selected Matching Model - Nearest Neighbor 1:1 Matching without Caliper or other options
  77: m.out.matched <- matchit(
  78:   treat ~ nonwhite + age + male + sesp_pell + engin_major,
  79:   data = example_dat,
  80:   distance = "logit",
  81:   method = "nearest" 
  82: )
  83: 
  84: ##### CODE IN THIS BLOCK NOT RUN DURING WEBINAR #######################################
  85: # This code runs the same matching process as above, but provides additional 
  86: # details related to different matching approaches - some details are provided below,
  87: # but see type ?matchit for package details on different options. 
  88: # m.out.matched <- matchit(
  89: #  treat ~ nonwhite + age + male + sesp_pell + engin_major,
  90: #  data = example_dat,
  91: #  #exact = "actcomp", # use this option to do exact matching
  92: #  distance = "logit",
  93: #  #caliper = 0.25, # use this to play with the caliper
  94: #  #replace = F, # without replacement is default "F" and with replacement is "T"
  95: #  #ratio = 1, # k:1 control to treatment units; default if 1:1
  96: #  method = "nearest" # other options include "optimal", "full", "subclass", and "exact"
  97: #)
  98: ######################################################################################
  99: 
 100: 
 101: # Extract and attach propensity scores to the original dataset
 102: example_dat$propensity_score <- m.out.matched$distance
 103: example_dat$matched <- m.out.matched$weights
 104: 
 105: # Display the initial rows showing the treatment and the corresponding propensity scores
 106: head(select(example_dat, treat, propensity_score))
 107: 
 108: # Display balance statistics after matching
 109: m.out.matched
 110: summary(m.out.matched)
 111: 
 112: # Graphically plot the Adjusted and Unadjusted Mean Differences for each Covariate. 
 113: plot_match <- love.plot(m.out.matched,colors = c("black", "gray70"))
 114: plot_match + geom_vline(xintercept = c(-0.25, 0.25), linetype = "dashed", linewidth = 0.5) +
 115:   theme(text = element_text(size = 14))
 116: 
 117: # -----------------------------------------------------------------------------
 118: # Group Exercise for Propensity Score (PS) Matching 
 119: # -----------------------------------------------------------------------------
 120: 
 121: # Experiment by adding a single covariate (e.g. "actcomp")
 122: 
 123: m.out.matched_r1 <- matchit(
 124:   treat ~ nonwhite + age + male + sesp_pell + engin_major + actcomp,
 125:   data = example_dat,
 126:   distance = "logit",
 127:   method = "nearest" 
 128: )
 129: 
 130: # Extract and attach propensity scores to the original dataset
 131: example_dat$propensity_score <- m.out.matched_r1$distance
 132: 
 133: # Display the initial rows showing the treatment and the corresponding propensity scores
 134: head(select(example_dat, treat, propensity_score))
 135: 
 136: # Display balance statistics after matching
 137: m.out.matched_r1
 138: summary(m.out.matched_r1)
 139: plot(summary(m.out.matched_r1))
 140: plot(m.out.matched_r1, type = "jitter", interactive = F)
 141: 
 142: ################## CODE IN THIS BLOCK NOT RUN DURING WEBINAR ########################################
 143: # ----------------------------------------------------------------------------------------------
 144: # Outcome Modeling - Simple Regression model with full dataset -Results not presented in webinar
 145: # ----------------------------------------------------------------------------------------------
 146: 
 147: # Simple regression for treatment effect
 148: naive.model <- lm(collgpa_yr1 ~ treat, data = example_dat)
 149: summary(naive.model)
 150: 
 151: # -----------------------------------------------------------------------------
 152: # Outcome Modeling after Matching - Results Not presented in webinar
 153: # -----------------------------------------------------------------------------
 154: 
 155: # Create matched dataset #
 156: example_dat_matched <- example_dat %>% filter(matched == 1)
 157: 
 158: # Run Model with no covariates with matched dataset
 159: matched.model <- lm(collgpa_yr1 ~ treat, data = example_dat_matched)
 160: summary(matched.model)
 161: 
 162: # -----------------------------------------------------------------------------
 163: # Outcome Modeling (Doubly Robust) after Matching - Results Not presented in webinar
 164: # -----------------------------------------------------------------------------
 165: 
 166: # Create matched dataset #
 167: example_dat_matched <- example_dat %>% filter(matched == 1)
 168: 
 169: # Run Model with all covariates - Doubly Robust Model
 170: matched.model <- lm(collgpa_yr1 ~ treat + nonwhite + age + male + sesp_pell + engin_major + actcomp,
 171:                     data = example_dat_matched)
 172: summary(matched.model)
 173: #############################################################################################
 174: 
 175: 
 176: 
 177: 
 178: 
 179: 
 180: #####################
 181: # Weighting Example #
 182: #####################
 183: 
 184: # -----------------------------------------------------------------------------
 185: # Initial Balance Assessment (No Matching or Weighting)
 186: # -----------------------------------------------------------------------------
 187: 
 188: # Here we assess the baseline balance between treatment and control groups
 189: # without any matching or weighting applied.
 190: 
 191: # Reload raw data
 192: example_dat <- read_excel("AAAS_pt2_ex_rdata.xlsx")
 193: 
 194: # Fit a basic model to assess baseline balance
 195: m.out <- matchit(
 196:   treat ~ nonwhite + age + male + sesp_pell + engin_major + actcomp,
 197:   data = example_dat,
 198:   method = NULL # Using NULL method means no matching
 199: )
 200: 
 201: # Generate summary statistics to assess initial balance between groups
 202: summary(m.out)
 203: 
 204: # Visualize the balance across all covariates before matching
 205: plot(summary(m.out))
 206: 
 207: # -----------------------------------------------------------------------------
 208: # Propensity Score (PS) Weighting
 209: # -----------------------------------------------------------------------------
 210: 
 211: # IPTW for Average Treatment Effect on the Treated (ATT)
 212: w.out <- weightit(
 213:   treat ~ nonwhite + age + male + sesp_pell + engin_major + actcomp,
 214:   data = example_dat,
 215:   method = "ps",
 216:   estimand = "ATT")
 217: 
 218: w.out
 219: summary(w.out)
 220: 
 221: bal.tab(w.out, stats = c("m", "v"), thresholds = c(m = .25))
 222: 
 223: ################## CODE IN THIS BLOCK NOT RUN DURING WEBINAR ############
 224: # IPTW for Average Treatment Effect on the Treated (ATE)
 225: #w.out_ATE <- weightit(
 226: #  treat ~ nonwhite + age + male + sesp_pell + engin_major + actcomp,
 227: #  data = example_dat,
 228: #  method = "ps",
 229: #  estimand = "ATE"
 230: #)
 231: #########################################################################
 232: 
 233: # -----------------------------------------------------------------------------
 234: # Outcome Modeling after Weighting
 235: # -----------------------------------------------------------------------------
 236: 
 237: # IPTW for Average Treatment Effect (ATT)
 238: w.out <- weightit(
 239:   treat ~ nonwhite + age + male + sesp_pell + engin_major + actcomp,
 240:   data = example_dat,
 241:   method = "ps",
 242:   estimand = "ATT"
 243: )
 244: 
 245: 
 246: # Fit a simple linear regression model on the weighted data without any covariates
 247: weighted.data <- example_dat
 248: weighted.data$weights <- w.out$weights
 249: 
 250: outcome.model.weighted <- lm(collgpa_yr1 ~ treat, 
 251:                              data = weighted.data, 
 252:                              weights = weighted.data$weights)  
 253: 
 254: summary(outcome.model.weighted)
 255: 
 256: # -----------------------------------------------------------------------------
 257: # Doubly Robust Estimation after Weighting
 258: # -----------------------------------------------------------------------------
 259: 
 260: # Incorporate covariates in the outcome model to make it doubly robust
 261: doubly.robust.model.weighted <- lm(collgpa_yr1 ~ treat + nonwhite + age + male + sesp_pell + engin_major + actcomp,
 262:                                    data = weighted.data, weights = weighted.data$weights)
 263: 
 264: summary(doubly.robust.model.weighted)