Modelo Pareto/NBD#

El modelo de Distribución Pareto/Binomial Negativa fue el primer modelo Buy-Till-You-Die (BTYD) para estimar la actividad de clientes no contractuales a lo largo de un periodo de tiempo continuo. Introducido por primera vez por Schmittlein, et. al. en 1987 y desarrollado posteriormente por Bruce Hardie y Peter Fader, se utiliza frecuentemente como un referente en la investigación de CLV debido a su sólido rendimiento y amplia gama de funcionalidades. Para derivaciones detalladas de este modelo, por favor consulte «A Note on Deriving the Pareto/NBD Model and Related Expressions.»

In this notebook we will use Bayesian inference to fit a Pareto/NBD model in PyMC-Marketing. We will also demonstrate the predictive functionality of this model, along with an example for time-invariant covariates.

Configurar el Notebook#

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sb
import xarray as xr
from fastprogress.fastprogress import progress_bar
from pymc_extras.prior import Prior

from pymc_marketing import clv

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

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

Cargar Datos#

En este cuaderno utilizaremos el conjunto de datos de muestra CDNOW, un conjunto de datos de referencia popular en la investigación de modelado de CLV y BTYD. Consulte aquí para obtener más información sobre el conjunto de datos.

url_cdnow = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/cdnow_transactions.csv"

raw_data = pd.read_csv(url_cdnow)

raw_data.info()
raw_data.head()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 6919 entries, 0 to 6918
Data columns (total 5 columns):
 #   Column      Non-Null Count  Dtype  
---  ------      --------------  -----  
 0   _id         6919 non-null   int64  
 1   id          6919 non-null   int64  
 2   date        6919 non-null   int64  
 3   cds_bought  6919 non-null   int64  
 4   spent       6919 non-null   float64
dtypes: float64(1), int64(4)
memory usage: 270.4 KB

	_id	id	date	cds_bought	spent
0	4	1	19970101	2	29.33
1	4	1	19970118	2	29.73
2	4	1	19970802	1	14.96
3	4	1	19971212	2	26.48
4	21	2	19970101	3	63.34

Los únicos requisitos para modelar el comportamiento de gasto con ParetoNBDModel son una columna de identificador de cliente y una columna de fecha y hora para cada compra. El número de CDs comprados y el dinero gastado por transacción también podrían ser covariables útiles, así que los tendremos en cuenta para más adelante.

Es común que las bases de datos de transacciones de clientes también contengan devoluciones, valores de descuento, etc., así que hagamos una rápida verificación de validación:

raw_data.describe()

	_id	id	date	cds_bought	spent
count	6919.000000	6919.000000	6.919000e+03	6919.000000	6919.000000
mean	11682.515826	1175.724816	1.997217e+07	2.381703	35.278500
std	6833.386793	679.426450	3.744182e+03	2.218380	34.074377
min	4.000000	1.000000	1.997010e+07	1.000000	0.000000
25%	5525.000000	570.500000	1.997022e+07	1.000000	14.490000
50%	11749.000000	1193.000000	1.997042e+07	2.000000	25.990000
75%	17717.000000	1766.000000	1.997103e+07	3.000000	42.970000
max	23569.000000	2357.000000	1.998063e+07	40.000000	506.970000

Tenga en cuenta que hubo algunas transacciones con valores de gasto de 0. ¡Quizás se trató de devoluciones o regalos promocionales! Casos como este no son actividades de compra reales y deben ser excluidos del modelado.

raw_data = raw_data[raw_data["spent"] > 0]

Utilice la utilidad rfm_summary para agregar datos para el modelado:

rfm_data = clv.rfm_summary(
    raw_data,
    customer_id_col="id",
    datetime_col="date",
    datetime_format="%Y%m%d",
    time_unit="W",
)

rfm_data.info()
rfm_data.head()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 2349 entries, 0 to 2348
Data columns (total 4 columns):
 #   Column       Non-Null Count  Dtype  
