R Notebook

Initial Data Preparation

raw_data <- read.csv("https://tinyurl.com/4736xjmw")

library(dplyr)

# Collapse the dataset to one observation per agreement (cas_id)
# We keep the row with the latest 'year' to capture the final outcome.
data <- raw_data %>%
  group_by(cas_id) %>%
  slice_max(order_by = year, n = 1) %>% # Selects the row with the highest year
  ungroup()

# ============ PREP VARIABLES IN ORIGINAL DATA =================
data$start_date_obj <- as.Date(data$strt_pa)
data$end_date_obj   <- as.Date(data$nd_pa)
data$start_date_num <- as.numeric(data$start_date_obj)
data$weeks_of_peace <- as.numeric(difftime(data$end_date_obj, data$start_date_obj, units = "weeks"))
data$days_of_peace  <- as.numeric(data$end_date_obj - data$start_date_obj)

# Check the new dimensions
str(data)

tibble [79 × 37] (S3: tbl_df/tbl/data.frame)
 $ PAID          : num [1:79] 7.0e+12 7.0e+12 7.0e+12 5.4e+12 5.4e+12 ...
 $ cas_id        : int [1:79] 1 2 3 4 5 6 7 8 9 10 ...
 $ new_con_id    : int [1:79] 333 333 333 387 327 327 322 398 389 287 ...
 $ iso3c         : chr [1:79] "AFG" "AFG" "AFG" "AGO" ...
 $ year          : int [1:79] 1993 1993 1996 2007 1998 2011 2011 2011 2011 2000 ...
 $ pa_name       : chr [1:79] "Islamabad accord" "Jalalabad agreement" "Mahipar agreement" "Memorandum of Understanding on Peace and National Reconciliation in Cabinda province" ...
 $ strt_pa       : chr [1:79] "1993-03-07" "1993-05-20" "1996-05-24" "2006-08-01" ...
 $ nd_pa         : chr [1:79] "1993-04-13" "1993-08-22" "1996-05-30" "2007-12-31" ...
 $ pc_nd         : int [1:79] 1 1 1 1 1 0 0 0 0 1 ...
 $ nd_pa_src     : chr [1:79] "GED" "GED" "GED" "ACD" ...
 $ strt_sgmnt    : chr [1:79] "1993-03-07" "1993-05-20" "1996-05-24" "2007-01-01" ...
 $ nd_sgmnt      : chr [1:79] "1993-04-13" "1993-08-22" "1996-05-30" "2007-12-31" ...
 $ fem_sig       : int [1:79] 0 0 0 0 0 0 0 0 0 0 ...
 $ con_dur       : int [1:79] 16 16 19 14 20 28 17 2 4 36 ...
 $ con_dur2      : int [1:79] 16 16 19 7 20 26 17 2 4 10 ...
 $ con_int       : num [1:79] 1 1 1 0 0.95 ...
 $ con_int2      : num [1:79] 1 1 1 0 0.95 ...
 $ con_int_bin   : int [1:79] 1 1 1 0 1 1 0 1 1 1 ...
 $ war_par       : int [1:79] 12 13 14 3 3 3 2 3 3 7 ...
 $ pow_sha       : int [1:79] 2 2 2 1 3 2 1 3 2 2 ...
 $ dem           : int [1:79] 0 0 0 0 0 0 1 0 0 0 ...
 $ gdp           : num [1:79] NA NA NA 3.31 2.7 ...
 $ hum_rig       : int [1:79] NA NA 3 5 3 4 7 6 6 5 ...
 $ cpa           : int [1:79] 0 0 0 0 1 0 1 0 1 0 ...
 $ con_iss       : int [1:79] 1 1 1 0 1 1 0 0 0 1 ...
 $ pea           : int [1:79] 0 0 0 0 1 0 0 0 0 0 ...
 $ fem_leg       : int [1:79] NA NA NA 2 2 3 2 2 2 1 ...
 $ fem_com       : int [1:79] 1 1 1 0 3 3 0 1 1 2 ...
 $ quo           : int [1:79] 0 0 0 0 0 0 1 0 0 1 ...
 $ lef_ide       : int [1:79] 0 0 0 0 0 0 0 0 0 0 ...
 $ civ_soc       : int [1:79] 0 0 0 0 0 0 1 0 1 1 ...
 $ un_pea_con    : num [1:79] 0 0 0 0 1216 ...
 $ start_date_obj: Date[1:79], format: "1993-03-07" "1993-05-20" ...
 $ end_date_obj  : Date[1:79], format: "1993-04-13" "1993-08-22" ...
 $ start_date_num: num [1:79] 8466 8540 9640 13361 9089 ...
 $ weeks_of_peace: num [1:79] 5.286 13.429 0.857 73.857 179.857 ...
 $ days_of_peace : num [1:79] 37 94 6 517 1259 ...

Hypothesis 1: Women participation impacts peace durability

# ===== CLEANING: Remove NAs from ALL necessary variables ====

# Define variables of interest

pretreatment_covariates <- c("con_dur", "con_int", "con_iss", "war_par", "gdp", "dem", "hum_rig", "fem_com", "lef_ide", "fem_leg", "start_date_num")

treatnent_var <- c("fem_sig")
outcome_vars <- c("days_of_peace", "pc_nd")

cleaning_vars <- c(treatnent_var, outcome_vars, pretreatment_covariates) # This list includes Y, Tr, and X. We use it ONLY for na.omit.

# Check which columns contain missing values
colSums(is.na(data[, cleaning_vars])) 

# Clean the data to enable genetic matching (cannot have NAs)
df_clean_H1 <- na.omit(data[, c("cas_id", cleaning_vars)])
df_clean_H1

cat("Number of control/treatment units in sample (before cleaning):")
table(data$fem_sig)

cat("Number of control/treatment units in sample (after cleaning):")
table(df_clean_H1$fem_sig)

library(rgenoud)
library(Matching)

set.seed(123)

# ====== DEFINE MATCHING OBJECTS ======

# Treatment 
Tr <- df_clean_H1$fem_sig

# Matching Matrix (X): ONLY contains pre-treatment covariates.
matching_vars <- pretreatment_covariates
X <- df_clean_H1[, matching_vars]
# Do not use outcome observations until matching is complete


# PREPARE CALIPER (for start_date_num)
# Calculate caliper value for start_date_num
desired_window_days <- 180
sd_date <- sd(df_clean_H1$start_date_num)
caliper_value <- desired_window_days / sd_date

# Create Caliper Vector
# We need a vector of Infs matching the number of columns in X
caliper_vec <- rep(Inf, ncol(X))

# Apply restriction ONLY to the start date column
date_index <- which(colnames(X) == "start_date_num")
caliper_vec[date_index] <- caliper_value

# ========== RUN GENMATCH ==============
genout_H1 <- GenMatch(Tr = Tr, X = X, BalanceMatrix = X, 
                   caliper = caliper_vec,
                   estimand = "ATT", pop.size = 5000, 
                   print.level = 0)

# ========== RUN MATCH ==============

