{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bayesian Interrupted Time Series\n",
"\n",
"Interrupted Time Series (ITS) analysis is a powerful approach for estimating the causal impact of an intervention or treatment when you have a single time series of observations. The key idea is to compare what actually happened after the intervention to what would have happened in the absence of the intervention (the \"counterfactual\"). To do this, we train a statistical model on the pre-intervention data (when no treatment has occurred) and then use this model to forecast the expected outcomes into the post-intervention period as-if treatment had not occurred. The difference between the observed outcomes and these model-based counterfactual predictions provides an estimate of the causal effect of the intervention, along with a measure of uncertainty if using a Bayesian approach."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What do we mean by _causal impact_ in Interrupted Time Series?\n",
"\n",
"In the context of interrupted time series analysis, the term **causal impact** refers to the estimated effect of an intervention or event on an outcome of interest. We can break this down into two components which tell us different aspects of the intervention's effect:\n",
"\n",
"- The **Instantaneous Bayesian Causal Effect** at each time point is the difference between the observed outcome and the model's posterior predictive distribution for the counterfactual (i.e., what would have happened without the intervention). This is not just a single number, but a full probability distribution that reflects our uncertainty.\n",
"- The **Cumulative Bayesian Causal Impact** is the sum of these instantaneous effects over the post-intervention period, again represented as a distribution.\n",
"\n",
"Let $y_t$ be the observed outcome at time $t$ (after the intervention), and $\\tilde{y}_t$ be the model's counterfactual prediction for the same time point. Then:\n",
"- **Instantaneous effect:** $\\Delta_t = y_t - \\tilde{y}_t$\n",
"- **Cumulative effect (up to time $T$):** $C_T = \\sum_{t=1}^T \\Delta_t$\n",
"\n",
"In Bayesian analysis, both $\\tilde{y}_t$ and $\\Delta_t$ are distributions, not just point estimates.\n",
"\n",
"### Why does this matter for decision making?\n",
"These metrics allow you to answer questions like:\n",
"- \"How much did the intervention change the outcome, compared to what would have happened otherwise?\"\n",
"- \"What is the total effect of the intervention over time, and how certain are we about it?\"\n",
"\n",
"By providing both instantaneous and cumulative effects, along with their uncertainty, you can make more informed decisions and better understand the impact of your interventions."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Interrupted Time Series example"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import arviz as az\n",
"import pandas as pd\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": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"image/png": {
"height": 811,
"width": 711
}
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = result.plot()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"==================================Pre-Post Fit==================================\n",
"Formula: y ~ 1 + t + C(month)\n",
"Model coefficients:\n",
" Intercept 23, 94% HDI [21, 24]\n",
" C(month)[T.2] 2.9, 94% HDI [0.95, 4.8]\n",
" C(month)[T.3] 1.2, 94% HDI [-0.82, 3.1]\n",
" C(month)[T.4] 7.1, 94% HDI [5.2, 9.2]\n",
" C(month)[T.5] 15, 94% HDI [13, 17]\n",
" C(month)[T.6] 25, 94% HDI [23, 27]\n",
" C(month)[T.7] 18, 94% HDI [16, 20]\n",
" C(month)[T.8] 33, 94% HDI [31, 36]\n",
" C(month)[T.9] 16, 94% HDI [14, 18]\n",
" C(month)[T.10] 9.2, 94% HDI [7.3, 11]\n",
" C(month)[T.11] 6.3, 94% HDI [4.4, 8.2]\n",
" C(month)[T.12] 0.58, 94% HDI [-1.3, 2.6]\n",
" t 0.21, 94% HDI [0.19, 0.23]\n",
" y_hat_sigma 2, 94% HDI [1.7, 2.3]\n"
]
}
],
"source": [
"result.summary()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"As well as the model coefficients, we might be interested in the estimated causal impact of the intervention over time - what we called the instantaneous Bayesian causal effect above. The post intervention causal impact estimates are contained in the `post_impact` attribute, which is of type `xarray.DataArray`. We can take a look at what this looks like, and we can see that it consists of 3 dimensions: `chain`, `draw`, and time (`obs_ind`). The `chain` and `draw` dimensions are used to store the samples from the posterior distribution, while the `obs_ind` dimension corresponds to the time points in the time series."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
""
],
"text/plain": [
" mean sd hdi_3% hdi_97%\n",
"x[unit_0] 65.711 21.026 27.114 106.002"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"az.summary(result.post_impact.sum(\"obs_ind\"), kind=\"stats\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, if we wanted to query the estimated impact over a specific time period, we can leverage the fact that this is an `xarray.DataArray` and we can use the `sel` method to select the time points we are interested in. We will leave this as an exercise for the reader.\n",
"\n",
"Moving on, it would also be possible to look at the mean post-intervention impact estimates, which would give us the average impact of the intervention over the post-intervention period. This can be done using the `mean` method and would equate to $\\bar{C_T} = \\Big[ \\sum_{t=1}^T \\Delta_t \\Big] / T$"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
mean
\n",
"
sd
\n",
"
hdi_3%
\n",
"
hdi_97%
\n",
"
\n",
" \n",
" \n",
"
\n",
"
x[unit_0]
\n",
"
1.825
\n",
"
0.584
\n",
"
0.753
\n",
"
2.944
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean sd hdi_3% hdi_97%\n",
"x[unit_0] 1.825 0.584 0.753 2.944"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"az.summary(result.post_impact.mean(\"obs_ind\"), kind=\"stats\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
":::{warning}\n",
"Care must be taken with the mean causal impact statistic $\\bar{C_T}$. It only makes sense to use this statistic if it looks like the intervention had a lasting (and roughly constant) effect on the outcome variable. If the effect is transient (like in the example here), then clearly there will be a lot of post-intervention period where the impact of the intervention has ‘worn off’. If so, then it will be hard to interpret the mean impacts real meaning.\n",
"\n",
"But if there was a roughly constant impact of the intervention over the post-intervention period, then the mean causal impact can be interpreted as the mean impact of the intervention (per time point) over the post-intervention period.\n",
":::"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Effect Summary Reporting\n",
"\n",
"For decision-making, you often need a concise summary of the causal effect with key statistics. The `effect_summary()` method provides a decision-ready report with:\n",
"\n",
"- Average and cumulative effects over a time window\n",
"- Highest Density Intervals (HDI) for uncertainty quantification\n",
"- Tail probabilities (e.g., P(effect > 0))\n",
"- Relative effects (% change vs counterfactual)\n",
"\n",
"This provides a comprehensive summary without manual post-processing.\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\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.825305
\n",
"
1.840825
\n",
"
0.656129
\n",
"
2.944493
\n",
"
0.99825
\n",
"
3.310732
\n",
"
1.161406
\n",
"
5.432045
\n",
"
\n",
"
\n",
"
cumulative
\n",
"
65.710976
\n",
"
66.269693
\n",
"
23.620628
\n",
"
106.001744
\n",
"
0.99825
\n",
"
3.310732
\n",
"
1.161406
\n",
"
5.432045
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean median hdi_lower hdi_upper p_gt_0 \\\n",
"average 1.825305 1.840825 0.656129 2.944493 0.99825 \n",
"cumulative 65.710976 66.269693 23.620628 106.001744 0.99825 \n",
"\n",
" relative_mean relative_hdi_lower relative_hdi_upper \n",
"average 3.310732 1.161406 5.432045 \n",
"cumulative 3.310732 1.161406 5.432045 "
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Generate effect summary for the full post-period\n",
"stats = result.effect_summary()\n",
"stats.table"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Post-period (2017-01-31 00:00:00 to 2019-12-31 00:00:00), the average effect was 1.83 (95% HDI [0.66, 2.94]), with a posterior probability of an increase of 0.998. The cumulative effect was 65.71 (95% HDI [23.62, 106.00]); probability of an increase 0.998. Relative to the counterfactual, this equals 3.31% on average (95% HDI [1.16%, 5.43%]).\n"
]
}
],
"source": [
"print(stats.text)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can customize the summary in several ways:\n",
"\n",
"- **Window**: Analyze a specific time period instead of the full post-period\n",
"- **Direction**: Specify whether you're testing for an increase, decrease, or two-sided effect\n",
"- **Alpha**: Set the HDI confidence level (default 95%)\n",
"- **Options**: Include/exclude cumulative or relative effects\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\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",
"
-1.582819
\n",
"
-1.570561
\n",
"
-2.641828
\n",
"
-0.511888
\n",
"
0.002
\n",
"
0.998
\n",
"
-3.191818
\n",
"
-5.324134
\n",
"
-1.155155
\n",
"
\n",
"
\n",
"
cumulative
\n",
"
-9.496912
\n",
"
-9.423368
\n",
"
-15.850971
\n",
"
-3.071330
\n",
"
0.002
\n",
"
0.998
\n",
"
-3.191818
\n",
"
-5.324134
\n",
"
-1.155155
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean median hdi_lower hdi_upper p_two_sided \\\n",
"average -1.582819 -1.570561 -2.641828 -0.511888 0.002 \n",
"cumulative -9.496912 -9.423368 -15.850971 -3.071330 0.002 \n",
"\n",
" prob_of_effect relative_mean relative_hdi_lower \\\n",
"average 0.998 -3.191818 -5.324134 \n",
"cumulative 0.998 -3.191818 -5.324134 \n",
"\n",
" relative_hdi_upper \n",
"average -1.155155 \n",
"cumulative -1.155155 "
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Example: Analyze first 6 months of post-period with two-sided test\n",
"post_dates = result.datapost.index\n",
"window_start = post_dates[0]\n",
"window_end = post_dates[5] # First 6 months\n",
"\n",
"stats_windowed = result.effect_summary(\n",
" window=(window_start, window_end),\n",
" direction=\"two-sided\",\n",
" alpha=0.05,\n",
" cumulative=True,\n",
" relative=True,\n",
")\n",
"stats_windowed.table"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Post-period (2017-01-31 00:00:00 to 2017-06-30 00:00:00), the average effect was -1.58 (95% HDI [-2.64, -0.51]), with a posterior probability of an effect of 0.998. The cumulative effect was -9.50 (95% HDI [-15.85, -3.07]); probability of an effect 0.998. Relative to the counterfactual, this equals -3.19% on average (95% HDI [-5.32%, -1.16%]).\n"
]
}
],
"source": [
"print(stats_windowed.text)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly, if we wanted, we would also be able to query the estimated cumulative impact of the intervention over the post-intervention period by instead looking at the `post_impact_cumulative` attribute, rather than the `post_impact` attribute."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "CausalPy",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.2"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "02f5385db19eab57520277c5168790c7855381ee953bdbb5c89c321e1f17586e"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}