---  ------       --------------  -----  
 0   customer_id  2349 non-null   int64  
 1   frequency    2349 non-null   float64
 2   recency      2349 non-null   float64
 3   T            2349 non-null   float64
dtypes: float64(3), int64(1)
memory usage: 73.5 KB

	customer_id	frequency	recency	T
0	1	3.0	49.0	78.0
1	2	1.0	2.0	78.0
2	3	0.0	0.0	78.0
3	4	0.0	0.0	78.0
4	5	0.0	0.0	78.0

Recuerde las definiciones de agregación de datos del CLV Quickstart:

customer_id es un índice de identificadores únicos para cada cliente.
frequency es el número de compras repetidas que un cliente ha realizado (es decir, el número total de compras menos uno).
recency indica el período de tiempo en el que un cliente realizó su compra más reciente. Si un cliente solo ha realizado 1 compra, la recency es 0.
T es la «edad» de un cliente, o el número de períodos de tiempo desde su primera compra.

Definición del Modelo#

El modelo Pareto/NBD se basa en las siguientes suposiciones para cada cliente:

Los clientes están activos durante un período de tiempo no observado, y luego se vuelven permanentemente inactivos.

Proceso de Compra#

Mientras está activo, el número de transacciones realizadas por un cliente sigue un proceso de Poisson con una tasa de transacción \(\lambda\):

\[P(X(t)=x|\lambda) = \frac{(\lambda t)^{x}e^{-\lambda t}}{x!}, x=0,1,2,...\]

Esto es equivalente a suponer que el tiempo entre transacciones se distribuye exponencialmente con una tasa de transacción \(\lambda\):

\[f(t_{j}-t_{j-1}| \lambda) = \lambda e^{-\lambda (t_{j} - t_{j - 1})}, \quad t_{j} \geq t_{j - 1} \geq 0\]

Donde \(t\) es el período de tiempo de la \(j\)-ésima compra.
La heterogeneidad en \(\lambda\) sigue una distribución Gamma con parámetro de forma \(r\) y parámetro de escala \(\alpha\):

\[g(\lambda|r, \alpha) = \frac{\alpha^{r}\lambda^{r - 1}e^{-\lambda \alpha}}{\Gamma(r)}\]

Proceso de abandono#

La duración de la vida activa no observada de un cliente se distribuye exponencialmente con una tasa de abandono \(\mu\).
La heterogeneidad en \(\mu\) también sigue una distribución Gamma con parámetro de forma \(s\) y parámetro de escala \(\beta\):

\[g(\mu|s, \beta) = \frac{\beta^{s}\mu^{s - 1}e^{-\mu \beta}}{\Gamma(s)}\]
La tasa de transacción \(\lambda\) y el tiempo hasta la deserción \(\mu\) varían de manera independiente para cada cliente.

Si tomamos la expectativa a través de las distribuciones de \(\lambda\) y \(\mu\), podemos derivar una función de verosimilitud para estimar los parámetros \(r\), \(\alpha\), \(s\) y \(\beta\) en la población de clientes. Para más detalles sobre la verosimilitud de ParetoNBD, por favor consulte la documentación.

Ajuste del modelo#

The Bayesian equivalent of Maximum Likelihood Estimation (MLE) is Maximum a Posteriori (MAP), where scalar values for fitted parameters are regularized with priors during estimation.

Un prior «plano» indica que el usuario es agnóstico, sin creencias o suposiciones previas sobre los datos. \(r\), \(\alpha\), \(s\) y \(\beta\) también deben ser valores positivos, así que configuremos nuestro ParetoNBDModel bayesiano con priors HalfFlat:

flat_config = {
    "r": Prior("HalfFlat"),
    "alpha": Prior("HalfFlat"),
    "s": Prior("HalfFlat"),
    "beta": Prior("HalfFlat"),
}

pnbd_pymc = clv.ParetoNBDModel(model_config=flat_config)

Construya el modelo para ver la elección de Priors utilizados para el modelado:

pnbd_pymc.build_model(data=rfm_data)  # optional step
pnbd_pymc

