Synthetic Panel: R fect vs Python pyfector

Author

Alan Huang

Published

December 31, 2025

This document generates synthetic panel data with known DGP parameters, estimates FE / IFE / MC in both R fect and Python pyfector (1000 bootstrap), and compares ATT, betas, SEs, and p-values side by side.

Two scenarios: (A) balanced panel, (B) 5 % randomly dropped observations.

DGP: 500 units, 12 periods, 40% eventually treated (staggered), true ATT = 1.5, 3 covariates, 1 interactive factor, two-way FE.

1 – Generate Data & Run R fect

Data is generated in R; R fect estimates are run there and exported to JSON.

Code 
import subprocess, json, textwrap

r_script = textwrap.dedent(r"""
library(fect)
library(jsonlite)

set.seed(42)

# ---- DGP parameters ----
N <- 500; TT <- 12; treat_frac <- 0.40
true_att <- 1.5; n_cov <- 3; beta_true <- c(0.8, -0.5, 0.3)

alpha_i <- rnorm(N, 0, 1)
xi_t    <- seq(-1, 1, length.out = TT)
factor_t  <- sin(2 * pi * (1:TT) / TT)
loading_i <- rnorm(N, 0, 0.5)
X_array <- array(rnorm(TT * N * n_cov), dim = c(TT, N, n_cov))

n_tr <- round(N * treat_frac)
treat_start <- sample(4:TT, n_tr, replace = TRUE)
D_mat <- matrix(0L, TT, N)
for (j in seq_len(n_tr)) D_mat[treat_start[j]:TT, j] <- 1L

eps <- matrix(rnorm(TT * N, 0, 0.5), TT, N)
Y_mat <- matrix(0, TT, N)
for (t in 1:TT) for (j in 1:N) {
  Y_mat[t, j] <- alpha_i[j] + xi_t[t] + factor_t[t] * loading_i[j] +
                  sum(X_array[t, j, ] * beta_true) + D_mat[t, j] * true_att + eps[t, j]
}

rows <- expand.grid(time = 1:TT, id = 1:N)
rows$Y <- as.vector(Y_mat); rows$D <- as.vector(D_mat)
rows$X1 <- as.vector(X_array[,,1]); rows$X2 <- as.vector(X_array[,,2])
rows$X3 <- as.vector(X_array[,,3])
df_full <- rows

set.seed(42)
drop_idx <- sort(sample(nrow(df_full), round(0.05 * nrow(df_full))))
df_drop  <- df_full[-drop_idx, ]

write.csv(df_full, "/tmp/synth_panel_full.csv", row.names = FALSE)
write.csv(df_drop, "/tmp/synth_panel_drop5.csv", row.names = FALSE)

cat(sprintf("Full: %d obs | Dropped: %d obs\n", nrow(df_full), nrow(df_drop)))

# ---- Estimate ----
run_fect <- function(df, nboots = 1000) {
  results <- list()
  for (meth in c("fe", "ife", "mc")) {
    cat(sprintf("  %s ... ", toupper(meth)))
    extra <- list()
    if (meth == "ife") extra <- list(r = c(0, 5), CV = TRUE)
    if (meth == "mc")  extra <- list(CV = TRUE)
    t0 <- proc.time()
    out <- do.call(fect, c(list(
      formula = Y ~ D + X1 + X2 + X3, data = df,
      index = c("id", "time"), force = "two-way", method = meth,
      se = TRUE, nboots = nboots, seed = 123,
      parallel = TRUE, cores = 16
    ), extra))
    elapsed <- (proc.time() - t0)["elapsed"]
    att <- out$est.avg[, "ATT.avg"]; se_ <- out$est.avg[, "S.E."]
    pv  <- out$est.avg[, "p.value"]
    beta_vals <- if (!is.null(out$beta)) as.list(setNames(as.numeric(out$beta), rownames(out$beta))) else list()
    r_cv <- if (!is.null(out$r.cv)) out$r.cv else NA
    lm_cv <- if (!is.null(out$lambda.cv)) out$lambda.cv else NA
    results[[meth]] <- list(method=meth, att=att, se=se_, p=pv,
                            beta=beta_vals, r_cv=r_cv, lambda_cv=lm_cv)
    sig <- ifelse(pv<0.01,"***",ifelse(pv<0.05,"**",ifelse(pv<0.1,"*","")))
    cat(sprintf("ATT=%7.4f SE=%7.4f p=%.4f%s (%.0fs)\n", att, se_, pv, sig, elapsed))
  }
  results
}

cat("\n=== Full Panel ===\n")
r_full <- run_fect(df_full)
cat("\n=== 5% Dropped ===\n")
r_drop <- run_fect(df_drop)

export <- list()
for (sc in c("full","drop")) {
  rr <- if (sc=="full") r_full else r_drop
  for (m in names(rr)) {
    e <- rr[[m]]
    export[[length(export)+1]] <- list(scenario=sc, method=e$method,
      att=e$att, se=e$se, p=e$p, beta=e$beta, r_cv=e$r_cv, lambda_cv=e$lambda_cv)
  }
}
write_json(export, "/tmp/synth_r_results.json", auto_unbox=TRUE)
cat("\nDone.\n")
""")

