Calibración de Prueba de Elevación#

Introducción#

Es posible que haya oído la frase «Todos los modelos son incorrectos, pero algunos son útiles.» Esto es cierto en muchas áreas, y es probable que después del primer intento, aún no haya creado un modelo que pueda determinar con precisión su verdadera atribución. Incluso si lo ha hecho, ¿cómo puede estar seguro? Creamos modelos para entender la atribución de nuestros canales de marketing, pero parece que incluso entonces, no siempre podemos confiar en lo que los modelos nos dicen.

Para garantizar que nuestros modelos se descomponen correctamente, podemos utilizar diversos métodos de prueba para recopilar datos del mundo real y compararlos con nuestros modelos. Esto nos ayudará a identificar cualquier discrepancia y mejorar la precisión de la descomposición de nuestros modelos.

Hoy, exploraremos una nueva forma de integrar experimentos en pymc-marketing. Esto nos acercará a representaciones precisas de los valores del mundo real y mejorará las estimaciones generadas por nuestros modelos.

¿Cuándo son útiles las pruebas de ascensores?#

En algunas situaciones, todos los datos de gasto relevantes de los experimentos se capturan dentro de los datos históricos que incluirá en su conjunto de entrenamiento. En estos casos, un enfoque válido podría ser simplemente entrenar su MMM con todos los datos disponibles (incluyendo períodos de experimentos y no experimentos) para informarle de manera óptima sobre los parámetros del modelo y, en última instancia, obtener mejores conocimientos. Sin embargo, hay una serie de situaciones en las que este enfoque no será apropiado.

  1. Imagina que tu plataforma de publicidad en línea te otorgara crédito adicional. Decides gastar este crédito aumentando la actividad en un canal de anuncios digitales y deseas utilizar esto como una forma de evaluar el aumento. En esta situación, habrás incrementado el nivel de publicidad que se lleva a cabo, pero este “gasto incrementado” no se carga a tu factura. Por lo tanto, a menos que tu canal de datos esté realmente en su punto, tus datos de gasto no reflejarán con precisión el nivel de gasto efectivo.

  2. Desea realizar una prueba en televisión. En ese caso, probablemente tenga que pagar meses antes de que se lleve a cabo la prueba o se emita el comercial. Por lo tanto, su gasto puede ser en febrero y la acción en marzo, pero durante el período del experimento en sí, no hay gasto. Aquí, el gasto puede estar contenido en sus datos de entrenamiento, pero el desfase temporal entre el gasto y la emisión de la publicidad es muy alto. Por lo tanto, la función de adstock que acomoda los efectos de desfase a corto plazo puede no capturar adecuadamente los desfases a largo plazo como este.

  3. Está experimentando con descuentos o promociones. Estos pueden ser cuantificables, pero no aparecer en sus canales de gasto en medios tradicionales, por lo que su experimentación puede no ser capturada por su MMM.

Estos son solo algunos ejemplos donde las pruebas de elevación pueden ser útiles. En estos casos, puede utilizar los resultados de la prueba de elevación para ajustar los parámetros del modelo y mejorar la precisión del modelo.

Requisitos#

Today, we won’t be discussing how to conduct lift tests, but instead, we will focus on their utilization. If you wish to acquire knowledge on how to generate results that are compatible with your MMM models, you can check out CausalPy for conducting experiments, such as using Interrupted Time Series for lift tests with no control groups.

Objetivo#

Después de leer este cuaderno, habrá adquirido la experiencia necesaria para incorporar los resultados (aumento detectado de sus experimentos) en nuestro modelo regresivo.

Este cuaderno mostrará utilizando el método add_lift_test_measurements de MMM y su flujo de trabajo:

  1. Construir modelo: mmm.build_model(X, y)

  2. Agregar mediciones de elevación: mmm.add_lift_test_measurements(df_lift_test)

  3. Muestra posterior: mmm.fit(X, y)

Este es un estudio de caso de dos canales correlacionados para ver cómo las pruebas de elevación ayudan a distinguir los efectos del canal.

import warnings
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 seaborn as sns
from xarray import DataArray

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM
from pymc_marketing.mmm.transformers import logistic_saturation

warnings.filterwarnings("ignore")

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, "Lift tests help distinguish channel effects"))
rng = np.random.default_rng(seed)

Generar Gastos Correlacionados y Modelar Objetivo#

Primero generaremos datos sintéticos para dos canales con gastos completamente correlacionados.

n_dates = 52
dates = pd.date_range(start="2024-01-01", periods=n_dates, freq="W-MON")
spend_rv = pm.Uniform.dist(lower=0.1, upper=1, size=n_dates)
spend = pm.draw(spend_rv, random_seed=rng)
spend1 = spend / spend.max()
spend2 = spend - 0.05
spend2 /= spend2.max()
df = pd.DataFrame(
    {
        "date": dates,
        "channel 1": spend1,
        "channel 2": spend2,
    }
)

ax = df.set_index("date").plot(ylabel="channel spend");
../../_images/4a4cf701486f5d8e1caa86d30be31ac7ab54a9de7a1e91c2aab3387b93a29648.png

Nota

Para este ejemplo, los datos de gasto sintético no se producen de una manera excesivamente realista. Sin embargo, hemos normalizado el gasto máximo a 1 para cada canal para imitar el preprocesamiento de datos que típicamente ocurre en un flujo de trabajo de pymc-marketing.

Podemos verificar que los dos canales están perfectamente correlacionados.

df.filter(regex="channel").corr()
channel 1 channel 2
channel 1 1.0 1.0
channel 2 1.0 1.0