Pareto/NBD
            alpha ~ HalfFlat()
             beta ~ HalfFlat()
                r ~ HalfFlat()
                s ~ HalfFlat()
recency_frequency ~ ParetoNBD(r, alpha, s, beta, <constant>)

Tenga en cuenta que no es necesario construir un modelo antes de la modelización.

Ahora ajustemos nuestro ParetoNBDModel con MAP.

idata_map = pnbd_pymc.fit(method="map")

Para el ajuste de MAP, pymc-marketing utiliza el optimizador L-BGFS-B de scipy.optimize, una alternativa más rápida y estable a Nelder-Mead.

flat_fit = pnbd_pymc.fit_summary()
print(flat_fit)

alpha    15.644
beta     13.796
r         0.614
s         0.447
Name: value, dtype: float64

The CDNOW sample we’re working with is quite small and comprises only 10% of the total CDNOW dataset, so it’s quite likely these estimates are overfitting if we attempt to run predictions on the full dataset.

Con distribuciones a priori, podemos informar el ajuste del modelo con nuestro propio conocimiento subjetivo del dominio e incluso mejorar la velocidad de los ajustes del modelo. La configuración a priori predeterminada para ParetoNBDModel funciona bien para una variedad de casos de uso:

pnbd_map = clv.ParetoNBDModel()
pnbd_map.build_model(data=rfm_data)  # required for prior predictive checks
pnbd_map

Pareto/NBD
            alpha ~ Weibull(2, 10)
             beta ~ Weibull(2, 10)
                r ~ Weibull(2, 1)
                s ~ Weibull(2, 1)
recency_frequency ~ ParetoNBD(r, alpha, s, beta, <constant>)

Comprobaciones Predictivas Anteriores y Posteriores#

Los PPCs nos permiten verificar la eficacia de nuestros priors y el rendimiento de los posteriors ajustados. Los PPCs generalmente no son una opción con modelos ajustados por MAP, pero aquí en realidad estamos muestreando de las distribuciones Gamma latentes \(\lambda\) y \(\mu\), por lo que los PPCs son posibles para ParetoNBDModel independientemente del método de ajuste.

Veamos cómo se desempeña el modelo en un chequeo predictivo previo, donde muestreamos de las distribuciones a priori predeterminadas antes de ajustar el modelo:

with pnbd_map.model:
    prior_idata = pm.sample_prior_predictive(random_seed=45, draws=1)

obs_freq = prior_idata.observed_data["recency_frequency"].sel(obs_var="frequency")
ppc_freq = prior_idata.prior_predictive["recency_frequency"].sel(obs_var="frequency")[
    0
][0]

# PPC histogram plot
clv.plot_expected_purchases_ppc(pnbd_map, ppc="prior");

Sampling: [alpha, beta, r, recency_frequency, s]
Sampling: [alpha, beta, r, recency_frequency, s]

../../_images/48299f4664b14b07d46fd27fb00c33f360f2dfb6564e2c77ffbccd26e6b1ca2a.png

Aquí el ParetoNBDModel está simulando las compras de los clientes a partir de las distribuciones previas para compararlas con los datos observados. Las distribuciones previas por defecto parecen razonables para los clientes no recurrentes, pero no modelan bien a los clientes recurrentes.

Ajustemos nuestro modelo y realicemos un chequeo predictivo posterior para comparación:

pnbd_map.fit()
map_fit = pnbd_map.fit_summary()  # save for plotting later

obs_freq = pnbd_map.idata.observed_data["recency_frequency"].sel(obs_var="frequency")
ppc_freq = pnbd_map.distribution_new_customer_recency_frequency(
    rfm_data,
    random_seed=42,
).sel(chain=0, draw=0, obs_var="frequency")

# PPC histogram plot
clv.plot_expected_purchases_ppc(pnbd_map, ppc="posterior");

Sampling: [recency_frequency]

Sampling: [recency_frequency]