import os, pathlib
env = os.environ.copy()
# Ensure R can find user-installed packages
r_user_lib = pathlib.Path.home() / "R" / "x86_64-pc-linux-gnu-library" / "4.3"
if r_user_lib.exists():
    env["R_LIBS_USER"] = str(r_user_lib)
result = subprocess.run(
    ["Rscript", "--vanilla", "-e", r_script],
    capture_output=True, text=True, timeout=1800, env=env,
)
print(result.stdout)
if result.returncode != 0:
    print("STDERR:", result.stderr[-2000:])
Full: 6000 obs | Dropped: 5700 obs

=== Full Panel ===
  FE ... ATT= 1.5173 SE= 0.0557 p=0.0000*** (18s)
  IFE ... ATT= 1.4900 SE= 0.0791 p=0.0000*** (72s)
  MC ... ATT= 1.5853 SE= 0.0574 p=0.0000*** (46s)

=== 5% Dropped ===
  FE ... ATT= 1.5070 SE= 0.0570 p=0.0000*** (18s)
  IFE ... ATT= 1.5098 SE= 0.0774 p=0.0000*** (70s)
  MC ... ATT= 1.5820 SE= 0.0627 p=0.0000*** (47s)

Done.

2 – Python pyfector Estimation

Code 
import json, time
import numpy as np
import polars as pl
import pyfector
from IPython.display import display, HTML

np.random.seed(123)

df_full = pl.read_csv("/tmp/synth_panel_full.csv")
df_drop = pl.read_csv("/tmp/synth_panel_drop5.csv")

with open("/tmp/synth_r_results.json") as f:
    r_res = {(e["scenario"], e["method"]): e for e in json.load(f)}

covariates = ["X1", "X2", "X3"]
py_res = {}

for scenario, df in [("full", df_full), ("drop", df_drop)]:
    label = "Full Panel" if scenario == "full" else "5% Dropped"
    print(f"\n{'='*60}")
    print(f"  {label}  ({len(df)} obs)")
    print(f"{'='*60}")

    for mname, mkwargs in [
        ("fe",  dict(method="fe",  force="two-way")),
        ("ife", dict(method="ife", force="two-way", CV=True, r=(0, 5))),
        ("mc",  dict(method="mc",  force="two-way", CV=True)),
    ]:
        t0 = time.time()
        result = pyfector.fect(
            data=df, Y="Y", D="D",
            index=("id", "time"), X=covariates,
            se=True, nboots=1000, n_jobs=16, seed=123,
            **mkwargs,
        )
        elapsed = time.time() - t0
        py_res[(scenario, mname)] = result
        inf = result.inference
        sig = "***" if inf.att_avg_pval < 0.01 else "**" if inf.att_avg_pval < 0.05 else "*" if inf.att_avg_pval < 0.1 else ""
        extra = ""
        if mname == "ife" and hasattr(result, "r_cv"):
            extra = f"  r_cv={result.r_cv}"
        if mname == "mc" and result.lambda_cv is not None:
            extra = f"  lam={result.lambda_cv:.6f}"
        print(f"  {mname.upper():4s}  ATT={result.att_avg:8.4f}  SE={inf.att_avg_se:8.4f}  "
              f"p={inf.att_avg_pval:.4f}{sig}{extra}  ({elapsed:.1f}s)")

