Modelado Predictivo Anterior

Esta guía proporciona una introducción a la modelización predictiva previa utilizando PyMC (y PyMC-Marketing) y la clase {class}Prior <pymc_extras.prior.Prior>` de PyMC-Marketing. Comenzamos examinando un ejemplo más simple y luego veremos cómo aplicarlo a escenarios de casos reales con modelos de mezcla de marketing en PyMC-Marketing.

Antes de profundizar en los detalles técnicos, comprendamos por qué los priors son cruciales en el análisis bayesiano y su importancia práctica en las aplicaciones industriales.

Comprendiendo la Inferencia Bayesiana

La inferencia bayesiana se basa en el teorema de Bayes, que proporciona una forma formal de actualizar nuestras creencias sobre los parámetros \(\theta\) (por ejemplo, la tasa de saturación o de decaimiento en modelos de mezcla de marketing) dado los datos observados \(y\):

\[p(\theta|y) = \frac{p(y|\theta)p(\theta)}{p(y)}\]

Dónde:

• \(p(\theta|y)\) es la probabilidad posterior (lo que queremos aprender)

• \(p(y|\theta)\) es la verosimilitud (cómo se generan los datos)

• \(p(\theta)\) es la probabilidad a priori (nuestras creencias iniciales)

• \(p(y)\) es la evidencia (una constante de normalización), que se puede escribir como \(p(y) = \displaystyle{ \int p(y|\theta)p(\theta)d\theta }\)

La distribución posterior combina nuestro conocimiento previo con los datos observados para proporcionarnos creencias actualizadas sobre los parámetros. En la práctica, a menudo trabajamos con la posterior no normalizada:

\[p(\theta|y) \propto p(y|\theta)p(\theta)\]

Esto se debe a que la constante de normalización \(p(y)\) a menudo es intractable de calcular directamente.

Por qué los priors son importantes en la industria

En aplicaciones industriales, los priors cumplen varios propósitos cruciales:

1. Integración del Conocimiento del Dominio:

   • Incorporando conocimiento experto en modelos

   • Aprovechando datos históricos de proyectos similares

   • Codificación de restricciones y requisitos empresariales

2. Gestión de Riesgos:

   • Prevención de predicciones poco realistas

   • Asegurando un comportamiento estable del modelo

   • Gestionar la incertidumbre en la toma de decisiones

3. Eficiencia de Datos:

   • Haciendo que los modelos funcionen con datos limitados

   • Convergencia más rápida hacia soluciones razonables

   • Predicciones robustas en nuevos escenarios

4. Regularización del Modelo:

   • Prevención del sobreajuste

   • Manejo de la multicolinealidad

   • Tratando con datos escasos

Escenarios Comunes de Especificación de Prioridades

En el análisis de marketing, a menudo se encontrará con estos escenarios:

1. Modelos de Mezcla de Marketing:

   • Efectividad del canal de medios (típicamente positiva)

   • Rendimientos decrecientes (restricciones de forma)

   • Calibración de pruebas de elevación

2. Valor del Tiempo de Vida del Cliente:

   • Tasas de compra (valores positivos)

   • Probabilidades de deserción (entre 0 y 1)

   • Distribuciones de valor monetario (positivas, a menudo logarítmicamente normales)

3. Pruebas A/B:

   • Tasas de conversión (limitadas entre 0 y 1)

   • Mediciones de elevación (centradas en efectos pequeños)

   • Impactos en los ingresos (potencialmente de cola pesada)

Las pruebas A/B son un campo fascinante donde los métodos bayesianos pueden ofrecer enormes beneficios. Sin embargo, omitir el análisis predictivo previo puede llevar a resultados negativos, como se describe en el excelente artículo «The Bet Test: Spotting Problems in Bayesian A/B Test Analysis» de Tyler Buffington.

¿Qué es el Modelado Predictivo Previo?

El modelado predictivo previo es un paso crucial en el flujo de trabajo bayesiano que nos ayuda a validar nuestras elecciones previas antes de ver los datos reales. El proceso implica:

1. Especificación:

   • Defina distribuciones a priori para los parámetros del modelo

   • Codificar el conocimiento del dominio y las restricciones

   • Documentar suposiciones y elecciones

2. Simulación:

   • Parámetros de muestra de distribuciones anteriores

   • Generar datos sintéticos utilizando la estructura del modelo

   • Cree múltiples escenarios de posibles resultados

3. Validación:

   • Verifique si los datos simulados coinciden con la experiencia del dominio.

   • Verifique que se excluyan los escenarios imposibles.

   • Asegure una cobertura razonable de los posibles resultados.

Beneficios en la Práctica

1. Detección Temprana de Problemas:

   • Identificar suposiciones poco realistas

   • Detectar problemas numéricos antes del ajuste del modelo

   • Validar la estructura del modelo

2. Comunicación con las partes interesadas:

   • Visualizar las implicaciones del modelo

   • Justificar las elecciones de modelado

   • Establecer expectativas realistas

3. Desarrollo de Modelos:

   • Iterar sobre elecciones anteriores de manera eficiente

   • Comparar especificaciones alternativas

   • Evolución del modelo de documento

4. Evaluación de Riesgos:

   • Entender las limitaciones del modelo

   • Identificar casos límite

   • Plan para modos de fallo

La distribución predictiva previa \(p(y)\) representa nuestras creencias sobre los datos antes de observarlos. Matemáticamente, es la distribución de los datos marginalizada sobre la previa:

\[p(y) = \int p(y|\theta)p(\theta)d\theta\]

En la práctica, podemos muestrear de esta distribución mediante:

1. Extrayendo parámetros de lo anterior: \(\theta^{(s)} \sim p(\theta)\)

2. Generando datos a partir de la verosimilitud: \(y^{(s)} \sim p(y|\theta^{(s)})\)

Este proceso nos ayuda a validar nuestro modelo de varias maneras:

1. Cobertura del Espacio de Parámetros: Las muestras \(\{\theta^{(s)}\}_{s=1}^S\) nos muestran qué valores de parámetros consideramos plausibles.

2. Cobertura del Espacio de Datos: Las muestras \(\{y^{(s)}\}_{s=1}^S\) nos muestran qué datos puede generar nuestro modelo.

3. Sensibilidad del Modelo: La relación entre \(\theta^{(s)}\) y \(y^{(s)}\) muestra cómo los parámetros influyen en las predicciones.

Exploramos estos conceptos a través de ejemplos prácticos utilizando la clase Prior de PyMC-Marketing.

Preparar el cuaderno

from itertools import product

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import preliz as pz
import pymc as pm
import seaborn as sns
from numpy.typing import NDArray
from pymc_extras.prior import Prior

from pymc_marketing.hsgp_kwargs import HSGPKwargs
from pymc_marketing.mmm import MMM, GeometricAdstock, LogisticSaturation
from pymc_marketing.paths import data_dir

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


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

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

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

Ejemplo simple: Distribución normal

Comencemos con un ejemplo simple utilizando una distribución normal. Vamos a:

1. Generar un conjunto de datos sintético

2. Estudie la distribución observada

3. Establezca un prior para la media y la desviación estándar.

4. Muestra de su distribución predictiva previa

Primero generamos un conjunto de datos sintético a partir de una distribución normal con una media de uno y una desviación estándar de dos.

true_mu: float = 1.0
true_sigma: float = 2.0
n_observations: int = 200

data: NDArray = rng.normal(loc=true_mu, scale=true_sigma, size=n_observations)

sample_mean: float = data.mean()
sample_std: float = data.std()

Visualicemos los datos observados.

fig, ax = plt.subplots()
sns.kdeplot(data, fill=True, color="C0", ax=ax)
ax.axvline(
    sample_mean, color="C1", linestyle="--", label=f"Sample Mean: {sample_mean:.2f}"
)
ax.axvline(true_mu, color="C2", linestyle="--", label=f"True Mean: {true_mu:.2f}")
ax.legend()
ax.set(title="Observed Data", xlabel="Value", ylabel="Density")
ax.set_title("Observed Data", fontsize=16, fontweight="bold");

Asumimos que no conocemos la media verdadera y la desviación estándar de los datos (como en casi todos los casos). Nuestra idea es ajustar un modelo bayesiano para intentar recuperar los parámetros verdaderos.

Consideramos la forma paramétrica de los datos:

\[::\]

¿Cuáles podrían ser priors sensatos para la media y la desviación estándar? Aquí es donde entra el modelado predictivo previo.

Media y Desviación Estándar Fijas

Primero, consideramos el caso simple en el que establecemos valores fijos para la media y la desviación estándar y muestreamos de la distribución predictiva previa (la distribución normal).

Considere los siguientes valores para la media y la desviación estándar:

mu_1, mu_2 = -2, 0
sigma_1, sigma_2 = 0.5, 3

Dado estos valores, podemos simplemente muestrear de una distribución normal con estos parámetros y compararlo con los datos observados.

fig, ax = plt.subplots()

sns.kdeplot(data, fill=True, color="C0", label="Observed Data", ax=ax)
for i, (mu, sigma) in enumerate(product((mu_1, mu_2), (sigma_1, sigma_2))):
    normal_prior = pz.Normal(mu=mu, sigma=sigma)
    normal_prior.plot_pdf(color=f"C{i + 1}", legend="legend", ax=ax)
ax.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
ax.set(xlabel="Value", ylabel="Density")
fig.suptitle(
    "Prior Predictive Samples (fixed mean and std)",
    fontsize=16,
    fontweight="bold",
    y=0.93,
);

Al observar los gráficos anteriores, podemos ver cómo diferentes combinaciones de los parámetros de media (μ) y desviación estándar (σ) afectan el ajuste a nuestros datos observados:

• \(\mu=-2, \sigma=0.5\): La distribución es demasiado estrecha y está centrada demasiado a la izquierda.

• \(\mu=-2, \sigma=3\): Mejor dispersión pero aún centrado demasiado a la izquierda

• \(\mu=0, \sigma=0.5\): Mejor centrado pero demasiado estrecho

• \(\mu=0, \sigma=3\): Esto parece ser el mejor ajuste, con un buen centro y dispersión que coinciden con los datos observados.

Esta comparación visual nos ayuda a entender que las distribuciones previas centradas alrededor de \(\mu=0\) con una desviación estándar más amplia (\(\sigma=3\)) pueden ser más apropiadas para nuestros datos.

Muestreo Predictivo Previo

Ahora podemos dar un paso adelante y, en lugar de fijar la media y la desviación estándar, podemos establecer un prior para ellas. Para comenzar, establecemos los siguientes priors para la media y la desviación estándar:

• \(\mu \sim \text{Normal}(0, 2)\)

• \(\sigma \sim \text{Exponential}(1/3)\)

Primero muestreamos de los parámetros del modelo.

# Set the number of prior samples
n_prior_samples: int = 1_000

# Sample from the prior distributions
mus = rng.normal(loc=0, scale=2, size=n_prior_samples)
sigmas = rng.exponential(scale=3, size=n_prior_samples)

fig, ax = plt.subplots(
    nrows=2,
    ncols=1,
    figsize=(10, 10),
    sharex=False,
    sharey=False,
    layout="constrained",
)

sns.kdeplot(mus, fill=True, color="C2", ax=ax[0])
ax[0].set(xlabel="Value", ylabel="Density")
ax[0].set_title(
    "Prior Predictive Samples (normal mean)", fontsize=16, fontweight="bold"
)

sns.kdeplot(sigmas, fill=True, color="C3", ax=ax[1])
ax[1].set(xlabel="Value", ylabel="Density")
ax[1].set_title(
    "Prior Predictive Samples (exponential std)", fontsize=16, fontweight="bold"
);

A continuación, pasamos estas muestras a través de una distribución normal para obtener las muestras predictivas previas.

prior_predictive_samples = rng.normal(loc=mus, scale=sigmas)

fig, ax = plt.subplots()
sns.kdeplot(
    prior_predictive_samples,
    fill=True,
    color="C1",
    label="Prior Predictive Samples",
    ax=ax,
)
sns.kdeplot(data, fill=True, color="C0", label="Observed Data", ax=ax)
ax.legend()
ax.set(xlabel="Value", ylabel="Density", xlim=(-30, 30))
ax.set_title("Prior Predictive Samples", fontsize=16, fontweight="bold");

La verificación predictiva previa muestra una buena concordancia entre nuestros datos simulados y las muestras predictivas previas. Las densidades superpuestas indican que nuestras distribuciones previas elegidas para la media (normal) y la desviación estándar (exponencial) son razonables y pueden generar datos similares a los que observamos. Esto sugiere que nuestras especificaciones previas son apropiadas para esta tarea de modelado.

Muestreo Predictivo Previo con PyMC

Comenzamos definiendo el modelo en PyMC. Tenga en cuenta que no necesitamos pasar los datos observados (aún) para muestrear de la distribución predictiva a priori.

with pm.Model(coords={"n_obs": range(data.shape[0])}) as model:
    # Define the prior distributions
    mu = pm.Normal("mu", mu=0, sigma=2)
    sigma = pm.Exponential("sigma", lam=1 / 3)

    # Define the likelihood
    y = pm.Normal("y", mu=mu, sigma=sigma, dims="n_obs")

model.to_graphviz()

PyMC ofrece una forma conveniente de muestrear de la distribución predictiva previa utilizando la función pymc.sample_prior_predictive <pymc.sample_prior_predictive>.

with model:
    # Sample from the prior predictive distribution
    idata = pm.sample_prior_predictive(draws=n_prior_samples)

Sampling: [mu, sigma, y]

fig, ax = plt.subplots()
sns.kdeplot(
    idata["prior"]["y"].to_numpy().flatten(),
    fill=True,
    color="C1",
    label="Prior Predictive Samples",
    ax=ax,
)
sns.kdeplot(data, fill=True, color="C0", label="Observed Data", ax=ax)
ax.legend()
ax.set(xlabel="Value", ylabel="Density", xlim=(-30, 30))
ax.set_title("Prior Predictive Samples", fontsize=16, fontweight="bold");

Las muestras predictivas anteriores son muy similares a las obtenidas a partir del muestreo manual.

Muestreo Predictivo Posterior

Ahora podemos ajustar el modelo a los datos observados y muestrear de la distribución predictiva posterior.

# Condition the model on the observed data
conditioned_model = pm.observe(model, {"y": data})

# Sample
with conditioned_model:
    # Sample from the posterior distribution
    idata.extend(pm.sample())
    # Sample from the posterior predictive distribution
    idata.extend(pm.sample_posterior_predictive(idata))

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mu, sigma]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 0 seconds.
Sampling: [y]

Vamos a trazar las distribuciones posteriores de los parámetros de la media y la desviación estándar.

axes = az.plot_posterior(idata, var_names=["mu", "sigma"], figsize=(12, 5))

axes[0].axvline(true_mu, color="C2", linestyle="--", label=f"True Mean: {true_mu:.1f}")
axes[1].axvline(
    true_sigma, color="C2", linestyle="--", label=f"True Std: {true_sigma:.1f}"
)
axes[0].legend()
axes[1].legend()

fig = plt.gcf()
fig.suptitle("Posterior Distributions", fontsize=16, fontweight="bold");

Obtenemos valores muy cercanos a los valores verdaderos (esto se conoce como recuperación de parámetros).

Podemos comparar las distribuciones previas y posteriores de los parámetros de la media y la desviación estándar. Estos gráficos son muy útiles para entender el impacto de los datos en los parámetros.

axes = az.plot_dist_comparison(idata, var_names=["mu"], figsize=(12, 7))
axes = axes.flatten()
for ax in axes:
    ax.axvline(true_mu, color="C2", linestyle="--", label=f"True Mean: {true_mu:.1f}")
    ax.legend()
fig = plt.gcf()
fig.suptitle(
    "Prior vs Posterior Distributions for the Mean", fontsize=16, fontweight="bold"
);

axes = az.plot_dist_comparison(idata, var_names=["sigma"], figsize=(12, 7))
axes = axes.flatten()
for ax in axes:
    ax.axvline(
        true_sigma, color="C2", linestyle="--", label=f"True Std: {true_sigma:.1f}"
    )
    ax.legend()
fig = plt.gcf()
fig.suptitle(
    "Prior vs Posterior Distributions for the Standard Deviation",
    fontsize=16,
    fontweight="bold",
);

Finalmente, podemos trazar la distribución predictiva posterior frente a los datos observados.

ax = az.plot_ppc(idata, var_names=["y"], figsize=(12, 7))
ax.set_title("Posterior Predictive Samples", fontsize=16, fontweight="bold");

Vemos que la distribución predictiva posterior es muy similar a los datos observados. Esto es una buena señal de que nuestro modelo está bien especificado.

Advertencia

Al trabajar con modelos lineales generalizados (GLMs) y modelos con transformaciones no lineales (por ejemplo, modelos de mezcla de medios), las distribuciones a priori normales sobre los coeficientes pueden, en ocasiones, conducir a comportamientos inesperados debido a las funciones de enlace no lineales. Asegúrese de revisar cuidadosamente la distribución predictiva a priori para estos modelos. Esto se explica de manera excelente en el (fantástico) libro «Statistical Rethinking» de Richard McElreath.

---

Muestreo Predictivo Previos para Modelos de Mezcla de Marketing

En esta sección nos enfocamos en cómo realizar muestreo predictivo previo para modelos de mezcla de marketing. Utilizaremos los mismos datos del cuaderno Cuaderno de ejemplo para MMM. En particular, nos centramos en el ROAS de los dos canales x1 y x2.

Leer datos

Comenzamos leyendo los datos.

data_path = data_dir / "mmm_example.csv"

data_df = pd.read_csv(data_path, parse_dates=["date_week"]).assign(
    y=lambda df: df["y"] / 1_000
)

data_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 179 entries, 0 to 178
Data columns (total 8 columns):
 #   Column     Non-Null Count  Dtype         
---  ------     --------------  -----         
 0   date_week  179 non-null    datetime64[ns]
 1   y          179 non-null    float64       
 2   x1         179 non-null    float64       
 3   x2         179 non-null    float64       
 4   event_1    179 non-null    float64       
 5   event_2    179 non-null    float64       
 6   dayofyear  179 non-null    int64         
 7   t          179 non-null    int64         
dtypes: datetime64[ns](1), float64(5), int64(2)
memory usage: 11.3 KB

Recuerde que queremos entender el ROAS y la contribución de los dos canales x1 y x2 a los datos de ventas y. Vamos a graficar los datos.

fig, ax = plt.subplots(
    nrows=3,
    ncols=1,
    figsize=(12, 9),
    sharex=True,
    sharey=False,
    layout="constrained",
)

sns.lineplot(data=data_df, x="date_week", y="y", color="black", ax=ax[0])
ax[0].set(title="Sales", xlabel="Date", ylabel="Sales")

sns.lineplot(data=data_df, x="date_week", y="x1", color="C0", ax=ax[1])
ax[1].set(title=r"$x_1$", xlabel="Date", ylabel="Spend")

sns.lineplot(data=data_df, x="date_week", y="x2", color="C1", ax=ax[2])
ax[2].set(title=r"$x_2$", xlabel="Date", ylabel="Spend")

fig.suptitle("Observed MMM Data", fontsize=16, fontweight="bold");

Internamente, durante el proceso de ajuste a través de la clase MMM, los canales y las variables objetivo se escalan dividiendo por el valor máximo (por variable). Visualicemos los datos escalados.

fig, ax = plt.subplots()
sns.lineplot(
    data=data_df.assign(y=lambda df: df["y"] / df["y"].max()),
    x="date_week",
    y="y",
    color="black",
    ax=ax,
)
ax.set(title="Sales", xlabel="Date", ylabel="Sales")
ax.set_title("Scaled Sales", fontsize=16, fontweight="bold");

Vemos que hay un componente de tendencia leve y algo de estacionalidad (por favor, consulte el cuaderno Cuaderno de ejemplo para MMM para más detalles).

Especificación de Priors sobre los Parámetros de Medios

Ahora describimos una forma de pensar sobre la especificación de los priors en los parámetros de los medios. En aplicaciones reales, el modelador generalmente tiene conocimientos y información valiosos adicionales del dominio que se pueden incorporar al modelo a través de priors (por ejemplo, pruebas de elevación, puntos de referencia, etc.). En este ejemplo, simplemente establecemos priors a partir de los datos de gasto de entrada.

Recuerde que la idea de estos tipos de modelos de mezcla de medios es pasar el gasto en medios a través de una transformación de adstock y luego a través de una función de saturación (ambas, funciones no lineales) antes de sumarlas al modelo lineal para añadir otras variables de control. Estas transformaciones no lineales están controladas por los parámetros \(\alpha\) (adstock), \(\lambda\) (saturación) y los coeficientes \(\beta\) (regresión). Debemos especificar priors para estos parámetros. Sin embargo, la clase MMM proporciona priors razonables por defecto.

Adstock

Comencemos examinando el parámetro de adstock. El parámetro de adstock controla la tasa de decaimiento del efecto mediático a lo largo del tiempo. Un valor más alto de \(\alpha\) indica un decaimiento más lento, lo que significa que el efecto mediático persiste por más tiempo. El valor de \(\alpha\) suele estar entre 0 y 1. Para muchos canales digitales, esperamos que este parámetro esté cerca de \(0\), mientras que para los canales offline (como la televisión) esperamos que este parámetro esté cerca de \(1\).

Asumiendo que \(x_1\) y \(x_2\) son canales digitales, esperamos que el parámetro de adstock esté cerca de \(0\). Podemos establecer un prior para el parámetro de adstock de la siguiente manera:

fig, ax = plt.subplots()
# The beta distribution has support in [0, 1]
adstock_alpha_prior = Prior("Beta", alpha=1, beta=3)
adstock_alpha_prior.preliz.plot_pdf(ax=ax)
ax.set(xlabel="Adstock Alpha", ylabel="Density")
ax.set_title("Adstock Alpha Prior", fontsize=16, fontweight="bold");

Podemos analizar las implicaciones de la tasa de decaimiento del adstock en el efecto de los medios.

# We usually set the maximum lags based on the context, but we could even try to learn it from the data.
# The challenge is the typical number of data points available for media mix models.
adstock = GeometricAdstock(priors={"alpha": adstock_alpha_prior}, l_max=8)
prior = adstock.sample_prior(samples=1_000, random_seed=rng)
curve = adstock.sample_curve(prior)
fig, axes = adstock.plot_curve(curve, sample_kwargs={"rng": rng})
axes[0].set_title("Geometric Adstock Curve", fontsize=16, fontweight="bold");

Sampling: [adstock_alpha]
Sampling: []

Saturación

A continuación, examinamos el parámetro de saturación. El parámetro de saturación controla el nivel de saturación del efecto de medios. Un valor más alto de \(\lambda\) indica un nivel de saturación más alto, lo que significa que el efecto de medios se satura más rápidamente. El valor de \(\lambda\) es siempre positivo. Un valor previo razonable para comenzar es algo entre \(0\) y \(3\), centrado alrededor de \(1\).

fig, ax = plt.subplots()
saturation_lam_prior = Prior("Gamma", alpha=4, beta=2)
saturation_lam_prior.preliz.plot_pdf(ax=ax)
ax.set(xlabel="Saturation Lambda", ylabel="Density")
ax.set_title("Saturation Lambda Prior", fontsize=16, fontweight="bold");

Para los coeficientes de regresión podemos utilizar las distribuciones de costos. Sin ver los datos, es razonable esperar que cuanto mayores sean los canales (en términos de gasto), mayor será el efecto (en términos de ROAS). Esto no es una restricción, sino al menos un punto de partida (también podemos usar los mismos priors para todos los canales si no queremos hacer ninguna suposición). Utilizamos la participación del gasto como un proxy para el tamaño del efecto a través del ancho del prior (distribución normal truncada, ya que queremos asegurarnos de tener valores positivos).

channel_columns = ["x1", "x2"]

n_channels = len(channel_columns)

total_spend_per_channel = data_df[channel_columns].sum(axis=0)

spend_share = total_spend_per_channel / total_spend_per_channel.sum()
# We scale by the number of channels since if we have a lot of them then some of these
# shares become very small. This is a heuristic and not a strict rule.
prior_sigma = n_channels * spend_share.to_numpy()

fig, ax = plt.subplots(
    nrows=n_channels,
    ncols=1,
    figsize=(12, 7),
    sharex=True,
    sharey=True,
    layout="constrained",
)

for i, sigma in enumerate(prior_sigma):
    saturation_beta_prior = Prior("HalfNormal", sigma=sigma)
    saturation_beta_prior.preliz.plot_pdf(ax=ax[i])
    ax[i].set(xlabel="Saturation Beta", ylabel="Density")
    ax[i].set_title(
        f"Saturation Beta Prior for Channel {i + 1}", fontsize=14, fontweight="bold"
    )

De manera similar a lo anterior, podemos analizar las implicaciones del parámetro de saturación en el efecto mediático.

saturation = LogisticSaturation(
    priors={
        "saturation_lam": Prior("Gamma", alpha=4, beta=2),
        "saturation_beta": Prior("HalfNormal", sigma=prior_sigma),
    },
    prefix="channel",
)
prior = saturation.sample_prior(samples=1_000, random_seed=rng)
curve = saturation.sample_curve(prior)
fig, axes = saturation.plot_curve(curve, sample_kwargs={"rng": rng})
axes[0].set_title("Logistic Saturation Curve", fontsize=16, fontweight="bold");

Sampling: [channel_beta, channel_lam]
Sampling: []

Estas visualizaciones, junto con el contexto empresarial, pueden ayudarnos a establecer las prioridades para el modelo de mezcla de medios. Asegúrate de incluirlas en tu flujo de trabajo.

Priors sobre el Intercepto Variable en el Tiempo

En contraste con el cuaderno Cuaderno de ejemplo para MMM, ahora incluimos un intercepto que varía con el tiempo. Esta es una forma de modelar el componente de tendencia de los datos sin asumir ninguna forma funcional específica. Utilizamos un Proceso Gaussiano en Espacio de Hilbert (HSGP) para modelar el intercepto que varía con el tiempo. Si no está familiarizado con los HSGP, consulte el cuaderno de ejemplo de PyMC «Procesos Gaussianos: Referencia HSGP y Primeros Pasos» o el video «Una introducción conceptual y práctica a los métodos de aproximación del Proceso Gaussiano en Espacio de Hilbert (HSGP)».

Dado que estableceremos priors sobre las escalas de longitud del componente temporal, debemos tener en cuenta el número de puntos de datos que tenemos:

data_df["date_week"].unique().size

El parámetro de escala de longitud representa la distancia entre dos puntos en la serie temporal donde esperamos un cambio en la intersección. En nuestro caso particular, dado que tenemos datos semanales y esperamos capturar una tendencia a largo plazo, queremos considerar escalas de longitud de alrededor de dos o tres semanas. El prior por defecto para la escala de longitud en PyMC-Marketing es una distribución gamma inversa. Esta distribución tiene muchas propiedades deseables: (1) es positiva, (2) tiende a cero muy rápidamente (por lo que no dividimos por cero) y (3) tiene una cola larga. Para traducir la intuición de la escala de longitud en los parámetros gamma inversos, podemos utilizar una función muy útil que maximiza la entropía (para una excelente explicación y motivación de la entropía máxima, consulte el Capítulo 10 en Statistical Rethinking) de la distribución dentro de un intervalo esperado dado:

ls_prior = Prior(
    "InverseGamma",
    # We set the mean to 4 as we expect the length scale to be around two or three weeks.
    mu=4,
).constrain(
    # We set the probability range we want to cover when trying to find the best parameters.
    mass=0.94,
    # We hand the length scale interval to have a lower limit of two weeks.
    lower=2,
    # We set an upper limit of six weeks.
    upper=6,
)


ax = ls_prior.preliz.plot_pdf()
ax.set(xlabel="Length Scale", ylabel="Density")
ax.set_title("Length Scale Prior", fontsize=16, fontweight="bold");

También necesitamos especificar la amplitud del proceso gaussiano. La distribución a priori predeterminada en PyMC-Marketing es una distribución exponencial (tenga en cuenta que debe ser positiva). Dado que estamos escalando los datos, un buen valor predeterminado es establecer esto en uno.

Ahora podemos codificar estos parámetros en un objeto HSGPKwargs.

hsgp_kwargs = HSGPKwargs(
    # m is the number of basis functions to use for the approximation.
    m=100,
    # This is the range we want to cover in the time series. As we have 179 points
    # it is always good to add more for the out of sample prediction.
    L=200,
    # This is the amplitude of the Gaussian process.
    eta_lam=1.0,
    # This is the mean of the length scale prior.
    # Note we extracted from the optimization result above
    ls_mu=ls_prior.parameters["mu"],
    # This is the standard deviation of the length scale prior.
    # Note we extracted from the optimization result above
    ls_sigma=ls_prior.parameters["sigma"],
)

Configuración previa del modelo de mezcla de medios

Podemos realizar un análisis similar para el resto de los componentes y pasar nuestros priors personalizados a la clase MMM como una configuración de diccionario.

my_model_config = {
    # Intercept: This we can see by looking into the (scaled!) data.
    # This defines the mean of the time-varying intercept.
    "intercept": Prior("Normal", mu=0.5, sigma=0.1),
    # Control variables: We expect a mild trend component and effect of events.
    "gamma_control": Prior("Normal", mu=0, sigma=0.3),
    # Fourier terms: We expect a yearly seasonality.
    "gamma_fourier": Prior("Laplace", mu=0, b=0.3),
    # Media parameters (see previous sections)
    "adstock_alpha": Prior("Beta", alpha=2, beta=3),
    "saturation_beta": Prior("HalfNormal", sigma=prior_sigma),
    "saturation_lam": Prior("Gamma", alpha=4, beta=2, dims="channel"),
    # Time-varying intercept: We set the HSGP configuration.
    "intercept_tvp_config": hsgp_kwargs,
    # Likelihood: We expect the data to be normally distributed and we set the scale
    # by looking into the (scaled!) data and expected deviations from the mean.
    # Nevertheless, we know sales can not be negative, so we set a truncated normal prior.
    "likelihood": Prior(
        "TruncatedNormal", lower=0, sigma=Prior("Exponential", lam=1 / 0.3)
    ),
}

Teniendo la configuración, ahora podemos inicializar el objeto del modelo:

mmm = MMM(
    model_config=my_model_config,
    date_column="date_week",
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
    channel_columns=channel_columns,
    control_columns=["event_1", "event_2"],
    time_varying_intercept=True,
    yearly_seasonality=2,
)

Muestreo Predictivo Previo

Con el modelo inicializado, ahora podemos muestrear de la distribución predictiva a priori.

X = data_df.drop("y", axis=1)
y = data_df["y"]

# Generate prior predictive samples
_ = mmm.sample_prior_predictive(X, y, samples=2_000)

Sampling: [adstock_alpha, gamma_control, gamma_fourier, intercept_baseline, intercept_temporal_latent_multiplier_raw_eta, intercept_temporal_latent_multiplier_raw_hsgp_coefs_offset, intercept_temporal_latent_multiplier_raw_ls, saturation_beta, saturation_lam, y, y_sigma]

Comencemos por analizar las ventas esperadas (recuerde que esto es antes de ajustar el modelo).

fig, ax = plt.subplots()
mmm.plot_prior_predictive(original_scale=True, ax=ax)
ax.legend(loc="lower center", bbox_to_anchor=(0.5, -0.2), ncol=4);

El rango abarca nuestros datos observados y está en el mismo orden de magnitud (décimas). La media predictiva previa está muy correlacionada con los datos observados. Esto puede explicarse por el hecho de que, en este ejemplo particular, estamos estableciendo priors atados para los parámetros de la media y las variables de la media son las que explicarán la mayor parte de la varianza en los datos.

Siempre es bueno observar los componentes del proceso gaussiano:

fig, ax = plt.subplots()


for draw in range(5):
    ax.plot(
        mmm.prior.coords["date"],
        mmm.prior["intercept_temporal_latent_multiplier_raw"].sel(chain=0, draw=draw),
        alpha=0.8,
    )
ax.set(xlabel="date", ylabel="Intercept")
ax.set_title("Intercept Prior Predictive", fontsize=16, fontweight="bold");

Confirmamos que el nivel de fluctuaciones está en el orden de meses, como se esperaba.

A continuación, podemos analizar el ROAS del canal, que es lo que queremos entender. Simplemente los calculamos dividiendo la contribución del canal por el gasto.

prior_channel_contribution_original_scale = (
    mmm.compute_channel_contribution_original_scale(prior=True)
)

spend_sum = X[channel_columns].sum().to_numpy()

prior_roas_samples = (
    prior_channel_contribution_original_scale.sum(dim="date")
    / spend_sum[np.newaxis, np.newaxis, :]
)

fig, axes = plt.subplots(
    nrows=2, ncols=1, figsize=(12, 7), sharex=True, sharey=False, layout="constrained"
)
az.plot_posterior(prior_roas_samples, ax=axes)
axes[0].set(title="Channel $x_{1}$")
axes[1].set(title="Channel $x_{2}$", xlabel="ROAS")
fig.suptitle("ROAS Prior Distributions", fontsize=18, fontweight="bold", y=1.06);

Como era de esperar, ¡los ROAS son positivos! Además, para nuestro caso particular, el conjunto de valores plausibles es muy razonable (ya que poseemos la simulación 😉). No podemos enfatizar lo suficiente cuán importante es este paso. Si los ROAS esperados no son razonables, antes de ver los datos, es muy probable que el modelo no se ajuste bien a los datos. Además, conceptualmente, el modelo estaría muy defectuoso. Así que, si hay algo que debes llevarte de este cuaderno, por favor, siempre realiza verificaciones predictivas previas sobre las cantidades que deseas inferir 🙂.

Advertencia

Hay un pequeño detalle respecto al cálculo del ROAS, ya que no estamos considerando las contribuciones del adstock y deberíamos ser más precisos sobre los intervalos de tiempo, ver fórmula (10) en «Métodos Bayesianos para la Modelización de Mezclas de Medios con Efectos de Continuidad y Forma». No queremos perdernos en detalles aquí porque nuestro objetivo es mostrar el muestreo predictivo previo y esperamos que esta pequeña corrección tenga poco impacto en este ejemplo.

Ajuste del modelo

Estamos listos para ajustar el modelo a los datos.

# Fit the model
_ = mmm.fit(
    X=X,
    y=y,
    target_accept=0.85,
    chains=4,
    draws=500,
    random_seed=rng,
)
# Sample from the posterior predictive distribution
_ = mmm.sample_posterior_predictive(X, extend_idata=True, combined=True)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept_baseline, intercept_temporal_latent_multiplier_raw_eta, intercept_temporal_latent_multiplier_raw_ls, intercept_temporal_latent_multiplier_raw_hsgp_coefs_offset, adstock_alpha, saturation_lam, saturation_beta, gamma_control, gamma_fourier, y_sigma]

Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 25 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

Sampling: [y]

# No divergences!
mmm.idata["sample_stats"]["diverging"].sum().item()

Para obtener más diagnósticos del modelo y análisis de errores, consulte el cuaderno Cuaderno de ejemplo para MMM.

Muestreo Predictivo Posterior

Finalmente, podemos analizar la distribución predictiva posterior.

mmm.plot_posterior_predictive(original_scale=True);

Los datos de ventas ajustados se ven muy bien.

No podemos examinar los componentes del modelo:

mmm.plot_components_contributions(original_scale=True);

Vemos cómo el modelo ha capturado una leve tendencia y un agradable componente de estacionalidad anual. Los resultados son muy consistentes con el cuaderno Cuaderno de ejemplo para MMM.

Analicemos el ROAS del canal.

channel_contribution_original_scale = mmm.compute_channel_contribution_original_scale(
    prior=False
)

roas_samples = (
    channel_contribution_original_scale.sum(dim="date")
    / spend_sum[np.newaxis, np.newaxis, :]
)

fig, axes = plt.subplots(
    nrows=2, ncols=1, figsize=(12, 7), sharex=True, sharey=False, layout="constrained"
)
az.plot_posterior(roas_samples, ax=axes)
axes[0].set(title="Channel $x_{1}$")
axes[1].set(title="Channel $x_{2}$", xlabel="ROAS")
fig.suptitle("ROAS Posterior Distributions", fontsize=18, fontweight="bold", y=1.06);

Estos son los mismos que los obtenidos en el cuaderno Cuaderno de ejemplo para MMM (y, por lo tanto, cercanos a los valores verdaderos).

Un consejo final es comparar las distribuciones de ROAS anteriores y posteriores. Esto puede ayudarnos a entender el impacto de los datos en los parámetros.

axes = az.plot_forest(
    data=[prior_roas_samples, roas_samples],
    model_names=["Prior", "Posterior"],
    combined=True,
    figsize=(10, 6),
)

fig = plt.gcf()
fig.suptitle(
    "ROAS Prior vs Posterior Distributions", fontsize=18, fontweight="bold", y=1.06
);

Vemos claramente cómo los datos han reducido la distribución previa para tener una estimación posterior más precisa.

Truco

Podemos utilizar pruebas de lift para calibrar nuestro modelo en PyMC-Marketing utilizando una probabilidad de calibración adicional como se explica en Calibración de Prueba de Elevación. De esta manera, no necesitamos «adivinar» los parámetros de los medios para restringir el modelo a que coincida con el rango esperado de las pruebas de lift.

Un enfoque alternativo es parametrizar el modelo utilizando directamente el ROAS, como se sugiere en «Calibración del Modelo de Mezcla de Medios con Priors Bayesianos». Por favor, encuentre más detalles sobre estos enfoques en el cuaderno de ejemplo Mitigating Unobserved Confounders in MMMs with Lift Test Likelihoods.

---

Conclusiones clave

1. La modelización predictiva previa nos ayuda a validar nuestras suposiciones del modelo antes de utilizar datos reales.

2. La clase Prior proporciona una interfaz conveniente para:

   • Creando y visualizando distribuciones a priori

   • Muestreo de distribuciones predictivas anteriores

3. Siempre visualice sus distribuciones predictivas anteriores para asegurarse de que se alineen con su conocimiento del dominio.

Errores Comunes y Mejores Prácticas

1. Siempre Verifique Sus Escalas

   • Asegúrese de que sus priors estén en la escala correcta para sus datos.

   • Utilice el conocimiento del dominio para establecer límites razonables.

2. Consideraciones sobre las restricciones

   • Utilice restricciones cuando tenga límites claros (por ejemplo, tasas entre 0 y 1)

   • Tenga cuidado de no sobreconstruir sus priors.

3. La visualización es clave

   • Siempre grafique sus distribuciones previas

   • Verifique las muestras predictivas anteriores en comparación con el conocimiento del dominio.

Truco

La conclusión principal de este cuaderno es que las verificaciones predictivas previas son un componente fundamental del flujo de trabajo de modelado y nunca deben ser omitidas. Observe que el análisis predictivo previo permite al modelador realizar modelados incluso sin utilizar datos. Este es un gran superpoder para usar simulaciones que desafíen las suposiciones del modelo.

Nota

¿Dónde aprender más? Hay muchos recursos excelentes en línea. Una gran referencia es la Guía de Recomendaciones de Elección de Priors de Stan. Aquí puedes encontrar excelentes consejos y referencias adicionales sobre cómo establecer priors.

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

Last updated: Tue Aug 19 2025

Python implementation: CPython
Python version       : 3.12.11
IPython version      : 9.4.0

pymc_marketing: 0.15.1
pytensor      : 2.31.7

pandas        : 2.3.1
pymc_marketing: 0.15.1
matplotlib    : 3.10.3
seaborn       : 0.13.2
pymc          : 5.25.1
arviz         : 0.22.0
pymc_extras   : 0.4.0
numpy         : 2.2.6
preliz        : 0.20.0

Watermark: 2.5.0