m_out_H1 <- Match(Tr = Tr, X = X, Weight.matrix = genout_H1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Report on matching summary ======

# Count how many times each treated unit appears in the match list
matches_per_treated <- table(m_out_H1$index.treated)

# Check the distribution (e.g., how many got 1 match, how many got 2, etc.)
table(matches_per_treated)

# Total treated units in original data
total_treated <- sum(Tr)

# Total treated units that successfully found a match
matched_treated <- length(unique(m_out_H1$index.treated))

cat("Original Treated:", total_treated, "\n")
cat("Matched Treated: ", matched_treated, "\n")
cat("Dropped:         ", total_treated - matched_treated, "\n")

# ===== Report Balance and Lowest P-Value ======
set.seed(123)
# treatment indicator is the dependent variable, check balance for all covariates
mb_formula_H1 <- fem_sig ~ con_dur + con_int + con_iss +  war_par + gdp + dem + hum_rig + fem_leg + fem_com + lef_ide + start_date_num

MatchBalance(mb_formula_H1, data = df_clean_H1, match.out = m_out_H1, nboots = 500)

#install.packages("cobalt")
library(cobalt)

love.plot(mb_formula_H1, 
          data = df_clean_H1, 
          weights = m_out_H1, # Use the match object as the source of weights
          binary = "std",
          thresholds = c(m = .25),
          alpha = 0.7,
          title = "Covariate Balance (SMD)",
          grid = TRUE)

love.plot(mb_formula_H1, 
          data = df_clean_H1, 
          weights = m_out_H1, # Use the match object as the source of weights
          stats = "ks.statistics", # Change to KS statistics
          binary = "std",
          thresholds = c(ks = 0.1), 
          alpha = 0.7,
          grid = TRUE,
          title = "Covariate Balance (KS Statistics)")

library(ggplot2)
library(gridExtra) # Required for arranging plots

# Panel A: Good Balance (Continuous - ECDF)
p1 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "start_date_num", type = "ecdf") + 
      ggtitle("A. Good Balance: Start Date")

# Panel B: Poor Balance (Continuous - ECDF)
p2 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "gdp", type = "ecdf") + 
      ggtitle("B. Poor Balance: GDP")

# Panel C: Good Balance (Discrete - Histogram)
p3 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "dem", type = "histogram") + 
      ggtitle("C. Good Balance: Democracy")

# Panel D: Poor Balance (Discrete - Histogram)
p4 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "hum_rig", type = "histogram") + 
      ggtitle("D. Poor Balance: Human Rights")

# Arrange plots into a 2x2 grid
grid.arrange(p1, p2, p3, p4, nrow = 2)

bal.tab(mb_formula_H1, 
        data = df_clean_H1, 
        weights = m_out_H1, # Use the match object as the source of weights
        estimand = "ATT",
        binary = "std",      # <--- Forces binary vars to use SMD
        continuous = "std",  # <--- Forces continuous vars (default) to use SMD
        un = TRUE,
        stats = c("m", "v", "ks"),   # Mean diffs, Variance ratios, KS stats
        thresholds = c(m = 0.25, ks = 0.1))

library(Matching)

get_balance_table <- function(formula, data, match_out, nboots = 500) {
  
  # 1. Run MatchBalance
  mb <- MatchBalance(formula, data = data, match.out = match_out, nboots = nboots, print.level = 0)
  
  # 2. Define the "Nested" Null Coalescing Operator 
  #    (This allows us to chain checks: "Try A, if NULL try B, if NULL try C")
  `%||%` <- function(a, b) if (is.null(a)) b else a

  # 3. Define the Extractors
  #    T-Test: Simple check
  get_tt <- function(x) x$tt$p.value %||% NA
  
  #    KS-Test: The "Nested" Approach
  #    1. Try nested bootstrap (x$ks$ks.boot.pvalue) <--- PRIORITY
  #    2. Try top-level bootstrap (x$ks.boot.pvalue)
  #    3. Try standard KS (x$ks$p.value)
  get_ks <- function(x) x$ks$ks.boot.pvalue %||% x$ks.boot.pvalue %||% x$ks$p.value %||% NA

  # 4. Extract efficiently with sapply
  idx <- seq_along(mb$AfterMatching)
  vars <- all.vars(formula)[-1]
  final_names <- if(length(vars) == length(idx)) vars else paste0("Var_", idx)

  df <- data.frame(
    Variable  = final_names,
    Pre_T_P   = round(sapply(idx, function(i) get_tt(mb$BeforeMatching[[i]])), 3),
    Post_T_P  = round(sapply(idx, function(i) get_tt(mb$AfterMatching[[i]])), 3),
    Pre_KS_P  = round(sapply(idx, function(i) get_ks(mb$BeforeMatching[[i]])), 3),
    Post_KS_P = round(sapply(idx, function(i) get_ks(mb$AfterMatching[[i]])), 3)
  )
  
  # 5. Status: Pass if (KS is NA or >0.05) AND (T-test > 0.05)
  df$Status <- ifelse((is.na(df$Post_KS_P) | df$Post_KS_P > 0.1) & df$Post_T_P > 0.1, "PASS", "FAIL")
  
  return(df)
}

set.seed(123)
# --- Call Function ---
balance_table_H1 <- get_balance_table(
  formula = mb_formula_H1, 
  data = df_clean_H1, 
  match_out = m_out_H1, 
  nboots = 500
)

balance_table_H1

Hypothesis 1A: fem-sig on days_of_peace (peace duration)

# Include outcome variable for effect estimate
Y  <- df_clean_H1$days_of_peace
m_out_H1A <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_H1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Estimate Treatment Effect, SE, and CI ======
summary(m_out_H1A)