Utilizamos la clase MMM para especificar nuestro modelo como de costumbre.

mmm = MMM(
    date_column="date",
    target_column="y",
    channel_columns=["channel 1", "channel 2"],
    adstock=GeometricAdstock(l_max=6).set_dims_for_all_priors("channel"),
    saturation=LogisticSaturation().set_dims_for_all_priors("channel"),
)

Para este ejemplo construido, estableceremos el parámetro del modelo con el operador pm.do y tomaremos una muestra aleatoria de la variable objetivo. Los parámetros fijos se encuentran a continuación, los cuales intentaremos recuperar.

En este momento, no se ha ajustado un modelo. Sin embargo, hemos creado nuestro conjunto de datos para ajustar nuestro modelo.

true_lam_c1 = 10
true_beta_c1 = 0.55

true_lam_c2 = 1.5
true_beta_c2 = 1.0

true_params = {
    "adstock_alpha": np.array([0.5, 0.5]),
    "saturation_lam": np.array([true_lam_c1, true_lam_c2]),
    "saturation_beta": np.array([true_beta_c1, true_beta_c2]),
    "intercept_contribution": 0.5,
    "y_sigma": 0.25,
}

df["y"] = np.ones(n_dates)

Hide code cell content

 
mmm.build_model(X=df.drop(columns=["y"]), y=df["y"])
fixed_model = pm.do(mmm.model, true_params)
df["y"] = pm.draw(fixed_model["y"], random_seed=rng)
del mmm.model

ax = df.set_index("date").plot(ylabel="", title="Channel spends and target variable")
ax.legend(bbox_to_anchor=(1.0, 1.05), title="column");
../../_images/76916752f891a46ac27000a11deb807519dfbf5c4f971c4f2983d7776634c5dc.png 
X = df.drop(columns=["y"])
y = df["y"]

mmm.build_model(X=X, y=y)
mmm.add_original_scale_contribution_variable(["channel_contribution"])

Estimaciones de Parámetros Indistinguibles#

Para mostrar el problema que presentan los canales completamente correlacionados, ajustemos el modelo.

fit_kwargs = {"nuts_sampler": "numpyro", "random_seed": rng}

idata_without = mmm.fit(X, y, **fit_kwargs)

There were 32 divergences after tuning. Increase `target_accept` or reparameterize.







Dado que los gastos están completamente correlacionados, no hay forma de distinguir los parámetros. No solo eso, sino que las estimaciones de los parámetros no se acercan a los valores reales.








Hide code cell content





def plot_true_value(value, channel: str, ax: plt.Axes, split: float = 0.42) -> plt.Axes:
    top = 2 * split
    ymin, ymax = (0, split) if channel == "channel 2" else (split, top)
    ax.axvline(value, ymin=ymin, ymax=ymax, color="black", linestyle="dashed")
    return ax















ax = az.plot_forest(idata_without, var_names=["saturation_lam"], combined=True)[0]

plot_true_value(true_lam_c1, "channel 1", ax=ax, split=0.4)
plot_true_value(true_lam_c2, "channel 2", ax=ax, split=0.4);








../../_images/b4c5550260dc92cd63941feca0be1451a892e673278f0bf029d4d95aa2a42f29.png










ax = az.plot_forest(idata_without, var_names=["saturation_beta"], combined=True)[0]

plot_true_value(true_beta_c1, "channel 1", ax=ax, split=0.4)
plot_true_value(true_beta_c2, "channel 2", ax=ax, split=0.4);








../../_images/107908dfccb61986f209002ed78c821f8b17c8e4e33f173f3c0d39e9ba9537b2.png






También podemos observar que, aunque las curvas de saturación reales son diferentes (porque especificamos los parámetros verdaderos para cada canal), el MMM no puede distinguir entre los dos canales.


Podemos mostrar esto a continuación trazando las curvas de saturación reales (líneas) y las curvas de respuesta directa (puntos).




Curvas de respuesta directa


Grafica las curvas de contribución directa para cada canal de marketing. El término «directo» se refiere al hecho de que graficamos costos frente a retornos inmediatos y no tomamos en cuenta los efectos rezagados de los canales, por ejemplo, las transformaciones de adstock. Por lo tanto, estas curvas consistirán en una serie de puntos, cada uno representando la contribución directa de un nivel de gasto.










Hide code cell source





def saturation_function(x, lam, beta):
    return (beta * logistic_saturation(x, lam)).eval()


step_size = 0.05
xx = DataArray(np.arange(0, spend.max() * 1.1, step_size), dims=("step",))

c1_curve_fn = partial(saturation_function, lam=true_lam_c1, beta=true_beta_c1)
c2_curve_fn = partial(saturation_function, lam=true_lam_c2, beta=true_beta_c2)

c1_curve = c1_curve_fn(xx)
c2_curve = c2_curve_fn(xx)















curve = mmm.saturation.sample_curve(
    mmm.idata.posterior[["saturation_beta", "saturation_lam"]], max_value=1.1
)
fig, axes = mmm.plot.saturation_curves(
    curve,
    original_scale=True,
    n_samples=10,
    hdi_probs=0.85,
    random_seed=rng,
    subplot_kwargs={"figsize": (12, 8), "ncols": 2},
    rc_params={
        "xtick.labelsize": 10,
        "ytick.labelsize": 10,
        "axes.labelsize": 10,
        "axes.titlesize": 10,
    },
)


