Custom Models with MMM components#

The underlying components used in the MMM class provide flexibility to build other, custom models. With a little knowledge of PyMC and how to customize these PyMC-Marketing components, a lot of different use-cases can be covered.

This notebook is not an introduction but rather an advance example for those trying to understand the PyMC-Marketing internals for flexibility for custom use-cases.

Overview#

This notebook will cover the currently exposed model components from the PyMC-Marketing API. At the moment, this includes:

media transformations
- adstock: how today’s media has an effect in the future
- saturation: the diminishing returns for media
recurring seasonality

For each of these, the flexibility and customization will be showcased and combined together in a toy model with with PyMC directly.

Setup#

from functools import partial

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import xarray as xr

from pymc_marketing import mmm
from pymc_marketing.mmm.plot import plot_curve
from pymc_marketing.prior import Prior

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

%config InlineBackend.figure_format = "retina"

seed = sum(map(ord, "PyMC-Marketing provides flexible model components"))
rng = np.random.default_rng(seed)

draw = partial(pm.draw, random_seed=rng)

Media Transformations#

There are classes for each of the adstock and saturation transformations. They can be imported from the pymc_marketing.mmm module.

saturation = mmm.MichaelisMentenSaturation()

Saturation curves can take many different forms. In this example, we will use the Michaelis Menten curve which we provide in the MichaelisMentenSaturation class.

This curve has two parameters, alpha and lam.

A characteristic of these curves are diminishing returns in order to indicate saturation of a media variable. This can be seen in the plateauing as x increases.

../../_images/5c8f7af436ec881bb17ec64fbb41558b1bb6ef0005f4effcc0771432bbaacd20.png

Adding Parameter Dimensions#

In most cases, a separate saturation function will be estimated for each media channel. A dimension needs to be added to the prior of the function parameters to account for this.

Let’s create some example data to work toward this transformation.

For this example, we will have 2 years of media spend for 4 channels

n_dates = 52 * 2
dates = pd.date_range("2023-01-01", periods=n_dates, freq="W-MON")

channels = ["C1", "C2", "C3", "C4"]

coords = {
    "date": dates,
    "channel": channels,
}

df_spends = random_spends(coords=coords).to_pandas()
df_spends.head()

channel	C1	C2	C3	C4
date
2023-01-02	0.0	0.992756	0.000000	0.000000
2023-01-09	0.0	0.902103	0.203395	0.000000
2023-01-16	0.0	0.559487	0.000000	0.000000
2023-01-23	0.0	0.990124	0.000000	2.013755
2023-01-30	0.0	0.753384	0.000000	0.377570

../../_images/09a0186537ac50de2cf92c306c55539313950574043cd5c7424e3b7a883f9fc1.png

As mentioned, the default priors do not have a channel dimension. In order to use with the in our model with “channel” dim, we have to add the dims to each of the function priors.

for dist in saturation.function_priors.values():
    dist.dims = "channel"

saturation.function_priors

{'alpha': Prior("Gamma", mu=2, sigma=1, dims="channel"),
 'lam': Prior("HalfNormal", sigma=1, dims="channel")}

The previous workflow can be used to understand our priors still. Just pass the coords to the sample_prior method in order to add dims to the appropriate variables.

Since each channel prior is the same, there will just be some noise between the HDI and curve samples.

curve = saturation.sample_curve(prior)
saturation.plot_curve(curve);

Sampling: []

../../_images/6256f8dc589b5451dd2e9b104ced11b2a0420ce377434d05941272b5e4b225b4.png

Using in PyMC Model#

When using the transformation in a larger PyMC model, the apply method will be used.

This method will:

create distributions based on prior specification of the instance
apply the transformation to the data

The dims parameter is the shape of parameters and not the data. The data has a different shape but will need to be broadcastable with the parameters!

with pm.Model(coords=coords) as model:
    saturated_spends = saturation.apply(df_spends, dims="channel")

Since independent alpha and lam were specified, we see that in the model graph below:

pm.model_to_graphviz(model)

