Causal Inference Mixtape

Causal Inference: The Mixtape — Code Skill

Practitioner-oriented causal inference skill built from Scott Cunningham's Causal Inference: The Mixtape repository. Covers 10 identification strategies with ready-to-run code templates in Python, R, and Stata.

---

Methods Covered

---

Core Workflow

Implement a Causal Method

1. Identify the method from the table above 2. Load the appropriate template from references/method-patterns.md 3. Adapt variable names, fixed effects, and clustering to the user's data 4. Add robustness checks (parallel trends for DiD, McCrary for RDD, first-stage F for IV)

Choose the Right Language

Cross-Language Equivalents

---

Key Python Patterns

DiD with Cluster-Robust SE

import statsmodels.formula.api as smf

model = smf.ols('y ~ C(treated)*C(post) + controls', data=df)
results = model.fit(cov_type='cluster', cov_kwds={'groups': df['firm_id']})

Event Study (Lead/Lag)

# Create relative time dummies
for k in range(-4, 5):
    col = f'rel_{k}' if k >= 0 else f'rel_m{abs(k)}'
    df[col] = (df['relative_time'] == k).astype(int)

# Drop t=-1 as reference
formula = 'y ~ ' + ' + '.join([c for c in rel_cols if c != 'rel_m1']) + ' + C(id) + C(year)'

IV / 2SLS

from linearmodels.iv import IV2SLS

model = IV2SLS.from_formula('y ~ 1 + exog + [endog ~ instrument]', data=df)
results = model.fit(cov_type='clustered', clusters=df['cluster_var'])

---

Robustness Check Patterns

---

Common Pitfalls

1. TWFE with staggered treatment — standard two-way FE is biased when treatment timing varies. Use Bacon decomposition or Sun & Abraham / Callaway & Sant'Anna estimators. 2. Synthetic Control with many treated units — the Synth package handles one treated unit. For multiple, use augmented synthetic control or stacked approach. 3. RDD without McCrary test — always test for manipulation at the cutoff before estimating. 4. IV weak instruments — report first-stage F-statistic. Below 10 indicates weak instrument bias. 5. Python Synth gap — no mature Python Synth package exists. Use rpy2 to call R's Synth from Python.

---

Additional Resources

Reference Files

• `references/method-patterns.md` — Detailed code templates for all 10 methods with full examples
• `references/r-stata-comparison.md` — Cross-language package comparison and method coverage gaps

Prompt Files

• `prompts/01-implement-method.md` — Copy-paste prompt for implementing any causal method
• `prompts/02-robustness-checks.md` — Copy-paste prompt for generating robustness check code

Prompt 1: Implement a Causal Inference Method

Copy and paste the prompt below into Claude with your details filled in.

---

You are an expert econometrician implementing causal inference methods.

Implement a complete [METHOD] analysis pipeline in [LANGUAGE: Python / R / Stata].

Requirements:
1. Data preparation (variable creation, sample restrictions)
2. Main estimation with correct standard errors
3. Key diagnostic / robustness check
4. Publication-ready output (coefficient table or plot)

Method-specific requirements:

- DiD: Include parallel trends event study plot. Cluster SE at [level]. Report DiD coefficient with baseline mean for economic magnitude.
- RDD: Include McCrary density test, bandwidth robustness (half/double), polynomial robustness. Report local linear estimate.
- IV: Report first-stage F-statistic. Defend exclusion restriction. Report Wu-Hausman test.
- Synthetic Control: Pre-treatment fit (RMSPE), placebo distribution, gaps plot.
- Matching/IPW: Covariate balance table before and after. Trimming at [0.1, 0.9].
- Event Study: Dynamic coefficients plot with 95% CI. Reference period = t-1.

