Introducción a la comprensión de las relaciones causales en el modelado de mezcla de medios#
La identificación causal consiste en determinar si podemos probar una relación de causa y efecto utilizando los datos que tenemos y nuestras suposiciones. Este proceso nos ayuda a establecer vínculos claros entre diferentes factores en lugar de centrarnos en sus correlaciones. Es esencial al trabajar con datos observacionales, que pueden verse afectados por factores ocultos y sesgos que dificultan la identificación de verdaderas relaciones causales.
Nota: Inspirado en el cuaderno de Bridging Causal Thinking and Media Mix Models: Carlos Trujillo & Ben Vincent.
¿Por qué es importante entender las relaciones causales en las regresiones?#
En los modelos de regresión, a menudo buscamos estimar cómo uno o más factores afectan un resultado. Sin embargo, si no consideramos la causalidad con cuidado, nuestras estimaciones pueden estar sesgadas debido a:
Sesgo confuso: Esto ocurre cuando un factor oculto influye tanto en el predictor como en el resultado, lo que lleva a conexiones engañosas.
Sesgo de Selección: Cuando utilizamos muestras no aleatorias, esto puede distorsionar las relaciones estimadas.
Sobrecontrol: Ajustar las variables que son influenciadas por el tratamiento puede llevar a estimaciones incorrectas de los efectos causales.
La identificación causal nos ayuda a ajustar las variables adecuadas en la regresión para aislar el verdadero efecto al abordar los factores de confusión y evitar ajustes innecesarios. Sin ella, podemos enfrentar varios problemas:
Sesgo en las estimaciones: Nuestra comprensión del tamaño y la dirección del efecto puede ser incorrecta, lo que puede llevar a decisiones equivocadas.
Correlaciones Espurias: Podríamos malinterpretar relaciones coincidentes como causales.
Modelos Ineficientes: Incluir variables irrelevantes o faltar variables de confusión importantes puede debilitar la precisión y claridad del modelo.
En la modelización de mezcla de medios (MMM), donde las empresas asignan presupuestos a diferentes canales para maximizar los retornos, estos problemas pueden llevar a una mala asignación de recursos, dinero desperdiciado o a atribuir erróneamente el éxito a estrategias ineficaces.
Conceptos clave en este cuaderno#
Este cuaderno abarca ideas clave en la inferencia causal, incluyendo:
Grafos Acíclicos Dirigidos Causales (DAGs): Estas son herramientas visuales que muestran las relaciones causales asumidas en los datos.
Criterio de Puerta Trasera: Esta regla ayuda a identificar qué variables bloquean los caminos que crean conexiones engañosas entre el tratamiento y el resultado.
Conjunto de Ajuste Mínimo: Este es el grupo más pequeño de factores necesarios para cumplir con el criterio de puerta trasera, asegurando estimaciones causales precisas.
Estas herramientas nos ayudan a comprender la inferencia causal de manera más clara y nos brindan pautas sólidas para elegir variables en modelos de mezcla de medios.
Al final de este cuaderno, aprenderá a utilizar pymc-marketing junto con los principios del pensamiento causal para un modelado efectivo de la mezcla de medios.
¡Comencemos importando las bibliotecas necesarias!
# avoid all warnings types
import warnings
warnings.filterwarnings("ignore")
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 graphviz import Digraph
from IPython.display import SVG, display
from pymc_extras.prior import Prior
from pymc_marketing.mmm import MMM, GeometricAdstock, MichaelisMentenSaturation
from pymc_marketing.mmm.transformers import geometric_adstock, michaelis_menten
az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 200
plt.rcParams.update({"figure.constrained_layout.use": True})
# We use a fixed random seed to ensure that our results are reproducible
seed = sum(map(ord, "Causal MMM"))
rng = np.random.default_rng(seed)
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = "retina"
Problema empresarial#
Imagina que diriges una empresa minorista que publicita sus productos. Has estado en el negocio durante un tiempo y has notado un crecimiento constante en las ventas. Durante las festividades, las ventas aumentan aún más. Sin embargo, no eres el único que publicita; un competidor también está promoviendo sus productos, probablemente con descuentos u ofertas especiales.
Exploremos cómo se relacionan entre sí diferentes factores:
Navidad (C): La temporada navideña incrementa el interés del consumidor, lo que lleva a más visualizaciones de anuncios (X1 y X2) y aumenta directamente sus ventas (T).
Canales de Marketing (X1, X2):
X1: Anuncios en redes sociales, como los de Facebook o TikTok.
X2: Anuncios en motores de búsqueda, que pueden aumentar el impacto de los anuncios en redes sociales (estos canales funcionan bien juntos).
Objetivo (T): Su meta es maximizar los ingresos por ventas.
Factores no visibles:
Ofertas de Competidores (I): Su competidor puede ofrecer descuentos agresivos durante las festividades, lo que podría disminuir el impacto de sus anuncios de búsqueda (X2) y reducir sus ventas (T).
Crecimiento del Mercado (G): El crecimiento económico durante la temporada navideña podría impulsar las ventas (T) independientemente de sus esfuerzos publicitarios.
Visualizando relaciones con un DAG#
Podemos ilustrar estas ideas utilizando un Grafo Acíclico Dirigido (DAG).
# Initialize a directed graph
dot = Digraph()
# Add nodes
dot.node("C", "Christmas", style="dashed")
dot.node("X1", "Marketing Impressions X1")
dot.node("X2", "Marketing Impressions X2")
dot.node("I", "Competitor Offers", style="dashed")
dot.node("G", "Market Growth", style="dashed")
dot.node("T", "Target")
# Add edges to represent the relationships
dot.edge("C", "X1", style="dashed")
dot.edge("C", "X2", style="dashed")
dot.edge("I", "X2", style="dashed")
dot.edge("X1", "X2")
## Variables that affect the target
dot.edge("C", "T", style="dashed")
dot.edge("X1", "T")
dot.edge("X2", "T")
dot.edge("I", "T", style="dashed")
dot.edge("G", "T", style="dashed")
# Render the graph to SVG and display it inline
svg_str = dot.pipe(format="svg")
display(SVG(svg_str))
Factores observados:#
Impresiones de Marketing X1 y X2: Dos canales (anuncios de búsqueda y anuncios en redes sociales) que pueden afectar las ventas (T). También pueden influirse mutuamente.
Objetivo (T): La métrica clave que desea mejorar, como los ingresos por ventas.
Factores no visibles (Elementos discontinuos):#
Navidad (C): Un evento estacional que influye tanto en las impresiones de marketing (X1, X2) como en las ventas (T).
Ofertas de Competidores (I): Acciones de los competidores, como descuentos, que pueden influir en los canales de marketing y en las ventas.
Crecimiento del Mercado (G): Tendencias económicas más amplias que afectan las ventas.
Nota: Los elementos en línea discontinua no se pueden intervenir, mientras que los elementos en línea continua sí se pueden intervenir.
Relaciones Clave:#
La Navidad (C) impulsa tanto las impresiones de marketing (X1, X2) como aumenta directamente las ventas (T) durante la temporada navideña.
Las impresiones de marketing (X1, X2) afectan las ventas (T) a través de sus esfuerzos publicitarios. X1 puede ayudar a X2 (por ejemplo, los anuncios en redes sociales apoyan a los anuncios de búsqueda).
Las Ofertas de Competidores (I) y el Crecimiento del Mercado (G) impactan los canales de marketing y las ventas, lo que puede confundir la comprensión de causa y efecto.
Ahora necesitamos definir cómo funcionan estas relaciones y crear un modelo para generar datos que reflejen estas suposiciones. De esta manera, podemos ver qué tan bien podemos identificar relaciones causales utilizando los datos que recopilamos.
# date range
min_date = pd.to_datetime("2022-01-01")
max_date = pd.to_datetime("2024-11-06")
date_range = pd.date_range(start=min_date, end=max_date, freq="D")
df = pd.DataFrame(data={"date_week": date_range}).assign(
year=lambda x: x["date_week"].dt.year,
month=lambda x: x["date_week"].dt.month,
dayofyear=lambda x: x["date_week"].dt.dayofyear,
)
n = df.shape[0]
print(f"Number of observations: {n}")
Number of observations: 1041
Establecimos un rango de fechas entre enero de 2022 y noviembre de 2024, con una frecuencia diaria.
Generando el Crecimiento del Mercado (Tendencia)#
A continuación, generamos los datos para el crecimiento del mercado. Asumimos que el crecimiento del mercado sigue una tendencia de ley de potencias, lo que significa que el crecimiento se acelera con el tiempo en lugar de permanecer constante. Esto se puede definir matemáticamente como:
Dónde:
\(t\): El índice de tiempo, que representa los días desde el inicio del rango de fechas.
\(baseline\): Una constante añadida a \(t\) para desplazar el punto de inicio de la tendencia. Este valor afecta el nivel inicial del crecimiento del mercado. El valor inicial de la función será \((baseline)^{exponent} - 1\), no 0.
\(exponent\): La potencia a la que se eleva el índice de tiempo, determinando la tasa a la que la tendencia se acelera con el tiempo.
df["market_growth"] = (np.linspace(start=0.0, stop=50, num=n) + 10) ** (1 / 4) - 1
fig, ax = plt.subplots()
sns.lineplot(
x="date_week", y="market_growth", color="C2", label="trend", data=df, ax=ax
)
ax.legend(loc="upper left")
ax.set(title="Trend & Seasonality Components", xlabel="date", ylabel=None);
Generando el Efecto Vacacional#
En esta sección, simulamos el efecto de las festividades en la tendencia de crecimiento del mercado. Ciertas festividades, como la Navidad, pueden tener un impacto significativo en el comportamiento del consumidor, lo que lleva a picos estacionales en las ventas. Para capturar estos efectos, introducimos una señal de festividad basada en distribuciones gaussianas (normales) centradas en fechas específicas de festividades.
La función utilizada para modelar el efecto de las vacaciones se define de la siguiente manera:
Dónde:
\(\Delta t\) es la diferencia de tiempo (en días) entre la fecha actual y la fecha de la festividad.
\(\sigma\) es la desviación estándar que controla la dispersión del efecto alrededor de la fecha de la festividad.
Para cada festividad, calculamos la señal de la festividad a lo largo del rango de fechas y añadimos una contribución de festividad escalando la señal con un coeficiente específico de la festividad. Este enfoque modela los picos estacionales de las festividades utilizando funciones gaussianas, que capturan el aumento transitorio en la actividad del mercado alrededor de las festividades y su respectivo decaimiento a lo largo del tiempo.
holiday_dates = ["24-12"] # List of holidays as month-day strings
std_devs = [25] # List of standard deviations for each holiday
holidays_coefficients = [2]
# Initialize the holiday effect array
holiday_signal = np.zeros(len(date_range))
holiday_contributions = np.zeros(len(date_range))
# Generate holiday signals
for holiday, std_dev, holiday_coef in zip(
holiday_dates, std_devs, holidays_coefficients, strict=False
):
# Find all occurrences of the holiday in the date range
holiday_occurrences = date_range[date_range.strftime("%d-%m") == holiday]
for occurrence in holiday_occurrences:
# Calculate the time difference in days
time_diff = (date_range - occurrence).days
# Generate the Gaussian basis for the holiday
_holiday_signal = np.exp(-0.5 * (time_diff / std_dev) ** 2)
# Add the holiday signal to the holiday effect
holiday_signal += _holiday_signal
holiday_contributions += _holiday_signal * holiday_coef
df["holiday_signal"] = holiday_signal
df["holiday_contributions"] = holiday_contributions
# Plot the holiday effect
fig, ax = plt.subplots()
sns.lineplot(x=date_range, y=holiday_signal, ax=ax)
ax.set(title="Holiday Effect Signal", xlabel="Date", ylabel="Signal Intensity")
plt.show()
Modelando Ofertas de Competidores#
Para simular el impacto de las ofertas de los competidores en el crecimiento del mercado.
Dónde:
\(A\): La amplitud de la oscilación, que representa la desviación máxima del centro. En este contexto, controla la intensidad de las ofertas de los competidores.
\(C\): El centro de oscilación, que establece el nivel promedio alrededor del cual fluctúan las ofertas de los competidores.
\(\omega\): La frecuencia angular, que determina la frecuencia de las oscilaciones. Se calcula como: \(\omega = \frac{\pi}{\frac{n}{2}}\) donde \(n\) es el número total de puntos de datos. Esto asegura un ciclo completo dentro del rango especificado.
A = 0.5 # Amplitude
C = 2.5 # Center of the oscillation
omega = np.pi / (n / 2)
df["competitor_offers"] = -A * np.cos(omega * df.index) + C
# plot the competition price and competitor_offers
fig, ax = plt.subplots(
nrows=1, ncols=1, figsize=(10, 7), sharex=True, layout="constrained"
)
sns.lineplot(x="date_week", y="competitor_offers", data=df, color="C3", ax=ax)
ax.set(title="Competitor Offers", xlabel="date", ylabel=None);
Con este formulario, simulamos un mundo donde las ofertas de los competidores aumentan durante los primeros días del año y luego disminuyen hacia el final del año.
Modelando Canales de Marketing#
En esta sección, simulamos dos canales de marketing, \(x1\) y \(x2\), que representan diferentes canales publicitarios (por ejemplo, anuncios de búsqueda y campañas en redes sociales). El comportamiento de cada canal está influenciado por la variabilidad aleatoria y los efectos confusos de las festividades estacionales y las acciones de los competidores. Así es como modelamos cada canal matemáticamente:
Canal \(x1\): Como se mencionó anteriormente, generamos \(x1\) que se ve afectado por la señal de vacaciones, podríamos definirlo como:
Canal \(x2\): Por otro lado, generamos \(x2\) que se ve afectado por la señal de vacaciones, la influencia de \(x1\) y las ofertas de los competidores. Podríamos definirlo como:
Interpretación#
Estas ecuaciones nos permiten capturar la dinámica compleja que influye en cada canal de marketing:
Efectos de Vacaciones aumentan la actividad del canal alrededor de fechas específicas, simulando picos estacionales.
Influencias entre canales introducen interdependencias, modelando cómo el éxito de un canal puede amplificar el de otro.
Efectos de la competencia introducen fluctuaciones debido a factores externos, mostrando cómo la presión competitiva puede disminuir el compromiso o la efectividad.
Este enfoque refleja las condiciones del mundo real, donde los canales de marketing no están aislados, sino que se ven afectados por tendencias estacionales, acciones de los competidores y sinergias entre diferentes canales.
x1 = pz.Normal(mu=5, sigma=3).rvs(n, random_state=rng)
cofounder_effect_holiday_x1 = 2.5
x1_conv = np.convolve(x1, np.ones(14) / 14, mode="same")
# Replace first and last 14 values with mean + noise
noise = pz.Normal(mu=0, sigma=0.1).rvs(28, random_state=rng)
x1_conv[:14] = x1_conv.mean() + noise[:14]
x1_conv[-14:] = x1_conv.mean() + noise[14:]
df["x1"] = x1_conv + (holiday_signal * cofounder_effect_holiday_x1)
x2 = pz.Normal(mu=5, sigma=2).rvs(n, random_state=rng)
cofounder_effect_holiday_x2 = 2.2
cofounder_effect_x1_x2 = 1.3
cofounder_effect_competitor_offers_x2 = -0.7
x2_conv = np.convolve(x2, np.ones(18) / 12, mode="same")
# Replace first and last 14 values with mean + noise
noise = pz.Normal(mu=0, sigma=0.1).rvs(28, random_state=rng)
x2_conv[:14] = x2_conv.mean() + noise[:14]
x2_conv[-14:] = x2_conv.mean() + noise[14:]
df["x2"] = (
x2_conv
+ (holiday_signal * cofounder_effect_holiday_x2)
+ (df["x1"] * cofounder_effect_x1_x2)
+ (df["competitor_offers"] * cofounder_effect_competitor_offers_x2)
)
fig, ax = plt.subplots(
nrows=2, ncols=1, figsize=(10, 7), sharex=True, layout="constrained"
)
sns.lineplot(x="date_week", y="x1", data=df, color="C0", ax=ax[0])
sns.lineplot(x="date_week", y="x2", data=df, color="C1", ax=ax[1])
ax[1].set(xlabel="date")
fig.suptitle("Media Costs Data", fontsize=16);
Las dos señales generadas representarán las «impresiones» de los dos canales de marketing. Básicamente, estamos simulando la exposición que observaríamos de una plataforma de medios, reportada después de ejecutar una campaña.
El nivel de exposición no impacta nuestras ventas de manera lineal, por lo que necesitamos aplicar una transformación para tenerlo en cuenta. Utilizaremos una transformación de adstock geométrico, que es una técnica común en la modelización de mezcla de marketing para considerar el efecto de adstock/retraso, y una saturación de michaelis-menten para tener en cuenta el efecto de rendimientos decrecientes.
Similar al ejemplo inicial en pymc-marketing.
# apply geometric adstock transformation
alpha1: float = 0.6
alpha2: float = 0.2
df["x1_adstock"] = geometric_adstock(
x=df["x1"].to_xarray(), alpha=alpha1, l_max=24, normalize=True, dim="index"
).eval()
df["x2_adstock"] = geometric_adstock(
x=df["x2"].to_xarray(), alpha=alpha2, l_max=24, normalize=True, dim="index"
).eval()
# apply saturation transformation
lam1: float = 5.0
lam2: float = 9.0
alpha_mm1: float = 6
alpha_mm2: float = 12
df["x1_adstock_saturated"] = michaelis_menten(
x=df["x1_adstock"].to_xarray(), lam=lam1, alpha=alpha_mm1
).eval()
df["x2_adstock_saturated"] = michaelis_menten(
x=df["x2_adstock"].to_xarray(), lam=lam2, alpha=alpha_mm2
).eval()
fig, ax = plt.subplots(
nrows=3, ncols=2, figsize=(16, 9), sharex=True, sharey=False, layout="constrained"
)
sns.lineplot(x="date_week", y="x1", data=df, color="C0", ax=ax[0, 0])
sns.lineplot(x="date_week", y="x2", data=df, color="C1", ax=ax[0, 1])
sns.lineplot(x="date_week", y="x1_adstock", data=df, color="C0", ax=ax[1, 0])
sns.lineplot(x="date_week", y="x2_adstock", data=df, color="C1", ax=ax[1, 1])
sns.lineplot(x="date_week", y="x1_adstock_saturated", data=df, color="C0", ax=ax[2, 0])
sns.lineplot(x="date_week", y="x2_adstock_saturated", data=df, color="C1", ax=ax[2, 1])
fig.suptitle("Media Costs Data - Transformed", fontsize=16)
# adjust size of X axis
ax[2, 0].tick_params(axis="x", labelsize=8)
ax[2, 1].tick_params(axis="x", labelsize=8);
Definiendo las Ventas (Objetivo)#
Finalmente, definimos la variable de ventas o ingresos (objetivo), que se ve afectada por el crecimiento del mercado, las ofertas de los competidores, las contribuciones de las festividades, los costos de medios saturados y un término de ruido.
La variable, \(y\), se define de la siguiente manera:
Dónde:
Intercepto: Un nivel base de ventas, establecido en 1.5, que representa el nivel base de ventas en ausencia de otros efectos.
Crecimiento del Mercado: Representa la tendencia subyacente en las ventas debido al crecimiento a largo plazo, con un coeficiente implícito de 1, añadiendo una influencia ascendente constante.
Ofertas de Competidores: Reduce directamente las ventas, con un coeficiente implícito de -1, simulando el impacto negativo de las acciones de los competidores en la variable objetivo.
Contribuciones Festivas: Añade picos de ventas durante los períodos festivos, capturando el aumento estacional en la demanda del consumidor.
\(m(Impressions_{x1_t})\) y \(m(Impressions_{x2_t})\): Representan los valores de adstock saturado para los canales de marketing \(x1\) y \(x2\).
Ruido \(\epsilon\): Un pequeño término de error aleatorio, extraído de una distribución normal con media 0 y desviación estándar 0.08, para tener en cuenta la variabilidad no explicada en las ventas.
df["intercept"] = 1.5
df["epsilon"] = rng.normal(loc=0.0, scale=0.08, size=n)
df["y"] = (
df["intercept"]
+ df["market_growth"] # implicit coef 1
- df["competitor_offers"] # explicit coef -1
+ df["holiday_contributions"]
+ df["x1_adstock_saturated"]
+ df["x2_adstock_saturated"]
+ df["epsilon"] # Noise
)
fig, ax = plt.subplots()
sns.lineplot(x="date_week", y="y", color="black", data=df, ax=ax)
ax.set(title="Sales (Target Variable)", xlabel="date", ylabel="y (thousands)");
columns_to_keep = [
"date_week",
"y",
"x1",
"x2",
]
data = df[columns_to_keep].copy()
Hagamos un resumen de lo que hemos hecho hasta ahora:
Creamos una hipótesis sobre cómo interactúan las cosas en el mundo y cómo están conectadas esas interacciones.
Convertimos nuestra hipótesis en un modelo matemático que refleja estas conexiones.
Generamos datos que se ajustan a nuestra hipótesis.
Ahora estamos listos para encontrar el modelo que mejor explique nuestros datos.
En términos bayesianos, llamamos a esto el proceso generador de datos, que corresponde al modelo causal estructural en términos causales.
Todos los modelos que escribimos en PyMC son modelos causales estructurales. Incluso si no buscamos una interpretación causal, nuestro modelo aún busca explicar un proceso específico de generación de datos.
Reconocer esto es crucial porque influye en nuestras decisiones al construir el modelo y elegir qué variables incluir. Después de nuestro trabajo anterior, deberíamos tener una buena idea de cuáles son las variables más importantes para explicar las ventas. Sin embargo, si olvidamos lo que discutimos anteriormente, no sabremos cómo empezar a construir el modelo.
El verdadero desafío#
En la vida real, no conocemos el proceso real que genera los datos; solo tenemos datos observados. A menudo, recibimos un conjunto de datos como el que se muestra a continuación.
data.head()
| date_week | y | x1 | x2 | |
|---|---|---|---|---|
| 0 | 2022-01-01 | 8.435918 | 4.901673 | 12.584369 |
| 1 | 2022-01-02 | 9.619534 | 5.097116 | 12.987150 |
| 2 | 2022-01-03 | 9.885896 | 4.930238 | 12.610042 |
| 3 | 2022-01-04 | 9.959450 | 4.921401 | 12.842278 |
| 4 | 2022-01-05 | 10.214027 | 4.894086 | 12.712858 |
Si tenemos PyMC instalado, el enfoque más simple sería incluir todas nuestras variables en un modelo y ejecutarlo. Esto es similar a utilizar cualquier modelo disponible en el mercado, verificando correlaciones y seleccionando las variables más relevantes. Este método significa esencialmente: «Controlaremos todo lo que pueda relacionarse con las ventas.»
Nota: Este enfoque no es malo, pero no utiliza información valiosa y conocimientos previos sobre el negocio y sus relaciones internas.
Podemos separar la variable objetivo de los factores potenciales que la influyen. Luego, podemos comenzar a entrenar nuestro modelo inicial basado en correlaciones.
X = data.drop("y", axis=1)
y = data["y"]
# sampling options for PyMC
sample_kwargs = {"draws": 500, "chains": 4, "nuts_sampler": "numpyro"}
correlational_mmm = MMM(
sampler_config=sample_kwargs,
date_column="date_week",
adstock=GeometricAdstock(l_max=24),
saturation=MichaelisMentenSaturation(),
channel_columns=["x1", "x2"],
)
Con las pequeñas líneas de código anteriores, nuestro modelo de mezcla de medios está configurado y ahora podemos comenzar a entrenarlo.
# sampling options for PyMC
correlational_mmm.fit(X=X, y=y, target_accept=0.90, random_seed=rng)
correlational_mmm.sample_posterior_predictive(
X, extend_idata=True, combined=True, random_seed=rng
)
There were 1 divergences after tuning. Increase `target_accept` or reparameterize.
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
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Sampling: [y]
<xarray.Dataset> Size: 17MB
Dimensions: (sample: 2000, date: 1041)
Coordinates:
* date (date) datetime64[ns] 8kB 2022-01-01 2022-01-02 ... 2024-11-06
* sample (sample) object 16kB MultiIndex
* chain (sample) int64 16kB 0 0 0 0 0 0 0 0 0 0 0 ... 3 3 3 3 3 3 3 3 3 3 3
* draw (sample) int64 16kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
Data variables:
y (date, sample) float64 17MB 7.058 6.55 6.831 ... 10.57 10.81 11.36
Attributes:
created_at: 2025-06-16T13:23:01.103002+00:00
arviz_version: 0.21.0
inference_library: pymc
inference_library_version: 5.23.0Dado que PyMC es estructural, podemos representar nuestro modelo como un grafo acíclico dirigido (DAG). Este grafo nos ayuda a ver cómo se relacionan los parámetros y las variables entre sí y a comprender nuestras decisiones en el proceso de modelado.
Podemos utilizar la función model_to_graphviz para crear esta representación visual.
pm.model_to_graphviz(model=correlational_mmm.model)
Nuestro DAG es bastante simple. Incluye nuestros canales, que tienen dos dimensiones: fecha y canal. También tiene una variable de intercepción que captura las ventas o ingresos base, que son factores fuera de nuestro control, y un término para reflejar la variación aleatoria.
En términos más simples, podemos expresar la función subyacente de este modelo de mezcla de medios como:
Aquí:
\(\mu\) representa la variable de resultado predicha (como ventas o ingresos) basada en los efectos combinados de nuestras variables de entrada.
\(\beta_0\) es el coeficiente para la variable base.
\(m(\text{x1}, \theta_{\text{x1}})\) y \(m(\text{x2}, \theta_{\text{x2}})\) son transformaciones no lineales de \(x1\) y \(x2\) que capturan sus impactos variables.
Basado en esto, creemos que esta función nos ayudará a entender el impacto de \(x1\) y \(x2\). Así, podemos esbozar el DAG causal que creemos que es correcto.
# Initialize a directed graph
naive_causal_mmm_graph = Digraph()
# Add nodes
naive_causal_mmm_graph.node("X1", "Marketing Impressions X1")
naive_causal_mmm_graph.node("X2", "Marketing Impressions X2")
naive_causal_mmm_graph.node("E", "Exogenous variables", style="dashed")
naive_causal_mmm_graph.node("T", "Target")
naive_causal_mmm_graph.edge("E", "T", style="dashed")
naive_causal_mmm_graph.edge("X1", "T")
naive_causal_mmm_graph.edge("X2", "T")
# Render the graph to SVG and display it inline
svg_str = naive_causal_mmm_graph.pipe(format="svg")
display(SVG(svg_str))
¿Cómo podemos inferir este DAG de nuestro modelo actual?
Durante el desarrollo de nuestro modelo, establecimos una estructura y un flujo específicos para nuestros datos. Concluimos que los impactos de nuestros canales operan de manera independiente entre sí. Además, determinamos que si falta algún componente de nuestro ecosistema, su influencia será contabilizada por el término base debido a esta ecuación. Como puede ver, incluso al adoptar este modelo básico, estamos haciendo suposiciones significativas.
Por un lado, usted está asumiendo que el impacto no es lineal al aplicar estas transformaciones, y está sugiriendo que el impacto es positivo y que puede haber un retraso máximo de un cierto número de días.
Incluso has definido la dirección de tus relaciones. Al definir estas relaciones y asumir que no hay conexiones causales directas entre nuestras variables, podemos concluir que, si la naturaleza de su relación está representada con precisión por la ecuación proporcionada, entonces al controlar los canales relevantes, podríamos descubrir sus verdaderos efectos.
Esto nos lleva a cuál DAG causal asumimos que es correcto, basado en nuestras suposiciones anteriores. Si reconoces este proceso, ¡felicitaciones! Has creado un modelo generativo o Modelo Causal Estructural, con una Ecuación Causal Estructural, utilizando PyMC-marketing.
Nota: Como puede ver, basándonos en nuestras suposiciones y las variables que incluimos, podemos retroceder e inferir el DAG causal que estamos asumiendo inherentemente. Este proceso de retroceso puede ser utilizado y formalizado en algoritmos de descubrimiento causal para inferir automáticamente los varios DAG causales en función de cómo nuestro modelo de regresión y sus estimaciones responden a los datos.
Sin embargo, este DAG causal no representa el verdadero DAG causal. Dado que nuestro modelo PyMC es estructural y causal, debemos preguntar: ¿Qué sucede si creo un modelo con una estructura causal diferente a la real?
# Number of diverging samples
correlational_mmm.idata["sample_stats"]["diverging"].sum().item()
Porque todos los modelos de PyMC son modelos generativos, son capaces de «hablar» sobre las relaciones causales que están presentes en el proceso de generación de datos. Si tenemos un modelo que no es consistente con las verdaderas relaciones causales, el modelo fallará durante la exploración del espacio de parámetros.
Por ejemplo, si asume que la variable A causa la variable B (es decir, A -> B), pero en realidad, B causa A o están confundidas por una variable no observada, la estructura de su modelo representa incorrectamente las relaciones reales. Este desajuste crea inconsistencias entre el modelo y los datos, lo que resulta en una distribución posterior que es difícil de explorar para el muestreador. La posterior puede exhibir características patológicas como multimodalidad, fuertes correlaciones entre parámetros o regiones de alta curvatura.
Un modelo como este presentará divergencias, porque estas ocurren cuando el integrador numérico utilizado en el muestreo de Monte Carlo Hamiltoniano (HMC) no puede simular con precisión la trayectoria a través del espacio de parámetros, a menudo debido a gradientes pronunciados o discontinuidades en el posterior. Si lo piensas, suposiciones causales incorrectas pueden forzar al muestreador a entrar en regiones incompatibles del espacio de parámetros, lo que provoca que la integración numérica falle.
Cuando evaluamos modelos generativos de series temporales, particularmente modelos bayesianos, siempre debemos evaluar los modelos con métricas adicionales además de las métricas de ajuste puro, como divergencias, rhat, ess, etc.
Observemos el ajuste de nuestro modelo actual 👀
correlational_mmm.plot_posterior_predictive(add_mean=False, original_scale=True)
r2 = az.r2_score(
y_true=df["y"].values,
y_pred=correlational_mmm.idata.posterior_predictive.stack(sample=("chain", "draw"))[
"y"
].values.T
* correlational_mmm.target_transformer["scaler"].scale_.item(),
).iloc[0]
plt.text(
0.05,
0.95,
f"R2: {r2:.2f}",
transform=plt.gca().transAxes,
fontsize=12,
verticalalignment="top",
bbox=dict(facecolor="white", alpha=0.8, edgecolor="none"),
)
plt.title("Estimated Target Variable Over Time")
plt.xlabel("Days")
plt.ylabel("Target Variable")
plt.legend()
plt.grid(True)
plt.show()
Genial, ¡nuestro modelo se ve increíble! A primera vista, la evaluación de nuestro modelo parece bastante satisfactoria. El valor de \(R^2\) es decente y la incertidumbre en los intervalos creíbles es manejable.
Sin embargo, la pregunta apremiante sigue siendo: ¿estamos realmente capturando las señales causales subyacentes?
initial_model_recover_effect = (
az.hdi(correlational_mmm.fit_result["channel_contribution"], hdi_prob=0.95)
* correlational_mmm.target_transformer["scaler"].scale_.item()
)
initial_model_mean_effect = (
correlational_mmm.fit_result.channel_contribution.mean(dim=["chain", "draw"])
* correlational_mmm.target_transformer["scaler"].scale_.item()
)
fig, ax = plt.subplots(2, 1, figsize=(14, 10), sharex=True)
# Social media
ax[0].plot(
date_range,
initial_model_mean_effect.sel(channel="x1"),
label="Mean Recover x1 Effect",
linestyle="--",
color="blue",
)
ax[0].fill_between(
date_range,
initial_model_recover_effect.channel_contribution.isel(hdi=0).sel(channel="x1"),
initial_model_recover_effect.channel_contribution.isel(hdi=1).sel(channel="x1"),
alpha=0.2,
label="95% Credible Interval",
color="blue",
)
ax[0].plot(
date_range, df["x1_adstock_saturated"], label="Real x1 Effect", color="black"
)
# google
ax[1].plot(
date_range,
initial_model_mean_effect.sel(channel="x2"),
label="Mean Recover x2 Effect",
linestyle="--",
color="orange",
)
ax[1].fill_between(
date_range,
initial_model_recover_effect.channel_contribution.isel(hdi=0).sel(channel="x2"),
initial_model_recover_effect.channel_contribution.isel(hdi=1).sel(channel="x2"),
alpha=0.2,
label="95% Credible Interval",
color="orange",
)
ax[1].plot(
date_range, df["x2_adstock_saturated"], label="Real x2 Effect", color="black"
)
# formatting
ax[0].legend()
ax[1].legend()
plt.grid(True)
ax[1].set(xlabel="date")
fig.suptitle("Media Contribution Recovery", fontsize=16)
plt.show()
Nuestro modelo tiene dificultades para capturar con precisión el verdadero efecto causal de las variables mediáticas, lo que sugiere que puede estar mal especificado.
En este contexto, aunque podemos utilizar el modelo verdadero para comparar los efectos reales con los efectos recuperados del modelo, las situaciones de la vida real a menudo carecen de esta claridad. Por lo tanto, necesitamos métodos para identificar y comparar los efectos recuperados con los verdaderos.
Este es el lugar donde CausalPy resulta útil. Permite la creación de cuasi-experimentos para determinar el impacto causal de las variables mediáticas, lo que permite la comparación entre los resultados de nuestro modelo y los de estos experimentos. Puede aprender más sobre CausalPy aquí.
Utilizando la Identificación Causal#
En lugar de agregar variables al azar a nuestra regresión, deberíamos utilizar un marco causal para construir nuestro modelo. Un lenguaje que nos ayude a traducir nuestras ideas e hipótesis en una representación gráfica que podamos interrogar. Al hacerlo, podemos entender si estamos pasando por alto una variable importante.
naive_causal_dag_str = """
digraph {
x1 -> y;
x2 -> y;
holiday_signal -> y;
}
"""
naive_causal_mmm = MMM(
sampler_config=sample_kwargs,
date_column="date_week",
adstock=GeometricAdstock(l_max=24),
saturation=MichaelisMentenSaturation(),
channel_columns=["x1", "x2"],
control_columns=["holiday_signal"],
# new columns
outcome_node="y",
dag=naive_causal_dag_str,
)
Hemos añadido una nueva columna llamada outcome_node para nuestra variable objetivo y definido el parámetro dag como una representación del Grafo Acíclico Dirigido (DAG) causal. Incluimos holiday_signal como una variable de control en nuestro modelo porque influye significativamente en la variable objetivo y puede haber sido pasada por alto anteriormente. Después de ejecutar eso, recibimos una advertencia de que las columnas de canal se están tratando como variables de tratamiento por defecto, lo que indica que son nuestro enfoque para estimar el verdadero impacto causal (comportamiento por defecto).
Ahora, podemos examinar el conjunto de variables de ajuste necesarias para revelar el verdadero efecto causal de nuestras variables de tratamiento.
(
naive_causal_mmm.causal_graphical_model.adjustment_set,
naive_causal_mmm.causal_graphical_model.minimal_adjustment_set,
)
([], ['x1', 'x2'])
Excelente, hemos obtenido una nueva pieza de información importante. Nos falta la variable holiday_signal en el conjunto de ajustes, lo que significa que, para encontrar el verdadero efecto causal de x1 y x2, no necesitamos controlar por holiday_signal. Entonces, ¿por qué no pudimos recuperar el verdadero efecto causal de x1 y x2 antes? Quizás el DAG especificado no era verdadero, lo que significa que es posible que la relación entre x1 y x2 esté confundida por la variable holiday_signal?
Veamos si este es el caso.
naive_causal_dag_str_v2 = """
digraph {
x1 -> y;
x2 -> y;
holiday_signal -> y;
holiday_signal -> x1;
holiday_signal -> x2;
}
"""
naive_causal_mmm_v2 = MMM(
sampler_config=sample_kwargs,
date_column="date_week",
adstock=GeometricAdstock(l_max=24),
saturation=MichaelisMentenSaturation(),
channel_columns=["x1", "x2"],
control_columns=["holiday_signal"],
# new columns
outcome_node="y",
dag=naive_causal_dag_str_v2,
)
(
naive_causal_mmm_v2.causal_graphical_model.adjustment_set,
naive_causal_mmm_v2.causal_graphical_model.minimal_adjustment_set,
)
(['holiday_signal'], ['holiday_signal', 'x1', 'x2'])
Ahora, con un nuevo conjunto de ajustes, vemos que si no controlamos por holiday_signal, no podemos identificar el verdadero efecto causal de x1 y x2, ya que están influenciados por esta variable. Esto se asemeja más a lo que ocurrió antes, después de comparar la verdadera contribución de x1 y x2 con la contribución recuperada.
Nota: Definimos este conjunto de ajuste utilizando la clase
CausalGraphModel. Dependemos dedowhypara ayudarnos a identificar el conjunto de ajuste. Puede aprender más sobredowhyaquí.
Utilizamos el criterio backdoor para elegir el conjunto de ajuste. Esto significa que buscamos variables que, al controlarlas, bloquean todos los caminos de retroceso entre x1 y x2. Estos caminos pueden crear vínculos falsos entre las variables de tratamiento y resultado. Si no abordamos estos caminos, nuestras estimaciones causales pueden estar equivocadas. Para entender el efecto causal con precisión, necesitamos bloquear estos caminos ajustando las variables de confusión. Puede leer más sobre el criterio backdoor en el libro de Aleksander Molak aquí.
Ahora, digamos que creamos una nueva variable llamada holiday_signal para controlar el efecto de las festividades. Agregar esta variable a nuestro conjunto de datos debería ayudarnos a encontrar el verdadero efecto causal de x1 y x2, ¿verdad?
data["holiday_signal"] = holiday_signal
data.head()
| date_week | y | x1 | x2 | holiday_signal | |
|---|---|---|---|---|---|
| 0 | 2022-01-01 | 8.435918 | 4.901673 | 12.584369 | 5.244234e-45 |
| 1 | 2022-01-02 | 9.619534 | 5.097116 | 12.987150 | 9.276916e-45 |
| 2 | 2022-01-03 | 9.885896 | 4.930238 | 12.610042 | 1.638439e-44 |
| 3 | 2022-01-04 | 9.959450 | 4.921401 | 12.842278 | 2.889097e-44 |
| 4 | 2022-01-05 | 10.214027 | 4.894086 | 12.712858 | 5.086267e-44 |
X = data.drop("y", axis=1)
y = data["y"]
naive_causal_mmm_v2.fit(X=X, y=y, target_accept=0.90, random_seed=rng)
naive_causal_mmm_v2.sample_posterior_predictive(
X, extend_idata=True, combined=True, random_seed=rng
)
There were 49 divergences after tuning. Increase `target_accept` or reparameterize.
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
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Sampling: [y]
<xarray.Dataset> Size: 17MB
Dimensions: (sample: 2000, date: 1041)
Coordinates:
* date (date) datetime64[ns] 8kB 2022-01-01 2022-01-02 ... 2024-11-06
* sample (sample) object 16kB MultiIndex
* chain (sample) int64 16kB 0 0 0 0 0 0 0 0 0 0 0 ... 3 3 3 3 3 3 3 3 3 3 3
* draw (sample) int64 16kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
Data variables:
y (date, sample) float64 17MB 7.605 7.466 7.616 ... 10.71 10.41 10.7
Attributes:
created_at: 2025-06-16T13:27:28.342091+00:00
arviz_version: 0.21.0
inference_library: pymc
inference_library_version: 5.23.0Tenemos aún más divergencias 🤪, ¿qué está sucediendo aquí? ¿Todo está peor? ¿Incluso si añadimos una nueva variable para controlar el holiday_signal? Esto se ve extraño, porque estamos añadiendo una variable que debería ayudarnos a identificar el verdadero efecto causal de x1 y x2.
Veamos el ajuste de nuestro modelo 👀
naive_causal_mmm_v2.plot_posterior_predictive(add_mean=False, original_scale=True)
r2 = az.r2_score(
y_true=df["y"].values,
y_pred=naive_causal_mmm_v2.idata.posterior_predictive.stack(
sample=("chain", "draw")
)["y"].values.T
* naive_causal_mmm_v2.target_transformer["scaler"].scale_.item(),
).iloc[0]
plt.text(
0.05,
0.95,
f"R2: {r2:.2f}",
transform=plt.gca().transAxes,
fontsize=12,
verticalalignment="top",
bbox=dict(facecolor="white", alpha=0.8, edgecolor="none"),
)
plt.title("Estimated Target Variable Over Time")
plt.xlabel("Days")
plt.ylabel("Target Variable")
plt.legend()
plt.grid(True)
plt.show()
De alguna manera, \(R^2\) es más alto que antes, ¿podría ser esto una buena señal? ¿Quizás estamos recuperando el verdadero efecto causal de x1 y x2?
initial_model_recover_effect = (
az.hdi(naive_causal_mmm_v2.fit_result["channel_contribution"], hdi_prob=0.95)
* naive_causal_mmm_v2.target_transformer["scaler"].scale_.item()
)
initial_model_mean_effect = (
naive_causal_mmm_v2.fit_result.channel_contribution.mean(dim=["chain", "draw"])
* naive_causal_mmm_v2.target_transformer["scaler"].scale_.item()
)
fig, ax = plt.subplots(2, 1, figsize=(14, 10), sharex=True)
# Social media
ax[0].plot(
date_range,
initial_model_mean_effect.sel(channel="x1"),
label="Mean Recover x1 Effect",
linestyle="--",
color="blue",
)
ax[0].fill_between(
date_range,
initial_model_recover_effect.channel_contribution.isel(hdi=0).sel(channel="x1"),
initial_model_recover_effect.channel_contribution.isel(hdi=1).sel(channel="x1"),
alpha=0.2,
label="95% Credible Interval",
color="blue",
)
ax[0].plot(
date_range, df["x1_adstock_saturated"], label="Real x1 Effect", color="black"
)
# google
ax[1].plot(
date_range,
initial_model_mean_effect.sel(channel="x2"),
label="Mean Recover x2 Effect",
linestyle="--",
color="orange",
)
ax[1].fill_between(
date_range,
initial_model_recover_effect.channel_contribution.isel(hdi=0).sel(channel="x2"),
initial_model_recover_effect.channel_contribution.isel(hdi=1).sel(channel="x2"),
alpha=0.2,
label="95% Credible Interval",
color="orange",
)
ax[1].plot(
date_range, df["x2_adstock_saturated"], label="Real x2 Effect", color="black"
)
# formatting
ax[0].legend()
ax[1].legend()
plt.grid(True)
ax[1].set(xlabel="date")
fig.suptitle("Media Contribution Recovery", fontsize=16)
plt.show()
En este punto, podríamos sentir que hemos perdido todo. Puede parecer que algún tipo de mala suerte está afectando nuestros datos. O, podríamos tener un factor desconocido que está influyendo en nuestros resultados de marketing sin que nos demos cuenta.
Después de una reflexión, podemos considerar que las acciones de nuestros competidores, como sus ofertas y descuentos, podrían estar impactando nuestros resultados de marketing. En otras palabras, cuando nuestros competidores gastan más en marketing, pueden llevarse algunas de nuestras posibles ventas incrementales. Incluso si gastamos la misma cantidad para obtener exposición, las ofertas agresivas de los competidores pueden desviar a los clientes de nosotros mientras compran otros productos en su lugar.
¿Afectará esto a todos nuestros canales de la misma manera? Probablemente no. Nuestros competidores pueden no impactar nuestra notoriedad de marca, pero probablemente influyen en la consideración de nuestros clientes. Podríamos pensar que \(x2\) está más afectado por nuestros competidores que \(x1\). Esto se debe a que \(x2\) está más cerca de la decisión de compra, mientras que \(x1\) se relaciona más con la exposición de la marca.
Esto sugiere que ‘competitor_offers’ podría ser un factor no visible que afecta a \(x2\) pero no a \(x1\). Formalmente, esto significa que podríamos expresar \(x2\) como una función de ‘competitor_offers’ mientras que \(x1\) no incluiría este factor.
Este es un buen enfoque para construir nuestro modelo. Sin embargo, necesitamos añadir ‘competitor_offers’ a nuestro modelo, y actualmente no tenemos datos sobre nuestros competidores, por lo que necesitamos pensar de manera creativa.
Uso de Procesos Gaussianos para Modelar Ofertas de Competidores#
Dado que no observamos directamente las actividades de los competidores, como ofertas o descuentos, podemos tratar ‘competitor_offers’ como una variable oculta. Esto significa que no la mediremos directamente, sino que inferiremos su comportamiento en función de su impacto en otras variables de nuestro modelo. Un Proceso Gaussiano (GP) puede ayudarnos a modelar este factor oculto, ya que es lo suficientemente flexible para reflejar los patrones subyacentes de los datos.
Un Proceso Gaussiano puede expresarse como:
Dónde:
\(f(t)\): La función oculta que representa ‘competitor_offers’ a lo largo del tiempo.
\(m(t)\): La función media, a menudo establecida en 0 a menos que tengamos conocimiento previo de una tendencia específica.
\(k(t, t')\): La función de covarianza (o núcleo), que determina la suavidad y la estructura de \(f(t)\).
PyMC-Marketing incluye Procesos Gaussianos, utilizando una aproximación de Proceso Gaussiano en Espacio de Hilbert (HSGP) para hacer que los cálculos sean más manejables. Esto nos permite aproximar la función oculta \(f(t)\) de manera eficiente.
¿Por qué utilizar un proceso gaussiano?#
Revele la estructura oculta de ‘competitor_offers’ y otras variables no observadas basadas en sus conexiones con \(x2\) y otras variables observadas.
Capture efectos no lineales y variables que pueden surgir de actividades de competidores no visibles.
Ahora veamos cómo podemos implementar esto en PyMC-Marketing.
causal_dag = """
digraph {
x1 -> y;
x2 -> y;
x1 -> x2;
holiday_signal -> y;
holiday_signal -> x1;
holiday_signal -> x2;
competitor_offers -> x2;
competitor_offers -> y;
market_growth -> y;
}
"""
causal_mmm = MMM(
sampler_config=sample_kwargs,
date_column="date_week",
adstock=GeometricAdstock(l_max=24),
saturation=MichaelisMentenSaturation(),
channel_columns=["x1", "x2"],
control_columns=["holiday_signal"],
# Define the outcome node and the causal DAG
outcome_node="y",
dag=causal_dag,
# Time varying intercept to account for the unobserved confounder
time_varying_intercept=True,
)
causal_mmm.model_config["intercept_tvp_config"].ls_mu = 180
causal_mmm.model_config["intercept"] = Prior("Normal", mu=1, sigma=2)
(
causal_mmm.causal_graphical_model.adjustment_set,
causal_mmm.causal_graphical_model.minimal_adjustment_set,
)
(['competitor_offers', 'holiday_signal'], ['holiday_signal', 'x1', 'x2'])
Profundicemos en la salida. Primero, es importante señalar que nuestro conjunto de ajustes indica claramente la necesidad de controlar tanto por competidores como por días festivos para estimar con precisión nuestras variables. Esto implica que las variables en nuestras regresiones deben incluir esos factores junto con nuestras variables de tratamiento. Dado que no contamos con ninguna variable observada en los datos para tener en cuenta a los competidores, establecemos el parámetro time_varying_intercept en True para que nuestro modelo pueda acomodar el confusor no observado.
Ahora, echemos un vistazo a qué tan bien se ajusta nuestro modelo.
X = data.drop("y", axis=1)
y = data["y"]
causal_mmm.fit(X=X, y=y, target_accept=0.95, random_seed=rng)
causal_mmm.sample_posterior_predictive(
X, extend_idata=True, combined=True, random_seed=rng
)
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
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Sampling: [y]
<xarray.Dataset> Size: 17MB
Dimensions: (sample: 2000, date: 1041)
Coordinates:
* date (date) datetime64[ns] 8kB 2022-01-01 2022-01-02 ... 2024-11-06
* sample (sample) object 16kB MultiIndex
* chain (sample) int64 16kB 0 0 0 0 0 0 0 0 0 0 0 ... 3 3 3 3 3 3 3 3 3 3 3
* draw (sample) int64 16kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
Data variables:
y (date, sample) float64 17MB 8.454 8.399 8.547 ... 11.26 11.35 11.28
Attributes:
created_at: 2025-06-16T13:40:28.655973+00:00
arviz_version: 0.21.0
inference_library: pymc
inference_library_version: 5.23.0# show divergencies
causal_mmm.idata["sample_stats"]["diverging"].sum().item()
Excelente en comparación con los modelos anteriores, no tenemos ninguna divergencia, como discutimos antes, eso es una buena señal 🔥.
Observemos el efecto recuperado de nuestras variables de tratamiento (\(x1\) y \(x2\)) y, como de costumbre, el ajuste del modelo que no es tan causal 👀
causal_mmm.plot_posterior_predictive(add_mean=False, original_scale=True)
r2 = az.r2_score(
y_true=df["y"].values,
y_pred=causal_mmm.idata.posterior_predictive.stack(sample=("chain", "draw"))[
"y"
].values.T
* causal_mmm.target_transformer["scaler"].scale_.item(),
).iloc[0]
plt.text(
0.05,
0.95,
f"R2: {r2:.2f}",
transform=plt.gca().transAxes,
fontsize=12,
verticalalignment="top",
bbox=dict(facecolor="white", alpha=0.8, edgecolor="none"),
)
plt.title("Estimated Target Variable Over Time")
plt.xlabel("Days")
plt.ylabel("Target Variable")
plt.legend()
plt.grid(True)
plt.show()
initial_model_recover_effect = (
az.hdi(causal_mmm.fit_result["channel_contribution"], hdi_prob=0.95)
* causal_mmm.target_transformer["scaler"].scale_.item()
)
fig, ax = plt.subplots(2, 1, figsize=(14, 10), sharex=True)
# Social media
ax[0].plot(
date_range,
causal_mmm.fit_result.channel_contribution.sel(channel="x1").mean(
dim=["chain", "draw"]
)
* causal_mmm.target_transformer["scaler"].scale_.item(),
label="Mean Recover x1 Effect",
linestyle="--",
color="blue",
)
ax[0].fill_between(
date_range,
initial_model_recover_effect.channel_contribution.isel(hdi=0).sel(channel="x1"),
initial_model_recover_effect.channel_contribution.isel(hdi=1).sel(channel="x1"),
alpha=0.2,
label="95% Credible Interval",
color="blue",
)
ax[0].plot(
date_range, df["x1_adstock_saturated"], label="Real x1 Effect", color="black"
)
# google
ax[1].plot(
date_range,
causal_mmm.fit_result.channel_contribution.sel(channel="x2").mean(
dim=["chain", "draw"]
)
* causal_mmm.target_transformer["scaler"].scale_.item(),
label="Mean Recover x2 Effect",
linestyle="--",
color="orange",
)
ax[1].fill_between(
date_range,
initial_model_recover_effect.channel_contribution.isel(hdi=0).sel(channel="x2"),
initial_model_recover_effect.channel_contribution.isel(hdi=1).sel(channel="x2"),
alpha=0.2,
label="95% Credible Interval",
color="orange",
)
ax[1].plot(
date_range, df["x2_adstock_saturated"], label="Real x2 Effect", color="black"
)
# formatting
ax[0].legend()
ax[1].legend()
plt.grid(True)
ax[1].set(xlabel="date")
fig.suptitle("Media Contribution Recovery", fontsize=16)
plt.show()
¡Grandes noticias! 🎉 Nuestro modelo está logrando un excelente desempeño al determinar el impacto real de nuestras variables de tratamiento. Hemos conseguido tener en cuenta tanto los factores de confusión no visibles como las cosas que podemos ver, como los días festivos.
Es posible que haya notado que tenemos variables como el crecimiento del mercado que no están en los datos—¿qué está sucediendo con eso? Antes de profundizar, observemos el intercepto variable en el tiempo.
# plot recover intercept
intercept_effect = (
az.hdi(causal_mmm.fit_result["intercept"], hdi_prob=0.95)
* causal_mmm.target_transformer["scaler"].scale_.item()
)
mean_intercept = (
causal_mmm.fit_result.intercept.mean(dim=["chain", "draw"])
* causal_mmm.target_transformer["scaler"].scale_.item()
)
fig, ax = plt.subplots()
(
df.set_index("date_week")["intercept"]
+ df.set_index("date_week")["market_growth"]
- df.set_index("date_week")["competitor_offers"]
).plot(ax=ax, label="f(Intercept) + f(Market Growth) - f(Competitor Offers)")
sns.lineplot(x=date_range, y=mean_intercept, label="Mean Varying Intercept")
ax.fill_between(
date_range,
intercept_effect.intercept.isel(hdi=0),
intercept_effect.intercept.isel(hdi=1),
alpha=0.2,
label="95% Credible Interval",
)
ax.set(title="Recovered Intercept", xlabel="date", ylabel="Intercept")
plt.show()
Este intercepto es realmente bueno para captar ese confusor no observado y el efecto del crecimiento del mercado. Esencialmente, el proceso gaussiano es lo suficientemente inteligente como para adivinar todas esas variables ocultas que afectan nuestro resultado.
Dado que tenemos una comprensión sólida de nuestro contexto, podemos identificar qué variables está determinando nuestro proceso gaussiano. Esta información nos ayuda a ajustar nuestro modelo. Por otro lado, si no entendemos bien nuestro mundo, el proceso gaussiano podría estar inferiendo las variables incorrectas, lo que llevaría a conclusiones erróneas. Por eso es crucial entender realmente nuestro ecosistema causal, el mundo que nos rodea y cómo están estructuradas nuestras relaciones causales. Cada elección y suposición que hacemos con nuestros modelos puede complicar las cosas, impactando la precisión de nuestros resultados.
No deberíamos complicar nuestro modelo sin una razón sólida. ¡Agregar nuevas variables o procesos gaussianos debe ser considerado con cuidado! Es esencial comprender las consecuencias de nuestras elecciones y las suposiciones que adoptamos una vez que las hacemos. Utilizar técnicas de identificación causal puede ayudarnos a determinar qué variables de control incluir, guiándonos en el viaje de descubrimiento causal.
Resumen#
En este cuaderno, exploramos cómo utilizar métodos de descubrimiento causal para identificar confusores no observados y cómo emplear conjuntos de ajuste para incluir variables de control. También analizamos cómo los procesos gaussianos pueden ayudarnos a modelar estos confusores ocultos, particularmente a través del intercepto variable en el tiempo.
Ahora es tu turno: ¡adelante y comienza a construir tus propios modelos!
Referencias#
%load_ext watermark
%watermark -n -u -v -iv -w -p pymc_marketing,pytensor
Last updated: Mon Jun 16 2025
Python implementation: CPython
Python version : 3.10.18
IPython version : 8.37.0
pymc_marketing: 0.14.0
pytensor : 2.31.3
seaborn : 0.13.2
pandas : 2.3.0
pytensor : 2.31.3
preliz : 0.11.0
IPython : 8.37.0
arviz : 0.21.0
matplotlib : 3.10.3
pymc : 5.23.0
pymc_marketing: 0.14.0
numpy : 1.26.4
graphviz : 0.20.3
Watermark: 2.5.0