def plot_actual_curves(
    axes_array: np.ndarray, linestyle: str | None = None
) -> np.ndarray:
    axes_array[0].plot(xx, c1_curve, label="channel 1", color="C0", linestyle=linestyle)
    axes_array[1].plot(xx, c2_curve, label="channel 2", color="C1", linestyle=linestyle)

    return axes_array


def plot_reference(ax: plt.Axes) -> plt.Axes:
    ax.plot(xx, xx, label="y=x", color="black", linestyle="--")

    return ax


axes = np.ravel(axes)
plot_actual_curves(axes, linestyle="dashed")
fig.suptitle(
    "Direct response curves (points) and actual saturation curves (dashed lines)"
);









Sampling: []







../../_images/72e076d17e9107c087ae209757ef0e6ccc87fd0c21366ca94cb764bac8c94cf3.png






Podemos ver que las verdaderas pero latentes curvas de saturación (líneas discontinuas) para cada canal son diferentes, pero que el MMM no puede distinguir que las curvas de saturación para cada canal (puntos) son diferentes. Y esto se debe a que los gastos están perfectamente correlacionados.


Veamos si podemos mejorar esto añadiendo mediciones de pruebas de elevación.




Nota


Para evitar confusiones con el gráfico anterior, las curvas de saturación reales están representadas por las líneas discontinuas, mientras que las curvas de contribución directa están representadas por los puntos. Si las pruebas de lift pueden ayudarnos en esta situación de alta correlación de gastos, entonces deberíamos poder detectar esto mediante: a) lograr mejores estimaciones de parámetros, y b) correspondientemente, obtener mejores estimaciones de las curvas de saturación.




También podemos visualizar las contribuciones del canal a lo largo del tiempo. Los amplios intervalos de HDI en el gráfico a continuación muestran que tenemos una incertidumbre considerable en las contribuciones del canal. También podemos ver que las medias posteriores de las contribuciones del canal están bastante alejadas de las contribuciones reales del canal. Nuevamente, estas verdaderas contribuciones del canal solo nos son conocidas en este ejemplo de datos sintéticos y no estarían disponibles para nosotros en un escenario del mundo real.








Hide code cell source





from pymc_marketing.mmm.transformers import geometric_adstock


def calculate_true_channel_contribitions(df, true_params):
    # apply geometric adstock transformation
    df["x1_adstock"] = geometric_adstock(
        x=df["channel 1"].to_xarray(),
        alpha=true_params["adstock_alpha"][0],
        l_max=8,
        normalize=True,
        dim="index",
    ).eval()
    df["x2_adstock"] = geometric_adstock(
        x=df["channel 2"].to_xarray(),
        alpha=true_params["adstock_alpha"][1],
        l_max=8,
        normalize=True,
        dim="index",
    ).eval()

    # apply saturation transformation
    df["x1_adstock_saturated"] = logistic_saturation(
        x=df["x1_adstock"].to_xarray(), lam=true_params["saturation_lam"][0]
    ).eval()

    df["x2_adstock_saturated"] = logistic_saturation(
        x=df["x2_adstock"].to_xarray(), lam=true_params["saturation_lam"][1]
    ).eval()
    return df


df = calculate_true_channel_contribitions(df, true_params)


def plot_channel_contributions(mmm, df):
    channels_contribution_original_scale = mmm.idata.posterior[
        "channel_contribution_original_scale"
    ]
    channels_contribution_original_scale_hdi = az.hdi(
        ary=channels_contribution_original_scale
    )

    get_mean_contributions_over_time_xr = channels_contribution_original_scale.mean(
        dim=["chain", "draw"]
    )

    _, ax = plt.subplots(
        nrows=2,
        figsize=(15, 8),
        ncols=1,
        sharex=True,
        sharey=False,
        layout="constrained",
    )

    for i, x in enumerate(["channel 1", "channel 2"]):
        # Estimate true contribution in the original scale from the data generating process
        sns.lineplot(
            x=df["date"],
            y=true_params["saturation_beta"][i]
            * df[f"x{i + 1}_adstock_saturated"],  # ??????
            color="black",
            label=f"{x} true contribution",
            ax=ax[i],
        )
        # HDI estimated contribution in the original scale
        ax[i].fill_between(
            x=df["date"],
            y1=channels_contribution_original_scale_hdi.sel(channel=x)[
                "channel_contribution_original_scale"
            ][:, 0],
            y2=channels_contribution_original_scale_hdi.sel(channel=x)[
                "channel_contribution_original_scale"
            ][:, 1],
            color=f"C{i}",
            label=rf"{x} $94\%$ HDI contribution",
            alpha=0.4,
        )
        # Mean estimated contribution in the original scale
        sns.lineplot(
            x=df["date"],
            y=get_mean_contributions_over_time_xr.sel(channel=x).to_numpy(),
            color=f"C{i}",
            label=f"{x} posterior mean contribution",
            alpha=0.8,
            ax=ax[i],
        )
        ax[i].legend(loc="center left", bbox_to_anchor=(1, 0.5))
        ax[i].set(title=x, ylim=(0, 1))















plot_channel_contributions(mmm, df)








../../_images/566aeb08c203613663535f9b6a696d10216764c3f05aed9382f2b4164fe26b70.png










Acerca de las Pruebas de Elevación#


En un estudio de elevación, se cambia temporalmente el presupuesto de un canal durante un período de tiempo fijo, y luego se utiliza algún método (por ejemplo, CausalPy) para inferir sobre el cambio en las ventas directamente causado por el ajuste.