../../_images/1f7956f649078431df61a88692564581d0ef6e9098ffc1b5d57c5d18d50b5470.svg

Note

Neither the df_spends nor saturated_spends show in the model. If needed, use pm.Data and pm.Deterministic to save off.

Our variable will be (date, channel) dims.

saturated_spends.type.shape

(104, 4)

We can manipulate this in anyway we’d like to connect it in with the larger model.

Changing Assumptions#

As hinted above, the priors for the function parameters are customizable which can lead to many different models. Change the priors, change the model.

The prior distributions just need to follow the distribution API here.

Instead of the defaults, we can use:

hierarchical parameter for lam parameter
common alpha parameter

hierarchical_lam = Prior(
    "HalfNormal",
    sigma=Prior("HalfNormal", sigma=1),
    dims="channel",
)
common_alpha = Prior("Gamma", mu=2, sigma=1)
priors = {
    "lam": hierarchical_lam,
    "alpha": common_alpha,
}

saturation = mmm.MichaelisMentenSaturation(priors=priors)

saturation.function_priors

{'alpha': Prior("Gamma", mu=2, sigma=1),
 'lam': Prior("HalfNormal", sigma=Prior("HalfNormal", sigma=1), dims="channel")}

Then this can be used in a new PyMC model which leads to a much different model graph than before!

with pm.Model(coords=coords) as model:
    saturated_spends = saturation.apply(df_spends, dims="channel")

pm.model_to_graphviz(model)

../../_images/1ffab4133abcb3e47dd3c74966f7302f070dc225bd90697e687af4382c508af4.svg

The shape of the output will still be (date, channel) even though some of the parameter’s dims has changed.

saturated_spends.type.shape

(104, 4)

The previous workflow still helps us understand the produced curves:

sample_prior
sample_curve
plot_curve

prior = saturation.sample_prior(coords=coords, random_seed=rng)

Sampling: [saturation_alpha, saturation_lam, saturation_lam_sigma]

Though they all look the same in the prior, the data generation process is indeed different as seen in the model graph.

curve = saturation.sample_curve(prior)
saturation.plot_curve(curve);

Sampling: [saturation_lam, saturation_lam_sigma]

../../_images/ba4c9042ccc055a07ae4a0ed0ecad0cd0abb5f5e74097456477f9dec0f202f67.png

Geo Hierarchical Model#

The dimensions of the parameters are not limited to 1D so additional hierarchies can be defined.

Below defines:

alpha which is hierarchical across channels
lam which is common across all geos but different channels

# For reference
mmm.MichaelisMentenSaturation.default_priors

{'alpha': Prior("Gamma", mu=2, sigma=1), 'lam': Prior("HalfNormal", sigma=1)}

hierarchical_alpha = Prior(
    "Gamma",
    mu=Prior("HalfNormal", sigma=1, dims="geo"),
    sigma=Prior("HalfNormal", sigma=1, dims="geo"),
    dims=("channel", "geo"),
)
common_lam = Prior("HalfNormal", sigma=1, dims="channel")
priors = {
    "alpha": hierarchical_alpha,
    "lam": common_lam,
}
saturation = mmm.MichaelisMentenSaturation(priors=priors)

Our new data set needs to have information for geo now. This is channel spends by date and geo. This is stored in an xarray.DataArray which can be converted to a 3D numpy.ndarray.

Displaying the data is easy with pandas.

geo_coords = {
    **coords,
    "geo": ["Region1", "Region2", "Region3"],
}

geo_spends = random_spends(coords=geo_coords)

geo_spends.to_series().unstack("channel").head(6)

	channel	C1	C2	C3	C4
date	geo
2023-01-02	Region1	0.000000	0.000000	0.541483	0.0
	Region2	0.213863	0.000000	0.000000	0.0
	Region3	1.100906	0.692670	1.493263	0.0
2023-01-09	Region1	0.309415	0.614882	0.818824	0.0
	Region2	2.264686	0.930828	0.028902	0.0
	Region3	0.433401	0.660699	1.318996	0.0

