Time-Slice-Cross-Validation and Parameter Stability#

In this notebook we will illustrate how to perform time-slice cross validation for a media mix model. This is an important step to evaluate the stability and quality of the model. We not only look into out of sample predictions but also the stability of the model parameters.

These imports and configurations form the fundamental setup necessary for the entire span of this notebook.

The expectation is that a model has already been trained using the functionalities provided in prior versions of the PyMC-Marketing library. Thus, the data generation and training processes will be replicated in a different notebook. Those unfamiliar with these procedures are advised to refer to the “MMM Example Notebook.”

Prepare Notebook#

import warnings

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from pymc_marketing.mmm.time_slice_cross_validation import TimeSliceCrossValidator
from pymc_marketing.paths import data_dir

warnings.simplefilter(action="ignore", category=FutureWarning)

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.facecolor"] = "white"


%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = "retina"

seed: int = sum(map(ord, "mmm"))
rng: np.random.Generator = np.random.default_rng(seed=seed)

Loading Data#

Here we will load our geo level dataset. This will then be used within our Time-Slice CV steps.

data_path = data_dir / "multidimensional_mock_data.csv"
data_df = pd.read_csv(data_path, parse_dates=["date"], index_col=0)
data_df.head()

	date	y	x1	dayofyear	t	geo
0	2018-04-02	3984.662237	159.290009	92	0	geo_a
1	2018-04-09	3762.871794	56.194238	99	1	geo_a
2	2018-04-16	4466.967388	146.200133	106	2	geo_a
3	2018-04-23	3864.219373	35.699276	113	3	geo_a
4	2018-04-30	4441.625278	193.372577	120	4	geo_a

X = data_df.drop(columns=["y"])
y = data_df["y"]

Specify Time-Slice-Cross-Validation Strategy#

The main idea of the time-slice cross validation process is to fit the model on a time slice of the data and then evaluate it on the next time slice. We repeat this process for each time slice of the data. As we want to simulate a production-like environment where we enlarge our training data over time, we make the time-slice size grow over time.

Data Leakage

It is very important to avoid data leakage when performing time-slice cross validation. This means that the model should not see any training data from the future. This also includes any data pre-processing steps!

For example, as mentioned above, we need to compute the costs share for each training time slice independently if we want to avoid data leakage. Other sources of data leakage include using a global feature for thr trend component. In our case, we simply use an increasing variable t so we are safe as we just increase it by one for each time slice.

Run Time-Slice-Cross-Validation Loop#

Depending on the business requirements, we need to decide the initial number of observations to use for fitting the model (n_init) and the forecast horizon (forecast_horizon). For this example, we use the first 342 observations to fit the model and then predict the next 12 observations (3 months).

# Initialize cross-validator
cv = TimeSliceCrossValidator(
    n_init=163,
    forecast_horizon=12,
    date_column="date",
    step_size=1,
)
cv.plot_suite = "new"

# We can check how many splits we will have
# As a reference, the number of splits is computed as:
# n_iterations = y.size - n_init - forecast_horizon + 1
n_splits = cv.get_n_splits(X, y)
print(f"Number of splits: {n_splits}")

Number of splits: 5

Let’s run it!

For more details on the build_mmm_from_yaml, consult the pymc-marketing documentation on Model Deployment.

Alternatively, load a model that has been saved to MLflow via pymc_marketing.mlflow.log_inference_data or has been autologged to MLflow via pymc_marketing.mlflow.autolog(log_mmm=True), from the PyMC-Marketing MLflow module.

results = cv.run(
    X,
    y,
    # You can also pass sampler_config here to speed things up
    sampler_config={
        "tune": 1_000,
        "draws": 1_000,
        "chains": 4,
        "random_seed": seed,
        "target_accept": 0.90,
        "nuts_sampler": "nutpie",
    },
    yaml_path=data_dir
    / "config_files"
    / "multi_dimensional_example_model_with_2_geos.yml",
)

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 11 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.11	63
2000	0.09	127
2000	0.11	63
2000	0.10	127

Sampling: [y]

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 11 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.10	63
2000	0.09	31
2000	0.11	127
2000	0.11	127