Una prueba de elevación se caracteriza por:


    

  • channel: el canal que fue probado
    • 

  • x: gasto del canal de pre-prueba
    • 

  • delta_x: cambio realizado en x
    • 

  • delta_y: cambio inferido en las ventas debido a delta_x
    • 

  • sigma: desviación estándar de delta_y
    • 



Un experimento caracterizado de esta manera puede ser visto como dos puntos en la curva de saturación para el canal. En consecuencia, la calibración de la prueba de elevación se implementa añadiendo un término a la verosimilitud del modelo, que hace que la curva de saturación del canal (la contribución como función del gasto) se alinee con la observación de la prueba.






def create_lift_test_from_actual_curve(
    channel: str, x: float, delta_x: float, sigma: float
) -> dict[str, float]:
    curve_fn = c1_curve_fn if channel == "channel 1" else c2_curve_fn
    delta_y = curve_fn(x + delta_x) - curve_fn(x)

    return {
        "channel": channel,
        "x": x,
        "delta_x": delta_x,
        "delta_y": delta_y,
        "sigma": sigma,
    }













df_lift_test = pd.DataFrame(
    [
        # Channel x1
        create_lift_test_from_actual_curve("channel 1", 0.05, 0.05, 0.05),
        create_lift_test_from_actual_curve("channel 1", 0.15, 0.05, 0.05),
        create_lift_test_from_actual_curve("channel 1", 0.3, -0.05, 0.05),
        # Channel x2
        create_lift_test_from_actual_curve("channel 2", 0.5, 0.05, 0.10),
    ]
)

df_lift_test













  
    

      
      channel
      x
      delta_x
      delta_y
      sigma
    

  
  
    

      0
      channel 1
      0.05
      0.05
      0.119459
      0.05
    

    

      1
      channel 1
      0.15
      0.05
      0.069545
      0.05
    

    

      2
      channel 1
      0.30
      -0.05
      -0.031276
      0.05
    

    

      3
      channel 2
      0.50
      0.05
      0.032236
      0.10
    

  







Visualicemos los resultados de la prueba de elevación para comprender cómo pueden informar mejor los parámetros del modelo, específicamente los parámetros de la curva de saturación.








Hide code cell source





def plot_triangle(
    x,
    delta_x,
    delta_y,
    color: str,
    ax: plt.Axes,
    offset: float = 0,
    label: str | None = None,
) -> plt.Axes:
    x_after = x + delta_x

    y = offset
    y_after = y + delta_y

    ax.plot([x, x_after], [y, y], color=color, label=label)
    ax.plot([x_after, x_after], [y, y_after], color=color)
    ax.plot([x, x_after], [y, y_after], color=color, linestyle="dashed")

    return ax


def plot_channel_triangles(
    df: pd.DataFrame, color: str, ax: plt.Axes, label: str
) -> plt.Axes:
    kwargs = {"label": label}
    for _, row in df.iterrows():
        plot_triangle(
            row["x"], row["delta_x"], row["delta_y"], ax=ax, color=color, **kwargs
        )
        if "label" in kwargs:
            kwargs.pop("label")
    return ax


def plot_lift_test_triangles(df: pd.DataFrame, ax: plt.Axes) -> plt.Axes:
    idx = df["channel"] == "channel 1"
    plot_channel_triangles(df.loc[idx], color="C0", ax=ax, label="channel 1")
    plot_channel_triangles(df.loc[~idx], color="C1", ax=ax, label="channel 2")
    return ax


_, ax = plt.subplots(1, 2, sharex=True, figsize=(12, 7))

# LEFT PLOT - Channel 1 curve with its triangles
ax[0].plot(xx, c1_curve, label="channel 1 curve", color="C0", linestyle="dashed")

# Add channel 1 triangles
ch1_data = df_lift_test[df_lift_test["channel"] == "channel 1"]
for i, (_, row) in enumerate(ch1_data.iterrows()):
    label = "channel 1 lift test" if i == 0 else None
    plot_triangle(
        row["x"],
        row["delta_x"],
        row["delta_y"],
        ax=ax[0],
        color="C0",
        label=label,
        offset=c1_curve_fn(row["x"]),
    )

ax[0].set(
    ylabel="Contribution",
    title="Channel 1: Curve and lift tests",
)
ax[0].legend()

# RIGHT PLOT - Channel 2 curve with its triangles
ax[1].plot(xx, c2_curve, label="channel 2 curve", color="C1", linestyle="dashed")

# Add channel 2 triangles
ch2_data = df_lift_test[df_lift_test["channel"] == "channel 2"]
for i, (_, row) in enumerate(ch2_data.iterrows()):
    label = "channel 2 lift test" if i == 0 else None
    plot_triangle(
        row["x"],
        row["delta_x"],
        row["delta_y"],
        ax=ax[1],
        color="C1",
        label=label,
        offset=c2_curve_fn(row["x"]),
    )

ax[1].set(
    xlabel="x",
    ylabel="Contribution",
    title="Channel 2: Curve and lift tests",
)
ax[1].legend()

plt.tight_layout();










../../_images/806f62f0a7979acb85f0982b59b6c04371cd683a7789ee7d9b2c4c30c6645057.png










fig, ax = plt.subplots(1, 1, figsize=(12, 5))
plot_lift_test_triangles(df_lift_test, ax)
ax.set_title("lift test showing different diminishing returns")
ax.set_ylabel("change in contribution")
ax.set_xlabel("x")
plt.tight_layout();








../../_images/d7afb56dd17432ec284531a741c0e903b79803418caee36da1ae04b57c409291.png