ATT_est <- m_out_H1A$est
ATT_se <- m_out_H1A$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 2)))
print(paste("Standard Error:", round(ATT_se, 2)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

Hypothesis 1B: fem-sig on pc_nd (peace ends and conflict resumes)

# Outcome
Y  <- df_clean_H1$pc_nd
m_out_H1B <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_H1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Estimate Treatment Effect, SE, and CI ======

summary(m_out_H1B)

ATT_est <- m_out_H1B$est
ATT_se <- m_out_H1B$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 2)))
print(paste("Standard Error:", round(ATT_se, 2)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

Model H1: Rosenbaum Bounds Sensitivity Analysis

# ===== CLEANING: Remove NAs from ALL necessary variables ====

# Define variables of interest

pretreatment_covariates <- c("con_dur", "con_int", "con_iss", "war_par", "gdp", "dem", "hum_rig", "fem_com", "lef_ide", "fem_leg", "start_date_num")

treatnent_var <- c("fem_sig")
outcome_vars <- c("days_of_peace", "pc_nd")

cleaning_vars <- c(treatnent_var, outcome_vars, pretreatment_covariates) # This list includes Y, Tr, and X. We use it ONLY for na.omit.

# Check which columns contain missing values
colSums(is.na(data[, cleaning_vars])) 

       fem_sig  days_of_peace          pc_nd        con_dur        con_int 
             0              0              0              0              0 
       con_iss        war_par            gdp            dem        hum_rig 
             0              0              9              0              9 
       fem_com        lef_ide        fem_leg start_date_num 
             9              2              9              0 

# Clean the data to enable genetic matching (cannot have NAs)
df_clean_H1 <- na.omit(data[, c("cas_id", cleaning_vars)])
df_clean_H1

cat("Number of control/treatment units in sample (before cleaning):")

Number of control/treatment units in sample (before cleaning):

table(data$fem_sig)

 0  1 
73  6 

cat("Number of control/treatment units in sample (after cleaning):")

Number of control/treatment units in sample (after cleaning):

table(df_clean_H1$fem_sig)

 0  1 
50  6 

library(rgenoud)
library(Matching)

set.seed(123)

# ====== DEFINE MATCHING OBJECTS ======

# Treatment 
Tr <- df_clean_H1$fem_sig

# Matching Matrix (X): ONLY contains pre-treatment covariates.
matching_vars <- pretreatment_covariates
X <- df_clean_H1[, matching_vars]
# Do not use outcome observations until matching is complete


# PREPARE CALIPER (for start_date_num)
# Calculate caliper value for start_date_num
desired_window_days <- 180
sd_date <- sd(df_clean_H1$start_date_num)
caliper_value <- desired_window_days / sd_date

# Create Caliper Vector
# We need a vector of Infs matching the number of columns in X
caliper_vec <- rep(Inf, ncol(X))

# Apply restriction ONLY to the start date column
date_index <- which(colnames(X) == "start_date_num")
caliper_vec[date_index] <- caliper_value

# ========== RUN GENMATCH ==============
genout_H1 <- GenMatch(Tr = Tr, X = X, BalanceMatrix = X, 
                   caliper = caliper_vec,
                   estimand = "ATT", pop.size = 5000, 
                   print.level = 0)

# ========== RUN MATCH ==============

m_out_H1 <- Match(Tr = Tr, X = X, Weight.matrix = genout_H1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Report on matching summary ======

# Count how many times each treated unit appears in the match list
matches_per_treated <- table(m_out_H1$index.treated)

# Check the distribution (e.g., how many got 1 match, how many got 2, etc.)
table(matches_per_treated)

matches_per_treated
1 
6 

# Total treated units in original data
total_treated <- sum(Tr)

# Total treated units that successfully found a match
matched_treated <- length(unique(m_out_H1$index.treated))

cat("Original Treated:", total_treated, "\n")

Original Treated: 6 

cat("Matched Treated: ", matched_treated, "\n")

Matched Treated:  6 

cat("Dropped:         ", total_treated - matched_treated, "\n")

Dropped:          0 

H1: Robustness check

# ============== Remove NAs from ALL necessary variables==================

# Specify variables of interest
pretreatment_covariates <- c("con_dur", "con_int", "con_iss", "war_par", "gdp", "dem", "hum_rig", "fem_com", "lef_ide", "fem_leg", "start_date_num")
treatnent_var <- c("fem_sig")
outcome_vars <- c("days_of_peace", "pc_nd")

# additional covariates for matching (used as controls by Krause et al.)
add_covariates <- c("pow_sha", "quo", "un_pea_con") 

cleaning_vars <- c(treatnent_var, outcome_vars, pretreatment_covariates, add_covariates)

# Check which columns contain missing values
colSums(is.na(data[, cleaning_vars])) 

# Clean the data to enable genetic matching (cannot have NAs)
df_clean_rb1 <- na.omit(data[, c("cas_id", cleaning_vars)])
df_clean_rb1

# Check for dropped units
cat("Number of control/treatment units in sample (before cleaning):")
table(data$fem_sig)

cat("Number of control/treatment units in sample (after cleaning):")
table(df_clean_rb1$fem_sig)

set.seed(123)
# ============ DEFINE MATCHING OBJECTS ================

# Treatment 
Tr <- df_clean_rb1$fem_sig

# Matching Matrix (X)
matching_vars <- c(pretreatment_covariates, add_covariates) #include additional covariates
X <- df_clean_rb1[, matching_vars]


# PREPARE CALIPER (For start_date_num)
desired_window_days <- 180 #90 causes treatment units to be dropped
sd_date <- sd(df_clean_rb1$start_date_num)
caliper_value <- desired_window_days / sd_date

# Create Caliper Vector
# We need a vector of Infs matching the number of columns in X
caliper_vec <- rep(Inf, ncol(X))

# Apply restriction ONLY to the date column
date_index <- which(colnames(X) == "start_date_num")
caliper_vec[date_index] <- caliper_value

# ================ RUN GENMATCH ===============
genout_rb1 <- GenMatch(Tr = Tr, X = X, BalanceMatrix = X,
                   caliper = caliper_vec,
                   estimand = "ATT", pop.size = 5000, 
                   print.level = 0)

# ================ RUN MATCH ===================
m_out_rb1 <- Match(Tr = Tr, X = X, Weight.matrix = genout_rb1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Report Summary of Matching Results ======
# Count how many times each treated unit appears in the match list
matches_per_treated <- table(m_out_rb1$index.treated)
# Check the distribution (e.g., how many got 1 match, how many got 2, etc.)
table(matches_per_treated)

# Total treated units in original data
total_treated <- sum(Tr)

# Total treated units that successfully found a match
matched_treated <- length(unique(m_out_rb1$index.treated))

cat("Original Treated:", total_treated, "\n")
cat("Matched Treated: ", matched_treated, "\n")
cat("Dropped:         ", total_treated - matched_treated, "\n")

# ===== Report Balance and Lowest P-Value ======
set.seed(123)
# treatment indicator is the dependent variable, check balance for all covariates
mb_formula_rb1 <- fem_sig ~ con_dur + con_int + con_iss +  war_par + gdp + dem + hum_rig + fem_leg + fem_com + lef_ide + start_date_num + pow_sha + quo + un_pea_con

MatchBalance(mb_formula_rb1, data = df_clean_rb1, match.out = m_out_rb1, nboots = 500)

#install.packages("cobalt")
library(cobalt)

love.plot(mb_formula_rb1, 
          data = df_clean_rb1, 
          weights = m_out_rb1, # Use the match object as the source of weights
          binary = "std",
          thresholds = c(m = .25),
          alpha = 0.7,
          title = "Covariate Balance (SMD)",
          grid = TRUE)

love.plot(mb_formula_rb1, 
          data = df_clean_rb1, 
          weights = m_out_rb1, # Use the match object as the source of weights
          stats = "ks.statistics", # Change to KS statistics
          binary = "std",
          thresholds = c(ks = 0.1), 
          alpha = 0.7,
          grid = TRUE,
          title = "Covariate Balance (KS Statistics)")

bal.tab(mb_formula_rb1, 
        data = df_clean_rb1, 
        weights = m_out_rb1, # Use the match object as the source of weights
        estimand = "ATT",
        binary = "std",      # <--- Forces binary vars to use SMD
        continuous = "std",  # <--- Forces continuous vars (default) to use SMD
        un = TRUE,
        stats = c("m", "v", "ks"),   # Mean diffs, Variance ratios, KS stats
        thresholds = c(m = 0.25, ks = 0.1))

set.seed(123)
# --- Call Function ---
balance_table_rb1 <- get_balance_table(
  formula = mb_formula_rb1, 
  data = df_clean_rb1, 
  match_out = m_out_rb1, 
  nboots = 500
)

balance_table_rb1

H1 vs R1: Comparison of Matching Results

matched_data_H1 <- rbind(
  df_clean_H1[m_out_H1$index.treated, ],
  df_clean_H1[m_out_H1$index.control, ]
)
n_pairs <- length(m_out_H1$index.treated)
matched_data_H1$pair_id <- c(1:n_pairs, 1:n_pairs)

matched_data_rb1 <- rbind(
  df_clean_rb1[m_out_rb1$index.treated, ],
  df_clean_rb1[m_out_rb1$index.control, ]
)
n_pairs <- length(m_out_rb1$index.treated)
matched_data_rb1$pair_id <- c(1:n_pairs, 1:n_pairs)

matches <- cbind(Pair_ID = matched_data_H1$pair_id, 
                 Treated = matched_data_H1$cas_id[matched_data_H1$fem_sig == 1], 
                 ModelH1_Matches = matched_data_H1$cas_id[matched_data_H1$fem_sig == 0], 
                 ModelR1_Matches = matched_data_rb1$cas_id[matched_data_rb1$fem_sig == 0])
matches[1:6, ]

# Compare outcome observations of matched units
matchcomp_days <- cbind(Pair_ID = matched_data_H1$pair_id,
                        cas_id = matched_data_H1$cas_id,
                        fem_sig = matched_data_H1$fem_sig,
                        days_of_peace_H1 = matched_data_H1$days_of_peace,
                        days_of_peace_rb1 = matched_data_rb1$days_of_peace)
matchcomp_pc <- cbind(Pair_ID = matched_data_H1$pair_id,
                      cas_id = matched_data_H1$cas_id,
                      fem_sig = matched_data_H1$fem_sig,
                      pc_nd_H1 = matched_data_H1$pc_nd,
                      pc_nd_rb1 = matched_data_rb1$pc_nd)
matchcomp_days
matchcomp_pc

H1 vs R1: Comparison of Treatment Effect Estimates

# ======== Run Match() with outcome variable =======

# Outcome
Y  <- df_clean_rb1$days_of_peace
m_out_rb1A <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_rb1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Estimate Treatment Effect, SE, and CI ======

summary(m_out_rb1A)

ATT_est <- m_out_rb1A$est
ATT_se <- m_out_rb1A$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 4)))
print(paste("Standard Error:", round(ATT_se, 4)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

# ======== Run Match() with outcome variable =======
# Outcome
Y  <- df_clean_rb1$pc_nd
m_out_rb1B <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_rb1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Estimate Treatment Effect, SE, and CI ====

summary(m_out_rb1B)

ATT_est <- m_out_rb1B$est
ATT_se <- m_out_rb1B$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 4)))
print(paste("Standard Error:", round(ATT_se, 4)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

print ("H1A vs. R1A")
summary(m_out_H1A)
summary(m_out_rb1A)
print ("H1B vs. R1B")
summary(m_out_H1B)
summary(m_out_rb1B)

Hypothesis 1: Survival Analysis

# ===== Report Balance and Lowest P-Value ======
set.seed(123)
# treatment indicator is the dependent variable, check balance for all covariates
mb_formula_H1 <- fem_sig ~ con_dur + con_int + con_iss +  war_par + gdp + dem + hum_rig + fem_leg + fem_com + lef_ide + start_date_num

MatchBalance(mb_formula_H1, data = df_clean_H1, match.out = m_out_H1, nboots = 500)

***** (V1) con_dur *****
                       Before Matching       After Matching
mean treatment........     27.333            27.333 
mean control..........      20.68              19.5 
std mean diff.........      47.54            55.971 

mean raw eQQ diff.....          6            7.8333 
med  raw eQQ diff.....        3.5               6.5 
max  raw eQQ diff.....         16                19 

mean eCDF diff........    0.14162           0.13636 
med  eCDF diff........    0.12667           0.16667 
max  eCDF diff........    0.35333           0.33333 

var ratio (Tr/Co).....    0.84093           0.76302 
T-test p-value........    0.31464           0.15974 
KS Bootstrap p-value..      0.382             0.826 
KS Naive p-value......    0.39891           0.93074 
KS Statistic..........    0.35333           0.33333 


***** (V2) con_int *****
                       Before Matching       After Matching
mean treatment........    0.11265           0.11265 
mean control..........    0.26814          0.082458 
std mean diff.........    -86.378             16.77 

mean raw eQQ diff.....    0.19184          0.040991 
med  raw eQQ diff.....    0.14263          0.016204 
max  raw eQQ diff.....       0.55           0.10729 

mean eCDF diff........    0.18784          0.071429 
med  eCDF diff........    0.18333                 0 
max  eCDF diff........    0.38667           0.16667 

var ratio (Tr/Co).....     0.3793            1.8253 
T-test p-value........   0.099899           0.72072 
KS Bootstrap p-value..      0.244             0.998 
KS Naive p-value......    0.30081                 1 
KS Statistic..........    0.38667           0.16667 


***** (V3) con_iss *****
                       Before Matching       After Matching
mean treatment........    0.66667           0.66667 
mean control..........       0.78               0.5 
std mean diff.........    -21.947            32.275 

mean raw eQQ diff.....          0           0.16667 
med  raw eQQ diff.....          0                 0 
max  raw eQQ diff.....          0                 1 

mean eCDF diff........   0.056667          0.083333 
med  eCDF diff........   0.056667          0.083333 
max  eCDF diff........    0.11333           0.16667 

var ratio (Tr/Co).....     1.5229           0.88889 
T-test p-value........    0.62384           0.57831 


***** (V4) war_par *****
                       Before Matching       After Matching
mean treatment........          5                 5 
mean control..........       6.86            6.1667 
std mean diff.........    -73.523           -46.117 

mean raw eQQ diff.....          3            2.1667 
med  raw eQQ diff.....        0.5               0.5 
max  raw eQQ diff.....         15                10 

mean eCDF diff........   0.086061          0.071429 
med  eCDF diff........        0.1                 0 
max  eCDF diff........       0.18           0.16667 

var ratio (Tr/Co).....    0.16592           0.17128 
T-test p-value........    0.19153           0.56147 
KS Bootstrap p-value..      0.858             0.994 
KS Naive p-value......    0.85647                 1 
KS Statistic..........       0.18           0.16667 


***** (V5) gdp *****
                       Before Matching       After Matching
mean treatment........     3.2898            3.2898 
mean control..........     2.7526            2.7601 
std mean diff.........     66.288            65.373 

mean raw eQQ diff.....    0.51179            0.5297 
med  raw eQQ diff.....    0.52796            0.5241 
max  raw eQQ diff.....    0.82098           0.93254 

mean eCDF diff........    0.20213           0.19444 
med  eCDF diff........    0.18667           0.16667 
max  eCDF diff........    0.46667               0.5 

var ratio (Tr/Co).....     2.2964            2.0773 
T-test p-value........    0.16866           0.11621 
KS Bootstrap p-value..      0.128             0.354 
KS Naive p-value......    0.13832           0.47403 
KS Statistic..........    0.46667               0.5 


***** (V6) dem *****
                       Before Matching       After Matching
mean treatment........    0.83333           0.83333 
mean control..........       0.32           0.66667 
std mean diff.........     68.192             22.14 

mean raw eQQ diff.....        0.5           0.16667 
med  raw eQQ diff.....        0.5                 0 
max  raw eQQ diff.....          1                 1 

mean eCDF diff........    0.17111          0.055556 
med  eCDF diff........    0.16667                 0 
max  eCDF diff........    0.34667           0.16667 

var ratio (Tr/Co).....     2.5521             2.125 
T-test p-value........    0.15838           0.57831 
KS Bootstrap p-value..      0.112             0.902 
KS Naive p-value......    0.17238                 1 
KS Statistic..........    0.34667           0.16667 


***** (V7) hum_rig *****
                       Before Matching       After Matching
mean treatment........     10.167            10.167 
mean control..........       6.32            8.8333 
std mean diff.........     106.82            37.028 

mean raw eQQ diff.....     3.8333                 2 
med  raw eQQ diff.....          4               1.5 
max  raw eQQ diff.....          6                 4 

mean eCDF diff........        0.3           0.16667 
med  eCDF diff........    0.27333           0.16667 
max  eCDF diff........    0.73333               0.5 

var ratio (Tr/Co).....     1.6682            1.4461 
T-test p-value........   0.046585           0.31494 
KS Bootstrap p-value.. < 2.22e-16              0.27 
KS Naive p-value......  0.0012302           0.35065 
KS Statistic..........    0.73333               0.5 


***** (V8) fem_leg *****
                       Before Matching       After Matching
mean treatment........     1.1667            1.1667 
mean control..........       1.46            1.6667 
std mean diff.........    -38.967           -66.421 

mean raw eQQ diff.....    0.33333               0.5 
med  raw eQQ diff.....          0               0.5 
max  raw eQQ diff.....          1                 1 

mean eCDF diff........       0.08             0.125 
med  eCDF diff........       0.08          0.083333 
max  eCDF diff........       0.16           0.33333 

var ratio (Tr/Co).....    0.59816           0.53125 
T-test p-value........    0.41188           0.45135 
KS Bootstrap p-value..      0.762             0.536 
KS Naive p-value......    0.87781           0.74026 
KS Statistic..........       0.16           0.33333 


***** (V9) fem_com *****
                       Before Matching       After Matching
mean treatment........     1.8333            1.8333 
mean control..........       1.12           0.83333 
std mean diff.........     61.018             85.54 

mean raw eQQ diff.....    0.83333                 1 
med  raw eQQ diff.....          1                 1 
max  raw eQQ diff.....          2                 2 

mean eCDF diff........    0.17833              0.25 
med  eCDF diff........    0.22333           0.33333 
max  eCDF diff........    0.26667           0.33333 

var ratio (Tr/Co).....     1.1691            1.4138 
T-test p-value........    0.20391           0.11626 
KS Bootstrap p-value..       0.39             0.648 
KS Naive p-value......    0.45736           0.84416 
KS Statistic..........    0.26667           0.33333 


***** (V10) lef_ide *****
                       Before Matching       After Matching
mean treatment........    0.33333           0.33333 
mean control..........       0.04                 0 
std mean diff.........     56.804             64.55 

mean raw eQQ diff.....    0.16667           0.33333 
med  raw eQQ diff.....          0                 0 
max  raw eQQ diff.....          1                 1 

mean eCDF diff........    0.14667           0.16667 
med  eCDF diff........    0.14667           0.16667 
max  eCDF diff........    0.29333           0.33333 

var ratio (Tr/Co).....     6.8056               Inf 
T-test p-value........    0.22441           0.14381 


***** (V11) start_date_num *****
                       Before Matching       After Matching
mean treatment........      10704             10704 
mean control..........      11066             10666 
std mean diff.........    -22.292            2.3179 

mean raw eQQ diff.....     647.83            83.333 
med  raw eQQ diff.....      432.5              92.5 
max  raw eQQ diff.....       1935               129 

mean eCDF diff........        0.1          0.083333 
med  eCDF diff........   0.076667          0.083333 
max  eCDF diff........       0.34           0.16667 

var ratio (Tr/Co).....    0.82872            1.0209 
T-test p-value........    0.62658           0.32108 
KS Bootstrap p-value..      0.436                 1 
KS Naive p-value......    0.46919                 1 
KS Statistic..........       0.34           0.16667 


Before Matching Minimum p.value: < 2.22e-16 
Variable Name(s): hum_rig  Number(s): 7 

After Matching Minimum p.value: 0.11621 
Variable Name(s): gdp  Number(s): 5 

#install.packages("cobalt")
library(cobalt)

love.plot(mb_formula_H1, 
          data = df_clean_H1, 
          weights = m_out_H1, # Use the match object as the source of weights
          binary = "std",
          thresholds = c(m = .25),
          alpha = 0.7,
          title = "Covariate Balance (SMD)",
          grid = TRUE)

love.plot(mb_formula_H1, 
          data = df_clean_H1, 
          weights = m_out_H1, # Use the match object as the source of weights
          stats = "ks.statistics", # Change to KS statistics
          binary = "std",
          thresholds = c(ks = 0.1), 
          alpha = 0.7,
          grid = TRUE,
          title = "Covariate Balance (KS Statistics)")

library(ggplot2)
library(gridExtra) # Required for arranging plots

# Panel A: Good Balance (Continuous - ECDF)
p1 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "start_date_num", type = "ecdf") + 
      ggtitle("A. Good Balance: Start Date")

# Panel B: Poor Balance (Continuous - ECDF)
p2 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "gdp", type = "ecdf") + 
      ggtitle("B. Poor Balance: GDP")