../../_images/518e32d0c5d414be6fd64a063b7bc7ee6db9a1c6407ae89c14515ef7a36d9d75.png

¡Nuestro modelo ajustado es capaz de simular de manera confiable el comportamiento del cliente!

Inferencia Bayesiana Completa#

MAP solo ajusta valores escalares para \(r\), \(\alpha\), \(s\) y \(\beta\), pero con un muestreo bayesiano completo podemos inferir las distribuciones de probabilidad posterior para estos parámetros, ilustrando la incertidumbre en nuestras estimaciones así como permitiendo intervalos de predicción.

NUTS es el muestreador predeterminado en pymc-marketing, que muestrea de la posterior explorando los gradientes del espacio de probabilidad. Sin embargo, el muestreo NUTS con ParetoNBDModel puede ser bastante lento debido a la complejidad de la expresión de verosimilitud. De hecho, la complejidad matemática de este modelo fue lo que motivó el desarrollo del BetaGeoModel en 2005. El modelo BG/NBD hace algunas suposiciones simplificadoras y sacrifica funcionalidad en la estimación de abandono de clientes para mejorar el rendimiento computacional.

Para ahorrar tiempo y costos computacionales, se recomienda utilizar el muestreador DEMetropolisZ sin gradiente. Esto a menudo requiere más muestras durante el ajuste, por lo que si se encuentran advertencias de estadística rhat, aumente el tamaño de los parámetros tune y draw hasta que la advertencia ya no aparezca.

pnbd_full.fit_summary()

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
alpha	15.769	1.114	13.717	17.876	0.038	0.027	882.0	1152.0	1.01
beta	12.565	3.630	6.598	19.738	0.124	0.088	821.0	1072.0	1.00
r	0.627	0.049	0.540	0.721	0.002	0.001	845.0	1100.0	1.00
s	0.428	0.060	0.321	0.541	0.002	0.001	838.0	1212.0	1.00

Utilice thin_fit_result si la gran cantidad de extracciones está causando problemas computacionales. Mantener cada segunda muestra reducirá el número de extracciones a la mitad:

axes = az.plot_trace(
    data=pnbd_full.idata,
    compact=True,
    kind="rank_bars",
    backend_kwargs={"figsize": (12, 7), "layout": "constrained"},
)
plt.gcf().suptitle("Pareto/NBD Model Trace", fontsize=18, fontweight="bold");

../../_images/e872c0de3196cdd02fabba39bfe6ac4eeade917cf14009fcded7fa76133ccfdc.png

Veamos cómo se comparan los posteriores de DEMZ con las estimaciones MAP:

_, axes = plt.subplots(
    nrows=2, ncols=2, figsize=(12, 7), sharex=False, sharey=False, layout="constrained"
)

axes = axes.flatten()

for i, var_name in enumerate(["r", "alpha", "s", "beta"]):
    ax = axes[i]
    az.plot_posterior(
        pnbd_full.idata.posterior[var_name].values.flatten(),
        color="C0",
        point_estimate="mean",
        ax=ax,
        label="DEMZ",
    )
    ax.axvline(x=map_fit[var_name], color="C1", linestyle="--", label="MAP")
    ax.legend(loc="upper right")
    ax.set_title(var_name)

plt.gcf().suptitle("Pareto/NBD Model Parameters", fontsize=18, fontweight="bold");

../../_images/b32401e1649f83edae5e66802620f88ceb7b5f992b846dbbec4f087eb4b63763.png

Después de ajustar, los modelos se pueden guardar con pnbd_pymc.save("pnbd.nc") y cargar con pnbd_pymc.load("pnbd.nc").

Métodos Predictivos#

El modelo Pareto/NBD admite una variedad de métodos predictivos:

compras_esperadas
probabilidad_esperada_viva
compras_esperadas_nuevo_cliente
probabilidad_de_compra_esperada

Tomemos una pequeña muestra de usuarios:

example_customer_ids = [1, 5, 10, 18, 46, 1413]

rfm_sample = rfm_data.query("customer_id.isin(@example_customer_ids)")