El gráfico superior muestra los resultados superpuestos sobre las curvas de saturación reales. En una situación real, no conoceríamos las curvas de saturación reales, pero el gráfico aquí se incluye para ayudar a comprender lo que está sucediendo.


El gráfico inferior muestra los mismos resultados pero en y=0, centrándose en el hecho de que solo estamos observando el cambio en la contribución. Esto es más cierto en relación con el conocimiento parcial que tendrías en una situación real.


Las pruebas de elevación se visualizan con los triángulos. La base del triángulo muestra el cambio en el gasto y la altura del triángulo representa el cambio en la contribución.


Aunque solo tenemos la información en el gráfico inferior, el gráfico superior muestra cómo la información de la prueba de elevación puede utilizarse para ayudar a estimar mejor las curvas de saturación.


Por ejemplo, solo a partir de las pruebas de elevación (triángulos inferiores) podemos ver que el canal 2 tarda más en saturarse que el canal 1 porque obtenemos un mayor cambio en la contribución a niveles de gasto más altos.




Nota


Las primeras 2 pruebas de lift en el canal 1 involucraron aumentar el gasto en ese canal. Sin embargo, la tercera prueba de lift en el canal 1 involucró una disminución en el gasto. Es útil poder considerar tanto aumentos como disminuciones en el gasto en las pruebas de lift. Y no hay nada que impida hacer que esta disminución sea aún más extrema al establecer temporalmente el gasto en cero.








Agregar Pruebas de Elevación al Modelo#


Habiendo creado una instancia del modelo MMM, mmm y construido utilizando el método build_model o ajustado con el método fit, podemos agregar los resultados de la prueba de elevación al modelo utilizando el método add_lift_test_measurments.


Primero, echemos un vistazo al gráfico del modelo antes de agregar las mediciones de la prueba de elevación.






mmm.graphviz()








../../_images/5ffb12c4f5ad5c9de1048fca38902959bfe6408b68e4de7684cf278c232eef24.svg





Y ahora añadiremos las mediciones de la prueba de elevación al modelo y veremos cómo ha cambiado el gráfico de nuestro modelo.






mmm.add_lift_test_measurements(df_lift_test)
mmm.graphviz()








../../_images/86525c348ca53b583733a95611d12b517e96982b50e8d8e811db662f3eb699fb.svg





Podemos ver que el gráfico del modelo se ha modificado con una nueva observación para nuestras mediciones de elevación. Se asume que la distribución de la observación es Gamma, ya que cada curva de saturación es monótonamente creciente dado un conjunto de parámetros.


En resumen, esto funciona al proporcionar un nuevo término de probabilidad de la forma



\[::\]


donde:


    

  • \(|\Delta_y|\) es el cambio absoluto observado en la contribución, estimado a partir de la prueba de elevación.
    • 

  • \(|\tilde{\text{lift}}|\) es el cambio absoluto en la contribución según lo estimado por el modelo, es decir, basado en la curva de saturación parametrizada.
    • 

  • \(\sigma\) es la desviación estándar del aumento en \(\Delta_y\) de la prueba de elevación. Es decir, tenemos incertidumbre en el resultado de la prueba de elevación, y \(\sigma\) representa la desviación estándar de esta incertidumbre.
    • 





How Lift Tests Are Implemented in PyMC


While the add_lift_test_measurements method handles all of this automatically in pymc-marketing, understanding the implementation mechanics is valuable. This section provides a step-by-step guide to how lift tests are integrated into the model.


Why the Gamma Distribution?


The choice of a Gamma distribution for the lift test likelihood is deliberate and mathematically motivated:


    

  • Saturation curves are monotonically increasing: Given a fixed set of parameters, as spend increases, contribution increases (or stays the same, but never decreases). This is a fundamental property of saturation functions used in MMM.
    • 

  • Lifts must be non-negative: The change in contribution from a change in spend must be non-negative (when taking absolute values).
    • 

  • Gamma is natural for positive quantities: The Gamma distribution is defined only for positive real numbers, making it a natural choice for modeling non-negative lift values.
    • 

  • Handling spend decreases: By taking absolute values of both the observed and predicted lift (\(|\Delta_y|\) and \(|\tilde{\text{lift}}|\)), we can handle both increases and decreases in spend. A decrease in spend from \(x\) to \(x - \delta\) produces the same (absolute) lift magnitude as an increase from \(x - \delta\) to \(x\).
    • 



While one could alternatively use a Normal distribution with absolute values, the Gamma distribution is more theoretically appropriate for strictly positive quantities and often provides better numerical stability in practice.


Implementation Algorithm


The following steps outline how lift test observations are added to a PyMC model:


Step 1: Scale the Lift Test Data


MMMs typically work with scaled/normalized data internally (e.g., using max scaling where each variable is divided by its maximum value). Lift tests, however, are specified in original units. Therefore, the first step is to transform the lift test measurements to match the model’s internal scale:


# Pseudo-code
# Scale channel-related measurements (x, delta_x)
x_scaled = channel_scaler.transform(x)
delta_x_scaled = channel_scaler.transform(x + delta_x) - x_scaled

# Scale target-related measurements (delta_y, sigma) 
delta_y_scaled = target_scaler.transform(delta_y)
sigma_scaled = target_scaler.transform(sigma)





This ensures that the lift test observations are compatible with the model’s internal representation of channels and target.


NOTE: This is done internally, users do not implement this step.


Step 2: Map DataFrame Coordinates to Model Indices


Each row in the lift test DataFrame corresponds to specific coordinates in the model. We need to map these coordinate values to their integer indices in the model:


