{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 4: Inference and Standard Errors\n\nThe marginal effects we computed in previous tutorials came with standard errors and confidence intervals. See {doc}`Mathematical motivation </math>` for details on how standard errors are computed. This tutorial shows alternative approaches: robust standard errors, the Krinsky-Robb simulation method, bootstrap resampling, and adjustments for multiple comparisons.\n\n## What you will learn\n\n- How to use robust (HC3) standard errors\n- How to use Krinsky-Robb simulation for inference\n- How to use bootstrap resampling for inference\n- How to adjust confidence intervals for multiple comparisons (Bonferroni, Sidak, sup-t)\n\n## Setup\n\nWe continue with the same fitted model and `Margins` object:"
   ],
   "id": "4ffd0840"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "import numpy as np\nimport pandas as pd\nimport statsmodels.formula.api as smf\nfrom smmargins import Margins\n\nrng = np.random.default_rng(7)\nN = 5_000\ndf = pd.DataFrame({\n    \"age\":    rng.normal(45, 12, N).clip(18, 90),\n    \"income\": rng.lognormal(10.5, 0.4, N),\n    \"educ\":   rng.choice([\"hs\", \"college\", \"grad\"], N, p=[0.4, 0.4, 0.2]),\n    \"female\": rng.integers(0, 2, N),\n})\neta = (-4.0 + 0.05 * df[\"age\"] + 0.00001 * df[\"income\"]\n       + 0.8 * (df[\"educ\"] == \"college\") + 1.4 * (df[\"educ\"] == \"grad\")\n       + 0.3 * df[\"female\"] - 0.0004 * df[\"age\"] * df[\"female\"])\ndf[\"voted\"] = (rng.uniform(0, 1, N) < 1 / (1 + np.exp(-eta))).astype(int)\n\nfit = smf.logit(\"voted ~ age + income + C(educ) + female + age:female\", data=df).fit(disp=False)\nM = Margins(fit)",
   "id": "3510b123"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## Robust standard errors (HC3)\n\nBy default, smmargins uses the variance-covariance matrix from the fitted model. If you want robust standard errors, you can pass `cov_type=\"HC3\"` when creating the `Margins` object. HC3 is a heteroskedasticity-consistent estimator that is widely recommended for nonlinear models.",
   "id": "b35d6e29"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "M_robust = Margins(fit, cov_type=\"HC3\")\nM_robust.dydx(\"age\").summary()",
   "id": "811ed268"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "Compare this to the default standard error from Tutorial 3 (0.0005). The HC3 standard error (0.0006) is slightly larger, reflecting additional uncertainty from heteroskedasticity.\n\nYou can also pass `cov_type=\"HC0\"`, `cov_type=\"HC1\"`, or `cov_type=\"HC2\"` for alternative robust estimators.\n\n## Krinsky-Robb simulation\n\nThe Krinsky-Robb (KR) method draws random parameter vectors from the estimated sampling distribution, recomputes the quantity of interest for each draw, and uses the distribution of those recomputed values as the basis for inference. This does not rely on the delta method.\n\nTo use KR simulation, pass `vce=\"simulation\"` along with the number of simulations and an optional seed:",
   "id": "07097af4"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "M.dydx(\"age\", vce=\"simulation\", n_sims=2000, sim_seed=42).summary()",
   "id": "59bd4427"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "The `n_sims` parameter controls the number of simulated parameter draws. More simulations yield more stable standard errors but take longer to compute. The `sim_seed` ensures reproducibility.\n\nYou can use KR simulation with any marginal effect or prediction:",
   "id": "14c0ef43"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "M.dydx(\"female\", vce=\"simulation\", n_sims=2000, sim_seed=42).summary()",
   "id": "d3b99fa2"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## Bootstrap resampling\n\nBootstrap inference resamples the data with replacement, re-fits the model, and recomputes the quantity of interest for each resample. The variation across bootstrap replicates provides the standard error.\n\nTo use bootstrap, pass `vce=\"bootstrap\"`:",
   "id": "20f6bc18"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "M.dydx(\"age\", vce=\"bootstrap\", n_boot=500, boot_seed=42).summary()",
   "id": "6b8ea56d"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "Here `n_boot=500` specifies 500 bootstrap replicates. The `boot_seed` ensures reproducibility. Bootstrap is more computationally intensive than KR simulation because it requires re-fitting the model for each resample, but it makes fewer assumptions about the sampling distribution.\n\n## Multiple comparison adjustments\n\nWhen you compute predictions at multiple points, you are making multiple simultaneous inferences. By default, confidence intervals are pointwise (each interval has its own coverage probability). You may want simultaneous coverage: the probability that all intervals contain their true values.\n\nsmmargins offers three adjustment methods. Let us first get predictions at several ages:",
   "id": "fa88c3e4"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "pred_grid = M.predict(atexog={\"age\": [25, 45, 65]})\npred_grid.summary()",
   "id": "b99a64e4"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "### Bonferroni adjustment\n\nThe Bonferroni method multiplies each p-value by the number of comparisons and widens the confidence intervals accordingly:",
   "id": "fade5ce2"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "M.predict(\n    atexog={\"age\": [25, 45, 65]},\n    vce=\"simulation\",\n    n_sims=4000,\n    ci_method=\"bonferroni\"\n).summary()",
   "id": "ceb8c090"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "### Sidak adjustment\n\nThe Sidak method is slightly less conservative than Bonferroni:",
   "id": "93324238"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "M.predict(\n    atexog={\"age\": [25, 45, 65]},\n    vce=\"simulation\",\n    n_sims=4000,\n    ci_method=\"sidak\"\n).summary()",
   "id": "44f94b26"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "### Sup-t adjustment\n\nThe sup-t (supremum t) method accounts for the correlation between estimates when constructing simultaneous bands. It typically produces tighter intervals than Bonferroni while maintaining correct simultaneous coverage:",
   "id": "94ecb603"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "M.predict(\n    atexog={\"age\": [25, 45, 65]},\n    vce=\"simulation\",\n    n_sims=4000,\n    ci_method=\"sup-t\"\n).summary()",
   "id": "423c4c61"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "Compare the three methods. All three widen the intervals relative to the unadjusted case, but Bonferroni is the most conservative (widest), sup-t is the least conservative (narrowest), and Sidak falls in between.\n\n## Summary of inference options\n\n| Method | When to use |\n|--------|-------------|\n| Default (delta method) | Fast, works well for most cases |\n| `cov_type=\"HC3\"` | You suspect heteroskedasticity |\n| `vce=\"simulation\"` | You want to avoid delta-method approximations |\n| `vce=\"bootstrap\"` | You want inference robust to distributional assumptions |\n| `ci_method=\"bonferroni\"` | Conservative multiple comparisons |\n| `ci_method=\"sidak\"` | Slightly less conservative than Bonferroni |\n| `ci_method=\"sup-t\"` | Correlated estimates, tighter simultaneous bands |\n\n## Recap\n\nIn this tutorial we covered four approaches to inference:\n\n1. **Robust standard errors** via `cov_type=\"HC3\"` on the `Margins` object\n2. **Krinsky-Robb simulation** via `vce=\"simulation\"`\n3. **Bootstrap resampling** via `vce=\"bootstrap\"`\n4. **Multiple comparison adjustments** via `ci_method` (Bonferroni, Sidak, sup-t)\n\nThese options can be combined. For example, you can use KR simulation with Bonferroni adjustment, or bootstrap with Sidak adjustment.\n\n## Next steps\n\n- Apply inference methods to difference-in-differences in {doc}`Tutorial 5: Difference-in-Differences </tutorials/did>`\n- See the reference for {doc}`inference options </api>`\n- Learn about robust and clustered standard errors in {doc}`How-To: Robust and Clustered Standard Errors </howto/robust_clustered_ses>`\n- Learn about KR simulation in depth in {doc}`How-To: KR Simulation VCE </howto/kr_simulation>`\n- Learn about bootstrap inference in depth in {doc}`How-To: Bootstrap VCE </howto/bootstrap>`\n- Learn about simultaneous confidence intervals in {doc}`How-To: Simultaneous Confidence Intervals </howto/simultaneous_ci>`\n- Compare inference methods in {doc}`Explanation: Inference Comparison </explanations/inference_comparison>`",
   "id": "6e4b1278"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}