============================================================
  Full Panel  (6000 obs)
============================================================
  FE    ATT=  1.5173  SE=  0.0569  p=0.0000***  (1.1s)
  IFE   ATT=  1.4064  SE=  0.1371  p=0.0000***  r_cv=1  (45.7s)
  MC    ATT=  1.5116  SE=  0.0538  p=0.0000***  lam=0.000379  (16.7s)

============================================================
  5% Dropped  (5700 obs)
============================================================
  FE    ATT=  1.5070  SE=  0.0580  p=0.0000***  (1.1s)
  IFE   ATT=  1.2803  SE=  0.1675  p=0.0000***  r_cv=1  (77.0s)
  MC    ATT=  1.5047  SE=  0.0549  p=0.0000***  lam=0.000865  (10.6s)

3 – Side-by-Side Comparison

ATT & Inference

Code 
TRUE_ATT = 1.5

rows = []
for scenario in ["full", "drop"]:
    for mname in ["fe", "ife", "mc"]:
        r = r_res.get((scenario, mname))
        py = py_res.get((scenario, mname))
        if not r or py is None:
            continue
        r_att, r_se, r_p = r["att"], r["se"], r["p"]
        py_att = py.att_avg
        py_se = py.inference.att_avg_se
        py_p = py.inference.att_avg_pval
        att_pct = abs(py_att - r_att) / max(abs(r_att), 1e-8) * 100
        def sig(p):
            if p < 0.01: return "***"
            if p < 0.05: return "**"
            if p < 0.1:  return "*"
            return ""
        rows.append({
            "Scenario": "Balanced" if scenario == "full" else "5% Drop",
            "Method": mname.upper(),
            "True": f"{TRUE_ATT:.1f}",
            "R ATT": f"{r_att:.4f}", "R SE": f"{r_se:.4f}",
            "R p": f"{r_p:.4f}{sig(r_p)}",
            "Py ATT": f"{py_att:.4f}", "Py SE": f"{py_se:.4f}",
            "Py p": f"{py_p:.4f}{sig(py_p)}",
            "% diff": f"{att_pct:.1f}%",
            "Sig": "Y" if sig(r_p) == sig(py_p) else "N",
        })

display(HTML(pl.DataFrame(rows).to_pandas().to_html(index=False)))
Scenario Method True R ATT R SE R p Py ATT Py SE Py p % diff Sig
Balanced FE 1.5 1.5173 0.0557 0.0000*** 1.5173 0.0569 0.0000*** 0.0% Y
Balanced IFE 1.5 1.4900 0.0791 0.0000*** 1.4064 0.1371 0.0000*** 5.6% Y
Balanced MC 1.5 1.5853 0.0574 0.0000*** 1.5116 0.0538 0.0000*** 4.6% Y
5% Drop FE 1.5 1.5070 0.0570 0.0000*** 1.5070 0.0580 0.0000*** 0.0% Y
5% Drop IFE 1.5 1.5098 0.0774 0.0000*** 1.2803 0.1675 0.0000*** 15.2% Y
5% Drop MC 1.5 1.5820 0.0627 0.0000*** 1.5047 0.0549 0.0000*** 4.9% Y

Covariate Coefficients