# Example: if df has channel values ["channel_1", "channel_2", "channel_1"]
# and the model coords are {"channel": ["channel_1", "channel_2", "channel_3"]}
# then we map to indices: [0, 1, 0]

indices = {}
for dim in required_dims:  # e.g., ["channel"] or ["channel", "geo"]
    lift_values = df_lift_test[dim].values
    model_coords = model.coords[dim]
    # Find index of each lift test value in model coordinates
    indices[dim] = [model_coords.index(val) for val in lift_values]





This coordinate mapping is essential for extracting the correct parameter values for each lift test.


Note: For a simple national-level MMM (with only a «channel» dimension), this step is straightforward—you’re just mapping channel names to indices. The complexity of this step increases for multi-dimensional models (e.g., channel × geo × product) where each lift test must be mapped across multiple coordinate dimensions simultaneously.


Step 3: Extract Parameter Values for Each Lift Test


For each lift test, we need to extract the specific parameter values that apply to that test’s coordinates. For example, if testing «channel_1», we need saturation_lam[0] and saturation_beta[0]:


# Create an indexer function that extracts parameters at specific coordinates
def get_parameter_at_indices(param_name):
    param = model[param_name]  # e.g., saturation_lam with dims=("channel",)
    dims = model.named_vars_to_dims[param_name]  # e.g., ("channel",)
    
    # Index into parameter using the coordinate indices
    idx = tuple([indices[dim] for dim in dims])
    return param[idx]

# Example: for lift tests on channels [0, 1, 0]
# get_parameter_at_indices("saturation_lam") returns [lam[0], lam[1], lam[0]]





Step 4: Evaluate the Saturation Curves


Now we compute the model’s prediction for each lift test by evaluating the saturation function at two points: before and after the spend change:


# Convert lift test data to PyTensor tensors
x_before = pt.as_tensor_variable(df_lift_test["x_scaled"])
x_after = x_before + pt.as_tensor_variable(df_lift_test["delta_x_scaled"])

# Define saturation curve evaluation with extracted parameters
def saturation_curve(x):
    return saturation_function(
        x,
        lam=get_parameter_at_indices("saturation_lam"),
        beta=get_parameter_at_indices("saturation_beta"),
    )

# Compute model-estimated lift: the key computation
model_estimated_lift = saturation_curve(x_after) - saturation_curve(x_before)





This model_estimated_lift is a PyTensor expression that depends on the model parameters, so it becomes part of the computational graph and will be different for each MCMC sample.


Step 5: Add the Likelihood Term


Finally, we add a Gamma observation to the PyMC model that links the model’s predicted lift to the empirically observed lift:


with model:
    pm.Gamma(
        name="lift_measurements",
        mu=pt.abs(model_estimated_lift),  # model's prediction
        sigma=df_lift_test["sigma_scaled"],  # measurement uncertainty
        observed=pt.abs(df_lift_test["delta_y_scaled"]),  # observed data
    )





This creates an additional likelihood term in the model. During MCMC sampling, the sampler will try to find parameter values that not only fit the main time series data but also produce saturation curves consistent with the lift test observations.


Summary


Lift tests act as additional observations of the saturation curve that constrain the model parameters during inference. Instead of only learning from the historical time series (which may have limited variation or correlated channels), the model also learns from controlled experiments that directly probe the saturation curve at specific points.


This is particularly valuable when:


    

  • Channels are highly correlated in historical data (as in this notebook’s example)
    • 

  • You want to validate that the model’s saturation curves match real-world behavior
    • 

  • Historical data has limited variation in spend levels for certain channels
    • 



By adding lift test measurements, you’re essentially saying: «I know that at spend level \(x\) for this channel, a change of \(\Delta x\) produces a contribution change of approximately \(\Delta y\).» This directly informs the saturation curve parameters and helps the model distinguish between otherwise confounded effects.




Podemos reacondicionar el modelo, pero con las pruebas de elevación incluidas.






idata_with = mmm.fit(X, y, **fit_kwargs)











Hide code cell output





There were 2 divergences after tuning. Increase `target_accept` or reparameterize.













El modelo se moldea a partir de las mediciones de la prueba de elevación y las curvas de respuesta comienzan a separarse.








Hide code cell content





def plot_channel_rug(
    df: pd.DataFrame, color: str, ax: plt.Axes, height: float
) -> plt.Axes:
    for x in df["x"].to_numpy():
        ax.axvline(x, ymin=0, ymax=height, color=color)
    return ax


def plot_lift_test_rug(df, ax, height: float = 0.05) -> plt.Axes:
    idx = df["channel"] == "channel 1"
    plot_channel_rug(df.loc[idx], color="C0", ax=ax, height=height)
    plot_channel_rug(df.loc[~idx], color="C1", ax=ax, height=height)

    return ax















curve = mmm.saturation.sample_curve(
    mmm.idata.posterior[["saturation_beta", "saturation_lam"]], max_value=1.1
)
fig, axes = mmm.plot.saturation_curves(
    curve,
    original_scale=True,
    n_samples=10,
    hdi_probs=0.85,
    random_seed=rng,
    subplot_kwargs={"figsize": (12, 8), "ncols": 2},
    rc_params={
        "xtick.labelsize": 10,
        "ytick.labelsize": 10,
        "axes.labelsize": 10,
        "axes.titlesize": 10,
    },
)