My details:
- Method: [e.g., Difference-in-Differences]
- Language: [e.g., Python]
- Outcome variable: [e.g., firm_investment]
- Treatment variable: [e.g., reform_exposure]
- Treatment timing: [e.g., 2014 for all treated units / staggered]
- Key controls: [e.g., firm size, leverage, ROA]
- Fixed effects: [e.g., firm + year]
- Clustering level: [e.g., firm]
- Data format: [e.g., panel, entity_id + year columns]
- Sample size: [e.g., ~50,000 firm-years]
- Key concern: [e.g., contemporaneous policies]

[PASTE SAMPLE OF YOUR DATA STRUCTURE OR DESCRIBE COLUMNS]

---

Method-Specific Variants

For Staggered DiD

Additional requirement: Treatment timing varies across units.
- Use [Callaway & Sant'Anna / Sun & Abraham / Bacon decomposition] to address TWFE bias.
- Show that standard TWFE is potentially biased.
- Report group-time ATTs and aggregated dynamic effects.

My staggered details:
- Treatment cohorts: [e.g., 2010, 2012, 2014, 2016]
- Never-treated group exists: [yes/no]
- Preferred estimator: [e.g., Callaway & Sant'Anna]

For Fuzzy RDD

Additional requirement: Treatment assignment is not sharp at the cutoff.
- Implement fuzzy RDD as IV where crossing the cutoff instruments for treatment.
- Report both reduced-form and 2SLS estimates.
- Show first-stage discontinuity in treatment probability.

My fuzzy RDD details:
- Running variable: [e.g., vote share]
- Cutoff: [e.g., 50%]
- Treatment: [e.g., policy implementation — not all units above cutoff comply]
- Compliance rate above cutoff: [e.g., ~75%]

Prompt 2: Generate Robustness Check Code

Copy and paste the relevant section below based on your identification strategy.

---

For DiD Papers

You are an expert econometrician. Generate robustness check code in [Python / R / Stata] for my DiD analysis.

Produce code for ALL of the following:

1. EVENT STUDY (parallel trends):
   - Dynamic specification with lead/lag dummies
   - Plot coefficients with 95% CI
   - Pre-period joint F-test

2. PLACEBO TREATMENT DATE:
   - Re-estimate using a fake treatment date [N] years before actual treatment
   - Expect null result

3. BACON DECOMPOSITION (if staggered):
   - Decompose TWFE into 2x2 comparisons
   - Identify problematic "already-treated vs later-treated" weight

4. ALTERNATIVE CONTROL GROUP:
   - Re-estimate dropping [specific units] from control group
   - Verify results hold

5. ENTROPY BALANCING / PSM-DID:
   - Re-weight sample to achieve covariate balance
   - Re-estimate on balanced sample

My details:
- Treatment: [e.g., 2014 SOE reform]
- Treated group: [e.g., listed SOEs]
- Control group: [e.g., non-SOE listed firms]
- Treatment timing: [e.g., 2014 for all / staggered]
- Outcome: [e.g., abnormal investment]
- Pre-period: [e.g., 2010-2013]
- Post-period: [e.g., 2015-2018]
- Covariates for balancing: [e.g., size, leverage, ROA, age]

---

For RDD Papers

You are an expert econometrician. Generate robustness check code in [Python / R / Stata] for my RDD analysis.

Produce code for ALL of the following:

1. McCRARY DENSITY TEST:
   - Test for manipulation at the cutoff
   - Report test statistic and p-value
   - Plot density

2. COVARIATE BALANCE:
   - Test each covariate for discontinuity at cutoff
   - Report coefficients and p-values in a table

3. BANDWIDTH ROBUSTNESS:
   - Re-estimate with bandwidths: h/2, h, 3h/2, 2h (h = IK optimal)
   - Table of estimates across bandwidths

4. POLYNOMIAL ROBUSTNESS:
   - Linear, quadratic, cubic specifications
   - Report all three estimates

5. PLACEBO CUTOFFS:
   - Re-estimate at false cutoffs (e.g., median of each side)
   - Expect null results

My details:
- Running variable: [e.g., vote share percentage]
- Cutoff: [e.g., 50% majority threshold]
- Outcome: [e.g., firm ESG score]
- IK optimal bandwidth: [e.g., 8.3 percentage points]
- Covariates for balance test: [e.g., firm size, age, leverage]

---

For IV Papers

You are an expert econometrician. Generate robustness check code in [Python / R / Stata] for my IV analysis.

Produce code for ALL of the following:

1. FIRST-STAGE DIAGNOSTICS:
   - First-stage regression with F-statistic
   - Report coefficient on instrument(s)
   - Cragg-Donald / Kleibergen-Paap F-stat

2. EXCLUSION RESTRICTION SUPPORT:
   - Placebo test: regress outcome on instrument controlling for endogenous variable
   - If coefficient on instrument ≈ 0, supports exclusion

3. OVER-IDENTIFICATION TEST (if multiple instruments):
   - Hansen J / Sargan test
   - Report test statistic and p-value

4. REDUCED FORM:
   - Regress outcome directly on instrument(s)
   - Should have same sign as 2SLS estimate

5. OLS vs IV COMPARISON:
   - Report OLS and IV side by side
   - Wu-Hausman endogeneity test

My details:
- Endogenous variable: [e.g., corruption level]
- Instrument(s): [e.g., ethnic fractionalization, distance to coast]
- Outcome: [e.g., patent count]
- Controls: [e.g., GDP per capita, education, population]
- First-stage F: [e.g., 23.4]

Causal Inference: The Mixtape — Claude Code Skill

A Claude Code skill providing ready-to-run code templates for causal inference methods, built from Scott Cunningham's Causal Inference: The Mixtape repository.

Languages: Python · R · Stata

---

What It Does

This skill helps you:

1. Implement causal inference methods — DiD, RDD, IV, Synthetic Control, Matching, and more 2. Choose the right language — cross-language equivalents and coverage gap analysis 3. Write robustness checks — parallel trends, McCrary tests, Bacon decomposition, bandwidth robustness 4. Avoid common pitfalls — staggered DiD bias, weak instruments, missing diagnostics

Methods Covered (10)

Trigger Phrases

Say any of the following to activate this skill:

• implement a DiD regression
• write a causal inference pipeline
• set up an event study
• implement instrumental variables
• run a regression discontinuity design
• build a synthetic control model
• implement propensity score matching
• implement Bacon decomposition

---

Installation

Copy the skill folder to your Claude Code skills directory:

cp -r causal-inference-mixtape ~/.claude/skills/

Or clone directly:

git clone https://github.com/Jill0099/causal-inference-mixtape.git ~/.claude/skills/causal-inference-mixtape

---

File Structure

causal-inference-mixtape/
├── SKILL.md                              # Core skill (auto-loaded when triggered)
├── references/
│   ├── method-patterns.md               # Full code templates for all 10 methods
│   └── r-stata-comparison.md            # Cross-language coverage gaps & packages
└── prompts/
    ├── 01-implement-method.md           # Copy-paste: implement any causal method
    └── 02-robustness-checks.md          # Copy-paste: DiD/RDD/IV robustness code

---

Key Features

Cross-Language Equivalents

Python Gaps Documented

Some methods lack mature Python implementations:

• Synthetic Control → use rpy2 to call R's Synth
• Bacon Decomposition → use R (bacondecomp) or Stata
• Coarsened Exact Matching → use Stata (cem) or R (MatchIt)
• McCrary Density Test → use R (rdd)

Robustness Check Patterns

---

Prompts (Copy-Paste Ready)

The prompts/ folder contains standalone prompts for use without Claude Code:

Each prompt has fill-in fields — replace with your paper's details and paste into any Claude chat.

---

Source

Built from systematic analysis of Scott Cunningham's Causal Inference: The Mixtape repository:

• 58 Python scripts
• ~56 R scripts
• ~60 Stata .do files
• Full course curriculum (9 sections)

---

License

MIT

Method Patterns — Full Code Templates

Detailed code templates extracted from 58 Python scripts, ~56 R scripts, and ~60 Stata .do files in the Mixtape repository.

---

§1 OLS / Regression

Python

import pandas as pd
import statsmodels.formula.api as smf

# Basic OLS with robust SE
model = smf.ols('outcome ~ treatment + control1 + control2', data=df)
results = model.fit(cov_type='HC1')
print(results.summary())

# WLS (weighted least squares)
model = smf.wls('outcome ~ treatment', data=df, weights=df['weight'])
results = model.fit()

R

library(estimatr)

# OLS with robust SE (HC1)
model <- lm_robust(outcome ~ treatment + control1 + control2, data = df, se_type = "HC1")
summary(model)

# Clustered SE
model <- lm_robust(outcome ~ treatment, data = df, clusters = firm_id, se_type = "stata")

Stata

* Basic OLS with robust SE
reg outcome treatment control1 control2, robust

* Cluster SE
reg outcome treatment control1 control2, cluster(firm_id)

* High-dimensional FE
reghdfe outcome treatment control1, absorb(firm_id year) cluster(firm_id)

---

§2 Difference-in-Differences

Python

import statsmodels.formula.api as smf

# Standard 2x2 DiD
model = smf.ols('y ~ C(treated)*C(post)', data=df)
results = model.fit(cov_type='cluster', cov_kwds={'groups': df['state']})

# The DiD coefficient is the interaction term: C(treated)[T.1]:C(post)[T.1]
did_coef = results.params['C(treated)[T.1]:C(post)[T.1]']
print(f"DiD estimate: {did_coef:.4f} (SE: {results.bse['C(treated)[T.1]:C(post)[T.1]']:.4f})")

# With controls and entity + time FE
model = smf.ols('y ~ C(treated)*C(post) + control1 + C(entity_id) + C(year)', data=df)
results = model.fit(cov_type='cluster', cov_kwds={'groups': df['state']})

R

library(lfe)

# DiD with two-way FE
model <- felm(y ~ treated:post + controls | entity_id + year | 0 | state, data = df)
summary(model)

# Alternative with fixest
library(fixest)
model <- feols(y ~ treated:post + controls | entity_id + year, data = df, cluster = ~state)

Stata

* Standard DiD
reg y treated##post, cluster(state)

* With two-way FE
reghdfe y treated_post controls, absorb(entity_id year) cluster(state)

* Triple difference
reghdfe y treated##post##group controls, absorb(entity_id year) cluster(state)

---

§3 Event Study (Dynamic DiD)

Python

import numpy as np
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

# Create relative time variable
df['rel_time'] = df['year'] - df['treatment_year']

# Create dummies (drop t=-1 as reference)
leads_lags = list(range(-4, 0)) + list(range(0, 5))  # exclude -1
for k in leads_lags:
    if k < 0:
        df[f'lead{abs(k)}'] = (df['rel_time'] == k).astype(int)
    else:
        df[f'lag{k}'] = (df['rel_time'] == k).astype(int)

# Regression
vars_str = ' + '.join([f'lead{abs(k)}' for k in range(-4, 0)] + [f'lag{k}' for k in range(0, 5)])
formula = f'y ~ {vars_str} + C(entity_id) + C(year)'
model = smf.ols(formula, data=df)
results = model.fit(cov_type='cluster', cov_kwds={'groups': df['state']})

# Plot coefficients
coefs = []
ses = []
periods = list(range(-4, 0)) + [0] + list(range(0, 5))
# Insert 0 for reference period t=-1
# ... extract from results.params and results.bse

fig, ax = plt.subplots(figsize=(10, 6))
ax.errorbar(periods, coefs, yerr=[1.96*s for s in ses], fmt='o-', capsize=3)
ax.axhline(y=0, color='red', linestyle='--')
ax.axvline(x=-0.5, color='grey', linestyle='--', alpha=0.5)
ax.set_xlabel('Periods Relative to Treatment')
ax.set_ylabel('Coefficient Estimate')
ax.set_title('Event Study Plot')
plt.tight_layout()

R

library(fixest)

# Sun & Abraham (2021) interaction-weighted estimator
model <- feols(y ~ sunab(treatment_year, year) | entity_id + year, data = df, cluster = ~state)
iplot(model, main = "Event Study")

Stata

* Event study with reghdfe
reghdfe y lead4 lead3 lead2 lag0 lag1 lag2 lag3 lag4, ///
    absorb(entity_id year) cluster(state)

* Plot
coefplot, keep(lead* lag*) vertical yline(0) xline(4.5, lpattern(dash))

---

§4 Staggered DiD / TWFE Issues

Bacon Decomposition (R)

library(bacondecomp)

# Goodman-Bacon (2021) decomposition
bacon_out <- bacon(y ~ treatment, data = df, id_var = "entity_id", time_var = "year")
print(bacon_out)

# Weighted sum = TWFE estimate
# Shows which 2x2 comparisons drive the estimate
# Flags problematic "already-treated vs later-treated" comparisons

Bacon Decomposition (Stata)

* Install: ssc install bacondecomp
bacondecomp y treatment, ddetail

Callaway & Sant'Anna (R)

library(did)

# Group-time ATT
att_gt <- att_gt(yname = "y", tname = "year", idname = "entity_id",
                 gname = "treatment_year", data = df)
summary(att_gt)
ggdid(att_gt)

# Aggregate to overall ATT
agg <- aggte(att_gt, type = "dynamic")
ggdid(agg)

---

§5 Regression Discontinuity Design

Python (Sharp RDD)

import statsmodels.formula.api as smf
import numpy as np

# Center running variable at cutoff
df['x_centered'] = df['running_var'] - cutoff
df['treated'] = (df['running_var'] >= cutoff).astype(int)

# Local linear regression (bandwidth h)
h = 10  # IK optimal or manually set
subset = df[abs(df['x_centered']) <= h].copy()

# Polynomial interaction
model = smf.ols('y ~ treated * x_centered', data=subset)
results = model.fit(cov_type='HC1')
rdd_effect = results.params['treated']

# Bandwidth robustness: repeat with h/2, 2h
for bw in [h/2, h, 2*h]:
    sub = df[abs(df['x_centered']) <= bw]
    m = smf.ols('y ~ treated * x_centered', data=sub).fit(cov_type='HC1')
    print(f"BW={bw:.1f}: effect={m.params['treated']:.3f} (SE={m.bse['treated']:.3f})")

R

library(rdrobust)

# Automatic bandwidth selection + local polynomial
rd_result <- rdrobust(y = df$y, x = df$running_var, c = cutoff)
summary(rd_result)

# Plot
rdplot(y = df$y, x = df$running_var, c = cutoff,
       title = "RD Plot", x.label = "Running Variable", y.label = "Outcome")

Stata

* RD plot
rdplot y running_var, c(cutoff) graph_options(title("RD Plot"))

* RD estimate with rdrobust
rdrobust y running_var, c(cutoff) kernel(triangular) bwselect(mserd)

McCrary Density Test (R)

library(rdd)
DCdensity(df$running_var, cutpoint = cutoff, plot = TRUE)

---

§6 Instrumental Variables / 2SLS

Python

from linearmodels.iv import IV2SLS

# 2SLS estimation
# Formula: dependent ~ exogenous + [endogenous ~ instruments]
model = IV2SLS.from_formula(
    'y ~ 1 + control1 + control2 + [endog_var ~ instrument1 + instrument2]',
    data=df
)
results = model.fit(cov_type='clustered', clusters=df['cluster_var'])
print(results.summary)

# First-stage F-statistic
print(f"First-stage F: {results.first_stage.diagnostics['f.stat']:.1f}")

# Manual first stage for inspection
first_stage = smf.ols('endog_var ~ instrument1 + instrument2 + control1 + control2', data=df)
fs_results = first_stage.fit(cov_type='HC1')

R

library(AER)

# 2SLS
model <- ivreg(y ~ control1 + control2 + endog_var | control1 + control2 + instrument1 + instrument2,
               data = df)
summary(model, diagnostics = TRUE)  # Includes weak instrument test, Wu-Hausman, Sargan

Stata

* 2SLS
ivregress 2sls y control1 control2 (endog_var = instrument1 instrument2), first robust

* Post-estimation diagnostics
estat firststage    /* First-stage F */
estat overid        /* Sargan-Hansen test */
estat endogenous    /* Wu-Hausman */

---

§7 Synthetic Control

Python (via rpy2)

import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
pandas2ri.activate()

# Transfer data to R
ro.globalenv['df'] = df

# Run Synth in R
ro.r('''
library(Synth)
dataprep.out <- dataprep(
    foo = df,
    predictors = c("predictor1", "predictor2"),
    predictors.op = "mean",
    dependent = "outcome",
    unit.variable = "unit_id",
    time.variable = "year",
    treatment.identifier = treated_unit,
    controls.identifier = control_units,
    time.predictors.prior = pre_period,
    time.optimize.ssr = pre_period,
    time.plot = full_period
)
synth.out <- synth(dataprep.out)
synth.tables <- synth.tab(dataprep.res = dataprep.out, synth.res = synth.out)
''')

R (Native)

library(Synth)

dataprep.out <- dataprep(
    foo = df,
    predictors = c("gdp", "trade", "infrate"),
    predictors.op = "mean",
    dependent = "outcome",
    unit.variable = "unit_id",
    time.variable = "year",
    treatment.identifier = 1,
    controls.identifier = c(2:10),
    time.predictors.prior = 1980:1990,
    time.optimize.ssr = 1980:1990,
    time.plot = 1980:2000
)

synth.out <- synth(dataprep.out)
path.plot(synth.res = synth.out, dataprep.res = dataprep.out)
gaps.plot(synth.res = synth.out, dataprep.res = dataprep.out)

# Placebo tests (permutation)
library(SCtools)
placebo <- generate.placebos(dataprep.out, synth.out, Sigf.ipop = 5)
plot_placebos(placebo)
mspe.plot(placebo, discard.extreme = TRUE, mspe.limit = 20)

Stata

* Synthetic control
synth outcome predictor1 predictor2 outcome(1980) outcome(1985), ///
    trunit(1) trperiod(1990) figure

---

§8 Matching / PSM / IPW / CEM

Python (Propensity Score + IPW)

import statsmodels.formula.api as smf
import numpy as np

# Step 1: Estimate propensity score
logit = smf.logit('treated ~ x1 + x2 + x3', data=df).fit()
df['pscore'] = logit.predict()

# Step 2: IPW weights
df['ipw'] = np.where(
    df['treated'] == 1,
    1 / df['pscore'],
    1 / (1 - df['pscore'])
)

# Step 3: Weighted regression (ATE)
model = smf.wls('y ~ treated', data=df, weights=df['ipw'])
results = model.fit(cov_type='HC1')

# Trimming extreme weights
df_trimmed = df[(df['pscore'] > 0.1) & (df['pscore'] < 0.9)]

R (MatchIt)

library(MatchIt)
library(Zelig)

# Nearest neighbor matching
m.out <- matchit(treated ~ x1 + x2 + x3, data = df, method = "nearest", ratio = 1)
summary(m.out)
plot(m.out, type = "jitter")

# Estimate treatment effect on matched data
m.data <- match.data(m.out)
model <- lm(y ~ treated + x1 + x2 + x3, data = m.data, weights = weights)

R (IPW)

library(ipw)

# IPW weights
temp <- ipwpoint(
    exposure = treated,
    family = "binomial",
    link = "logit",
    numerator = ~ 1,
    denominator = ~ x1 + x2 + x3,
    data = df
)
df$ipw <- temp$ipw.weights

# Weighted model
library(survey)
design <- svydesign(ids = ~1, weights = ~ipw, data = df)
model <- svyglm(y ~ treated, design = design)

Stata (CEM + teffects)

* Coarsened Exact Matching
cem x1 (#5) x2 (#3) x3, treatment(treated)
reg y treated [iweight = cem_weights]

* Propensity Score Matching via teffects
teffects psmatch (y) (treated x1 x2 x3), atet

* IPW
teffects ipw (y) (treated x1 x2 x3), atet

---

§9 DAGs and Collider Bias

Conceptual Framework

# DAGs are primarily conceptual tools
# Use dagitty.net for interactive DAG drawing

# Key rules from Mixtape:
# 1. Condition on confounders (common causes of treatment and outcome)
# 2. Never condition on colliders (common effects of treatment and outcome)
# 3. Never condition on mediators (if estimating total effect)
# 4. Backdoor criterion: block all backdoor paths from treatment to outcome

R (ggdag)

library(ggdag)
library(dagitty)

# Define DAG
dag <- dagitty('dag {
    X -> Y
    Z -> X
    Z -> Y
    M -> X
    M -> Y
}')

# Identify adjustment sets
adjustmentSets(dag, exposure = "X", outcome = "Y")

# Plot
ggdag(dag) + theme_dag()

---

§10 Randomization Inference

Python

import numpy as np

def permutation_test(treatment, outcome, n_permutations=1000):
    """Sharp null hypothesis test via randomization inference."""
    observed_diff = outcome[treatment == 1].mean() - outcome[treatment == 0].mean()

    null_diffs = []
    for _ in range(n_permutations):
        perm_treatment = np.random.permutation(treatment)
        diff = outcome[perm_treatment == 1].mean() - outcome[perm_treatment == 0].mean()
        null_diffs.append(diff)

    p_value = np.mean(np.abs(null_diffs) >= np.abs(observed_diff))
    return observed_diff, p_value

obs_diff, p_val = permutation_test(df['treated'].values, df['y'].values)
print(f"Observed difference: {obs_diff:.4f}, RI p-value: {p_val:.4f}")

R

library(ri2)

# Declare randomization procedure
declaration <- declare_ra(N = nrow(df), m = sum(df$treated))

# Conduct randomization inference
ri_out <- conduct_ri(
    y ~ treated,
    declaration = declaration,
    sharp_hypothesis = 0,
    data = df
)
summary(ri_out)
plot(ri_out)

Stata

* Randomization inference
ritest treated _b[treated], reps(1000): reg y treated

R & Stata Comparison — Methods Not Available in Python

Cross-language coverage gaps and package recommendations.

---

Method Coverage Matrix

---

Python Gaps — Recommended Workarounds

Synthetic Control

No mature Python package. Two options:

1. rpy2 bridge (recommended for production):

import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
pandas2ri.activate()
ro.globalenv['df'] = df
ro.r('library(Synth); ...')

2. SparseSC (experimental): pip install SparseSC — limited functionality

Bacon Decomposition

No Python implementation. Use R:

library(bacondecomp)
bacon(y ~ treatment, data = df, id_var = "id", time_var = "year")

Coarsened Exact Matching (CEM)

Stata's cem command is the gold standard. R alternative:

library(MatchIt)
m.out <- matchit(treated ~ x1 + x2, data = df, method = "cem")

McCrary Density Test

R implementation:

library(rdd)
DCdensity(running_var, cutpoint = cutoff, plot = TRUE)

---

Package Quick Reference

R Packages

Stata Packages

Python Packages

---

When to Switch Languages