Code 
TRUE_BETA = {"X1": 0.8, "X2": -0.5, "X3": 0.3}
rows = []
for scenario in ["full", "drop"]:
    for mname in ["fe", "ife", "mc"]:
        r = r_res.get((scenario, mname))
        py = py_res.get((scenario, mname))
        if not r or py is None: continue
        r_beta_raw = r.get("beta", {})
        # Handle both dict {"X1": 0.8} and list [0.8, -0.5, 0.3] formats
        if isinstance(r_beta_raw, list):
            r_beta = dict(zip(covariates, r_beta_raw))
        else:
            r_beta = r_beta_raw
        py_beta = dict(zip(py.covariate_names, py.beta)) if py.beta is not None else {}
        for cov in covariates:
            r_b, py_b = r_beta.get(cov), py_beta.get(cov)
            if r_b is None or py_b is None: continue
            rows.append({
                "Scenario": "Balanced" if scenario == "full" else "5% Drop",
                "Method": mname.upper(), "Cov": cov,
                "True": f"{TRUE_BETA[cov]:.1f}",
                "R": f"{r_b:.4f}", "Py": f"{py_b:.4f}",
                "|diff|": f"{abs(py_b - r_b):.4f}",
            })
display(HTML(pl.DataFrame(rows).to_pandas().to_html(index=False)))
Scenario Method Cov True R Py |diff|
Balanced FE X1 0.8 0.7980 0.7980 0.0000
Balanced FE X2 -0.5 -0.4918 -0.4918 0.0000
Balanced FE X3 0.3 0.2851 0.2851 0.0000
Balanced IFE X1 0.8 0.7975 0.7973 0.0002
Balanced IFE X2 -0.5 -0.4908 -0.4894 0.0014
Balanced IFE X3 0.3 0.2884 0.2892 0.0008
Balanced MC X1 0.8 0.7956 0.7964 0.0008
Balanced MC X2 -0.5 -0.4921 -0.4936 0.0015
Balanced MC X3 0.3 0.2861 0.2850 0.0011
5% Drop FE X1 0.8 0.7977 0.7977 0.0000
5% Drop FE X2 -0.5 -0.4895 -0.4895 0.0000
5% Drop FE X3 0.3 0.2895 0.2895 0.0000
5% Drop IFE X1 0.8 0.7965 0.7956 0.0009
5% Drop IFE X2 -0.5 -0.4919 -0.4899 0.0020
5% Drop IFE X3 0.3 0.2940 0.2941 0.0001
5% Drop MC X1 0.8 0.7956 0.7952 0.0004
5% Drop MC X2 -0.5 -0.4922 -0.4938 0.0016
5% Drop MC X3 0.3 0.2911 0.2899 0.0012

CV Selection

Code 
rows = []
for scenario in ["full", "drop"]:
    label = "Balanced" if scenario == "full" else "5% Drop"
    for mname, param in [("ife", "r"), ("mc", "lambda")]:
        r = r_res.get((scenario, mname))
        py = py_res.get((scenario, mname))
        if not r or py is None: continue
        r_val = r.get("r_cv") if param == "r" else r.get("lambda_cv")
        py_val = getattr(py, "r_cv", None) if param == "r" else getattr(py, "lambda_cv", None)
        r_str = f"{r_val:.6f}" if isinstance(r_val, float) else str(r_val)
        py_str = f"{py_val:.6f}" if isinstance(py_val, float) else str(py_val)
        rows.append({"Scenario": label, "Method": mname.upper(), "Param": param,
                      "R": r_str, "Python": py_str})
display(HTML(pl.DataFrame(rows).to_pandas().to_html(index=False)))
Scenario Method Param R Python
Balanced IFE r 1 1
Balanced MC lambda 0.000900 0.000379
5% Drop IFE r 1 1
5% Drop MC lambda 0.000900 0.000865

4 – Event Study Plots

Code 
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 3, figsize=(18, 12))
method_labels = {"fe": "FE", "ife": "IFE", "mc": "MC"}