rfm_sample.sort_values(by="frequency")

	customer_id	frequency	recency	T
4	5	0.0	0.0	78.0
9	10	0.0	0.0	78.0
17	18	1.0	5.0	78.0
0	1	3.0	49.0	78.0
1405	1413	19.0	54.0	71.0
45	46	21.0	73.0	78.0

Observe que los clientes 5 y 10 son compradores no recurrentes, mientras que 1413 y 46 son compradores frecuentes.

Número Esperado de Compras#

Vamos a trazar el número esperado de compras de cada cliente durante los próximos \(90\) períodos de tiempo:

time_periods = 90

expected_purchases_over_time = xr.concat(
    objs=[
        pnbd_full.expected_purchases(
            data=rfm_sample,
            future_t=t,
        )
        for t in progress_bar(range(time_periods))
    ],
    dim="t",
).transpose(..., "t")

100.00% [90/90 00:24<00:00]

_, axes = plt.subplots(
    nrows=len(example_customer_ids),
    ncols=1,
    figsize=(12, 14),
    sharex=True,
    sharey=True,
    layout="constrained",
)

axes = axes.flatten()

for i, customer_id in enumerate(example_customer_ids):
    ax = axes[i]
    customer_expected_purchases = expected_purchases_over_time.sel(
        customer_id=customer_id
    )
    az.plot_hdi(
        range(time_periods),
        customer_expected_purchases,
        hdi_prob=0.94,
        color="C0",
        fill_kwargs={"alpha": 0.3, "label": "$94 \\%$ HDI"},
        ax=ax,
    )
    az.plot_hdi(
        range(time_periods),
        customer_expected_purchases,
        hdi_prob=0.5,
        color="C0",
        fill_kwargs={"alpha": 0.5, "label": "$50 \\%$ HDI"},
        ax=ax,
    )
    ax.plot(
        range(time_periods),
        customer_expected_purchases.mean(dim=("chain", "draw")),
        color="C0",
        label="posterior mean",
    )
    ax.legend(loc="upper left")
    ax.set(title=f"Customer {customer_id}", xlabel="t", ylabel="# of purchases")

axes[-1].set(xlabel="t")
plt.gcf().suptitle("Expected Number of Purchases", fontsize=18, fontweight="bold");

../../_images/c13325b12af2e0e0e877ad244f404564384f0d45bb2b4d928fe59c4f3007242a.png

Tenga en cuenta que los intervalos de predicción del IDH solo están disponibles si el modelo se ajusta con posteriors completos.

Observe el gran número de compras que se esperan de los compradores frecuentes (Clientes 46 y 1413), mientras que se espera poca o ninguna actividad futura de los demás clientes.

Probabilidad Esperada de Estar Vivo#

Verifiquemos la probabilidad de que nuestros clientes sigan vivos y comparemos con los resultados del MAP.

_, axes = plt.subplots(
    nrows=3,
    ncols=2,
    figsize=(12, 14),
    layout="constrained",
)

axes = axes.flatten()

for i, customer_id in enumerate(example_customer_ids):
    ax = axes[i]
    demz_alive = pnbd_full.expected_probability_alive(
        rfm_sample,
        future_t=0,
    ).sel(customer_id=customer_id)
    map_alive = pnbd_map.expected_probability_alive(rfm_sample, future_t=0).sel(
        customer_id=customer_id
    )

    az.plot_density(demz_alive, hdi_prob=1, colors="C1", shade=0.3, bw=0.005, ax=ax)
    ax.axvline(x=map_alive, color="C3", linestyle="--", label="MAP")
    ax.legend(loc="upper right")
    ax.set(
        title=f"Customer {customer_id}",
        xlabel="Probability Alive",
        ylabel="Posterior Density",
    )

plt.gcf().suptitle("Expected Probability Alive", fontsize=18, fontweight="bold");

../../_images/5536872a91d4792200a299c00ae0a3a6d4ad1dd8e2b890c17fc7641f5a470e4c.png