# Panel C: Good Balance (Discrete - Histogram)
p3 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "dem", type = "histogram") + 
      ggtitle("C. Good Balance: Democracy")

# Panel D: Poor Balance (Discrete - Histogram)
p4 <- bal.plot(m_out_H1, treat = Tr, covs = X, 
               var.name = "hum_rig", type = "histogram") + 
      ggtitle("D. Poor Balance: Human Rights")

# Arrange plots into a 2x2 grid
grid.arrange(p1, p2, p3, p4, nrow = 2)

bal.tab(mb_formula_H1, 
        data = df_clean_H1, 
        weights = m_out_H1, # Use the match object as the source of weights
        estimand = "ATT",
        binary = "std",      # <--- Forces binary vars to use SMD
        continuous = "std",  # <--- Forces continuous vars (default) to use SMD
        un = TRUE,
        stats = c("m", "v", "ks"),   # Mean diffs, Variance ratios, KS stats
        thresholds = c(m = 0.25, ks = 0.1))

Balance Measures
                  Type Diff.Un V.Ratio.Un  KS.Un Diff.Adj         M.Threshold V.Ratio.Adj KS.Adj
con_dur        Contin.  0.4754     0.8409 0.3533   0.5597 Not Balanced, >0.25      0.7630 0.3333
con_int        Contin. -0.8638     0.3793 0.3867   0.1677     Balanced, <0.25      1.8253 0.1667
con_iss         Binary -0.2404          . 0.1133   0.3536 Not Balanced, >0.25           . 0.1667
war_par        Contin. -0.7352     0.1659 0.1800  -0.4612 Not Balanced, >0.25      0.1713 0.1667
gdp            Contin.  0.6629     2.2964 0.4667   0.6537 Not Balanced, >0.25      2.0773 0.5000
dem            Contin.  0.6819     2.5521 0.3467   0.2214     Balanced, <0.25      2.1250 0.1667
hum_rig        Contin.  1.0682     1.6682 0.7333   0.3703 Not Balanced, >0.25      1.4461 0.5000
fem_leg        Contin. -0.3897     0.5982 0.1600  -0.6642 Not Balanced, >0.25      0.5312 0.3333
fem_com        Contin.  0.6102     1.1691 0.2667   0.8554 Not Balanced, >0.25      1.4138 0.3333
lef_ide         Binary  0.6223          . 0.2933   0.7071 Not Balanced, >0.25           . 0.3333
start_date_num Contin. -0.2229     0.8287 0.3400   0.0232     Balanced, <0.25      1.0209 0.1667
                     KS.Threshold