for row_idx, (scenario, slab) in enumerate([("full", "Balanced"), ("drop", "5% Dropped")]):
    for col_idx, mname in enumerate(["fe", "ife", "mc"]):
        ax = axes[row_idx, col_idx]
        res = py_res.get((scenario, mname))
        if res is None:
            ax.text(0.5, 0.5, "N/A", ha="center", va="center")
            continue
        time_on = res.time_on; att_on = res.att_on
        ci_lo = res.inference.att_on_ci_lower if res.inference else None
        ci_hi = res.inference.att_on_ci_upper if res.inference else None
        ax.axhline(0, color="gray", ls="--", lw=0.8)
        ax.axhline(1.5, color="red", ls=":", lw=1, alpha=0.6, label="True ATT=1.5")
        ax.axvline(-0.5, color="gray", ls=":", lw=0.8, alpha=0.5)
        if ci_lo is not None:
            ax.fill_between(time_on, ci_lo, ci_hi, alpha=0.15, color="steelblue")
        ax.plot(time_on, att_on, color="steelblue", lw=1.5)
        if row_idx == 0: ax.set_title(f"{method_labels[mname]}", fontsize=14, fontweight="bold")
        if col_idx == 0: ax.set_ylabel(slab, fontsize=12, fontweight="bold")
        ax.set_xlabel("Periods Rel. to Treatment" if row_idx == 1 else "")
        ax.tick_params(labelsize=9); ax.legend(fontsize=8, loc="upper left")

plt.suptitle("Dynamic Treatment Effects (pyfector) -- True ATT = 1.5",
             fontsize=16, fontweight="bold", y=1.01)
plt.tight_layout(); plt.show()

5 – Diagnostic Tests

Code 
for scenario, slab in [("full", "Balanced"), ("drop", "5% Dropped")]:
    print(f"\n{'='*60}")
    print(f"  {slab}")
    print(f"{'='*60}")
    for mname in ["fe", "ife", "mc"]:
        res = py_res.get((scenario, mname))
        if res is None or res.inference is None:
            continue
        diag = res.diagnose(
            f_threshold=0.5,
            tost_threshold=0.36,
            placebo_period=(-3, -1),
            loo=True,
        )
        print(f"\n--- {mname.upper()} ---")
        print(diag.summary())

============================================================
  Balanced
============================================================

--- FE ---
Diagnostic Tests
==================================================
Pre-trend F-test:
  F(11,989) = 1.2616, p = 0.2418
Equivalence F-test:
  p = 0.0000 (threshold=0.5)
TOST per pre-period:
  t=-11: p=0.0002 [PASS]
  t=-10: p=0.0000 [PASS]
  t=-9: p=0.0000 [PASS]
  t=-8: p=0.0000 [PASS]
  t=-7: p=0.0000 [PASS]
  t=-6: p=0.0000 [PASS]
  t=-5: p=0.0000 [PASS]
  t=-4: p=0.0000 [PASS]
  t=-3: p=0.0000 [PASS]
  t=-2: p=0.0000 [PASS]
  t=-1: p=0.0000 [PASS]
Placebo test:
  ATT=-0.0071, p=0.7740
Leave-one-out:
  Max ATT change = 0.022509

--- IFE ---
Diagnostic Tests
==================================================
Pre-trend F-test:
  F(11,989) = 0.9870, p = 0.4562
Equivalence F-test:
  p = 0.0000 (threshold=0.5)
TOST per pre-period:
  t=-11: p=0.0007 [PASS]
  t=-10: p=0.0000 [PASS]
  t=-9: p=0.0000 [PASS]
  t=-8: p=0.0000 [PASS]
  t=-7: p=0.0000 [PASS]
  t=-6: p=0.0000 [PASS]
  t=-5: p=0.0000 [PASS]
  t=-4: p=0.0000 [PASS]
  t=-3: p=0.0000 [PASS]
  t=-2: p=0.0000 [PASS]
  t=-1: p=0.0000 [PASS]
Placebo test:
  ATT=-0.0033, p=0.6940