El cliente 1413 tiene una probabilidad de permanencia bastante baja a pesar de ser un comprador frecuente. Este sería un buen ejemplo de un cliente al que dirigirse con una oferta especial para la retención.

Estas probabilidades se estiman en el período de tiempo 0, pero también podemos estimar las probabilidades de que los clientes sigan activos en el futuro. Calculemos las densidades posteriores 90 períodos de tiempo a partir de ahora y comparemos con el MAP en el período de tiempo 0:

_, axes = plt.subplots(
    nrows=3,
    ncols=2,
    figsize=(12, 14),
    layout="constrained",
)

axes = axes.flatten()

for i, customer_id in enumerate(example_customer_ids):
    ax = axes[i]
    demz_alive = pnbd_full.expected_probability_alive(
        rfm_sample,
        future_t=90,
    ).sel(customer_id=customer_id)
    map_alive = pnbd_map.expected_probability_alive(rfm_sample, future_t=0).sel(
        customer_id=customer_id
    )

    az.plot_density(demz_alive, hdi_prob=1, colors="C1", shade=0.3, bw=0.005, ax=ax)
    ax.axvline(x=map_alive, color="C3", linestyle="--", label="MAP at Time 0")
    ax.legend(loc="upper left")
    ax.set(
        title=f"Customer {customer_id}",
        xlabel="Probability Alive",
        ylabel="Posterior Density at Time 90",
    )

plt.gcf().suptitle("Expected Probability Alive Delta", fontsize=18, fontweight="bold");

../../_images/adb50af6b15d510c514776d1898260d5c103f73748b2226312ddfb01603c914e.png

Preste atención a los ejes x para cada cliente - Las probabilidades apenas cambiaron para los clientes no recurrentes, pero hay un delta significativo para los compradores frecuentes.

Una buena regla general es que una probabilidad de actividad de .25-.30 suele indicar un cliente en riesgo o inactivo. Las proyecciones futuras pueden proporcionar información adicional sobre el riesgo de deserción de clientes.

Probabilidad de \(n\) Compras a lo largo del Tiempo \(t\)#

El cliente 46 es nuestro mejor cliente en este pequeño conjunto de muestras y se espera que realice al menos \(15\) compras en los próximos \(90\) períodos de tiempo. ¿Cuál es la probabilidad de que se realicen tantas compras y cómo cambiará con el tiempo?

Vamos a trazar un mapa de calor para pintar el panorama completo:

# create arrays of parameter combinations
n_purchases = np.repeat([0, 3, 6, 9, 12, 15], 6)
time_periods = np.tile([15, 30, 45, 60, 75, 90], 6)

expected_purchase_prob_heatmap = xr.concat(
    objs=[
        pnbd_map.expected_purchase_probability(
            rfm_sample,
            n_purchases=params[0],
            future_t=params[1],
        ).sel(customer_id=46)
        for params in zip(n_purchases, time_periods, strict=False)
    ],
    dim="customer_id",
).transpose(..., "customer_id")

heatmap_reshape = expected_purchase_prob_heatmap.values.reshape(6, 6)

sb.heatmap(heatmap_reshape, annot=True)

plt.xlabel("Time Periods")
plt.xticks(np.arange(6) + 0.5, [15, 30, 45, 60, 75, 90], rotation=0)
plt.ylabel("Number of Purchases")
plt.yticks(np.arange(6) + 0.5, [0, 3, 6, 9, 12, 15], rotation=0)
plt.gcf().suptitle(
    "Expected Purchase Probabilities for Customer 46", fontsize=18, fontweight="bold"
);

../../_images/140e942462628b8b6bc825ebb02f96e5dbea1d9cb7a7dc313d8163fb8b23f786.png

Este mapa de calor destaca cómo se espera que el Cliente 46 realice al menos 15 compras hasta el período de tiempo 90, pero las probabilidades de que se realicen 15 compras antes del período de tiempo 75 o incluso del período de tiempo 60 son ligeramente más altas. También tenga en cuenta que estas probabilidades asumen expectativas exactas (es decir, hay un 6.2% de probabilidad de que la 15ª compra se realice precisamente durante el período de tiempo 60).

