{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Synthetic control with `pymc` models"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n",
"/Users/benjamv/mambaforge/envs/CausalPy/lib/python3.13/site-packages/pymc_extras/model/marginal/graph_analysis.py:10: FutureWarning: `pytensor.graph.basic.io_toposort` was moved to `pytensor.graph.traversal.io_toposort`. Calling it from the old location will fail in a future release.\n",
" from pytensor.graph.basic import io_toposort\n"
]
}
],
"source": [
"import arviz as az\n",
"\n",
"import causalpy as cp"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2\n",
"%config InlineBackend.figure_format = 'retina'\n",
"seed = 42"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"df = cp.load_data(\"sc\")\n",
"treatment_time = 70"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Run the analysis\n",
"\n",
":::{note}\n",
"The `random_seed` keyword argument for the PyMC sampler is not necessary. We use it here so that the results are reproducible.\n",
":::"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" a \n",
" b \n",
" c \n",
" d \n",
" e \n",
" f \n",
" g \n",
" counterfactual \n",
" causal effect \n",
" actual \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0.793234 \n",
" 1.277264 \n",
" -0.055407 \n",
" -0.791535 \n",
" 1.075170 \n",
" 0.817384 \n",
" -2.607528 \n",
" 0.144888 \n",
" -0.0 \n",
" 0.398287 \n",
" \n",
" \n",
" 1 \n",
" 1.841898 \n",
" 1.185068 \n",
" -0.221424 \n",
" -1.430772 \n",
" 1.078303 \n",
" 0.890110 \n",
" -3.108099 \n",
" 0.601862 \n",
" -0.0 \n",
" 0.491644 \n",
" \n",
" \n",
" 2 \n",
" 2.867102 \n",
" 1.922957 \n",
" -0.153303 \n",
" -1.429027 \n",
" 1.432057 \n",
" 1.455499 \n",
" -3.149104 \n",
" 1.060285 \n",
" -0.0 \n",
" 1.232330 \n",
" \n",
" \n",
" 3 \n",
" 2.816255 \n",
" 2.424558 \n",
" 0.252894 \n",
" -1.260527 \n",
" 1.938960 \n",
" 2.088586 \n",
" -3.563201 \n",
" 1.520801 \n",
" -0.0 \n",
" 1.672995 \n",
" \n",
" \n",
" 4 \n",
" 3.865208 \n",
" 2.358650 \n",
" 0.311572 \n",
" -2.393438 \n",
" 1.977716 \n",
" 2.752152 \n",
" -3.515991 \n",
" 1.983661 \n",
" -0.0 \n",
" 1.775940 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {
"image/png": {
"height": 811,
"width": 711
}
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = result.plot(plot_predictors=True)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"================================SyntheticControl================================\n",
"Control units: ['a', 'b', 'c', 'd', 'e', 'f', 'g']\n",
"Treated unit: actual\n",
"Model coefficients:\n",
" a 0.34, 94% HDI [0.3, 0.38]\n",
" b 0.049, 94% HDI [0.01, 0.089]\n",
" c 0.3, 94% HDI [0.26, 0.35]\n",
" d 0.054, 94% HDI [0.0092, 0.099]\n",
" e 0.026, 94% HDI [0.0015, 0.068]\n",
" f 0.19, 94% HDI [0.1, 0.26]\n",
" g 0.039, 94% HDI [0.0028, 0.09]\n",
" y_hat_sigma 0.26, 94% HDI [0.22, 0.31]\n"
]
}
],
"source": [
"result.summary()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"As well as the model coefficients, we might be interested in the average causal impact and average cumulative causal impact."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" sd \n",
" hdi_3% \n",
" hdi_97% \n",
" mcse_mean \n",
" mcse_sd \n",
" ess_bulk \n",
" ess_tail \n",
" r_hat \n",
" \n",
" \n",
" \n",
" \n",
" x[actual] \n",
" -1.717 \n",
" 0.216 \n",
" -2.133 \n",
" -1.338 \n",
" 0.006 \n",
" 0.003 \n",
" 1260.0 \n",
" 1812.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" sd \n",
" hdi_3% \n",
" hdi_97% \n",
" mcse_mean \n",
" mcse_sd \n",
" ess_bulk \n",
" ess_tail \n",
" r_hat \n",
" \n",
" \n",
" \n",
" \n",
" x[actual] \n",
" -51.498 \n",
" 6.479 \n",
" -63.982 \n",
" -40.131 \n",
" 0.182 \n",
" 0.104 \n",
" 1260.0 \n",
" 1812.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"