Skip to content

Commit 215fff5

Browse files
igerberigerber
authored
Merge pull request #220 from igerber/power-analysis-notebook
Update power analysis tutorial with simulation-based features
2 parents fe3d79a + e26c093 commit 215fff5

1 file changed

Lines changed: 198 additions & 65 deletions

File tree

‎docs/tutorials/06_power_analysis.ipynb‎

Lines changed: 198 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -4,40 +4,15 @@
44
"cell_type": "markdown",
55
"id": "cell-0",
66
"metadata": {},
7-
"source": [
8-
"# Power Analysis for Difference-in-Differences\n",
9-
"\n",
10-
"This notebook demonstrates how to use the power analysis tools in `diff-diff` for study design. We'll cover:\n",
11-
"\n",
12-
"1. Computing minimum detectable effects (MDE)\n",
13-
"2. Calculating required sample sizes\n",
14-
"3. Estimating statistical power\n",
15-
"4. Creating power curves for visualization\n",
16-
"5. Simulation-based power analysis for complex designs\n",
17-
"6. Panel data considerations (ICC, multiple periods)"
18-
]
7+
"source": "# Power Analysis for Difference-in-Differences\n\nThis notebook demonstrates how to use the power analysis tools in `diff-diff` for study design. We'll cover:\n\n1. Computing minimum detectable effects (MDE)\n2. Calculating required sample sizes\n3. Estimating statistical power\n4. Creating power curves for visualization\n5. Panel data considerations (ICC, multiple periods)\n6. Simulation-based power analysis for complex designs\n7. Power analysis for any estimator (staggered, synthetic DiD, triple difference)\n8. Finding MDE via simulation (bisection search)\n9. Finding required sample size via simulation (bisection search)\n10. Custom data generators\n11. Convenience functions\n12. Practical recommendations"
198
},
209
{
2110
"cell_type": "code",
2211
"execution_count": null,
2312
"id": "cell-1",
2413
"metadata": {},
2514
"outputs": [],
26-
"source": [
27-
"import numpy as np\n",
28-
"import pandas as pd\n",
29-
"import matplotlib.pyplot as plt\n",
30-
"\n",
31-
"from diff_diff import (\n",
32-
" PowerAnalysis,\n",
33-
" DifferenceInDifferences,\n",
34-
" simulate_power,\n",
35-
" compute_mde,\n",
36-
" compute_power,\n",
37-
" compute_sample_size,\n",
38-
" plot_power_curve,\n",
39-
")"
40-
]
15+
"source": "import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom diff_diff import (\n PowerAnalysis,\n DifferenceInDifferences,\n CallawaySantAnna,\n SyntheticDiD,\n TripleDifference,\n simulate_power,\n simulate_mde,\n simulate_sample_size,\n compute_mde,\n compute_power,\n compute_sample_size,\n plot_power_curve,\n)"
4116
},
4217
{
4318
"cell_type": "markdown",
@@ -475,15 +450,197 @@
475450
")"
476451
]
477452
},
453+
{
454+
"cell_type": "markdown",
455+
"id": "cjpvh2ze7lh",
456+
"source": "## 8. Power Analysis for Any Estimator\n\nThe simulation-based approach works with **all 12 supported estimators** — not just basic DiD. An internal registry automatically selects the appropriate data-generating process (DGP) and fit signature for each registered estimator. Just swap in the estimator object and everything else is handled. See the support table below for the full list, and Section 11 for using custom DGPs with unsupported estimators.\n\n### Staggered Adoption Estimators",
457+
"metadata": {}
458+
},
459+
{
460+
"cell_type": "code",
461+
"id": "la7ps3nxufq",
462+
"source": "# Power analysis with Callaway-Sant'Anna — the registry auto-selects\n# generate_staggered_data as the DGP and the correct fit kwargs\ncs = CallawaySantAnna()\n\ncs_results = simulate_power(\n estimator=cs,\n n_units=100,\n n_periods=6,\n treatment_effect=5.0,\n treatment_fraction=0.5,\n treatment_period=3,\n sigma=5.0,\n n_simulations=100,\n seed=42,\n progress=False,\n)\n\nprint(cs_results.summary())",
463+
"metadata": {},
464+
"execution_count": null,
465+
"outputs": []
466+
},
467+
{
468+
"cell_type": "markdown",
469+
"id": "q61cjchqjrd",
470+
"source": "### Factor Model Estimators (Synthetic DiD)\n\nFor `SyntheticDiD` with the default placebo variance method, the DGP must generate **more control than treated units** (`treatment_fraction < 0.5`). The registry uses `generate_factor_data` automatically.",
471+
"metadata": {}
472+
},
473+
{
474+
"cell_type": "code",
475+
"id": "x068rpe24gf",
476+
"source": "# Synthetic DiD — note treatment_fraction=0.3 (placebo variance requires\n# more control units than treated units)\nsdid = SyntheticDiD()\n\nsdid_results = simulate_power(\n estimator=sdid,\n n_units=60,\n n_periods=6,\n treatment_effect=5.0,\n treatment_fraction=0.3,\n treatment_period=3,\n sigma=3.0,\n n_simulations=100,\n seed=42,\n progress=False,\n)\n\nprint(sdid_results.summary())",
477+
"metadata": {},
478+
"execution_count": null,
479+
"outputs": []
480+
},
481+
{
482+
"cell_type": "markdown",
483+
"id": "6qpu05hi18s",
484+
"source": "### Triple Difference\n\n`TripleDifference` uses a fixed 2×2×2 factorial design (group × partition × time). Sample sizes are **rounded via `n_per_cell = max(2, n_units // 8)`**, so the minimum effective N is 16 (2 units per cell × 8 cells). The `effective_n_units` field in results tracks any rounding. Note that `simulate_sample_size()` uses a higher search floor of 64 from the registry.",
485+
"metadata": {}
486+
},
487+
{
488+
"cell_type": "code",
489+
"id": "91uackfwqp",
490+
"source": "# Triple Difference — n_units snaps to multiples of 8\nddd = TripleDifference()\n\nddd_results = simulate_power(\n estimator=ddd,\n n_units=64,\n treatment_effect=3.0,\n sigma=2.0,\n n_simulations=100,\n seed=42,\n progress=False,\n)\n\nprint(ddd_results.summary())\nif ddd_results.effective_n_units is not None:\n print(f\"\\nEffective N (after grid rounding): {ddd_results.effective_n_units}\")",
491+
"metadata": {},
492+
"execution_count": null,
493+
"outputs": []
494+
},
495+
{
496+
"cell_type": "markdown",
497+
"id": "6kb8ovmue4m",
498+
"source": "### Supported Estimators\n\nThe following 12 estimators are supported by the simulation power analysis registry. Each is automatically paired with the correct data-generating process:\n\n| DGP Family | Estimators | Min N |\n|---|---|---|\n| **Basic DiD** (`generate_did_data`) | DifferenceInDifferences, TwoWayFixedEffects, MultiPeriodDiD | 20 |\n| **Staggered** (`generate_staggered_data`) | CallawaySantAnna, SunAbraham, ImputationDiD, TwoStageDiD, StackedDiD, EfficientDiD | 40 |\n| **Factor Model** (`generate_factor_data`) | TROP, SyntheticDiD | 30 |\n| **Triple Difference** (`generate_ddd_data`) | TripleDifference | 16* |\n\n\\* DDD effective N rounds to `max(2, n_units // 8) * 8` with minimum 16. `simulate_sample_size()` uses a higher search floor of 64.\n\n> **Note:** `ContinuousDiD` is not in the registry because continuous/dose-response treatments require a different DGP structure. `BaconDecomposition` and `HonestDiD` are diagnostic/sensitivity tools rather than treatment effect estimators. For unsupported estimators, you can pass a custom `data_generator` and `result_extractor` (see Section 11).",
499+
"metadata": {}
500+
},
501+
{
502+
"cell_type": "markdown",
503+
"id": "c4erqwll1af",
504+
"source": "### Power Curve for a Staggered Estimator",
505+
"metadata": {}
506+
},
507+
{
508+
"cell_type": "code",
509+
"id": "ox3uab7h5bj",
510+
"source": "# Power curve across effect sizes for Callaway-Sant'Anna\ncs_curve = simulate_power(\n estimator=CallawaySantAnna(),\n n_units=100,\n n_periods=6,\n effect_sizes=[1.0, 2.0, 3.0, 5.0, 7.0],\n treatment_period=3,\n sigma=5.0,\n n_simulations=100,\n seed=42,\n progress=False,\n)\n\nplot_power_curve(\n cs_curve.power_curve_df(),\n target_power=0.80,\n title=\"CS Power Curve (100 units, 6 periods, SD=5)\",\n figsize=(10, 6),\n)",
511+
"metadata": {},
512+
"execution_count": null,
513+
"outputs": []
514+
},
515+
{
516+
"cell_type": "markdown",
517+
"id": "kqw6y4du5u",
518+
"source": "## 9. Finding MDE via Simulation\n\nThe analytical `PowerAnalysis.mde()` works for basic DiD, but for complex estimators there is no closed-form formula. `simulate_mde()` uses **bisection search** to find the minimum detectable effect: it repeatedly calls `simulate_power()` at different effect sizes, narrowing the bracket until it finds the smallest effect that achieves the target power.",
519+
"metadata": {}
520+
},
521+
{
522+
"cell_type": "code",
523+
"id": "p9a03aycu2",
524+
"source": "# Find MDE for basic DiD via simulation\nmde_result = simulate_mde(\n DifferenceInDifferences(),\n n_units=100,\n n_periods=4,\n sigma=5.0,\n n_simulations=100,\n seed=42,\n progress=False,\n)\n\nprint(mde_result.summary())",
525+
"metadata": {},
526+
"execution_count": null,
527+
"outputs": []
528+
},
529+
{
530+
"cell_type": "markdown",
531+
"id": "1mgrw7qmish",
532+
"source": "### Inspecting the Search Path\n\nThe `search_path` attribute records the effect size and power at each bisection step, which is useful for diagnosing convergence:",
533+
"metadata": {}
534+
},
535+
{
536+
"cell_type": "code",
537+
"id": "3p2tmxivi2g",
538+
"source": "# search_path is a List[Dict] — convert to DataFrame for display\nsearch_df = pd.DataFrame(mde_result.search_path)\nprint(search_df.to_string(index=False))",
539+
"metadata": {},
540+
"execution_count": null,
541+
"outputs": []
542+
},
543+
{
544+
"cell_type": "markdown",
545+
"id": "o4se5ofngjh",
546+
"source": "### MDE for a Staggered Estimator\n\nThe same function works with any registered estimator:",
547+
"metadata": {}
548+
},
549+
{
550+
"cell_type": "code",
551+
"id": "gwou7kv0ht9",
552+
"source": "# MDE for Callaway-Sant'Anna\ncs_mde = simulate_mde(\n CallawaySantAnna(),\n n_units=100,\n n_periods=6,\n treatment_period=3,\n sigma=5.0,\n n_simulations=100,\n seed=42,\n progress=False,\n)\n\nprint(cs_mde.summary())",
553+
"metadata": {},
554+
"execution_count": null,
555+
"outputs": []
556+
},
557+
{
558+
"cell_type": "markdown",
559+
"id": "ljyojxybdhc",
560+
"source": "**Key parameters for `simulate_mde()`:**\n- `effect_range=(lo, hi)` — custom search bracket (auto-detected if omitted)\n- `tol` — convergence tolerance on power (default 0.02)\n- `max_steps` — maximum bisection steps (default 15)\n- `n_simulations` — simulations per step (use 500+ for production analyses)",
561+
"metadata": {}
562+
},
563+
{
564+
"cell_type": "markdown",
565+
"id": "qrfhvtizg58",
566+
"source": "## 10. Finding Required Sample Size via Simulation\n\n`simulate_sample_size()` uses bisection search over `n_units` to find the smallest sample size that achieves the target power for a given effect size.",
567+
"metadata": {}
568+
},
569+
{
570+
"cell_type": "code",
571+
"id": "8o018e2dnl6",
572+
"source": "# Find required sample size for basic DiD\nn_result = simulate_sample_size(\n DifferenceInDifferences(),\n treatment_effect=5.0,\n sigma=5.0,\n n_simulations=100,\n seed=42,\n progress=False,\n)\n\nprint(n_result.summary())",
573+
"metadata": {},
574+
"execution_count": null,
575+
"outputs": []
576+
},
577+
{
578+
"cell_type": "markdown",
579+
"id": "gwtzwm3oedl",
580+
"source": "### Inspecting the Search Path",
581+
"metadata": {}
582+
},
583+
{
584+
"cell_type": "code",
585+
"id": "ubg3hqe9ypa",
586+
"source": "# View the bisection steps\nn_search_df = pd.DataFrame(n_result.search_path)\nprint(n_search_df.to_string(index=False))",
587+
"metadata": {},
588+
"execution_count": null,
589+
"outputs": []
590+
},
591+
{
592+
"cell_type": "markdown",
593+
"id": "dblnwe076kf",
594+
"source": "### Comparing Analytical and Simulation Results\n\nFor basic DiD, we can compare the simulation result against the analytical formula. With only 100 simulations per bisection step there will be Monte Carlo noise, so we expect **approximate** — not exact — agreement:",
595+
"metadata": {}
596+
},
597+
{
598+
"cell_type": "code",
599+
"id": "bgpin5xmsud",
600+
"source": "# Analytical sample size\nanalytical = pa.sample_size(effect_size=5.0, sigma=5.0)\n\nprint(f\"{'Method':<25} {'Required N':>12}\")\nprint(\"-\" * 40)\nprint(f\"{'Analytical:':<25} {analytical.required_n:>12}\")\nprint(f\"{'Simulation:':<25} {n_result.required_n:>12}\")\nprint(f\"\\nSimulation power at N: {n_result.power_at_n:.1%}\")",
601+
"metadata": {},
602+
"execution_count": null,
603+
"outputs": []
604+
},
605+
{
606+
"cell_type": "markdown",
607+
"id": "rnte1h09hra",
608+
"source": "**Key parameters for `simulate_sample_size()`:**\n- `n_range=(lo, hi)` — custom search bracket for sample size (auto-detected if omitted)\n- `max_steps` — maximum bisection steps (default 15)\n- `n_simulations` — simulations per step (use 500+ for production analyses)",
609+
"metadata": {}
610+
},
611+
{
612+
"cell_type": "markdown",
613+
"id": "usp5iwyacop",
614+
"source": "## 11. Custom Data Generators\n\nThe default DGPs cover common designs, but you can customize them in two ways:\n1. **Tweak the default DGP** with `data_generator_kwargs` (e.g., add multiple treatment cohorts)\n2. **Supply a fully custom DGP** with `data_generator`\n\n### Tweaking the Default DGP\n\nPass additional keyword arguments to the registry's DGP via `data_generator_kwargs`. For example, the default staggered DGP generates a single treatment cohort — here we create a multi-cohort design:\n\n> **Note:** Some keys are *protected* and cannot be overridden via `data_generator_kwargs` because they are controlled by the simulation function itself: `treatment_effect`, `noise_sd`, `n_units`, `n_periods`, `treatment_fraction`, `treatment_period`, `n_pre`, `n_post`.",
615+
"metadata": {}
616+
},
617+
{
618+
"cell_type": "code",
619+
"id": "igft0epkiic",
620+
"source": "# Multi-cohort staggered design: treatment starts at periods 2 and 4\ncs_multi = simulate_power(\n estimator=CallawaySantAnna(),\n n_units=120,\n n_periods=6,\n treatment_effect=5.0,\n sigma=5.0,\n n_simulations=100,\n seed=42,\n progress=False,\n data_generator_kwargs={\n \"cohort_periods\": [2, 4],\n \"never_treated_frac\": 0.3,\n },\n)\n\nprint(cs_multi.summary())",
621+
"metadata": {},
622+
"execution_count": null,
623+
"outputs": []
624+
},
625+
{
626+
"cell_type": "markdown",
627+
"id": "hznafqrhzoq",
628+
"source": "### Fully Custom Data Generator\n\nFor designs not covered by the built-in DGPs, supply your own `data_generator` function. It receives the standard simulation parameters and must return a DataFrame. You may also need a custom `result_extractor` if your estimator returns non-standard results, and `estimator_kwargs` to pass the right column names to `fit()`.",
629+
"metadata": {}
630+
},
631+
{
632+
"cell_type": "code",
633+
"id": "v06p7ubbj9p",
634+
"source": "def my_dgp(n_units, n_periods, treatment_effect, treatment_fraction,\n treatment_period, noise_sd, seed=None):\n \"\"\"Custom DGP with heterogeneous unit effects.\"\"\"\n rng = np.random.default_rng(seed)\n n_treat = int(n_units * treatment_fraction)\n\n rows = []\n for i in range(n_units):\n unit_fe = rng.normal(0, 3) # heterogeneous unit effect\n treated_unit = i < n_treat\n for t in range(n_periods):\n post = int(t >= treatment_period)\n effect = treatment_effect * post if treated_unit else 0.0\n y = unit_fe + 2.0 * t + effect + rng.normal(0, noise_sd)\n rows.append({\n \"unit\": i, \"period\": t, \"outcome\": y,\n \"ever_treated\": int(treated_unit), \"post\": post,\n })\n return pd.DataFrame(rows)\n\n# Use the custom DGP with simulate_power\ncustom_results = simulate_power(\n estimator=DifferenceInDifferences(),\n n_units=80,\n n_periods=4,\n treatment_effect=4.0,\n sigma=3.0,\n n_simulations=100,\n seed=42,\n progress=False,\n data_generator=my_dgp,\n estimator_kwargs={\"outcome\": \"outcome\", \"treatment\": \"ever_treated\", \"time\": \"post\"},\n)\n\nprint(custom_results.summary())",
635+
"metadata": {},
636+
"execution_count": null,
637+
"outputs": []
638+
},
478639
{
479640
"cell_type": "markdown",
480641
"id": "cell-28",
481642
"metadata": {},
482-
"source": [
483-
"## 8. Convenience Functions\n",
484-
"\n",
485-
"For quick calculations, use the convenience functions:"
486-
]
643+
"source": "## 12. Convenience Functions\n\nFor quick calculations, use the convenience functions:"
487644
},
488645
{
489646
"cell_type": "code",
@@ -509,18 +666,7 @@
509666
"cell_type": "markdown",
510667
"id": "cell-30",
511668
"metadata": {},
512-
"source": [
513-
"## 9. Practical Recommendations\n",
514-
"\n",
515-
"### Estimating Sigma (Residual SD)\n",
516-
"\n",
517-
"The residual standard deviation is crucial for power calculations. Options:\n",
518-
"\n",
519-
"1. **Pilot data**: Fit a model on historical data and get residual SD\n",
520-
"2. **Literature**: Find similar studies and use their reported SDs\n",
521-
"3. **Domain knowledge**: Expert judgment about outcome variability\n",
522-
"4. **Sensitivity analysis**: Calculate power for a range of sigma values"
523-
]
669+
"source": "## 13. Practical Recommendations\n\n### Estimating Sigma (Residual SD)\n\nThe residual standard deviation is crucial for power calculations. Options:\n\n1. **Pilot data**: Fit a model on historical data and get residual SD\n2. **Literature**: Find similar studies and use their reported SDs\n3. **Domain knowledge**: Expert judgment about outcome variability\n4. **Sensitivity analysis**: Calculate power for a range of sigma values"
524670
},
525671
{
526672
"cell_type": "code",
@@ -569,30 +715,17 @@
569715
" print(f\"{power:>10.0%} {result.required_n:>15}\")"
570716
]
571717
},
718+
{
719+
"cell_type": "markdown",
720+
"id": "eub4cew045u",
721+
"source": "### Analytical vs. Simulation: When to Use Each\n\n| Approach | Best for | Advantages |\n|---|---|---|\n| **Analytical** (`PowerAnalysis`) | Basic 2×2 DiD, panel DiD | Fast, exact, closed-form |\n| **Simulation** (`simulate_power/mde/sample_size`) | Staggered, SDID, TROP, DDD, custom designs | Works with any estimator, reports bias/RMSE/coverage |\n\n**Rule of thumb:** Start with analytical power analysis for basic designs. Move to simulation when using specialized estimators or non-standard DGPs.",
722+
"metadata": {}
723+
},
572724
{
573725
"cell_type": "markdown",
574726
"id": "cell-34",
575727
"metadata": {},
576-
"source": [
577-
"## Summary\n",
578-
"\n",
579-
"Key takeaways for DiD power analysis:\n",
580-
"\n",
581-
"1. **Always do a power analysis** before running a study\n",
582-
"2. **MDE decreases** with sample size, more periods, and lower variance\n",
583-
"3. **ICC matters** for panel data - high autocorrelation reduces effective sample size\n",
584-
"4. **Use simulation** for complex designs (staggered, synthetic DiD)\n",
585-
"5. **Be realistic about sigma** - err on the side of larger values\n",
586-
"6. **Consider your smallest meaningful effect** - don't just target statistical significance\n",
587-
"\n",
588-
"For more on DiD estimation, see the other tutorials:\n",
589-
"- `01_basic_did.ipynb` - Basic DiD estimation\n",
590-
"- `02_staggered_did.ipynb` - Staggered adoption designs\n",
591-
"- `03_synthetic_did.ipynb` - Synthetic DiD\n",
592-
"- `04_parallel_trends.ipynb` - Testing assumptions\n",
593-
"- `05_honest_did.ipynb` - Sensitivity analysis\n",
594-
"- `07_pretrends_power.ipynb` - Pre-trends power analysis (Roth 2022)"
595-
]
728+
"source": "## Summary\n\nKey takeaways for DiD power analysis:\n\n1. **Always do a power analysis** before running a study\n2. **MDE decreases** with sample size, more periods, and lower variance\n3. **ICC matters** for panel data — high autocorrelation reduces effective sample size\n4. **Use simulation** for complex designs (staggered, synthetic DiD, triple difference)\n5. **12 estimators are supported** out of the box via the auto-registry — just swap in the estimator\n6. **`simulate_mde()` and `simulate_sample_size()`** extend MDE and sample size calculations to any estimator via bisection search\n7. **Custom DGPs** let you model non-standard designs with `data_generator` and `data_generator_kwargs`\n8. **Be realistic about sigma** — err on the side of larger values\n9. **Consider your smallest meaningful effect** — don't just target statistical significance\n\nFor more on DiD estimation, see the other tutorials:\n- `01_basic_did.ipynb` — Basic DiD estimation\n- `02_staggered_did.ipynb` — Staggered adoption designs\n- `03_synthetic_did.ipynb` — Synthetic DiD\n- `04_parallel_trends.ipynb` — Testing assumptions\n- `05_honest_did.ipynb` — Sensitivity analysis\n- `07_pretrends_power.ipynb` — Pre-trends power analysis (Roth 2022)"
596729
}
597730
],
598731
"metadata": {
@@ -602,4 +735,4 @@
602735
},
603736
"nbformat": 4,
604737
"nbformat_minor": 5
605-
}
738+
}

0 commit comments

Comments
 (0)