Número Esperado de Compras para Nuevos Clientes#

Hasta ahora, solo hemos estado realizando predicciones para clientes existentes, pero también podemos estimar el número esperado de transacciones a lo largo del tiempo para un nuevo cliente:

expected_purchases_over_time_new_customer = xr.concat(
    objs=[
        pnbd_full.expected_purchases_new_customer(
            data=rfm_sample,
            t=t,
        ).sel(customer_id=1)  # customer_id is arbitrary here
        for t in range(90)
    ],
    dim="t",
).transpose(..., "t")


# plot results
ax = plt.axes()

az.plot_hdi(
    range(90),
    expected_purchases_over_time_new_customer,
    hdi_prob=0.94,
    color="C2",
    fill_kwargs={"alpha": 0.3, "label": "$94 \\%$ HDI"},
    ax=ax,
)
az.plot_hdi(
    range(90),
    expected_purchases_over_time_new_customer,
    hdi_prob=0.5,
    color="C2",
    fill_kwargs={"alpha": 0.5, "label": "$50 \\%$ HDI"},
    ax=ax,
)
ax.plot(
    range(90),
    expected_purchases_over_time_new_customer.mean(dim=("chain", "draw")),
    color="C2",
    label="posterior mean",
)
ax.legend(loc="upper left")
ax.set(
    title="Expected Number of Purchases by New Customers",
    ylabel="# of purchases",
    xlabel="t",
);

../../_images/c1a0cc49f06ab649896e930314137066c41026546f7586995e44bd83bd12bbc2.png

Veamos cómo cambian estas estimaciones cuando añadimos covariables al modelo.

Covariables Invariantes en el Tiempo#

Recuerde que \(\alpha\) y \(\beta\) representan los parámetros de escala para las distribuciones de la tasa de compra y de abandono, respectivamente. Para modelar covariables invariables en el tiempo, simplemente modificamos estos parámetros de la siguiente manera:

