library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(MatchIt)
library(fixest)
library(rollmatch)
library(tidyr)

Define the Simulation Parameters

# Set simulation parameters
num_simulations <- 100
N <- 1000  # Number of units (projects)
T <- 24    # Number of time periods (months)
true_beta <- 2  # True treatment effect
t_start_treatment <- 6  # No units receive treatment before this time

# Storage for results
estimates_psm <- numeric(num_simulations)
estimates_rem <- numeric(num_simulations)

Loop Over Simulations

for (sim in 1:num_simulations) {
  # Set seed for reproducibility
  set.seed(sim)
  
  # Generate unit fixed effects (alpha_i)
  alpha_i <- rnorm(N, mean = 0, sd = 1)
  
  # Generate time fixed effects (lambda_t)
  lambda_t <- rnorm(T, mean = 0, sd = 1)
  
  # Initialize matrices
  X_it <- matrix(0, nrow = N, ncol = T)
  Treatment_it <- matrix(0, nrow = N, ncol = T)
  Y_it <- matrix(0, nrow = N, ncol = T)
  
  # Parameters for the DGP
  intercept <- -7
  strong_coefficient <- 5.0
  strong_outcome_coefficient <- 10.0
  
  # Generate X_it as a non-stationary process (random walk)
  for (t in 1:T) {
    if (t == 1) {
      X_it[, t] <- rnorm(N, mean = 0, sd = 1)
    } else {
      # Random walk for X_it
      X_it[, t] <- X_it[, t - 1] + rnorm(N, mean = 0, sd = 1)
    }
  }
  
  # Initialize Treatment_it for t = 1
  Treatment_it[, 1] <- 0  # No units are treated at time t = 1
  
  for (t in 1:(T - 1)) {
    if (t >= t_start_treatment) {
      # Compute probability of treatment for time t + 1 based on X_it at time t
      logit_p <- intercept + strong_coefficient * X_it[, t]
      p_treatment <- exp(logit_p) / (1 + exp(logit_p))
      
      # Constrain probabilities between 0.001 and 0.8 to prevent too many units from being treated
      # p_treatment <- pmin(pmax(p_treatment, 0.001), 0.999)
      
      # Identify units not yet treated at time t
      untreated_units <- which(Treatment_it[, t] == 0)
      
      # Assign treatment status for time t + 1
      Treatment_it[untreated_units, t + 1] <- rbinom(length(untreated_units), 1, p_treatment[untreated_units])
      
      # Ensure "once treated, always treated"
      Treatment_it[, t + 1] <- pmax(Treatment_it[, t + 1], Treatment_it[, t])
    } else {
      # For t < t_start_treatment, carry over treatment status
      Treatment_it[, t + 1] <- Treatment_it[, t]
    }
    
    # Generate outcome at time t based on treatment status at time t
    epsilon_it <- rnorm(N, mean = 0, sd = 1)
    Y_it[, t] <- true_beta * Treatment_it[, t] + strong_outcome_coefficient * X_it[, t] + alpha_i + lambda_t[t] + epsilon_it
  }
  
  # For the last period T, generate the outcome
  epsilon_it <- rnorm(N, mean = 0, sd = 1)
  Y_it[, T] <- true_beta * Treatment_it[, T] + strong_outcome_coefficient * X_it[, T] + alpha_i + lambda_t[T] + epsilon_it
  
  # Reshape Data into Long Format
  data_list <- list()
  for (t in 1:T) {
    data_t <- data.frame(
      unit = 1:N,
      time = t,
      X = X_it[, t],
      Y = Y_it[, t],
      Treatment = Treatment_it[, t],
      alpha_i = alpha_i,
      lambda_t = lambda_t[t]
    )
    data_list[[t]] <- data_t
  }
  
  data_long <- do.call(rbind, data_list)
  
  # Identify never-treated units
  never_treated_units <- which(rowSums(Treatment_it) == 0)
  num_never_treated <- length(never_treated_units)
  prop_never_treated <- num_never_treated / N
  
  cat("Number of never-treated units:", num_never_treated, "\n")
  cat("Proportion of never-treated units:", prop_never_treated, "\n")
  
  # Apply Traditional PSM with TWFE DiD
  t_fixed <- t_start_treatment - 1  # Fixed period before any treatments occur

  # Data at fixed time
  data_fixed <- data_long %>% filter(time == t_fixed)

  # Create future_treatment variable by checking if unit is treated after t_fixed
  future_treatment_df <- data_long %>%
    filter(time >= t_start_treatment) %>%
    group_by(unit) %>%
    summarise(future_treatment = as.integer(any(Treatment == 1))) %>%
    ungroup()

  # Merge future_treatment back into data_fixed
  data_fixed <- data_fixed %>%
    left_join(future_treatment_df, by = "unit") %>%
    mutate(future_treatment = ifelse(is.na(future_treatment), 0, future_treatment))

  # Check if there is variation in future_treatment
  if (length(unique(data_fixed$future_treatment)) < 2) {
    warning("Not enough variation in future_treatment in simulation ", sim)
    next  # Skip to next simulation
  }

  # Subset treated and control units
  treated_units <- data_fixed %>% filter(future_treatment == 1)
  control_units <- data_fixed %>% filter(future_treatment == 0)

  # Combine treated and control units
  data_psm <- bind_rows(treated_units, control_units)

  # Estimate propensity scores using X at t_fixed
  formula_psm <- as.formula("future_treatment ~ X")

  psm_model <- matchit(formula_psm, data = data_psm, method = "nearest", distance = "logit")

  # Check if matching was successful
  matched_data <- match.data(psm_model)
  if (nrow(matched_data) == 0) {
    warning("No matches found in PSM for simulation ", sim)
    next  # Skip to next simulation
  }

  # Get matched unit IDs
  matched_units <- unique(matched_data$unit)
  print(paste("PSM, Number of matched pairs: ", length(matched_units)/2 ))

  # Subset the full data_long for these units
  data_matched <- data_long %>% filter(unit %in% matched_units)

  # Estimate TWFE DiD model
  model_psm <- feols(Y ~ Treatment + X | unit + time, data = data_matched)

  # Extract treatment effect estimate
  treatment_effect_psm <- coef(model_psm)["Treatment"]

  estimates_psm[sim] <- treatment_effect_psm
  print(paste("PSM estimated coefficient:", treatment_effect_psm))
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  # 
  # Prepare the data for rollmatch
  data_long <- data_long %>%
    rename(
      indiv_id = unit,
      quarter = time,
      treat_status = Treatment  # Rename to avoid confusion
    ) %>%
    arrange(indiv_id, quarter) %>%
    group_by(indiv_id) %>%
    mutate(
      # Identify the entry quarter when treatment starts
      entry_q = ifelse(treat_status == 1 & lag(treat_status, default = 0) == 0, quarter, NA),
      # Determine if the unit ever receives treatment
      ever_treated = ifelse(any(treat_status == 1), 1, 0)
    ) %>%
    # Fill entry_q for all periods within each unit
    mutate(
      entry_q = first(entry_q[!is.na(entry_q)], default = NA_real_)
    ) %>%
    ungroup()

  # Determine the maximum time period
  max_time <- max(data_long$quarter)

  # Assign entry_q to control units and ensure it's integer
  data_long <- data_long %>%
    mutate(
      entry_q = ifelse(is.na(entry_q), max_time + 1, entry_q),
      entry_q = as.integer(entry_q)
    )

  # Check for missing entry_q values
  if (any(is.na(data_long$entry_q))) {
    warning("Missing entry_q values in simulation ", sim)
    next  # Skip to next simulation
  }

  # Set lookback period (e.g., 1 period)
  lookback_period <- 1

  reduced_data <- reduce_data(
    data = data_long,
    treat = "ever_treated",
    tm = "quarter",
    entry = "entry_q",
    id = "indiv_id",
    lookback = lookback_period
  )
  
  # Check if reduced_data is empty
  if (nrow(reduced_data) == 0) {
    warning("Reduced data is empty in simulation ", sim)
    next  # Skip to next simulation
  }
  
  # Define formula for propensity score model
  fm <- as.formula("ever_treated ~ X")
  vars <- all.vars(fm)
  
  scored_data <- score_data(
    reduced_data = reduced_data,
    model_type = "logistic",
    match_on = "logit",
    fm = fm,
    treat = "ever_treated",
    tm = "quarter",
    entry = "entry_q",
    id = "indiv_id"
  )
  
  # Perform rollmatch
  output <- rollmatch(
    scored_data,
    data = data_long,
    treat = "ever_treated",
    tm = "quarter",
    entry = "entry_q",
    id = "indiv_id",
    vars = vars,
    lookback = lookback_period,
    alpha = 0,            # Adjust as needed
    standard_deviation = "average",
    num_matches = 1,       # One-to-one matching
    replacement = FALSE    # Without replacement
  )
  
  # Extract matched pairs
  matched_pairs <- output$matched_data
  
  # Check if matches were found
  if (is.null(matched_pairs) || nrow(matched_pairs) == 0) {
    warning("No matches found in REM for simulation ", sim)
    next  # Skip to next simulation
  } else {
    print(paste("REM, Number of matched pairs:", nrow(matched_pairs)))
  }
  
  # Get list of matched units
  matched_units_rem <- unique(c(matched_pairs$treat_id, matched_pairs$control_id))
  
  # Subset data_long for these units
  data_matched_rem <- data_long %>% filter(indiv_id %in% matched_units_rem)
  
  # Estimate TWFE DiD model
  model_rem <- feols(Y ~ treat_status + X | indiv_id + quarter, data = data_matched_rem)
  
  # Extract treatment effect estimate
  treatment_effect_rem <- coef(model_rem)["treat_status"]
  
  estimates_rem[sim] <- treatment_effect_rem
  print(paste("REM estimated coefficient:", treatment_effect_rem))
}
## Number of never-treated units: 348 
## Proportion of never-treated units: 0.348
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  348"
## [1] "PSM estimated coefficient: 1.99751707417471"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: `group_by_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `group_by()` instead.
## ℹ See vignette('programming') for more help
## ℹ The deprecated feature was likely used in the rollmatch package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [1] "REM, Number of matched pairs: 348"
## [1] "REM estimated coefficient: 2.02626934388433"
## Number of never-treated units: 348 
## Proportion of never-treated units: 0.348
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  348"
## [1] "PSM estimated coefficient: 2.01420553043874"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 348"
## [1] "REM estimated coefficient: 2.01551282519787"
## Number of never-treated units: 379 
## Proportion of never-treated units: 0.379
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  379"
## [1] "PSM estimated coefficient: 1.95509730816307"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 379"
## [1] "REM estimated coefficient: 1.9837269110692"
## Number of never-treated units: 374 
## Proportion of never-treated units: 0.374
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  374"
## [1] "PSM estimated coefficient: 2.00299759631667"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 374"
## [1] "REM estimated coefficient: 1.98931563001385"
## Number of never-treated units: 348 
## Proportion of never-treated units: 0.348
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  348"
## [1] "PSM estimated coefficient: 1.92699086314758"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 348"
## [1] "REM estimated coefficient: 1.9745875740136"
## Number of never-treated units: 331 
## Proportion of never-treated units: 0.331
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  331"
## [1] "PSM estimated coefficient: 2.03895603000338"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 331"
## [1] "REM estimated coefficient: 2.00526115560238"
## Number of never-treated units: 346 
## Proportion of never-treated units: 0.346
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  346"
## [1] "PSM estimated coefficient: 1.92579480572275"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 346"
## [1] "REM estimated coefficient: 1.93624882929778"
## Number of never-treated units: 347 
## Proportion of never-treated units: 0.347
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  347"
## [1] "PSM estimated coefficient: 1.99229671510224"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 347"
## [1] "REM estimated coefficient: 1.9810316441063"
## Number of never-treated units: 351 
## Proportion of never-treated units: 0.351
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  351"
## [1] "PSM estimated coefficient: 1.97160807892128"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 351"
## [1] "REM estimated coefficient: 2.00630929134979"
## Number of never-treated units: 365 
## Proportion of never-treated units: 0.365
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  365"
## [1] "PSM estimated coefficient: 2.04303311491245"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 365"
## [1] "REM estimated coefficient: 2.07051173582406"
## Number of never-treated units: 358 
## Proportion of never-treated units: 0.358
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  358"
## [1] "PSM estimated coefficient: 2.02272696004372"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 358"
## [1] "REM estimated coefficient: 2.00794146509412"
## Number of never-treated units: 345 
## Proportion of never-treated units: 0.345
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  345"
## [1] "PSM estimated coefficient: 2.00080728943662"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 345"
## [1] "REM estimated coefficient: 1.95804046802234"
## Number of never-treated units: 379 
## Proportion of never-treated units: 0.379
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  379"
## [1] "PSM estimated coefficient: 1.96049497101514"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 379"
## [1] "REM estimated coefficient: 1.99582021936134"
## Number of never-treated units: 352 
## Proportion of never-treated units: 0.352
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  352"
## [1] "PSM estimated coefficient: 2.0301808245299"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 352"
## [1] "REM estimated coefficient: 2.04481986163211"
## Number of never-treated units: 373 
## Proportion of never-treated units: 0.373
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  373"
## [1] "PSM estimated coefficient: 1.99713049370312"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 373"
## [1] "REM estimated coefficient: 1.97144480251889"
## Number of never-treated units: 358 
## Proportion of never-treated units: 0.358
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  358"
## [1] "PSM estimated coefficient: 2.05181369003828"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 358"
## [1] "REM estimated coefficient: 2.01872668750583"
## Number of never-treated units: 329 
## Proportion of never-treated units: 0.329
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  329"
## [1] "PSM estimated coefficient: 2.01913837290052"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 329"
## [1] "REM estimated coefficient: 1.95183885192952"
## Number of never-treated units: 343 
## Proportion of never-treated units: 0.343
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  343"
## [1] "PSM estimated coefficient: 2.0085775345575"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 343"
## [1] "REM estimated coefficient: 2.03426834319808"
## Number of never-treated units: 358 
## Proportion of never-treated units: 0.358
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  358"
## [1] "PSM estimated coefficient: 1.94950052895801"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 358"
## [1] "REM estimated coefficient: 1.97891571418355"
## Number of never-treated units: 355 
## Proportion of never-treated units: 0.355
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  355"
## [1] "PSM estimated coefficient: 2.0302314969296"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 355"
## [1] "REM estimated coefficient: 1.99247063589155"
## Number of never-treated units: 355 
## Proportion of never-treated units: 0.355
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  355"
## [1] "PSM estimated coefficient: 2.03287086242435"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 355"
## [1] "REM estimated coefficient: 2.00391950672673"
## Number of never-treated units: 359 
## Proportion of never-treated units: 0.359
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  359"
## [1] "PSM estimated coefficient: 1.97424250875686"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 359"
## [1] "REM estimated coefficient: 1.99421185735658"
## Number of never-treated units: 353 
## Proportion of never-treated units: 0.353
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  353"
## [1] "PSM estimated coefficient: 2.04623022071966"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 353"
## [1] "REM estimated coefficient: 2.01925548452993"
## Number of never-treated units: 343 
## Proportion of never-treated units: 0.343
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  343"
## [1] "PSM estimated coefficient: 1.969098314061"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 343"
## [1] "REM estimated coefficient: 1.94303769748043"
## Number of never-treated units: 363 
## Proportion of never-treated units: 0.363
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  363"
## [1] "PSM estimated coefficient: 2.02167693543782"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 363"
## [1] "REM estimated coefficient: 2.03078466515209"
## Number of never-treated units: 355 
## Proportion of never-treated units: 0.355
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  355"
## [1] "PSM estimated coefficient: 1.95647143030997"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 355"
## [1] "REM estimated coefficient: 2.00121266352379"
## Number of never-treated units: 400 
## Proportion of never-treated units: 0.4
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  400"
## [1] "PSM estimated coefficient: 2.04449344708229"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 400"
## [1] "REM estimated coefficient: 2.03786240104665"
## Number of never-treated units: 331 
## Proportion of never-treated units: 0.331
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  331"
## [1] "PSM estimated coefficient: 2.03987213185687"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 331"
## [1] "REM estimated coefficient: 2.00795408121672"
## Number of never-treated units: 372 
## Proportion of never-treated units: 0.372
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  372"
## [1] "PSM estimated coefficient: 2.05127410979254"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 372"
## [1] "REM estimated coefficient: 2.07129282613197"
## Number of never-treated units: 348 
## Proportion of never-treated units: 0.348
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  348"
## [1] "PSM estimated coefficient: 1.96120774633278"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 348"
## [1] "REM estimated coefficient: 2.0297421573697"
## Number of never-treated units: 370 
## Proportion of never-treated units: 0.37
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  370"
## [1] "PSM estimated coefficient: 2.0291659825013"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 370"
## [1] "REM estimated coefficient: 1.98377330430819"
## Number of never-treated units: 312 
## Proportion of never-treated units: 0.312
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  312"
## [1] "PSM estimated coefficient: 2.02201496382493"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 312"
## [1] "REM estimated coefficient: 1.99353746877186"
## Number of never-treated units: 359 
## Proportion of never-treated units: 0.359
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  359"
## [1] "PSM estimated coefficient: 2.02342144066708"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 359"
## [1] "REM estimated coefficient: 2.06978937218258"
## Number of never-treated units: 348 
## Proportion of never-treated units: 0.348
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  348"
## [1] "PSM estimated coefficient: 2.01082619784776"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 348"
## [1] "REM estimated coefficient: 1.99454192021336"
## Number of never-treated units: 332 
## Proportion of never-treated units: 0.332
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  332"
## [1] "PSM estimated coefficient: 2.04306169582328"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 332"
## [1] "REM estimated coefficient: 1.99686239693622"
## Number of never-treated units: 350 
## Proportion of never-treated units: 0.35
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  350"
## [1] "PSM estimated coefficient: 2.06287432001322"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 350"
## [1] "REM estimated coefficient: 2.00641242454099"
## Number of never-treated units: 342 
## Proportion of never-treated units: 0.342
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  342"
## [1] "PSM estimated coefficient: 2.05896564898554"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 342"
## [1] "REM estimated coefficient: 2.04790131557654"
## Number of never-treated units: 351 
## Proportion of never-treated units: 0.351
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  351"
## [1] "PSM estimated coefficient: 2.01755053163174"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 351"
## [1] "REM estimated coefficient: 1.93742361269572"
## Number of never-treated units: 353 
## Proportion of never-treated units: 0.353
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  353"
## [1] "PSM estimated coefficient: 1.96506370281556"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 353"
## [1] "REM estimated coefficient: 1.99395115936389"
## Number of never-treated units: 344 
## Proportion of never-treated units: 0.344
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  344"
## [1] "PSM estimated coefficient: 2.00420857830645"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 344"
## [1] "REM estimated coefficient: 1.99934754670728"
## Number of never-treated units: 355 
## Proportion of never-treated units: 0.355
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  355"
## [1] "PSM estimated coefficient: 1.96804290061315"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 355"
## [1] "REM estimated coefficient: 1.99513560524944"
## Number of never-treated units: 367 
## Proportion of never-treated units: 0.367
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  367"
## [1] "PSM estimated coefficient: 2.05914697917639"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 367"
## [1] "REM estimated coefficient: 2.01425204845393"
## Number of never-treated units: 360 
## Proportion of never-treated units: 0.36
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  360"
## [1] "PSM estimated coefficient: 2.01622096823627"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 360"
## [1] "REM estimated coefficient: 2.0084853425474"
## Number of never-treated units: 353 
## Proportion of never-treated units: 0.353
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  353"
## [1] "PSM estimated coefficient: 2.05364403877027"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 353"
## [1] "REM estimated coefficient: 2.03516960530752"
## Number of never-treated units: 348 
## Proportion of never-treated units: 0.348
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  348"
## [1] "PSM estimated coefficient: 1.99604766154266"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 348"
## [1] "REM estimated coefficient: 2.03051192563419"
## Number of never-treated units: 337 
## Proportion of never-treated units: 0.337
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  337"
## [1] "PSM estimated coefficient: 2.07578332002677"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 337"
## [1] "REM estimated coefficient: 1.99014362042466"
## Number of never-treated units: 374 
## Proportion of never-treated units: 0.374
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  374"
## [1] "PSM estimated coefficient: 1.97891941564136"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 374"
## [1] "REM estimated coefficient: 1.99279493553792"
## Number of never-treated units: 347 
## Proportion of never-treated units: 0.347
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  347"
## [1] "PSM estimated coefficient: 2.02079328036034"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 347"
## [1] "REM estimated coefficient: 2.03894176051121"
## Number of never-treated units: 354 
## Proportion of never-treated units: 0.354
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  354"
## [1] "PSM estimated coefficient: 2.03041445799021"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 354"
## [1] "REM estimated coefficient: 1.97493404161447"
## Number of never-treated units: 336 
## Proportion of never-treated units: 0.336
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  336"
## [1] "PSM estimated coefficient: 2.01725229239022"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 336"
## [1] "REM estimated coefficient: 2.00388254037402"
## Number of never-treated units: 379 
## Proportion of never-treated units: 0.379
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  379"
## [1] "PSM estimated coefficient: 2.04732102917368"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 379"
## [1] "REM estimated coefficient: 2.04239258245831"
## Number of never-treated units: 355 
## Proportion of never-treated units: 0.355
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  355"
## [1] "PSM estimated coefficient: 2.06680158917762"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 355"
## [1] "REM estimated coefficient: 2.04986403056073"
## Number of never-treated units: 363 
## Proportion of never-treated units: 0.363
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  363"
## [1] "PSM estimated coefficient: 1.96629400891326"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 363"
## [1] "REM estimated coefficient: 2.00323439743502"
## Number of never-treated units: 366 
## Proportion of never-treated units: 0.366
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  366"
## [1] "PSM estimated coefficient: 2.01535759566935"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 366"
## [1] "REM estimated coefficient: 2.01008705238611"
## Number of never-treated units: 370 
## Proportion of never-treated units: 0.37
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  370"
## [1] "PSM estimated coefficient: 2.01620171361813"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 370"
## [1] "REM estimated coefficient: 2.03386582435213"
## Number of never-treated units: 345 
## Proportion of never-treated units: 0.345
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  345"
## [1] "PSM estimated coefficient: 2.02188762542769"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 345"
## [1] "REM estimated coefficient: 1.98258261243481"
## Number of never-treated units: 366 
## Proportion of never-treated units: 0.366
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  366"
## [1] "PSM estimated coefficient: 2.00115143509135"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 366"
## [1] "REM estimated coefficient: 1.98075445477594"
## Number of never-treated units: 361 
## Proportion of never-treated units: 0.361
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  361"
## [1] "PSM estimated coefficient: 2.04730736910322"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 361"
## [1] "REM estimated coefficient: 2.00036812399465"
## Number of never-treated units: 340 
## Proportion of never-treated units: 0.34
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  340"
## [1] "PSM estimated coefficient: 1.97033699739703"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 340"
## [1] "REM estimated coefficient: 1.97344045564758"
## Number of never-treated units: 340 
## Proportion of never-treated units: 0.34
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  340"
## [1] "PSM estimated coefficient: 2.11627032480757"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 340"
## [1] "REM estimated coefficient: 2.02493227987766"
## Number of never-treated units: 368 
## Proportion of never-treated units: 0.368
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  368"
## [1] "PSM estimated coefficient: 2.039120783796"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 368"
## [1] "REM estimated coefficient: 1.97676027402399"
## Number of never-treated units: 363 
## Proportion of never-treated units: 0.363
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  363"
## [1] "PSM estimated coefficient: 1.99620407768397"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 363"
## [1] "REM estimated coefficient: 2.03423200827351"
## Number of never-treated units: 364 
## Proportion of never-treated units: 0.364
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  364"
## [1] "PSM estimated coefficient: 2.033396910359"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 364"
## [1] "REM estimated coefficient: 2.00799529435587"
## Number of never-treated units: 362 
## Proportion of never-treated units: 0.362
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  362"
## [1] "PSM estimated coefficient: 1.97885535282717"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 362"
## [1] "REM estimated coefficient: 1.98060432785217"
## Number of never-treated units: 346 
## Proportion of never-treated units: 0.346
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  346"
## [1] "PSM estimated coefficient: 2.02318886393617"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 346"
## [1] "REM estimated coefficient: 2.01918611845527"
## Number of never-treated units: 375 
## Proportion of never-treated units: 0.375
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  375"
## [1] "PSM estimated coefficient: 1.95209905348364"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 375"
## [1] "REM estimated coefficient: 1.95054056485495"
## Number of never-treated units: 352 
## Proportion of never-treated units: 0.352
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  352"
## [1] "PSM estimated coefficient: 1.99125234171698"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 352"
## [1] "REM estimated coefficient: 2.02066017706775"
## Number of never-treated units: 368 
## Proportion of never-treated units: 0.368
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  368"
## [1] "PSM estimated coefficient: 1.93587197796197"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 368"
## [1] "REM estimated coefficient: 1.98649188965585"
## Number of never-treated units: 372 
## Proportion of never-treated units: 0.372
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  372"
## [1] "PSM estimated coefficient: 2.1077151372635"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 372"
## [1] "REM estimated coefficient: 2.08369227803152"
## Number of never-treated units: 334 
## Proportion of never-treated units: 0.334
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  334"
## [1] "PSM estimated coefficient: 2.02795845602078"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 334"
## [1] "REM estimated coefficient: 2.04942073463064"
## Number of never-treated units: 333 
## Proportion of never-treated units: 0.333
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  333"
## [1] "PSM estimated coefficient: 2.02273134297677"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 333"
## [1] "REM estimated coefficient: 2.01108427413985"
## Number of never-treated units: 391 
## Proportion of never-treated units: 0.391
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  391"
## [1] "PSM estimated coefficient: 2.04937126894637"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 391"
## [1] "REM estimated coefficient: 2.03405286094833"
## Number of never-treated units: 361 
## Proportion of never-treated units: 0.361
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  361"
## [1] "PSM estimated coefficient: 1.99155157050949"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 361"
## [1] "REM estimated coefficient: 1.94359547183308"
## Number of never-treated units: 348 
## Proportion of never-treated units: 0.348
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  348"
## [1] "PSM estimated coefficient: 2.00252740189863"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 348"
## [1] "REM estimated coefficient: 1.99988055357196"
## Number of never-treated units: 376 
## Proportion of never-treated units: 0.376
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  376"
## [1] "PSM estimated coefficient: 2.01127312301118"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 376"
## [1] "REM estimated coefficient: 2.02939726116345"
## Number of never-treated units: 346 
## Proportion of never-treated units: 0.346
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  346"
## [1] "PSM estimated coefficient: 2.08393597431658"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 346"
## [1] "REM estimated coefficient: 2.02426827324982"
## Number of never-treated units: 355 
## Proportion of never-treated units: 0.355
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  355"
## [1] "PSM estimated coefficient: 1.97123671458737"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 355"
## [1] "REM estimated coefficient: 1.99417822037237"
## Number of never-treated units: 347 
## Proportion of never-treated units: 0.347
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  347"
## [1] "PSM estimated coefficient: 2.04659160801783"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 347"
## [1] "REM estimated coefficient: 2.00972620933297"
## Number of never-treated units: 370 
## Proportion of never-treated units: 0.37
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  370"
## [1] "PSM estimated coefficient: 2.01408502016363"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 370"
## [1] "REM estimated coefficient: 2.0325568057719"
## Number of never-treated units: 346 
## Proportion of never-treated units: 0.346
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  346"
## [1] "PSM estimated coefficient: 2.0209211873248"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 346"
## [1] "REM estimated coefficient: 1.97412085345398"
## Number of never-treated units: 362 
## Proportion of never-treated units: 0.362
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  362"
## [1] "PSM estimated coefficient: 2.02790264164904"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 362"
## [1] "REM estimated coefficient: 2.01668988863577"
## Number of never-treated units: 362 
## Proportion of never-treated units: 0.362
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  362"
## [1] "PSM estimated coefficient: 2.03029240894137"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 362"
## [1] "REM estimated coefficient: 1.99312477938118"
## Number of never-treated units: 358 
## Proportion of never-treated units: 0.358
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  358"
## [1] "PSM estimated coefficient: 2.01031436210072"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 358"
## [1] "REM estimated coefficient: 1.98489384971573"
## Number of never-treated units: 354 
## Proportion of never-treated units: 0.354
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  354"
## [1] "PSM estimated coefficient: 1.95428700606808"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 354"
## [1] "REM estimated coefficient: 2.02115891523368"
## Number of never-treated units: 344 
## Proportion of never-treated units: 0.344
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  344"
## [1] "PSM estimated coefficient: 2.01096135065153"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 344"
## [1] "REM estimated coefficient: 1.97507122713095"
## Number of never-treated units: 359 
## Proportion of never-treated units: 0.359
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  359"
## [1] "PSM estimated coefficient: 1.97811784363017"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 359"
## [1] "REM estimated coefficient: 2.0121466319968"
## Number of never-treated units: 350 
## Proportion of never-treated units: 0.35
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  350"
## [1] "PSM estimated coefficient: 2.0448157903397"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 350"
## [1] "REM estimated coefficient: 2.00135101112025"
## Number of never-treated units: 345 
## Proportion of never-treated units: 0.345
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  345"
## [1] "PSM estimated coefficient: 2.01878959090384"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 345"
## [1] "REM estimated coefficient: 2.05207735672259"
## Number of never-treated units: 362 
## Proportion of never-treated units: 0.362
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  362"
## [1] "PSM estimated coefficient: 1.98298000325393"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 362"
## [1] "REM estimated coefficient: 1.99908299471614"
## Number of never-treated units: 344 
## Proportion of never-treated units: 0.344
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  344"
## [1] "PSM estimated coefficient: 2.07813658595388"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 344"
## [1] "REM estimated coefficient: 2.06280215742179"
## Number of never-treated units: 350 
## Proportion of never-treated units: 0.35
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  350"
## [1] "PSM estimated coefficient: 2.03221778081676"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 350"
## [1] "REM estimated coefficient: 2.01072138451725"
## Number of never-treated units: 342 
## Proportion of never-treated units: 0.342
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  342"
## [1] "PSM estimated coefficient: 2.00975544950927"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 342"
## [1] "REM estimated coefficient: 2.02239805580212"
## Number of never-treated units: 352 
## Proportion of never-treated units: 0.352
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  352"
## [1] "PSM estimated coefficient: 2.03284582390631"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 352"
## [1] "REM estimated coefficient: 1.97863579521842"
## Number of never-treated units: 356 
## Proportion of never-treated units: 0.356
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  356"
## [1] "PSM estimated coefficient: 2.0107891696498"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 356"
## [1] "REM estimated coefficient: 1.96777527714988"
## Number of never-treated units: 356 
## Proportion of never-treated units: 0.356
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  356"
## [1] "PSM estimated coefficient: 2.02710804589433"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 356"
## [1] "REM estimated coefficient: 1.97141691603048"
## Number of never-treated units: 364 
## Proportion of never-treated units: 0.364
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  364"
## [1] "PSM estimated coefficient: 2.03593581286479"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 364"
## [1] "REM estimated coefficient: 1.981229878234"
## Number of never-treated units: 358 
## Proportion of never-treated units: 0.358
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  358"
## [1] "PSM estimated coefficient: 1.98622093709243"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 358"
## [1] "REM estimated coefficient: 1.94358021262224"
## Number of never-treated units: 344 
## Proportion of never-treated units: 0.344
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  344"
## [1] "PSM estimated coefficient: 1.93712296298286"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 344"
## [1] "REM estimated coefficient: 1.99786595175811"
## Number of never-treated units: 370 
## Proportion of never-treated units: 0.37
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  370"
## [1] "PSM estimated coefficient: 1.98735506590986"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 370"
## [1] "REM estimated coefficient: 1.9942296290705"
## Number of never-treated units: 355 
## Proportion of never-treated units: 0.355
## Warning: Fewer control units than treated units; not all treated units will get
## a match.
## [1] "PSM, Number of matched pairs:  355"
## [1] "PSM estimated coefficient: 1.98815383622"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## [1] "REM, Number of matched pairs: 355"
## [1] "REM estimated coefficient: 1.96764620286405"
# Remove NAs from estimates (from simulations that were skipped)
estimates_psm <- estimates_psm[!is.na(estimates_psm)]
estimates_rem <- estimates_rem[!is.na(estimates_rem)]