Leave-one-out:
  Max ATT change = 0.037610

--- MC ---
Diagnostic Tests
==================================================
Pre-trend F-test:
  F(11,989) = 1.2532, p = 0.2470
Equivalence F-test:
  p = 0.0000 (threshold=0.5)
TOST per pre-period:
  t=-11: p=0.0000 [PASS]
  t=-10: p=0.0000 [PASS]
  t=-9: p=0.0000 [PASS]
  t=-8: p=0.0000 [PASS]
  t=-7: p=0.0000 [PASS]
  t=-6: p=0.0000 [PASS]
  t=-5: p=0.0000 [PASS]
  t=-4: p=0.0000 [PASS]
  t=-3: p=0.0000 [PASS]
  t=-2: p=0.0000 [PASS]
  t=-1: p=0.0000 [PASS]
Placebo test:
  ATT=-0.0018, p=0.5540
Leave-one-out:
  Max ATT change = 0.020828

============================================================
  5% Dropped
============================================================

--- FE ---
Diagnostic Tests
==================================================
Pre-trend F-test:
  F(11,989) = 0.9560, p = 0.4853
Equivalence F-test:
  p = 0.0000 (threshold=0.5)
TOST per pre-period:
  t=-11: p=0.0002 [PASS]
  t=-10: p=0.0004 [PASS]
  t=-9: p=0.0000 [PASS]
  t=-8: p=0.0001 [PASS]
  t=-7: p=0.0000 [PASS]
  t=-6: p=0.0000 [PASS]
  t=-5: p=0.0000 [PASS]
  t=-4: p=0.0000 [PASS]
  t=-3: p=0.0000 [PASS]
  t=-2: p=0.0000 [PASS]
  t=-1: p=0.0000 [PASS]
Placebo test:
  ATT=0.0116, p=0.6020
Leave-one-out:
  Max ATT change = 0.059113

--- IFE ---
Diagnostic Tests
==================================================
Pre-trend F-test:
  F(11,989) = 0.8201, p = 0.6200
Equivalence F-test:
  p = 0.0000 (threshold=0.5)
TOST per pre-period:
  t=-11: p=0.0003 [PASS]
  t=-10: p=0.0000 [PASS]
  t=-9: p=0.0000 [PASS]
  t=-8: p=0.0000 [PASS]
  t=-7: p=0.0000 [PASS]
  t=-6: p=0.0000 [PASS]
  t=-5: p=0.0000 [PASS]
  t=-4: p=0.0000 [PASS]
  t=-3: p=0.0000 [PASS]
  t=-2: p=0.0000 [PASS]
  t=-1: p=0.0000 [PASS]
Placebo test:
  ATT=-0.0006, p=0.9780
Leave-one-out:
  Max ATT change = 0.068158

--- MC ---
Diagnostic Tests
==================================================
Pre-trend F-test:
  F(11,989) = 1.0002, p = 0.4441
Equivalence F-test:
  p = 0.0000 (threshold=0.5)
TOST per pre-period:
  t=-11: p=0.0000 [PASS]
  t=-10: p=0.0000 [PASS]
  t=-9: p=0.0000 [PASS]
  t=-8: p=0.0000 [PASS]
  t=-7: p=0.0000 [PASS]
  t=-6: p=0.0000 [PASS]
  t=-5: p=0.0000 [PASS]
  t=-4: p=0.0000 [PASS]
  t=-3: p=0.0000 [PASS]
  t=-2: p=0.0000 [PASS]
  t=-1: p=0.0000 [PASS]
Placebo test:
  ATT=0.0010, p=0.8620
Leave-one-out:
  Max ATT change = 0.056556

6 – Built-in Plots

Gap Plots (Event Study via pyfector.plot)

Code 
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
method_labels = {"fe": "FE", "ife": "IFE", "mc": "MC"}