\[\alpha = \alpha_0e^{-\gamma_1'z_1}\]

\[\beta = \beta_0e^{-\gamma_2'z_2}\]

Donde \(\gamma_1\) y \(\gamma_2\) son coeficientes que capturan el impacto de las covariables, y \(z_1\) y \(z_2\) son los arreglos de covariables para cada cliente.

Echemos un vistazo a las covariables disponibles en los datos en bruto:

# aggregate raw data by customer id
covar_df = raw_data[["id", "cds_bought", "spent"]].groupby("id").mean().reset_index()

# plot covariate histograms
_, axes = plt.subplots(
    nrows=1,
    ncols=2,
    figsize=(10, 5),
    layout="constrained",
)

axes = axes.flatten()

covars = ["cds_bought", "spent"]
colors = ["C0", "C2"]

for ax in zip(axes, covars, colors, strict=False):
    ax[0].hist(
        x=covar_df[ax[1]],
        bins=5,
        color=ax[2],
    )
    ax[0].set(title=f"{ax[1]}", xlabel="value", ylabel="count")

../../_images/c476cf90709c295a505f27997b10211f07562e90396460f9ba526d2ea444fac3.png

Las distribuciones de una cola con valores grandes como este complicarán el ajuste del modelo, así que registremos y estandaricemos nuestras covariables:

for covar in ["cds_bought", "spent"]:
    covar_df[f"log_std_{covar}"] = np.log(covar_df[covar]).copy()
    covar_df[f"log_std_{covar}"] -= np.nanmean(covar_df[f"log_std_{covar}"])
    covar_df[f"log_std_{covar}"] /= np.nanstd(covar_df[f"log_std_{covar}"])
    covar_df[f"log_std_{covar}"] = covar_df[f"log_std_{covar}"].fillna(0)

_, axes = plt.subplots(
    nrows=1,
    ncols=2,
    figsize=(10, 5),
    layout="constrained",
)

axes = axes.flatten()

for ax in zip(axes, covars, colors, strict=False):
    ax[0].hist(
        x=covar_df[f"log_std_{ax[1]}"],
        bins=5,
        color=ax[2],
    )
    ax[0].set(title=f"transformed {ax[1]}", xlabel="value", ylabel="count")

../../_images/d74af3034b2c70407f6a4c218d648a42776d070ea90674bf94d6301ffce5abca.png

Para parametrizar ParetoNBDModel con covariables, une las covariables a los datos RFM existentes y especifica los nombres de las columnas en la model_config. Las covariables para las tasas de compra y abandono pueden especificarse por separado, lo que permite experimentar con varias combinaciones para encontrar lo que mejor funciona:

rfm_covar = rfm_data.merge(covar_df, left_on="customer_id", right_on="id", how="inner")

pnbd_covar = clv.ParetoNBDModel(
    model_config={
        "purchase_covariate_cols": ["log_std_cds_bought", "log_std_spent"],
        "dropout_covariate_cols": ["log_std_cds_bought", "log_std_spent"],
    },
)

pnbd_covar.build_model(rfm_covar)
pnbd_covar

Pareto/NBD
purchase_coefficient ~ Normal(0, 1)
         alpha_scale ~ Weibull(2, 10)
 dropout_coefficient ~ Normal(0, 1)
          beta_scale ~ Weibull(2, 10)
                   r ~ Weibull(2, 1)
                   s ~ Weibull(2, 1)
               alpha ~ Deterministic(f(alpha_scale, purchase_coefficient))
                beta ~ Deterministic(f(beta_scale, dropout_coefficient))
   recency_frequency ~ ParetoNBD(r, alpha, s, beta, <constant>)

Los parámetros adicionales se crean automáticamente cuando se añaden covariables.

Hagamos un ajuste rápido del MAP y verifiquemos los resultados:

pnbd_covar.fit(method="map")

print("Fitted Model Parameters")
summary = pnbd_covar.fit_summary(
    var_names=[
        "r",
        "alpha_scale",
        "s",
        "beta_scale",
        "purchase_coefficient",
        "dropout_coefficient",
    ]
)
print(summary)

Fitted Model Parameters
r                                            0.668
alpha_scale                                 17.245
s                                            0.469
beta_scale                                  13.446
purchase_coefficient[log_std_cds_bought]     0.062
purchase_coefficient[log_std_spent]          0.189
dropout_coefficient[log_std_cds_bought]     -0.954
dropout_coefficient[log_std_spent]           0.337
Name: value, dtype: float64

Los parámetros purchase_coefficient y dropout_coefficient indican los respectivos impactos de cada covariable; un signo negativo puede interpretarse como «menos probable que realice una compra o que abandone».

Podemos ver que el número promedio de CDs por compra solo tiene un pequeño impacto en el tiempo entre compras, pero un impacto bastante grande en la tasa de abandono. Los clientes que compran múltiples CDs con frecuencia son los menos propensos a abandonar.

El gasto promedio por compra es significativo tanto para la compra como para el tiempo hasta la deserción, pero tenga en cuenta que si utiliza el modelo Gamma-Gamma para estimar el valor de vida del cliente según el Inicio Rápido, entonces el gasto promedio no puede ser utilizado como una covariable porque una suposición importante del modelo Gamma-Gamma es que el gasto y la frecuencia son no correlacionados.

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

Last updated: Sat Nov 23 2024

Python implementation: CPython
Python version       : 3.10.14
IPython version      : 8.22.2

pymc    : 5.15.1
pytensor: 2.22.1

matplotlib    : 3.8.4
arviz         : 0.18.0
pymc          : 5.15.1
seaborn       : 0.13.2
numpy         : 1.26.4
xarray        : 2024.10.0
pymc_marketing: 0.10.0
pytensor      : 2.22.1
pandas        : 2.2.2

Watermark: 2.4.3