con_dur        Not Balanced, >0.1
con_int        Not Balanced, >0.1
con_iss        Not Balanced, >0.1
war_par        Not Balanced, >0.1
gdp            Not Balanced, >0.1
dem            Not Balanced, >0.1
hum_rig        Not Balanced, >0.1
fem_leg        Not Balanced, >0.1
fem_com        Not Balanced, >0.1
lef_ide        Not Balanced, >0.1
start_date_num Not Balanced, >0.1

Balance tally for mean differences
                    count
Balanced, <0.25         3
Not Balanced, >0.25     8

Variable with the greatest mean difference
 Variable Diff.Adj         M.Threshold
  fem_com   0.8554 Not Balanced, >0.25

Balance tally for KS statistics
                   count
Balanced, <0.1         0
Not Balanced, >0.1    11

Variable with the greatest KS statistic
 Variable KS.Adj       KS.Threshold
      gdp    0.5 Not Balanced, >0.1

Effective sample sizes
           Control Treated
Unadjusted      50       6
Adjusted         6       6

library(Matching)

get_balance_table <- function(formula, data, match_out, nboots = 500) {
  
  # 1. Run MatchBalance
  mb <- MatchBalance(formula, data = data, match.out = match_out, nboots = nboots, print.level = 0)
  
  # 2. Define the "Nested" Null Coalescing Operator 
  #    (This allows us to chain checks: "Try A, if NULL try B, if NULL try C")
  `%||%` <- function(a, b) if (is.null(a)) b else a

  # 3. Define the Extractors
  #    T-Test: Simple check
  get_tt <- function(x) x$tt$p.value %||% NA
  
  #    KS-Test: The "Nested" Approach
  #    1. Try nested bootstrap (x$ks$ks.boot.pvalue) <--- PRIORITY
  #    2. Try top-level bootstrap (x$ks.boot.pvalue)
  #    3. Try standard KS (x$ks$p.value)
  get_ks <- function(x) x$ks$ks.boot.pvalue %||% x$ks.boot.pvalue %||% x$ks$p.value %||% NA

  # 4. Extract efficiently with sapply
  idx <- seq_along(mb$AfterMatching)
  vars <- all.vars(formula)[-1]
  final_names <- if(length(vars) == length(idx)) vars else paste0("Var_", idx)

  df <- data.frame(
    Variable  = final_names,
    Pre_T_P   = round(sapply(idx, function(i) get_tt(mb$BeforeMatching[[i]])), 3),
    Post_T_P  = round(sapply(idx, function(i) get_tt(mb$AfterMatching[[i]])), 3),
    Pre_KS_P  = round(sapply(idx, function(i) get_ks(mb$BeforeMatching[[i]])), 3),
    Post_KS_P = round(sapply(idx, function(i) get_ks(mb$AfterMatching[[i]])), 3)
  )
  
  # 5. Status: Pass if (KS is NA or >0.05) AND (T-test > 0.05)
  df$Status <- ifelse((is.na(df$Post_KS_P) | df$Post_KS_P > 0.1) & df$Post_T_P > 0.1, "PASS", "FAIL")
  
  return(df)
}