Sampling: [y]

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 11 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.09	63
2000	0.12	63
2000	0.10	127
2000	0.10	127

Sampling: [y]

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 11 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.11	31
2000	0.11	127
2000	0.11	63
2000	0.10	63

Sampling: [y]

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 11 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.09	127
2000	0.11	127
2000	0.10	127
2000	0.08	63

Sampling: [y]

# We can view the cross-validation results!
# The CV object is an instance of ArviZ InferenceData
results

arviz.InferenceData

posterior

posterior_predictive

sample_stats

prior

prior_predictive
observed_data

constant_data

fit_data

posterior_predictive_constant_data

cv_metadata

Model Diagnostics#

First, we evaluate whether we have any divergences in the model (we can extend the analysis more more model diagnostics).

# Let's check if there are any divergences
diverging_count = int(results.sample_stats["diverging"].values.sum())
print("Diverging transitions:", diverging_count)

Diverging transitions: 0

We have no divergences in the model 😃!

Evaluate Parameter Stability#

Next, we look at the stability of the model parameters. For a good model, these should not change abruptly over time.

Adstock Alpha

cv.plot.param_stability(
    var_names=["adstock_alpha"],
    # dims={"geo": ["geo_b"]} # to plot specific dimensions only
    figsize=(16, 12),
);

../../_images/c1e9d5efe38f3804e849a7dccf427c1190ee9bc9870ec1fcf67be035f578e301.png

Saturation Beta

cv.plot.param_stability(results, var_names=["saturation_beta"], figsize=(16, 12));

../../_images/ef948a33bc43a826d6633ae9805059c215540f6be1428e6c3be16d0e93bd844e.png

Saturation Lambda

cv.plot.param_stability(results, var_names=["saturation_lam"], figsize=(16, 12));

../../_images/54cd5a942d48e93e793de9b7906852176a83a4c1dda49e6f7137cdc1699fd2b6.png

The parameters seem to be stable over time. This implies that the estimates ROAS will not change abruptly over time.

Evaluate Out of Sample Predictions#

Finally, we evaluate the out of sample predictions. To begin with, we can simply plot the posterior predictive distributions for each iteration for both the training and test data.

# Plot model predictions across time slices
cv.plot.predictions(figsize=(20, 40));

../../_images/61028a105b1198153ab5f0c63ea8402bd857fe0e3d79f35a11696ff684b1a42a.png

Overall, the out of sample predictions look very good 🚀!

We can quantify the model performance using the Continuous Ranked Probability Score (CRPS).

“The CRPS — Continuous Ranked Probability Score — is a score function that compares a single ground truth value to a Cumulative Distribution Function. It can be used as a metric to evaluate a model’s performance when the target variable is continuous and the model predicts the target’s distribution; Examples include Bayesian Regression or Bayesian Time Series models.”

For a nice explanation of the CRPS, check out this blog post.

In PyMC-Marketing, we provide the function crps to compute this metric. We can use it to compute the CRPS score for each iteration.

# Compute the CRPS score for each iteration and plot!
cv.plot.crps(
    # dims={"geo": ["geo_b"]} # to plot specific dimensions only
);

../../_images/4f8351b80d2ded6e3d2208cf269e262944ef9274f7958d8513cdbede204e1756.png

Event though the visual results look great, we see that the CRPS mildly decreases for the training data while it increases for the test data as we increase the size of the training data. This is a sign that we are overfitting the model to the training data. Some strategies to overcome this issue include using regularization techniques and re-evaluate the model specification. This should be an iterative process.

%load_ext watermark
%watermark -n -u -v -iv -w -p pymc_marketing,pytensor,numpyro

Last updated: Mon, 11 May 2026

Python implementation: CPython
Python version       : 3.12.13
IPython version      : 9.12.0

pymc_marketing: 0.19.3
pytensor      : 2.38.2
numpyro       : 0.20.1

arviz         : 0.23.4
matplotlib    : 3.10.8
numpy         : 2.4.3
pandas        : 2.3.3
platform      : 1.0.8
pymc_marketing: 0.19.3

Watermark: 2.6.0