def plot_lift_test_rugs_on_axes(
    axes_array: np.ndarray, df: pd.DataFrame, height: float = 0.05
) -> np.ndarray:
    # Channel 1 lift tests on first axis
    idx_c1 = df["channel"] == "channel 1"
    if idx_c1.any():
        plot_channel_rug(df.loc[idx_c1], color="C0", ax=axes_array[0], height=height)

    # Channel 2 lift tests on second axis
    idx_c2 = df["channel"] == "channel 2"
    if idx_c2.any():
        plot_channel_rug(df.loc[idx_c2], color="C1", ax=axes_array[1], height=height)

    return axes_array


axes = np.ravel(axes)
plot_actual_curves(axes, linestyle="dashed")
plot_lift_test_rugs_on_axes(axes, df_lift_test)
fig.suptitle(
    "Direct response curves (points) and actual saturation curves (dashed lines)"
);









Sampling: []







../../_images/c134040ef99a03efcbeef50178037f173d3ccab9612eca7dea710c97b1f66d99.png






Las marcas de “alfombra” en el gráfico anterior muestran el gasto inicial del canal antes de la prueba de elevación.


A continuación, mostramos los cambios (actualmente modestos) en las estimaciones del parámetro de saturación. La idea es que a medida que añadimos más pruebas de elevación (ver más adelante), podemos estimar mejor las curvas de saturación.








Hide code cell content





def plot_comparison(data, model_names, var_name: str) -> plt.Axes:
    return az.plot_forest(
        data,
        model_names=model_names,
        var_names=[var_name],
        combined=True,
        figsize=(8, 4),
    )[0]















data = [idata_without, idata_with]
model_names = ["without lift tests", "with lift tests"]

ax = plot_comparison(data, model_names, "saturation_lam")
plot_true_value(true_lam_c1, "channel 1", ax)
plot_true_value(true_lam_c2, "channel 2", ax);








../../_images/30d177789d041c1a1edf35245da9ca38b3de3321ccb92d9d04bffa92709f47fa.png










ax = plot_comparison(data, model_names, "saturation_beta")
plot_true_value(true_beta_c1, "channel 1", ax)
plot_true_value(true_beta_c2, "channel 2", ax);








../../_images/eed7ff82eb64098318f65436abb41d99f6cae3f42b212dfdfea019d9dde336d6.png






Un examen cuidadoso de estos 2 gráficos muestra que la media posterior se ha desplazado más cerca del valor verdadero del parámetro y/o el HDI se está reduciendo, lo que indica una mejora en la precisión de nuestras estimaciones.


Volvamos a nuestras contribuciones estimadas de canal a lo largo del tiempo. Los HDI son más estrechos que antes, lo que indica que hemos reducido la incertidumbre en nuestras estimaciones de contribución de canal. Aún más impresionante, las estimaciones de la media posterior ahora están mucho más cerca de las contribuciones reales del canal.






plot_channel_contributions(mmm, df)








../../_images/3d81917c5b79494d4d02ee6ef3f7d88dc04deec6f4feb7b48ddd1b7db8bc5c2b.png










Agregar Pruebas de Elevación Adicionales#


Podemos agregar incluso más pruebas de elevación.


Pueden ser añadidos todos a la vez o por separado.






df_additional_lift_test = pd.DataFrame(
    [
        # More for Channel x1
        create_lift_test_from_actual_curve("channel 1", 0.1, 0.05, sigma=0.01),
        create_lift_test_from_actual_curve("channel 1", 0.5, 0.05, sigma=0.01),
        # More for channel x2
        create_lift_test_from_actual_curve("channel 2", 0.3, 0.05, sigma=0.01),
    ]
)

df_additional_lift_test













  
    

      
      channel
      x
      delta_x
      delta_y
      sigma
    

  
  
    

      0
      channel 1
      0.1
      0.05
      0.095167
      0.01
    

    

      1
      channel 1
      0.5
      0.05
      0.002885
      0.01
    

    

      2
      channel 2
      0.3
      0.05
      0.035354
      0.01
    

  







Utilice el parámetro name para separar los dos conjuntos de observaciones en el gráfico del modelo.






mmm.add_lift_test_measurements(df_additional_lift_test, name="more_lift_measurements")

mmm.graphviz()








../../_images/37277787401075641cfa47970cc7c7fb32a77a219dd1f9fbedd0d678acaa79fb.svg









idata_with_more = mmm.fit(X, y, **fit_kwargs)















La curva de respuesta se está desplazando cada vez más.






df_all_lift_tests = pd.concat([df_lift_test, df_additional_lift_test])

curve = mmm.saturation.sample_curve(
    mmm.idata.posterior[["saturation_beta", "saturation_lam"]], max_value=1.1
)
fig, axes = mmm.plot.saturation_curves(
    curve,
    original_scale=True,
    n_samples=10,
    hdi_probs=0.85,
    random_seed=rng,
    subplot_kwargs={"figsize": (12, 8), "ncols": 2},
    rc_params={
        "xtick.labelsize": 10,
        "ytick.labelsize": 10,
        "axes.labelsize": 10,
        "axes.titlesize": 10,
    },
)

axes = np.ravel(axes)
plot_actual_curves(axes, linestyle="dashed")
plot_lift_test_rugs_on_axes(axes, df_all_lift_tests)
fig.suptitle(
    "Direct response curves (points) and actual saturation curves (dashed lines)"
);









Sampling: []







../../_images/9f0dbd201e14b9c4d576f71ee7be1648e6046ace9e6fac698dda15a516245e97.png










data = [idata_without, idata_with, idata_with_more]
model_names = ["without lift tests", "with lift tests", "with more lift tests"]