As long as the dims argument of apply can broadcast with the data going in, then the media transformations can be used!

Here, the data is in the shape (date, channel, geo) so it can broadcast with the parameters in shape (channel, geo) to create the saturated spends.

with pm.Model(coords=geo_coords) as geo_model:
    geo_data = pm.Data(
        "geo_data",
        geo_spends.to_numpy(),
        dims=("date", "channel", "geo"),
    )
    saturated_geo_spends = pm.Deterministic(
        "saturated_geo_spends",
        saturation.apply(geo_data, dims=("channel", "geo")),
        dims=("date", "channel", "geo"),
    )

The saturation assumptions can be seen in the model graph:

pm.model_to_graphviz(geo_model)

../../_images/27c557e2afc443f146e11dd4db29a93b38022123ed59d62da083854f709c2ceb.svg

Tip

The PyMC model context will stay the same but changing model assumptions will happen with input data and prior configuration!

Seasonality#

Recurring seasonality can be modeled with either a MonthlyFourier or YearlyFourier instance.

yearly = mmm.YearlyFourier(n_order=2)

There is a similar workflow to understand these priors as before:

sample_prior: Sample all the priors
sample_curve: Sample the curve across the whole period
plot_curve: Plot the HDI and few samples

prior = yearly.sample_prior()
curve = yearly.sample_curve(prior)
yearly.plot_curve(curve);

Sampling: [fourier_beta]
Sampling: []

../../_images/67b3ce211948d44252e0e7ced34388bde491cee4f22d2ba9ede876e3d3c56cc5.png

This also supports arbitrary hierarchies that can be defined with the Prior class. Pass these in with the prior parameters.

Note

A dimension associated with the prefix will be required! By default it is fourier

prior = Prior(
    "Normal",
    mu=[0, 0, -1, 0],
    sigma=Prior("Gamma", mu=0.15, sigma=0.1, dims="fourier"),
    dims=("geo", "fourier"),
)
yearly = mmm.YearlyFourier(n_order=2, prior=prior)

The above workflow works here as well! The coords just need to be passed like in pm.Model.

coords = {
    "geo": ["A", "B"],
}
prior = yearly.sample_prior(coords=coords)
curve = yearly.sample_curve(prior)

Sampling: [fourier_beta, fourier_beta_sigma]

Sampling: []

Based on the hierarchical priors, we can see similar seasonality betweens geos. However, they are not exactly the same!

subplot_kwargs = {"ncols": 1}
sample_kwargs = {"n": 3}
fig, _ = yearly.plot_curve(
    curve, subplot_kwargs=subplot_kwargs, sample_kwargs=sample_kwargs
)
fig.suptitle("Prior seasonality");

../../_images/a82899787dea7e68b197eb114a0c0d5fd210892e494c6bbcb8735222506490c5.png

Summary#

Custom models are possible using the components that build up the MMM class and PyMC distributions themselves. With some prior distribution configuration and the components that PyMC-Marketing provides, novel models can be built up to fit various use-cases and various model assumptions.

Much of the flexibility will come from the prior distribution configuration rather then the transformation themselves. This is meant to keep a standard interface while working with them regardless what their role is.

If there is any suggestions or feedback on how to make better custom models with the package, create a GitHub Issue or chime into the various discussions.

Though models can be built up like this, the prebuilt structures provide many benefits as well. For instance, the MMM class provides:

scaling of input and output data
plotting methods for parameters, predictive data, contributions, etc
customized adstock and saturation transformations
out of sample predictions
lift test integration
budget optimization

Our recommendation is to start with the prebuilt models and work up from there.

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

Last updated: Fri Aug 02 2024

Python implementation: CPython
Python version       : 3.12.4
IPython version      : 8.25.0

pymc_marketing: 0.8.0
pytensor      : 2.22.1

matplotlib    : 3.9.0
numpy         : 1.26.4
arviz         : 0.18.0
pymc          : 5.15.1
pandas        : 2.2.2
xarray        : 2024.6.0
pymc_marketing: 0.8.0

Watermark: 2.4.3