sensmaker Sensitivity Test for Survival Analysis

set.seed(123)
# --- Call Function ---
balance_table_H1 <- get_balance_table(
  formula = mb_formula_H1, 
  data = df_clean_H1, 
  match_out = m_out_H1, 
  nboots = 500
)

balance_table_H1

# Include outcome variable for effect estimate
Y  <- df_clean_H1$days_of_peace
m_out_H1A <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_H1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Estimate Treatment Effect, SE, and CI ======
summary(m_out_H1A)

Estimate...  2624.2 
AI SE......  1031.8 
T-stat.....  2.5433 
p.val......  0.01098 

Original number of observations..............  56 
Original number of treated obs...............  6 
Matched number of observations...............  6 
Matched number of observations  (unweighted).  6 

Caliper (SDs)........................................   Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf 0.1023685 
Number of obs dropped by 'exact' or 'caliper'  0 

ATT_est <- m_out_H1A$est
ATT_se <- m_out_H1A$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 2)))

[1] "Estimated ATT: 2624.17"

print(paste("Standard Error:", round(ATT_se, 2)))

[1] "Standard Error: 1031.78"

print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))

[1] "95% CI: [ 601.876 , 4646.458 ]"

print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

[1] "99% CI: [ -33.701 , 5282.035 ]"

Hypothesis 2: Women participation and quality of peace agreement (content)

# ============ Remove NAs from ALL necessary variables =============

treatnent_var <- c("fem_sig")
outcome_vars <- c("cpa", "civ_soc")
pretreatment_covariates <- c("con_dur", "con_int", "con_iss", "war_par", "gdp", "dem", "hum_rig", "fem_com", "lef_ide", "fem_leg", "start_date_num") 

cleaning_vars <- c(treatnent_var, outcome_vars, pretreatment_covariates)

# Check which columns contain missing values
colSums(is.na(data[, cleaning_vars])) 

# Clean the data to enable genetic matching (cannot have NAs)
df_clean_H2 <- na.omit(data[, c("cas_id", cleaning_vars)])
df_clean_H2

#Check how many rows have been dropped
cat("Number of control/treatment units in sample (before cleaning):")
table(data$fem_sig)

cat("Number of control/treatment units in sample (after cleaning):")
table(df_clean_H2$fem_sig)

set.seed(123)

# ============== DEFINE MATCHING OBJECTS =============
# Treatment 
Tr <- df_clean_H2$fem_sig

# Matching Matrix (X)
matching_vars <- pretreatment_covariates
X <- df_clean_H2[, matching_vars]
# do not include outcomes

# ============= RUN GENMATCH ===============

genout_H2 <- GenMatch(Tr = Tr, X = X, BalanceMatrix = X, 
                   estimand = "ATT", pop.size = 5000, 
                   replace=FALSE,
                   print.level = 0)