# Analyze Simulation Results

# Bias
bias_psm <- mean(estimates_psm) - true_beta
bias_rem <- mean(estimates_rem) - true_beta

# Variance
variance_psm <- var(estimates_psm)
variance_rem <- var(estimates_rem)

# Mean Squared Error (MSE)
mse_psm <- bias_psm^2 + variance_psm
mse_rem <- bias_rem^2 + variance_rem

# Print Results
cat("PSM Bias:", bias_psm, "\n")
## PSM Bias: 0.01242876
cat("REM Bias:", bias_rem, "\n")
## REM Bias: 0.004079917
cat("PSM MSE:", mse_psm, "\n")
## PSM MSE: 0.001553947
cat("REM MSE:", mse_rem, "\n")
## REM MSE: 0.0009605021
# Visualize Results
library(ggplot2)

# Combine estimates into a data frame
results_df <- data.frame(
  Estimate = c(estimates_psm, estimates_rem),
  Method = rep(c("PSM", "REM"), times = c(length(estimates_psm), length(estimates_rem)))
)

# Density plot
ggplot(results_df, aes(x = Estimate, fill = Method)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = true_beta, linetype = "dashed") +
  labs(title = "Distribution of Estimated Treatment Effects",
       x = "Estimated Treatment Effect",
       y = "Density") +
  theme_minimal()