ax = plot_comparison(data, model_names, "saturation_lam")
plot_true_value(true_lam_c1, "channel 1", ax, split=0.435)
plot_true_value(true_lam_c2, "channel 2", ax, split=0.435);








../../_images/5f5ac3b32d5bd45179ef80b6798445378dad59d14ccc4e31bfe8d6bbd077d57f.png










ax = plot_comparison(data, model_names, "saturation_beta")
plot_true_value(true_beta_c1, "channel 1", ax, split=0.435)
plot_true_value(true_beta_c2, "channel 2", ax, split=0.435);








../../_images/ec4f1eba2a646e30f6a8699663fe0db33f6692700088695ed5f1e091191f2b23.png






Así podemos ver en los 2 gráficos anteriores que cuando añadimos aún más datos de prueba de elevación a nuestro MMM, las estimaciones de los parámetros relacionadas con las curvas de saturación se están acercando a los valores verdaderos.


Ahora, por última vez, revisemos nuestras contribuciones de canal a lo largo del tiempo. Los HDI son mucho más estrechos cuando comparamos con el modelo original; añadir una serie de resultados de pruebas de elevación ha reducido nuestra incertidumbre en nuestras estimaciones de contribución de canal. Y las medias posteriores ahora están mucho más cerca de las contribuciones reales del canal en comparación con antes de que añadimos los datos de la prueba de elevación al MMM.






plot_channel_contributions(mmm, df)








../../_images/5f05f2e0cc7316e43abdca9c33d452ddef972a43b829a6e254083b2f84677994.png










Restricción basada en el Costo Por Objetivo#


Si no dispone de información experimental en esa forma, puede proporcionar el costo estimado por objetivo para cada canal (por ejemplo: costo por instalación, costo por registro, costo por acción, costo por venta unitario) y calibrar la respuesta del modelo en función de ello.






mmm.add_cost_per_target_calibration?









Signature:
mmm.add_cost_per_target_calibration(
    data: 'pd.DataFrame',
    calibration_data: 'pd.DataFrame',
    name_prefix: 'str' = 'cpt_calibration',
) -> 'None'
Docstring:
Calibrate cost-per-target using an observed Normal likelihood.

This computes cost-per-target as
``channel_data_spend / channel_contribution_original_scale`` and adds
an observed ``Normal`` likelihood for each calibration row:

``Normal(mu=cpt_mean, sigma=sigma, observed=target)``

.. versionchanged:: 0.14.0
    Replaced the manual quadratic ``Potential`` with an observed
    ``Normal`` likelihood for compatibility with ``pymc.dims``.

Parameters
----------
data : pd.DataFrame
    Feature-like DataFrame with columns matching training ``X`` but with
    channel values representing spend (original units). Must include the
    same ``date`` and any model ``dims`` columns.
calibration_data : pd.DataFrame
    DataFrame with rows specifying calibration targets. Must include:
      - ``channel``: channel name in ``self.channel_columns``
      - ``cost_per_target``: desired CPT value
      - ``sigma``: accepted deviation; larger => weaker penalty
    and one column per dimension in ``self.dims``.
cpt_variable_name : str
    Name for the cost-per-target Deterministic in the model.
name_prefix : str
    Prefix to use for generated potential names.

Examples
--------
Build a model and calibrate CPT for selected (dims, channel):

.. code-block:: python

    # spend data in original scale with the same structure as X
    spend_df = X.copy()
    # e.g., if X contains impressions, replace with monetary spend
    # spend_df[channels] = ...

    calibration_df = pd.DataFrame(
        {
            "channel": ["C1", "C2"],
            "geo": ["US", "US"],  # dims columns as needed
            "cost_per_target": [30.0, 45.0],
            "sigma": [2.0, 3.0],
        }
    )

    mmm.add_cost_per_target_calibration(
        data=spend_df,
        calibration_data=calibration_df,
        name_prefix="cpt_calibration",
    )
File:      ~/Documents/GitHub/pymc-marketing/pymc_marketing/mmm/multidimensional.py
Type:      method













avg_cost_per_target_df = pd.DataFrame(
    {
        "channel": ["channel 1"],
        "cost_per_target": [1.03],
        "sigma": [0.1],
    }
)
mmm.add_cost_per_target_calibration(
    data=df[["date", "channel 1", "channel 2"]],
    calibration_data=avg_cost_per_target_df,
)
mmm.graphviz()








../../_images/024a6e9679fc7b6b7b913888868347fbbe037c63dde322b2610845d65e0a977c.svg









idata_with_more_and_penalty = mmm.fit(X, y, **fit_kwargs)



















plot_channel_contributions(mmm, df)








../../_images/a7460f82c56cb84590bf2f4771a05356ff86ba159438668f8f5bb0d8b583bf31.png










Conclusión#


El método add_lift_test_measurements se puede utilizar para incorporar experimentos en la verosimilitud de nuestro modelo y ajustar los parámetros del modelo más cerca de los valores reales en este ejemplo.


Realizar diversos experimentos para cada canal con diferentes gastos traerá los mejores resultados.






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









Last updated: Mon Mar 16 2026

Python implementation: CPython
Python version       : 3.12.11
IPython version      : 9.6.0

pymc_marketing: 0.18.2
pytensor      : 2.38.2

pymc_marketing: 0.18.2
numpy         : 2.3.3
matplotlib    : 3.10.6
pandas        : 2.3.3
pymc          : 5.28.1
xarray        : 2025.10.1
seaborn       : 0.13.2
arviz         : 0.22.0

Watermark: 2.5.0