diff --git a/causalpy/pymc_experiments.py b/causalpy/pymc_experiments.py index d7e006fc..31e6922d 100644 --- a/causalpy/pymc_experiments.py +++ b/causalpy/pymc_experiments.py @@ -12,7 +12,7 @@ """ import warnings # noqa: I001 -from typing import Union +from typing import Union, Dict import arviz as az import matplotlib.pyplot as plt @@ -32,7 +32,10 @@ IVDataValidator, ) from causalpy.plot_utils import plot_xY -from causalpy.utils import round_num +from causalpy.utils import ( + round_num, + compute_bayesian_tail_probability, +) LEGEND_FONT_SIZE = 12 az.style.use("arviz-darkgrid") @@ -321,7 +324,9 @@ def plot(self, counterfactual_label="Counterfactual", round_to=None, **kwargs): return fig, ax - def summary(self, round_to=None) -> None: + def summary( + self, round_to=None, version: str = "coefficients", **kwargs + ) -> Union[None, pd.DataFrame]: """ Print text output summarising the results @@ -329,10 +334,343 @@ def summary(self, round_to=None) -> None: Number of decimals used to round results. Defaults to 2. Use "None" to return raw numbers. """ - print(f"{self.expt_type:=^80}") - print(f"Formula: {self.formula}") - # TODO: extra experiment specific outputs here - self.print_coefficients(round_to) + if version == "coefficients": + print(f"{self.expt_type:=^80}") + print(f"Formula: {self.formula}") + # TODO: extra experiment specific outputs here + self.print_coefficients() + elif version == "intervention": + return self._summary_intervention(**kwargs) + + def _power_estimation(self, alpha: float = 0.05, correction: bool = False) -> Dict: + """ + Estimate the statistical power of an intervention based on cumulative and mean results. + This function calculates posterior estimates, systematic differences, credible intervals, and + minimum detectable effects (MDE) for both cumulative and mean measures. It can apply corrections to + account for systematic differences in the data if the mean pre-intervention is consider far from the + real mean. + + Parameters + ---------- + - alpha (float, optional): The significance level for confidence interval calculations. + Should be a fraction between 0 and 1, not a percentage between 0 and 100. Default is 0.05. + - correction (bool, optional): If True, applies corrections to account for systematic differences in + cumulative and mean calculations. Default is False. + + Returns + ------- + - Dict: A dictionary containing key statistical measures such as posterior estimation, + systematic differences, credible intervals, and posterior MDE for both cumulative and mean results. + """ + assert 0 <= alpha <= 1, "Alpha must be in the range [0, 1]." + + if not isinstance(correction, bool): + raise ValueError("Correction must be either True or False.") + + results = {} + ci = (alpha * 100) / 2 + # Cumulative calculations + cumulative_results = self.post_y.sum() + _mu_samples_cumulative = ( + self.post_pred["posterior_predictive"] + .mu.stack(sample=("chain", "draw")) + .sum("obs_ind") + ) + # Mean calculations + mean_results = self.post_y.mean() + _mu_samples_mean = ( + self.post_pred["posterior_predictive"] + .mu.stack(sample=("chain", "draw")) + .mean("obs_ind") + ) + # Posterior Mean + results["posterior_estimation"] = { + "cumulative": np.mean(_mu_samples_cumulative.values), + "mean": np.mean(_mu_samples_mean.values), + } + results["results"] = { + "cumulative": cumulative_results, + "mean": mean_results, + } + results["_systematic_differences"] = { + "cumulative": results["results"]["cumulative"] + - results["posterior_estimation"]["cumulative"], + "mean": results["results"]["mean"] + - results["posterior_estimation"]["mean"], + } + if correction: + _mu_samples_cumulative += results["_systematic_differences"]["cumulative"] + _mu_samples_mean += results["_systematic_differences"]["mean"] + results["ci"] = { + "cumulative": [ + np.percentile(_mu_samples_cumulative, ci), + np.percentile(_mu_samples_cumulative, 100 - ci), + ], + "mean": [ + np.percentile(_mu_samples_mean, ci), + np.percentile(_mu_samples_mean, 100 - ci), + ], + } + cumulative_upper_mde = ( + results["ci"]["cumulative"][1] + - results["posterior_estimation"]["cumulative"] + ) + cumulative_lower_mde = ( + results["posterior_estimation"]["cumulative"] + - results["ci"]["cumulative"][0] + ) + mean_upper_mde = ( + results["ci"]["mean"][1] - results["posterior_estimation"]["mean"] + ) + mean_lower_mde = ( + results["posterior_estimation"]["mean"] - results["ci"]["mean"][0] + ) + results["posterior_mde"] = { + "cumulative": (cumulative_upper_mde + cumulative_lower_mde) / 2, + "mean": (mean_upper_mde + mean_lower_mde) / 2, + } + return results + + def _summary_intervention(self, alpha: float = 0.05, **kwargs) -> pd.DataFrame: + """ + Calculate and summarize the intervention analysis results in a DataFrame format. + + This function performs cumulative and mean calculations on the posterior predictive distributions, + computes Bayesian tail probabilities, posterior estimations, causal effects, and credible intervals. + It optionally applies corrections to the cumulative and mean calculations. + + Parameters + ---------- + alpha (float, optional): The significance level for confidence interval calculations. Default is 0.05. + kwargs (Dict[str, Any], optional): Additional keyword arguments. + - "correction" (bool | Dict[str, float] | pd.Series): If True, applies predefined corrections to cumulative and mean results. + If a dictionary, the corrections for 'cumulative' and 'mean' should be provided. Default is False. + + Returns + ------- + - pd.DataFrame: A DataFrame where each row represents different statistical measures such as + Bayesian tail probability, posterior estimation, causal effect, and credible intervals for cumulative and mean results. + """ + + correction = kwargs.get("correction", False) + + results = {} + ci = (alpha * 100) / 2 + + # Cumulative calculations + cumulative_results = self.post_y.sum() + _mu_samples_cumulative = ( + self.post_pred["posterior_predictive"] + .mu.stack(sample=("chain", "draw")) + .sum("obs_ind") + ) + + # Mean calculations + mean_results = self.post_y.mean() + _mu_samples_mean = ( + self.post_pred["posterior_predictive"] + .mu.stack(sample=("chain", "draw")) + .mean("obs_ind") + ) + + if not isinstance(correction, bool): + if not isinstance(correction, (pd.Series, Dict)): + raise ValueError( + "Correction must be a Pandas series (`pd.Series`), Dictionary." + ) + elif isinstance(correction, pd.Series): + if set(correction.index) != {"cumulative", "mean"}: + raise ValueError( + "Correction index must have ['cumulative', 'mean'] values." + ) + if not all( + isinstance(value, (float, int)) and not isinstance(value, bool) + for value in correction.values + ): + raise ValueError( + "All values in the correction Pandas Series must be integers or floats, and not boolean." + ) + + elif isinstance(correction, Dict): + if set(correction.keys()) != {"cumulative", "mean"}: + raise ValueError( + "Correction dictionary must have keys ['cumulative', 'mean']." + ) + if not all( + isinstance(value, (float, int)) and not isinstance(value, bool) + for value in correction.values() + ): + raise ValueError( + "All values in the correction dictionary must be integers or floats, and not boolean." + ) + + _mu_samples_cumulative += correction["cumulative"] + _mu_samples_mean += correction["mean"] + + # Bayesian Tail Probability + results["bayesian_tail_probability"] = { + "cumulative": compute_bayesian_tail_probability( + posterior=_mu_samples_cumulative, x=cumulative_results + ), + "mean": compute_bayesian_tail_probability( + posterior=_mu_samples_mean, x=mean_results + ), + } + + # Posterior Mean + results["posterior_estimation"] = { + "cumulative": np.mean(_mu_samples_cumulative.values), + "mean": np.mean(_mu_samples_mean.values), + } + + results["results"] = { + "cumulative": cumulative_results, + "mean": mean_results, + } + + # Causal Effect + results["causal_effect"] = { + "cumulative": cumulative_results + - results["posterior_estimation"]["cumulative"], + "mean": mean_results - results["posterior_estimation"]["mean"], + } + + # credible intervals + results["ci"] = { + "cumulative": [ + np.percentile(_mu_samples_cumulative, ci), + np.percentile(_mu_samples_cumulative, 100 - ci), + ], + "mean": [ + np.percentile(_mu_samples_mean, ci), + np.percentile(_mu_samples_mean, 100 - ci), + ], + } + + # Convert to DataFrame + results_df = pd.DataFrame(results) + + return results_df + + def power_summary( + self, alpha: float = 0.05, correction: bool = False + ) -> pd.DataFrame: + """ + Summarize the power estimation results in a DataFrame format. + + This function perform a power estimation and then + converts the resulting dictionary into a pandas DataFrame for easier analysis and visualization. + + Parameters + ---------- + - alpha (float, optional): The significance level for confidence interval calculations used in power estimation. Default is 0.05. + - correction (bool, optional): Indicates whether to apply corrections in the power estimation process. Default is False. + + Returns + ------- + - pd.DataFrame: A DataFrame representing the power estimation results, including posterior estimations, + systematic differences, credible intervals, and posterior MDE for cumulative and mean results. + + References + ---------- + https://causalpy.readthedocs.io/en/latest/notebooks/sc_power_analysis.html + """ + warnings.warn( + "The power function is experimental and the API may change in the future." + ) + return pd.DataFrame(self._power_estimation(alpha=alpha, correction=correction)) + + def plot_power(self, alpha: float = 0.05, correction: bool = False) -> plt.Figure: + """ + Generate and return a figure containing plots that visualize power estimation results. + + This function creates a two-panel plot (for mean and cumulative measures) to visualize the posterior distributions + along with the credible intervals, real mean, and posterior mean values. It allows for adjustments based on + systematic differences if the correction is applied. + + Parameters + ---------- + - alpha (float, optional): The significance level for confidence interval calculations used in power estimation. Default is 0.05. + - correction (bool, optional): Indicates whether to apply corrections for systematic differences in the plotting process. Default is False. + + Returns + ------- + - plt.Figure: A matplotlib figure object containing the plots. + + References + ---------- + https://causalpy.readthedocs.io/en/latest/notebooks/sc_power_analysis.html + """ + if not isinstance(correction, bool): + raise ValueError("Correction must be either True or False.") + + _estimates = self._power_estimation(alpha=alpha, correction=correction) + + fig, axs = plt.subplots(1, 2, figsize=(20, 6)) # Two subplots side by side + + # Adjustments for Mean and Cumulative plots + for i, key in enumerate(["mean", "cumulative"]): + _mu_samples = self.post_pred["posterior_predictive"].mu.stack( + sample=("chain", "draw") + ) + if key == "mean": + _mu_samples = _mu_samples.mean("obs_ind") + elif key == "cumulative": + _mu_samples = _mu_samples.sum("obs_ind") + + if correction: + _mu_samples += _estimates["_systematic_differences"][key] + + # Histogram and KDE + sns.histplot( + _mu_samples, + bins=30, + kde=True, + ax=axs[i], + color="C0", + stat="density", + alpha=0.6, + ) + kde_x, kde_y = ( + sns.kdeplot(_mu_samples, color="C1", fill=True, ax=axs[i]) + .get_lines()[0] + .get_data() + ) + + # Adjust y-limits based on max density + max_density = max(kde_y) + axs[i].set_ylim(0, max_density + 0.05 * max_density) # Adding 5% buffer + + # Fill between for the percentile interval + axs[i].fill_betweenx( + y=np.linspace(0, max_density + 0.05 * max_density, 100), + x1=_estimates["ci"][key][0], + x2=_estimates["ci"][key][1], + color="C0", + alpha=0.3, + label="C.I", + ) + + # Vertical lines for the means + axs[i].axvline( + _estimates["results"][key], color="C3", linestyle="-", label="Real Mean" + ) + if not correction: + axs[i].axvline( + _estimates["posterior_estimation"][key], + color="C4", + linestyle="--", + label="Posterior Mean", + ) + + axs[i].set_title(f"Posterior of mu ({key.capitalize()})") + axs[i].set_xlabel("mu") + axs[i].set_ylabel("Density") + axs[i].legend() + warnings.warn( + "The power function is experimental and the API may change in the future." + ) + return fig, axs class InterruptedTimeSeries(PrePostFit): diff --git a/causalpy/tests/test_utils.py b/causalpy/tests/test_utils.py index 34af4252..1756442a 100644 --- a/causalpy/tests/test_utils.py +++ b/causalpy/tests/test_utils.py @@ -2,9 +2,15 @@ Tests for utility functions """ +import numpy as np import pandas as pd -from causalpy.utils import _is_variable_dummy_coded, _series_has_2_levels, round_num +from causalpy.utils import ( + _is_variable_dummy_coded, + _series_has_2_levels, + compute_bayesian_tail_probability, + round_num, +) def test_dummy_coding(): @@ -44,3 +50,45 @@ def test_round_num(): assert round_num(123.456, 5) == "123.46" assert round_num(123.456, 6) == "123.456" assert round_num(123.456, 7) == "123.456" + + +def test_compute_bayesian_tail_probability(): + """ + Re-running all tests for the compute_bayesian_tail_probability function with the corrected understanding + and expectations for various scenarios. + """ + # Test 1: Posterior is a standard normal distribution, x = mean = 0 + posterior_standard_normal = np.random.normal(0, 1, 10000) + x_at_mean = 0 + prob_at_mean = compute_bayesian_tail_probability( + posterior_standard_normal, x_at_mean + ) + assert np.isclose(prob_at_mean, 1, atol=0.05), f"Expected 1, got {prob_at_mean}" + + # Test 2: Posterior is a standard normal distribution, x = 1 + x_one_std_above = 1 + prob_one_std_above = compute_bayesian_tail_probability( + posterior_standard_normal, x_one_std_above + ) + assert ( + 0 < prob_one_std_above < 1 + ), "Probability should decrease from 1 as x moves away from mean" + + # Test 3: Posterior is a standard normal distribution, x well outside the distribution + x_far_out = 5 + prob_far_out = compute_bayesian_tail_probability( + posterior_standard_normal, x_far_out + ) + # Expect a very low probability for a value far outside the distribution + assert prob_far_out < 0.01, f"Expected a value < 0.01, got {prob_far_out}" + + # Test 4: Posterior is a normal distribution with mean=5, std=2, x = mean + posterior_shifted = np.random.normal(5, 2, 10000) + x_at_shifted_mean = 5 + prob_at_shifted_mean = compute_bayesian_tail_probability( + posterior_shifted, x_at_shifted_mean + ) + # Expect the probability at the mean of a shifted distribution to be close to 1 + assert np.isclose( + prob_at_shifted_mean, 1, atol=0.05 + ), f"Expected 1, got {prob_at_shifted_mean}" diff --git a/causalpy/utils.py b/causalpy/utils.py index 2b794503..48c43172 100644 --- a/causalpy/utils.py +++ b/causalpy/utils.py @@ -3,6 +3,7 @@ """ import numpy as np import pandas as pd +from scipy.stats import norm def _is_variable_dummy_coded(series: pd.Series) -> bool: @@ -48,3 +49,29 @@ def _format_sig_figs(value, default=None): if value == 0: return 1 return max(int(np.log10(np.abs(value))) + 1, default) + + +def compute_bayesian_tail_probability(posterior, x) -> float: + """ + Calculate the probability of a given value being in a distribution defined by the posterior, + + Args: + - data: a list or array-like object containing the data to define the distribution + - x: a numeric value for which to calculate the probability of being in the distribution + + Returns: + - prob: a numeric value representing the probability of x being in the distribution + """ + lower_bound, upper_bound = min(posterior), max(posterior) + mean, std = np.mean(posterior), np.std(posterior) + + cdf_lower = norm.cdf(lower_bound, mean, std) + cdf_upper = 1 - norm.cdf(upper_bound, mean, std) + cdf_x = norm.cdf(x, mean, std) + + if cdf_x <= 0.5: + probability = 2 * (cdf_x - cdf_lower) / (1 - cdf_lower - cdf_upper) + else: + probability = 2 * (1 - cdf_x + cdf_lower) / (1 - cdf_lower - cdf_upper) + + return abs(probability) diff --git a/docs/source/_static/classes.png b/docs/source/_static/classes.png index 5842bfff..f437e614 100644 Binary files a/docs/source/_static/classes.png and b/docs/source/_static/classes.png differ diff --git a/docs/source/_static/interrogate_badge.svg b/docs/source/_static/interrogate_badge.svg index 9a7d6c00..427eebeb 100644 --- a/docs/source/_static/interrogate_badge.svg +++ b/docs/source/_static/interrogate_badge.svg @@ -1,58 +1,58 @@ - - interrogate: 96.6% - - - - - - - - - - - interrogate - interrogate - 96.6% - 96.6% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + interrogate: 96.7% + + + + + + + + + + + interrogate + interrogate + 96.7% + 96.7% + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/_static/packages.png b/docs/source/_static/packages.png index 9713e0b7..5dbd4d91 100644 Binary files a/docs/source/_static/packages.png and b/docs/source/_static/packages.png differ diff --git a/docs/source/examples.rst b/docs/source/examples.rst index e4850d58..577dde2d 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -17,6 +17,7 @@ Synthetic Control notebooks/sc_skl.ipynb notebooks/sc_pymc_brexit.ipynb notebooks/geolift1.ipynb + notebooks/power_analysis.ipynb Difference in Differences diff --git a/docs/source/notebooks/sc_power_analysis.ipynb b/docs/source/notebooks/sc_power_analysis.ipynb new file mode 100644 index 00000000..f9663428 --- /dev/null +++ b/docs/source/notebooks/sc_power_analysis.ipynb @@ -0,0 +1,22043 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e57b43b7", + "metadata": {}, + "source": [ + "# From Sensitivity to Power: A Methodological Framework to Optimize experimental design in Synthetic Control\n", + "\n", + "Our addition to `causalpy` is aimed at improving the tool's usefulness for experimenters in the planning and validation stages of quasi-experiments. We have introduced a sensitivity analysis method in a Bayesian framework that provides a systematic approach to determine the power of a model in detecting expected effects. We will explore practical applications, comprehend the underlying principles, and discuss how this approach can improve our decision-making process in selecting models and designing experiments.\n", + "\n", + "This notebook is dedicated to exploring and implementing Bayesian Power Analysis in the context of `causalpy`.\n", + "\n", + "## How it Works\n", + "The method involves creating a null model that does not capture any effect during a period where no effect is present. By analyzing the posterior distribution derived from this null model, experimenters can estimate the magnitude of effect necessary for it to be considered significant. This estimation allows for an assessment of the model's sensitivity to changes and the experiment's feasibility.\n", + "\n", + "**Application during the pre-experimentation phase!**\n", + "By applying this method before the experiment period we will be able to determine what is the setup of our most optimal model to reduce our MDE and increase the power. Using this methon we can answer questions like:\n", + "1. Should we generate more samples? \n", + "2. Should we increase the number of chains? \n", + "3. What are the best set of regressors?\n", + "4. What is the best number of observations to train?\n", + "\n", + "Our proposed power (sensitivity) analysis method is designed to be universally applicable across different regression models in `causalpy` such as **synthetic controls** and **interrupted time series**. By assessing the null model's posterior, we can validate that our regression does not inadvertently capture effects during a control period. This validation process not only aids in determining the required effect size for significance but also helps in evaluating the natural bias of the model, thus ensuring more reliable and accurate experimental planning and analysis.\n", + "\n", + "\n", + "### Similarities and Differences to Frequentist Methods\n", + "\n", + "Both Bayesian and frequentist methods aim to provide insights into the effectiveness of interventions or treatments. Our power analysis method in the Bayesian context differs from traditional frequentist approaches by focusing on the probability distributions of outcomes instead of point estimates and p-values. Frequentist methods rely on p-values and confidence intervals to reject or fail to reject a null hypothesis. In contrast, Bayesian approaches use the posterior distribution to estimate the probability of various outcomes, providing more nuanced insights into the expected effects and giving the observer the possibility to determine their significance based on the risk of false positives." + ] + }, + { + "cell_type": "markdown", + "id": "0d26c308", + "metadata": {}, + "source": [ + "## Getting Started\n", + "\n", + "Before diving into the implementation, let's set up our working environment. This setup ensures we have all the necessary tools at our disposal for conducting Bayesian Power Analysis using `causalpy`. \n", + "\n", + "1. **Installation**: First, we'll install the `causalpy` package, which is essential for our analysis. Simply run `%pip install causalpy` in your Jupyter notebook.\n", + "\n", + "2. **Importing Libraries**: Next, we import key libraries:\n", + " - `causalpy` for statistical modeling and analysis.\n", + "\n", + "3. **Environment Setup**:\n", + " - We use `%load_ext autoreload` and `%autoreload 2` for automatic reloading of modules before executing a new line of code, which is helpful during development.\n", + " - The `seed` variable is set to ensure reproducibility of our results.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "id": "6adee44a", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "import causalpy as cp" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "id": "bf91ea16", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "%config InlineBackend.figure_format = 'svg'\n", + "\n", + "seed = 42\n", + "pd.set_option('display.precision', 2)" + ] + }, + { + "cell_type": "markdown", + "id": "d3151407", + "metadata": {}, + "source": [ + "## Loading Data\n", + "\n", + "The dataset is loaded using `causalpy`'s built-in data loading functionality, which provides a suitable dataset for our quasi-experimental design.\n", + "\n", + "Once we've loaded the data, our first step is to identify the point in time where the treatment or intervention is expected to occur. In this case, we've set a specific time (`actual_treatment_time`) to represent when the intervention is expected to take place. Following this, we filter the dataset to exclude any data beyond this treatment time to reflect the pre-intervention period, which aligns with the realistic scenario of not having post-intervention data available during the planning phase of an experiment.\n", + "\n", + "By focusing on the pre-intervention period, we aim to understand the dataset's characteristics and establish a baseline to measure the expected effects of our intervention. This baseline is crucial for determining the magnitude of effect required to be detectable in our framework." + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "id": "a4508a65", + "metadata": {}, + "outputs": [], + "source": [ + "df = cp.load_data(\"sc\")\n", + "actual_treatment_time = 70\n", + "fake_treatment_time = 60 #Period before interventation to check the power of your model" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "id": "5a687a4c", + "metadata": {}, + "outputs": [], + "source": [ + "power_df = df[:actual_treatment_time].copy()" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "cc5be440", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "RangeIndex: 70 entries, 0 to 69\n", + "Data columns (total 10 columns):\n", + " # Column Non-Null Count Dtype \n", + "--- ------ -------------- ----- \n", + " 0 a 70 non-null float64\n", + " 1 b 70 non-null float64\n", + " 2 c 70 non-null float64\n", + " 3 d 70 non-null float64\n", + " 4 e 70 non-null float64\n", + " 5 f 70 non-null float64\n", + " 6 g 70 non-null float64\n", + " 7 counterfactual 70 non-null float64\n", + " 8 causal effect 70 non-null float64\n", + " 9 actual 70 non-null float64\n", + "dtypes: float64(10)\n", + "memory usage: 5.6 KB\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
abcdefgcounterfactualcausal effectactual
00.791.28-0.06-0.791.080.82-2.610.14-0.00.40
11.841.19-0.22-1.431.080.89-3.110.60-0.00.49
22.871.92-0.15-1.431.431.46-3.151.06-0.01.23
32.822.420.25-1.261.942.09-3.561.52-0.01.67
43.872.360.31-2.391.982.75-3.521.98-0.01.78
\n", + "
" + ], + "text/plain": [ + " a b c d e f g counterfactual causal effect \\\n", + "0 0.79 1.28 -0.06 -0.79 1.08 0.82 -2.61 0.14 -0.0 \n", + "1 1.84 1.19 -0.22 -1.43 1.08 0.89 -3.11 0.60 -0.0 \n", + "2 2.87 1.92 -0.15 -1.43 1.43 1.46 -3.15 1.06 -0.0 \n", + "3 2.82 2.42 0.25 -1.26 1.94 2.09 -3.56 1.52 -0.0 \n", + "4 3.87 2.36 0.31 -2.39 1.98 2.75 -3.52 1.98 -0.0 \n", + "\n", + " actual \n", + "0 0.40 \n", + "1 0.49 \n", + "2 1.23 \n", + "3 1.67 \n", + "4 1.78 " + ] + }, + "execution_count": 50, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "power_df.info()\n", + "power_df.head()" + ] + }, + { + "cell_type": "markdown", + "id": "9983a6a8", + "metadata": {}, + "source": [ + "## Creating the Model for Power Analysis\n", + "\n", + "Let's move on to creating our model using the `causalpy` library. Our goal is to predict the period right before the intervention occurs, and to do so, we use a specific time point as a reference for our predictions. The idea is to capture the true nature of the data before any intervention effects come into play.\n", + "\n", + "For our purposes, the Synthetic Control method from `causalpy` fits the bill perfectly. If our model is well-calibrated, we should see the difference between the real data and the estimated inference approach zero, indicating that we have accurately captured the with our method the pre-intervention scenario." + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "id": "58bf5b47", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (4 chains in 4 jobs)\n", + "NUTS: [beta, sigma]\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + " 100.00% [8000/8000 00:01<00:00 Sampling 4 chains, 0 divergences]\n", + "
\n", + " " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.\n", + "Sampling: [beta, sigma, y_hat]\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n" + ] + } + ], + "source": [ + "result = cp.pymc_experiments.SyntheticControl(\n", + " power_df,\n", + " fake_treatment_time,\n", + " formula=\"actual ~ 0 + a + b + c\",\n", + " model=cp.pymc_models.WeightedSumFitter(\n", + " sample_kwargs={\"target_accept\": 0.95, \"random_seed\": seed}\n", + " ),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "95aa16e7", + "metadata": {}, + "source": [ + "## Evaluating the model\n", + "\n", + "Essentially, if the model accurately represents our reality, and we'll base our decision on the model intervals. By grasping these intervals, we can ascertain how significant an effect must be. In turn, this will allow us to determine the size of the effect that will be considered significant.\n", + "\n", + "In order to ensure that our model isn't overestimating the positive or negative effects, it's important to expect that the errors will be symmetrically distributed around zero. This means that the true and predicted means are aligned. If the model is failing systematically towards either the positive or negative sides, it indicates that it's not accurately capturing the full range of the data." + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "id": "df967213", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-03-10T22:37:34.653005\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = result.plot(plot_predictors=True)" + ] + }, + { + "cell_type": "markdown", + "id": "f19a0ae1", + "metadata": {}, + "source": [ + "In this case, we can see that the model based on these regressors is systematically failing and although there is still no intervention, it is capturing an effect. But how big is this effect? It seems relevant from the perspective of the images but how to quantify it?" + ] + }, + { + "cell_type": "markdown", + "id": "a655a236", + "metadata": {}, + "source": [ + "## Power Analysis\n", + "\n", + "We can employ the `power_summary` function to observe the actual cumulative value during the test period and compare it with the projected cumulative outcomes for that same period. Ideally, we should expect errors to be near zero; if they are significantly negative or positive, our model will tend to err consistently in one direction.\n", + "\n", + "Furthermore, we can view the quantiles and determine the minimum detectable effect (MDE) required for the effect to be within those quantiles. Essentially, if your action produces an incremental value that is equal to or greater than the MDE, then your effect will be outside the quantile range and could be considered substantial.\n", + "\n", + "The alpha parameter plays a crucial role in deciding the cumulative confidence interval (CCI). In simple words, the alpha value determines the width of the confidence interval. For instance, if you set an alpha value of `0.1`, it implies that the interval is **90%**. Similarly, an alpha value of `0.05` corresponds to a **95%** interval, and so forth." + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "id": "18692ec0", + "metadata": {}, + "outputs": [], + "source": [ + "alpha=0.1 #" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "id": "5fd7b620", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/Documents/GitHub/CausalPy/causalpy/pymc_experiments.py:573: UserWarning: The power function is experimental and the API may change in the future.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
posterior_estimationresults_systematic_differencesciposterior_mde
cumulative114.31111.43-2.87[112.430854649129, 116.1599388110477]1.86
mean11.4311.14-0.29[11.2430854649129, 11.61599388110477]0.19
\n", + "
" + ], + "text/plain": [ + " posterior_estimation results _systematic_differences \\\n", + "cumulative 114.31 111.43 -2.87 \n", + "mean 11.43 11.14 -0.29 \n", + "\n", + " ci posterior_mde \n", + "cumulative [112.430854649129, 116.1599388110477] 1.86 \n", + "mean [11.2430854649129, 11.61599388110477] 0.19 " + ] + }, + "execution_count": 54, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result.power_summary(alpha=alpha)" + ] + }, + { + "cell_type": "markdown", + "id": "f044d6b9", + "metadata": {}, + "source": [ + "We could appreciate this better in an image, and for that we can use the `power_plot` function. This would visually display the information presented by `power_summary`." + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "db112259", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/Documents/GitHub/CausalPy/causalpy/pymc_experiments.py:661: UserWarning: The power function is experimental and the API may change in the future.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-03-10T22:37:35.175942\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = result.plot_power(alpha=alpha)" + ] + }, + { + "cell_type": "markdown", + "id": "068bb8d2", + "metadata": {}, + "source": [ + "We can observe that the real mean (solid line) deviates quite a bit from the estimated average and cumulative mean (dashed line) of the intervention period. \n", + "\n", + "It's evident that the model fails to accurately represent reality. Even during periods without direct actions affecting the target, the actual values significantly deviate from the model's predicted distribution.\n", + "\n", + ":::{note}\n", + "It is important to remember that in this methodology, the **power analysis must be executed in a period of the same duration before the intervention** and we have the knowledge that no other action is affecting our target variable (`mu`).\n", + "\n", + "If we do not have pre-intervention data, we cannot make this comparison of the model performance prior to the intervention, losing insight into its power.\n", + "\n", + "Additionally, if during the period prior to the intervention our assumption about other activities affecting the target variable is not maintained, we will not be able to trust the results.\n", + ":::\n", + "\n", + "By using the `summary` function and specifying the version, we can generate an output that treats the data as if it were subjected to an intervention. This approach allows for a more comprehensive analysis of the results. In this scenario, the likelihood of our observed value falling within the model's posterior distribution is extremely low. These findings indicate that our model erroneously attributes significance to this particular time frame.\n", + "\n", + "We can observe this after observe the low value of the *Bayesian tail probability*. **How the bayesian tail probability works**? Think of our posterior distribution as a range of possible values that we might see, with the mean value representing the most probable outcome. In this way, we can evaluate the probability of a new value being part of this distribution by measuring how far it deviates from the mean value.\n", + "\n", + "If a value is precisely at the mean, it has a probability of 1 to fall in our posterior. As the value moves away from the average towards both extremes of the distribution, the probability decreases and approaches zero. This process allows us to determine how 'typical' or 'atypical' a new observation is, based on our model estimated posterior. It is an effective tool for interpreting and comprehending the model's inference.\n", + "\n", + ":::{note}\n", + "The function shows the values on their regular scale. We are seeing the possible values by our target (´mu´) in every observation or during the whole period.\n", + ":::" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "1307b896", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
bayesian_tail_probabilityposterior_estimationresultscausal_effectci
cumulative0.01114.31111.43-2.87[112.430854649129, 116.1599388110477]
mean0.0111.4311.14-0.29[11.2430854649129, 11.61599388110477]
\n", + "
" + ], + "text/plain": [ + " bayesian_tail_probability posterior_estimation results \\\n", + "cumulative 0.01 114.31 111.43 \n", + "mean 0.01 11.43 11.14 \n", + "\n", + " causal_effect ci \n", + "cumulative -2.87 [112.430854649129, 116.1599388110477] \n", + "mean -0.29 [11.2430854649129, 11.61599388110477] " + ] + }, + "execution_count": 56, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result.summary(version=\"intervention\", alpha=alpha)" + ] + }, + { + "cell_type": "markdown", + "id": "215fd3bc", + "metadata": {}, + "source": [ + ":::{note}\n", + "The `summary` function provides insights into both cumulative and average results. Under normal circumstances, it helps in understanding the absolute cumulative and average estimations and effects, along with their respective confidence intervals.\n", + ":::\n", + "\n", + "## Update the model\n", + "\n", + "Following the shortcomings of our previous model, we have incorporated new regressors to improve its efficacy. These additional features are instrumental in assessing whether the updated model achieves greater accuracy and precision in mirroring reality. This enhancement is crucial for more accurately estimating experimental outcomes.\n", + "\n", + "During this same process we can iterate to determine the model parameters and configurations that give us a better representation of reality during said period." + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "a7ff64e4", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (4 chains in 4 jobs)\n", + "NUTS: [beta, sigma]\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + " 100.00% [8000/8000 00:10<00:00 Sampling 4 chains, 0 divergences]\n", + "
\n", + " " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 11 seconds.\n", + "Sampling: [beta, sigma, y_hat]\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n" + ] + } + ], + "source": [ + "result = cp.pymc_experiments.SyntheticControl(\n", + " power_df,\n", + " fake_treatment_time,\n", + " formula=\"actual ~ 0 + a + b + c + d + e + f + g\",\n", + " model=cp.pymc_models.WeightedSumFitter(\n", + " sample_kwargs={\"target_accept\": 0.95, \"random_seed\": seed}\n", + " ),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "a71d97a0", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-03-10T22:37:56.609163\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = result.plot(plot_predictors=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "id": "c0911b46", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/Documents/GitHub/CausalPy/causalpy/pymc_experiments.py:573: UserWarning: The power function is experimental and the API may change in the future.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
posterior_estimationresults_systematic_differencesciposterior_mde
cumulative111.37111.435.86e-02[109.88865914515354, 112.89282712094747]1.50
mean11.1411.145.86e-03[10.988865914515353, 11.289282712094748]0.15
\n", + "
" + ], + "text/plain": [ + " posterior_estimation results _systematic_differences \\\n", + "cumulative 111.37 111.43 5.86e-02 \n", + "mean 11.14 11.14 5.86e-03 \n", + "\n", + " ci posterior_mde \n", + "cumulative [109.88865914515354, 112.89282712094747] 1.50 \n", + "mean [10.988865914515353, 11.289282712094748] 0.15 " + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result.power_summary(alpha=alpha)" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "id": "0a82fb59", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/Documents/GitHub/CausalPy/causalpy/pymc_experiments.py:661: UserWarning: The power function is experimental and the API may change in the future.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-03-10T22:37:57.143833\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = result.plot_power(alpha=alpha)" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "id": "293aab40", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
bayesian_tail_probabilityposterior_estimationresultscausal_effectci
cumulative0.95111.37111.435.86e-02[109.88865914515354, 112.89282712094747]
mean0.9511.1411.145.86e-03[10.988865914515353, 11.289282712094748]
\n", + "
" + ], + "text/plain": [ + " bayesian_tail_probability posterior_estimation results \\\n", + "cumulative 0.95 111.37 111.43 \n", + "mean 0.95 11.14 11.14 \n", + "\n", + " causal_effect ci \n", + "cumulative 5.86e-02 [109.88865914515354, 112.89282712094747] \n", + "mean 5.86e-03 [10.988865914515353, 11.289282712094748] " + ] + }, + "execution_count": 61, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result.summary(version=\"intervention\", alpha=alpha)" + ] + }, + { + "cell_type": "markdown", + "id": "5124e0d3", + "metadata": {}, + "source": [ + "This model is outperforming its predecessor significantly. We now observe that the actual value falls within the range estimated by our model. However, there's still a noticeable discrepancy between the estimated total value and the actual one. As previously mentioned, this suggests that we might be overestimating or underestimating certain effects on one side of the distributions.\n", + "\n", + "This becomes apparent when we examine the Bayesian Tail Probability using the `summary` function once more. We now observe a significantly higher value, which is logical. If the actual value is equal to the mean, the probability should be 1. This implies that obtaining that specific value is the most probable outcome according to the posterior distribution.\n", + "\n", + "### Estimation Bias\n", + "\n", + "**Consider this scenario**: If our model exhibits a positive or negative bias, we can recenter our data within the distribution to gauge what the Minimum Detectable Effect (MDE) would be in the absence of this bias. You can apply this bias correction using either the `power_summary` or `power_plot` functions. This adjustment will allow you to observe how the interval of estimation changes.\n", + "\n", + "Essentially, we're implementing a correction to shift our posterior distribution, envisioning a scenario where our model predicts reality with perfect accuracy. This **approach assumes that the same estimation bias present during the model's development will persist throughout the intervention**. While not guaranteed, this is a reasonable assumption, especially considering that the **power test period should closely mirror the actual intervention**, ideally occurring just a few days prior." + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "id": "1f58a2c4", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/Documents/GitHub/CausalPy/causalpy/pymc_experiments.py:573: UserWarning: The power function is experimental and the API may change in the future.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
posterior_estimationresults_systematic_differencesciposterior_mde
cumulative111.37111.435.86e-02[109.94724463521162, 112.95141261100555]1.50
mean11.1411.145.86e-03[10.994724463521159, 11.295141261100554]0.15
\n", + "
" + ], + "text/plain": [ + " posterior_estimation results _systematic_differences \\\n", + "cumulative 111.37 111.43 5.86e-02 \n", + "mean 11.14 11.14 5.86e-03 \n", + "\n", + " ci posterior_mde \n", + "cumulative [109.94724463521162, 112.95141261100555] 1.50 \n", + "mean [10.994724463521159, 11.295141261100554] 0.15 " + ] + }, + "execution_count": 62, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "power_correction=result.power_summary(alpha=0.1, correction=True)\n", + "power_correction" + ] + }, + { + "cell_type": "code", + "execution_count": 63, + "id": "122ed5f5", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/Documents/GitHub/CausalPy/causalpy/pymc_experiments.py:661: UserWarning: The power function is experimental and the API may change in the future.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-03-10T22:38:05.430025\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = result.plot_power(alpha=0.1, correction=True)" + ] + }, + { + "cell_type": "markdown", + "id": "89547ff2", + "metadata": {}, + "source": [ + "With the recent adjustments, our intervals have shifted to the right, aligning the mean estimation with the actual value. This alignment is why the estimation is not visible on the plot. However, it's crucial to acknowledge the real uncertainty that we should account for. The intervals have moved from a range of (100 - 103) to (101 - 104) at the cumulative distribution, indicating the need for an additional unit on each side.\n", + "\n", + "The MDE remains unchanged because the sigma (standard deviation) is the same for our posterior. However, since the absolute value has shifted, we understand the uncertanty on a model without biases. We should apply this correction during the experimentation period. This adjustment will enable us to observe the probable effect more accurately, devoid of any biases from the model's estimation.\n", + "\n", + ":::{note}\n", + "The `power_summary` function by default use ´correction´ equal ´False´.\n", + ":::\n", + "\n", + "## Running our experiment\n", + "\n", + "Having identified which model is most appropriate, we can now proceed with our experiment using the best-suited model. Therefore, we will conduct the experiment during the treatment period, utilizing the complete dataset, now that we have access to all the necessary data." + ] + }, + { + "cell_type": "code", + "execution_count": 64, + "id": "9efc93e9", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (4 chains in 4 jobs)\n", + "NUTS: [beta, sigma]\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + " 100.00% [8000/8000 00:08<00:00 Sampling 4 chains, 2 divergences]\n", + "
\n", + " " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 9 seconds.\n", + "There were 2 divergences after tuning. Increase `target_accept` or reparameterize.\n", + "Sampling: [beta, sigma, y_hat]\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n", + "/Users/carlos.trujillo/anaconda3/envs/CausalPy/lib/python3.11/site-packages/pymc/model/core.py:560: FutureWarning: Model.model property is deprecated. Just use Model.\n", + " warnings.warn(\"Model.model property is deprecated. Just use Model.\", FutureWarning)\n", + "Sampling: [y_hat]\n" + ] + } + ], + "source": [ + "result = cp.pymc_experiments.SyntheticControl(\n", + " df,\n", + " actual_treatment_time,\n", + " formula=\"actual ~ 0 + a + b + c + d + e + f + g\",\n", + " model=cp.pymc_models.WeightedSumFitter(\n", + " sample_kwargs={\"target_accept\": 0.95, \"random_seed\": seed}\n", + " ),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 65, + "id": "d4742fa1", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-03-10T22:38:17.436198\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.3, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = result.plot(plot_predictors=True)" + ] + }, + { + "cell_type": "markdown", + "id": "a078d6a5", + "metadata": {}, + "source": [ + "Excellent, the effect is clearly visible, almost unmistakable to the naked eye. While we might not necessarily need statistical analysis to confirm it, it's still prudent to examine the mean and cumulative effect for thoroughness.\n", + "\n", + "We can utilize the summary function once again to gauge the extent to which our actions have influenced the behavior of the target variable." + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "id": "81fe464d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
bayesian_tail_probabilityposterior_estimationresultscausal_effectci
cumulative5.84e-03563.49512.18-51.30[553.6898523368234, 574.2563747948708]
mean5.84e-0318.7817.07-1.71[18.456328411227446, 19.141879159829028]
\n", + "
" + ], + "text/plain": [ + " bayesian_tail_probability posterior_estimation results \\\n", + "cumulative 5.84e-03 563.49 512.18 \n", + "mean 5.84e-03 18.78 17.07 \n", + "\n", + " causal_effect ci \n", + "cumulative -51.30 [553.6898523368234, 574.2563747948708] \n", + "mean -1.71 [18.456328411227446, 19.141879159829028] " + ] + }, + "execution_count": 66, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result.summary(version=\"intervention\", alpha=0.1)" + ] + }, + { + "cell_type": "markdown", + "id": "8b0e9422", + "metadata": {}, + "source": [ + "As anticipated, the total causal effect of our actions during the intervention is approximately `-52.20` units. This indicates a significant impact on our target variable, with a decrease of about `50` units over the entire period.\n", + "\n", + "When broken down to a daily scale, which matches the granularity of our data, it becomes evident that we're experiencing a loss of approximately `-1.74` units daily on our target variable on average due to this action. This result is significant, especially considering how it deviates from the intervals of our posterior distribution.\n", + "\n", + "The minimum expected value was `564.38` cumulatively and `18.81` on average, and we get `512` and `17` respectively. This discrepancy is why our tail probability is very low, suggesting that these results are unlikely to be a part of our posterior distribution.\n", + "\n", + "**Wait a moment!** This could be a biased result, based on our previous power analysis. If the posterior was shifted to the left relative to the real mean in our power analysis, it implies that **we might be UNDERESTIMATING the decrease**. If our actions are indeed reducing the variable, the actual impact could be even more pronounced than we currently estimate.\n", + "\n", + "To address this, let's reapply the parameter correction, this time incorporating the value needed to adjust for the systematic differences observed in the power analysis." + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "id": "d76c99ed", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
bayesian_tail_probabilityposterior_estimationresultscausal_effectci
cumulative5.84e-03563.55512.18-51.36[553.7484378268815, 574.3149602849289]
mean5.84e-0318.7917.07-1.72[18.46218696023325, 19.147737708834832]
\n", + "
" + ], + "text/plain": [ + " bayesian_tail_probability posterior_estimation results \\\n", + "cumulative 5.84e-03 563.55 512.18 \n", + "mean 5.84e-03 18.79 17.07 \n", + "\n", + " causal_effect ci \n", + "cumulative -51.36 [553.7484378268815, 574.3149602849289] \n", + "mean -1.72 [18.46218696023325, 19.147737708834832] " + ] + }, + "execution_count": 67, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result.summary(\n", + " version=\"intervention\", \n", + " alpha=0.1, \n", + " correction=power_correction[\"_systematic_differences\"]\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "a3fa10a1", + "metadata": {}, + "source": [ + "It appears that the potential causal effect was `52` as we initially thought. We also notice differences in the intervals, yet this additional value doesn't significantly alter our Bayesian Tail Probability (BTP).\n", + "\n", + "## Conclusion\n", + "Power analysis proves to be an invaluable tool in determining the feasibility of an experiment. As we model reality, it guides us in selecting appropriate models and experimenting with those that better interpret reality. Consequently, this enhances our ability to measure the impact of our actions.\n", + "\n", + "However, the application of power analysis still relies on several assumptions, such as:\n", + "\n", + "1. The model's accuracy before the intervention is presumed to remain consistent afterward. Essentially, we assume stable model behavior both pre- and post-intervention.\n", + "2. When applying corrections, we assume that systematic errors will persist throughout the intervention period.\n", + "\n", + "Netherless, if we can uphold these assumptions, then we are positioned to:\n", + "\n", + "1. Determine the likely impact of our actions and potentially discard those that may not produce significant changes, addressing them with other types of tests.\n", + "2. Choose models based on factors other than R-squared ($R^2$) or Mean Absolute Percentage Error (MAPE).\n", + "3. Correct inherent biases in the model, leading to more accurate estimations of effects." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "528386df", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}