# ============= RUN MATCH ===============

m_out_H2 <- Match(Tr = Tr, X = X, Weight.matrix = genout_H2, 
                  estimand = "ATT",
                  replace=FALSE)


# ===== Report on matching summary ======

# Count how many times each treated unit appears in the match list
matches_per_treated <- table(m_out_H2$index.treated)
# Check the distribution (e.g., how many got 1 match, how many got 2, etc.)
table(matches_per_treated)

# Total treated units in original data
total_treated <- sum(Tr)

# Total treated units that successfully found a match
matched_treated <- length(unique(m_out_H2$index.treated)) #all treated units have been matched

cat("Original Treated:", total_treated, "\n")
cat("Matched Treated: ", matched_treated, "\n")
cat("Dropped:         ", total_treated - matched_treated, "\n")

# ===== Report Balance and Lowest P-Value ======
set.seed(123)
# treatment indicator is the dependent variable, check balance for all covariates
mb_formula_H2 <- fem_sig ~ con_dur + con_int + con_iss +  war_par + gdp + dem + hum_rig + fem_leg + fem_com + lef_ide + start_date_num

MatchBalance(mb_formula_H2, data = df_clean_H2, match.out = m_out_H2, nboots = 500)

#install.packages("cobalt")
library(cobalt)

love.plot(mb_formula_H2, 
          data = df_clean_H2, 
          weights = m_out_H2, # Use the match object as the source of weights
          binary = "std",
          thresholds = c(m = .25),
          alpha = 0.7,
          title = "Covariate Balance (SMD)",
          grid = TRUE)

love.plot(mb_formula_H2, 
          data = df_clean_H2, 
          weights = m_out_H2, # Use the match object as the source of weights
          stats = "ks.statistics", # Change to KS statistics
          binary = "std",
          thresholds = c(ks = 0.1), 
          alpha = 0.7,
          grid = TRUE,
          title = "Covariate Balance (KS Statistics)")

bal.tab(mb_formula_H2, 
        data = df_clean_H2, 
        weights = m_out_H2, # Use the match object as the source of weights
        estimand = "ATT",
        binary = "std",      # <--- Forces binary vars to use SMD
        continuous = "std",  # <--- Forces continuous vars (default) to use SMD
        un = TRUE,
        stats = c("m", "v", "ks"),   # Mean diffs, Variance ratios, KS stats
        thresholds = c(m = 0.25, ks = 0.1))

set.seed(123)
# --- Call Function ---
balance_table_H2 <- get_balance_table(
  formula = mb_formula_H2, 
  data = df_clean_H2, 
  match_out = m_out_H2, 
  nboots = 500
)

balance_table_H2

Hypothesis 2A: fem_sig on cpa

set.seed(123)

# ========== Run Match() with outcome variable =======

Y  <- df_clean_H2$cpa # Outcome
m_out_H2A <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_H2, estimand = "ATT")

# ===== Estimate Treatment Effect, SE, and CI ======
summary(m_out_H2A)

ATT_est <- m_out_H2A$est
ATT_se <- m_out_H2A$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 4)))
print(paste("Standard Error:", round(ATT_se, 4)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

Hypothesis 2B: fem_sig on civ_soc

set.seed(123)
# ========== Run Match() with outcome variable =======

Y  <- df_clean_H2$civ_soc # Outcome
m_out_H2B <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_H2, estimand = "ATT")


# ===== Estimate Treatment Effect, SE, and CI ======
summary(m_out_H2B)

ATT_est <- m_out_H2B$est
ATT_se <- m_out_H2B$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 2)))
print(paste("Standard Error:", round(ATT_se, 2)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

Model H2: Rosenbaum Bounds Sensitivity Analysis

# Outcome
Y  <- df_clean_H1$pc_nd
m_out_H1B <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_H1, 
               caliper = caliper_vec, 
               estimand = "ATT")

# ===== Estimate Treatment Effect, SE, and CI ======

summary(m_out_H1B)

Estimate...  -0.5 
AI SE......  0.20412 
T-stat.....  -2.4495 
p.val......  0.014306 

Original number of observations..............  56 
Original number of treated obs...............  6 
Matched number of observations...............  6 
Matched number of observations  (unweighted).  6 

Caliper (SDs)........................................   Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf 0.1023685 
Number of obs dropped by 'exact' or 'caliper'  0 

ATT_est <- m_out_H1B$est
ATT_se <- m_out_H1B$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 2)))

[1] "Estimated ATT: -0.5"

print(paste("Standard Error:", round(ATT_se, 2)))

[1] "Standard Error: 0.2"

print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))

[1] "95% CI: [ -0.9 , -0.1 ]"

print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

[1] "99% CI: [ -1.026 , 0.026 ]"

# Run Sensitivity Analysis on Durability (Days of Peace)
#install.packages("rbounds")
library(rbounds)

# Obtain indexes for matched treatment and control units
treated_indices <- m_out_H1$index.treated
control_indices <- m_out_H1$index.control

# Extract 'Days of Peace' for both groups
days_Tr <- as.numeric(df_clean_H1$days_of_peace[m_out_H1$index.treated])
days_Co <- as.numeric(df_clean_H1$days_of_peace[m_out_H1$index.control])

# P-Value
psens_H1A <- psens(x = days_Tr, 
                   y = days_Co, 
                   Gamma = 3, # We check up to a "tripling" of the odds of hidden bias
                   GammaInc = 0.5) # We move in steps of 0.5
# The probability of assignment to treatment under hidden bias Gamma
# For Upper Bound (worst case), we assume bias helps the observed effect

# Hodges-Lehmann Estimate (non-parametric)
hlsens_H1A <- hlsens(x = days_Tr, 
                   y = days_Co, 
                   pr = 0.5,
                   Gamma = 10, 
                   GammaInc = 1) 
psens_H1A

 Rosenbaum Sensitivity Test for Wilcoxon Signed Rank P-Value 
 
Unconfounded estimate ....  0.0579 

 Note: Gamma is Odds of Differential Assignment To
 Treatment Due to Unobserved Factors 
 

# result is only "marginally significant" (p=0.06) assuming no hidden bias
# statistical significance disappears with even a tiny amount of hidden bias
hlsens_H1A

 Rosenbaum Sensitivity Test for Hodges-Lehmann Point Estimate 
 
Unconfounded estimate ....  2661 

 Note: Gamma is Odds of Differential Assignment To
 Treatment Due to Unobserved Factors 
 

# magnitude of the effect point estimate is quite robust
#it would take a massive amount of unobserved bias (Gamma = 7) to reduce the estimated effect to zero
# even if unobserved factors caused massive selection bias, the underlying signal still suggests a positive effect on duration

Robustness Check

# ============ Remove NAs from ALL necessary variables =============

# Check which columns contain missing values
colSums(is.na(data[, cleaning_vars])) 

# Clean the data to enable genetic matching (cannot have NAs)
pretreatment_covariates <- c("con_dur", "con_int", "con_iss", "war_par", "gdp", "dem", "hum_rig", "fem_com", "lef_ide", "fem_leg", "start_date_num")
treatnent_var <- c("fem_sig")
outcome_vars <- c("cpa", "civ_soc")

