{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Synthetic control with `pymc` models"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"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.33, 94% HDI [0.3, 0.38]\n",
" b 0.05, 94% HDI [0.011, 0.089]\n",
" c 0.3, 94% HDI [0.26, 0.35]\n",
" d 0.054, 94% HDI [0.01, 0.1]\n",
" e 0.024, 94% HDI [0.001, 0.065]\n",
" f 0.2, 94% HDI [0.11, 0.26]\n",
" g 0.038, 94% HDI [0.0028, 0.087]\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.72 \n",
" 0.212 \n",
" -2.117 \n",
" -1.35 \n",
" 0.006 \n",
" 0.004 \n",
" 1152.0 \n",
" 1649.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.61 \n",
" 6.35 \n",
" -63.52 \n",
" -40.51 \n",
" 0.186 \n",
" 0.105 \n",
" 1152.0 \n",
" 1649.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" median \n",
" hdi_lower \n",
" hdi_upper \n",
" p_gt_0 \n",
" relative_mean \n",
" relative_hdi_lower \n",
" relative_hdi_upper \n",
" \n",
" \n",
" \n",
" \n",
" average \n",
" -1.720330 \n",
" -1.705730 \n",
" -2.141285 \n",
" -1.344371 \n",
" 0.0 \n",
" -9.142522 \n",
" -11.144332 \n",
" -7.299538 \n",
" \n",
" \n",
" cumulative \n",
" -51.609893 \n",
" -51.171893 \n",
" -64.238552 \n",
" -40.331122 \n",
" 0.0 \n",
" -9.142522 \n",
" -11.144332 \n",
" -7.299538 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" median \n",
" hdi_lower \n",
" hdi_upper \n",
" p_two_sided \n",
" prob_of_effect \n",
" relative_mean \n",
" relative_hdi_lower \n",
" relative_hdi_upper \n",
" \n",
" \n",
" \n",
" \n",
" average \n",
" -3.057512 \n",
" -3.045187 \n",
" -3.430665 \n",
" -2.706604 \n",
" 0.0 \n",
" 1.0 \n",
" -18.635032 \n",
" -20.456322 \n",
" -16.867122 \n",
" \n",
" \n",
" cumulative \n",
" -48.920194 \n",
" -48.722991 \n",
" -54.890645 \n",
" -43.305657 \n",
" 0.0 \n",
" 1.0 \n",
" -18.635032 \n",
" -20.456322 \n",
" -16.867122 \n",
" \n",
" \n",
"
\n",
"