for row_idx, (scenario, slab) in enumerate([("full", "Balanced"), ("drop", "5% Dropped")]):
    for col_idx, mname in enumerate(["fe", "ife", "mc"]):
        ax = axes[row_idx, col_idx]
        res = py_res.get((scenario, mname))
        if res is None:
            ax.text(0.5, 0.5, "N/A", ha="center", va="center")
            continue
        res.plot(kind="gap", ax=ax, show_ci=True)
        ax.axhline(1.5, color="red", ls=":", lw=1, alpha=0.6, label="True ATT")
        if row_idx == 0:
            ax.set_title(f"{method_labels[mname]}", fontsize=14, fontweight="bold")
        if col_idx == 0:
            ax.set_ylabel(slab, fontsize=12, fontweight="bold")
        ax.legend(fontsize=8, loc="upper left")

plt.suptitle("Gap Plots (pyfector built-in) -- True ATT = 1.5",
             fontsize=16, fontweight="bold", y=1.01)
plt.tight_layout(); plt.show()

Equivalence Plots (TOST)

Code 
fig, axes = plt.subplots(2, 3, figsize=(18, 10))

for row_idx, (scenario, slab) in enumerate([("full", "Balanced"), ("drop", "5% Dropped")]):
    for col_idx, mname in enumerate(["fe", "ife", "mc"]):
        ax = axes[row_idx, col_idx]
        res = py_res.get((scenario, mname))
        if res is None:
            ax.text(0.5, 0.5, "N/A", ha="center", va="center")
            continue
        res.plot(kind="equiv", ax=ax, threshold=0.36)
        if row_idx == 0:
            ax.set_title(f"{method_labels[mname]}", fontsize=14, fontweight="bold")
        if col_idx == 0:
            ax.set_ylabel(slab, fontsize=12, fontweight="bold")

plt.suptitle("Pre-treatment Equivalence (TOST bounds = 0.36)",
             fontsize=16, fontweight="bold", y=1.01)
plt.tight_layout(); plt.show()

Treatment Status Heatmap (FE, Balanced)

Code 
fig, ax = plt.subplots(1, 1, figsize=(14, 6))
res = py_res.get(("full", "fe"))
if res is not None:
    res.plot(kind="status", ax=ax)
    ax.set_title("Treatment Status Heatmap (Balanced, FE)", fontsize=14, fontweight="bold")
plt.tight_layout(); plt.show()

7 – Summary

Code 
import sys
print(f"True ATT  : 1.5")
print(f"True beta : X1=0.8, X2=-0.5, X3=0.3\n")
for scenario, slab in [("full", "Balanced"), ("drop", "5% Dropped")]:
    n_total = n_sig = n_close = 0
    for mname in ["fe", "ife", "mc"]:
        r = r_res.get((scenario, mname)); py = py_res.get((scenario, mname))
        if not r or py is None: continue
        n_total += 1
        def sig(p):
            if p < 0.01: return "***"
            if p < 0.05: return "**"
            if p < 0.1:  return "*"
            return ""
        if sig(r["p"]) == sig(py.inference.att_avg_pval): n_sig += 1
        if abs(r["att"]) > 0.01 and abs(py.att_avg - r["att"]) / abs(r["att"]) < 0.10: n_close += 1
    print(f"--- {slab} ---")
    print(f"  Models compared      : {n_total}")
    print(f"  Significance match   : {n_sig}/{n_total}")
    print(f"  ATT within 10% of R  : {n_close}/{n_total}\n")
print(f"Python: {sys.version.split()[0]} | pyfector: {pyfector.__version__} | numpy: {np.__version__}")
True ATT  : 1.5
True beta : X1=0.8, X2=-0.5, X3=0.3

--- Balanced ---
  Models compared      : 3
  Significance match   : 3/3
  ATT within 10% of R  : 3/3

--- 5% Dropped ---
  Models compared      : 3
  Significance match   : 3/3
  ATT within 10% of R  : 2/3

Python: 3.12.3 | pyfector: 0.1.0 | numpy: 2.4.4