add_covariates <- c("pow_sha", "quo", "un_pea_con") # additional covariates

# This list includes Y, Tr, and X. We use it ONLY for na.omit.
cleaning_vars <- c(treatnent_var, outcome_vars, pretreatment_covariates, add_covariates)

df_clean_rb2 <- na.omit(data[, c("cas_id", cleaning_vars)])
df_clean_rb2

# Check how many rows have been dropped
cat("Number of control/treatment units in sample (before cleaning):")
table(data$fem_sig)
cat("Number of control/treatment units in sample (after cleaning):")
table(df_clean_rb2$fem_sig)

# ========== DEFINE MATCHING OBJECTS ===============

set.seed(123)

# Treatment 
Tr <- df_clean_rb2$fem_sig

# Matching Matrix (X)
matching_vars <- c(pretreatment_covariates, add_covariates) # include additional covariates
X <- df_clean_rb2[, matching_vars]

# =========== RUN GENMATCH =================

genout_rb2 <- GenMatch(Tr = Tr, X = X, BalanceMatrix = X,
                   replace=FALSE, 
                   estimand = "ATT", pop.size = 5000, 
                   print.level = 0)

# =============== RUN MATCH ===============

m_out_rb2 <- Match(Tr = Tr, X = X, Weight.matrix = genout_rb2, estimand = "ATT", replace=FALSE)


# ========== Report on matching summary =========
# Count how many times each treated unit appears in the match list
matches_per_treated <- table(m_out_rb2 $index.treated)
# Check the distribution (e.g., how many got 1 match, how many got 2, etc.)
table(matches_per_treated)

# Total treated units in original data
total_treated <- sum(Tr)

# Total treated units that successfully found a match
matched_treated <- length(unique(m_out_rb2 $index.treated)) #all treated units have been matched

cat("Original Treated:", total_treated, "\n")
cat("Matched Treated: ", matched_treated, "\n")
cat("Dropped:         ", total_treated - matched_treated, "\n")

set.seed(123)
# ===== Report Balance and Lowest P-Value ======

# treatment indicator is the dependent variable, check balance for all covariates
mb_formula_rb2 <- fem_sig ~ con_dur + con_int + con_iss +  war_par + gdp + dem + hum_rig + fem_leg + fem_com + lef_ide + pow_sha + quo + un_pea_con + start_date_num

MatchBalance(mb_formula_rb2, data = df_clean_rb2, match.out = m_out_rb2, nboots = 500)

H2 vs R2: Comparison of Matching Results

bal.tab(mb_formula_rb2, 
        data = df_clean_rb2, 
        weights = m_out_rb2, # Use the match object as the source of weights
        estimand = "ATT",
        binary = "std",      # <--- Forces binary vars to use SMD
        continuous = "std",  # <--- Forces continuous vars (default) to use SMD
        un = TRUE,
        stats = c("m", "v", "ks"),   # Mean diffs, Variance ratios, KS stats
        thresholds = c(m = 0.25, ks = 0.1))

set.seed(123)
# --- Call Function ---
balance_table_rb2 <- get_balance_table(
  formula = mb_formula_rb2, 
  data = df_clean_rb2, 
  match_out = m_out_rb2, 
  nboots = 500
)

balance_table_rb2

matched_data_H2 <- rbind(
  df_clean_H2[m_out_H2$index.treated, ],
  df_clean_H2[m_out_H2$index.control, ]
)
n_pairs <- length(m_out_H2$index.treated)
matched_data_H2$pair_id <- c(1:n_pairs, 1:n_pairs)

matched_data_rb2 <- rbind(
  df_clean_rb2[m_out_rb2$index.treated, ],
  df_clean_rb2[m_out_rb2$index.control, ]
)
n_pairs <- length(m_out_rb2$index.treated)
matched_data_rb2$pair_id <- c(1:n_pairs, 1:n_pairs)

matches <- cbind(Pair_ID = matched_data_H2$pair_id, 
                 Treated = matched_data_H2$cas_id[matched_data_H2$fem_sig == 1], 
                 ModelH2_Matches = matched_data_H2$cas_id[matched_data_H2$fem_sig == 0], 
                 ModelR2_Matches = matched_data_rb2$cas_id[matched_data_rb2$fem_sig == 0])
matches[1:6, ]

# Compare outcome observations of matched units
matchcomp_cpa <- cbind(Pair_ID = matched_data_H2$pair_id,
                        cas_id = matched_data_H2$cas_id,
                        fem_sig = matched_data_H2$fem_sig,
                        cpa_H2 = matched_data_H2$cpa,
                        cpa_rb2 = matched_data_rb2$cpa)
matchcomp_civ_soc <- cbind(Pair_ID = matched_data_H2$pair_id,
                      cas_id = matched_data_H2$cas_id,
                      fem_sig = matched_data_H2$fem_sig,
                      civ_soc_H2 = matched_data_H2$civ_soc,
                      civ_soc_rb2 = matched_data_rb2$civ_soc)
matchcomp_cpa
matchcomp_civ_soc

H2 vs R2: Comparison of Treatment Effect Estimates

# ========== Run Match() with outome variable (cpa) ==============
set.seed(123)
# Outcome
Y  <- df_clean_rb2$cpa
m_out_rb2A <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_rb2, estimand = "ATT", replace = FALSE)

# ===== Estimate Treatment Effect, SE, and CI ======

summary(m_out_rb2A)

ATT_est <- m_out_rb2A$est
ATT_se <- 0.20412#m_out_rb2A$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 2)))
print(paste("Standard Error:", round(ATT_se, 2)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

set.seed(123)
# ========== Run Match() with outome variable (civ_soc) ==============
# Outcome
Y  <- df_clean_rb2$civ_soc
m_out_rb2B <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = genout_rb2, estimand = "ATT", replace = FALSE)

# ===== Estimate Treatment Effect, SE, and CI ======
summary(m_out_rb2B)

ATT_est <- m_out_rb2B$est
ATT_se <- 0.15215 #m_out_rb2B$se

#----- Get Confidence Intervals ------
# CI bounds given by point estimate +/- margin of error
# CV for 95% interval is 1.96
ci_lower95 <- ATT_est - 1.96 * ATT_se
ci_upper95 <- ATT_est + 1.96 * ATT_se
# CV for 99% interval is 2.576
ci_lower99 <- ATT_est - 2.576 * ATT_se
ci_upper99 <- ATT_est + 2.576 * ATT_se

print(paste("Estimated ATT:", round(ATT_est, 3)))
print(paste("Standard Error:", round(ATT_se, 3)))
print(paste("95% CI: [", round(ci_lower95, 3), ",", round(ci_upper95, 3), "]"))
print(paste("99% CI: [", round(ci_lower99, 3), ",", round(ci_upper99, 3), "]"))

print ("H2A vs. R2A")
summary(m_out_H2A)
summary(m_out_rb2A)
print ("H2B vs. R2B")
summary(m_out_H2B)
summary(m_out_rb2B)