Uplift modeling: ¿a quién contactar?

"""#3 Uplift modeling: ¿a quién vale la pena contactar?

Datos: Kevin Hillstrom MineThatData E-Mail Analytics (64,000 clientes; tercio sin correo = control).
Se estima el efecto incremental del correo sobre la visita al sitio con un S-learner / T-learner y se evalúa
con la curva Qini, comparándolo con el enfoque tradicional de "propensión".
"""
import numpy as np
import pandas as pd
from lightgbm import LGBMClassifier, LGBMRegressor
from sklearn.model_selection import StratifiedKFold, train_test_split

from common import SEED, bars, download, export, line, table
from validacion import model_card

URL = "https://hillstorm1.s3.us-east-2.amazonaws.com/hillstorm_no_indices.csv.gz"


def load() -> pd.DataFrame:
    df = pd.read_csv(download(URL, "hillstrom.csv.gz"))
    df["tratamiento"] = (df["segment"] != "No E-Mail").astype(int)
    for c in ["history_segment", "zip_code", "channel"]:
        df[c] = df[c].astype("category")
    return df


def qini_curve(y, t, uplift, points=40):
    """Qini: incremento acumulado de conversiones atribuibles al tratamiento al contactar top-k."""
    o = np.argsort(-uplift)
    y, t = y[o], t[o]
    nt, nc = np.cumsum(t), np.cumsum(1 - t)
    yt, yc = np.cumsum(y * t), np.cumsum(y * (1 - t))
    q = yt - yc * np.divide(nt, nc, out=np.zeros_like(nt, dtype=float), where=nc > 0)
    idx = np.linspace(0, len(y) - 1, points).astype(int)
    pts = [[0.0, 0.0]] + [[round((i + 1) / len(y) * 100, 1), round(float(q[i]), 1)] for i in idx]
    return pts, float(q[-1])


def qini_coef(curve, total):
    """Área entre la curva y la diagonal aleatoria, normalizada por el número de puntos."""
    xs = np.array([c[0] for c in curve]) / 100
    ys = np.array([c[1] for c in curve])
    return float(np.trapezoid(ys - xs * total, xs))


FEATS = ["recency", "history_segment", "history", "mens", "womens", "zip_code", "newbie", "channel"]
# Árboles poco profundos y hojas grandes: el efecto incremental es una señal débil y ruidosa
PARAMS = dict(n_estimators=200, learning_rate=0.03, num_leaves=4, min_child_samples=500,
              random_state=SEED, verbose=-1)


def learners(X_tr, y_tr, t_tr, X_te) -> dict[str, np.ndarray]:
    """Scores de uplift de varios meta-learners (y la propensión como referencia)."""
    m_t = LGBMClassifier(**PARAMS).fit(X_tr[t_tr == 1], y_tr[t_tr == 1])
    m_c = LGBMClassifier(**PARAMS).fit(X_tr[t_tr == 0], y_tr[t_tr == 0])
    m_s = LGBMClassifier(**PARAMS).fit(X_tr.assign(t=t_tr), y_tr)
    # Transformed outcome: Z = Y·(T/e − (1−T)/(1−e)) tiene esperanza igual al uplift
    e = t_tr.mean()
    z = y_tr * (t_tr / e - (1 - t_tr) / (1 - e))
    m_to = LGBMRegressor(**PARAMS).fit(X_tr, z)
    return {
        "S-learner": m_s.predict_proba(X_te.assign(t=1))[:, 1] - m_s.predict_proba(X_te.assign(t=0))[:, 1],
        "T-learner": m_t.predict_proba(X_te)[:, 1] - m_c.predict_proba(X_te)[:, 1],
        "Transformed outcome": m_to.predict(X_te),
        "Propensión": m_t.predict_proba(X_te)[:, 1],
    }


def smd(df: pd.DataFrame, t: np.ndarray) -> pd.Series:
    """Diferencia de medias estandarizada tratados vs. control (balance de la aleatorización)."""
    X = pd.get_dummies(df, dtype=float)
    a, b = X[t == 1], X[t == 0]
    return ((a.mean() - b.mean()) / np.sqrt((a.var() + b.var()) / 2)).fillna(0)


def dif_ic(y, t):
    """Diferencia de proporciones tratados − control con IC 95 % (aproximación normal)."""
    p1, p0 = y[t == 1].mean(), y[t == 0].mean()
    se = np.sqrt(p1 * (1 - p1) / (t == 1).sum() + p0 * (1 - p0) / (t == 0).sum())
    return p1 - p0, p1 - p0 - 1.96 * se, p1 - p0 + 1.96 * se


def run():
    df = load()
    y, t = df["visit"].values, df["tratamiento"].values

    # 1) Validación del experimento: ¿la aleatorización balanceó las covariables?
    balance = smd(df[FEATS], t)
    max_smd = float(balance.abs().max())

    # 2) Efectos promedio por brazo con IC 95 % (visita, conversión y gasto)
    brazos = {}
    for brazo in ["Mens E-Mail", "Womens E-Mail"]:
        m = df["segment"].isin([brazo, "No E-Mail"]).values
        tt = (df["segment"].values[m] == brazo).astype(int)
        brazos[brazo] = {k: dif_ic(df[k].values[m], tt) for k in ["visit", "conversion"]}
        g = df["spend"].values[m]
        brazos[brazo]["spend"] = (g[tt == 1].mean() - g[tt == 0].mean(),)

    # 3) Validación cruzada 5-fold de los learners (Qini en cada fold de validación)
    X = df[FEATS]
    cv = []
    for k, (tr, va) in enumerate(StratifiedKFold(5, shuffle=True, random_state=SEED).split(X, t * 2 + y)):
        sc = learners(X.iloc[tr], y[tr], t[tr], X.iloc[va])
        for nombre, u in sc.items():
            curva, total = qini_curve(y[va], t[va], u)
            cv.append({"learner": nombre, "fold": k + 1, "qini": qini_coef(curva, total)})
    cv = pd.DataFrame(cv)
    cv_res = cv.groupby("learner", sort=False).qini.agg(["mean", "std"])

    # 4) Ajuste final con split 60/40 y curvas en prueba
    X_tr, X_te, y_tr, y_te, t_tr, t_te = train_test_split(X, y, t, test_size=0.4, stratify=t * 2 + y,
                                                          random_state=SEED)
    sc = learners(X_tr, y_tr, t_tr, X_te)
    uplift = sc["S-learner"]
    rng = np.random.default_rng(SEED)
    curvas = {k: qini_curve(y_te, t_te, v)[0] for k, v in sc.items() if k != "Transformed outcome"}
    curvas["Al azar"], _ = qini_curve(y_te, t_te, rng.random(len(y_te)))
    ate, ate_lo, ate_hi = dif_ic(y_te, t_te)

    d = pd.DataFrame({"y": y_te, "t": t_te, "u": uplift})
    d["decil"] = pd.qcut(-d["u"].rank(method="first"), 10, labels=[f"D{i}" for i in range(1, 11)])
    dec = d.groupby("decil", observed=True).apply(
        lambda g: pd.Series({"pred": g.u.mean() * 100, "obs": (g.y[g.t == 1].mean() - g.y[g.t == 0].mean()) * 100}),
        include_groups=False).reset_index()
    top30 = d.nlargest(int(len(d) * 0.3), "u")
    up_top30, lo30, hi30 = dif_ic(top30.y.values, top30.t.values)
    # Correlación de rangos entre uplift predicho y observado por decil (¿el orden se sostiene?)
    rho = float(dec[["pred", "obs"]].corr(method="spearman").iloc[0, 1])

    # 5) Explicabilidad honesta: efecto observado (con IC) por segmento, sin modelo de por medio
    seg_rows = []
    for col, nombre in [("history_segment", "Historial"), ("channel", "Canal"), ("zip_code", "Zona"),
                        ("newbie", "Cliente nuevo")]:
        for v, g in df.groupby(col, observed=True):
            e, lo, hi = dif_ic(g.visit.values, g.tratamiento.values)
            seg_rows.append([f"{nombre}: {v}", len(g), round(e * 100, 2), round(lo * 100, 2), round(hi * 100, 2)])
    seg = pd.DataFrame(seg_rows, columns=["Segmento", "Clientes", "Efecto (pp)", "IC 2.5 %", "IC 97.5 %"])
    seg = seg.sort_values("Efecto (pp)", ascending=False)

    # 6) Cambio de enfoque: política de 3 brazos (¿correo de hombres, de mujeres o ninguno?)
    #    S-learner multi-tratamiento y valor de la política estimado con IPW en el set de prueba.
    arm_tr = df.loc[X_tr.index, "segment"].values
    arm_te = df.loc[X_te.index, "segment"].values
    ms = LGBMClassifier(**PARAMS).fit(X_tr.assign(arm=pd.Categorical(arm_tr)), y_tr)
    cats = pd.Categorical(arm_tr).categories
    pred_arm = np.column_stack([ms.predict_proba(X_te.assign(arm=pd.Categorical([a] * len(X_te), categories=cats)))[:, 1]
                                for a in cats])
    politica = cats[np.argmax(pred_arm, axis=1)]
    prob = pd.Series(arm_tr).value_counts(normalize=True)
    def valor(asignacion):  # IPW: promedio de y en quienes recibieron lo que la política indica
        m = arm_te == np.asarray(asignacion)
        return float((y_te * m / prob.reindex(arm_te).values).mean() * 100)
    politicas = {"Nadie": valor(["No E-Mail"] * len(y_te)), "Todos: correo hombres": valor(["Mens E-Mail"] * len(y_te)),
                 "Todos: correo mujeres": valor(["Womens E-Mail"] * len(y_te)), "Política del modelo": valor(politica)}
    mix = pd.Series(politica).value_counts(normalize=True)

    cv_tbl = cv_res.reset_index()
    cv_tbl["Qini (media ± desv.)"] = cv_tbl.apply(lambda r: f"{r['mean']:.1f} ± {r['std']:.1f}", axis=1)
    cv_tbl = cv_tbl[["learner", "Qini (media ± desv.)"]].rename(columns={"learner": "Learner"})
    brazo_tbl = pd.DataFrame([[b, f"{v['visit'][0] * 100:+.2f} pp [{v['visit'][1] * 100:.2f}, {v['visit'][2] * 100:.2f}]",
                               f"{v['conversion'][0] * 100:+.2f} pp [{v['conversion'][1] * 100:.2f}, {v['conversion'][2] * 100:.2f}]",
                               f"${v['spend'][0]:+.2f}"] for b, v in brazos.items()],
                             columns=["Brazo vs. sin correo", "Visita (IC 95 %)", "Conversión (IC 95 %)", "Gasto por cliente"])
    dec_tbl = dec.round(2)
    dec_tbl.columns = ["Decil", "Uplift predicho (pp)", "Uplift observado (pp)"]
    mejor = cv_res["mean"].idxmax()

    result = {
        "title": "Uplift modeling: ¿a quién contactar?",
        "summary": [
            f"Experimento válido: el mayor desbalance entre grupos es |SMD| = {max_smd:.3f} (< 0.1). Efecto promedio "
            f"del correo sobre la visita: +{ate * 100:.1f} pp (IC 95 % {ate_lo * 100:.1f}–{ate_hi * 100:.1f}).",
            f"El correo para hombres rinde más: +{brazos['Mens E-Mail']['visit'][0] * 100:.1f} pp de visita y "
            f"+${brazos['Mens E-Mail']['spend'][0]:.2f} de gasto por cliente vs. +{brazos['Womens E-Mail']['visit'][0] * 100:.1f} pp "
            f"y +${brazos['Womens E-Mail']['spend'][0]:.2f} del de mujeres.",
            f"CV 5-fold del Qini: {', '.join(f'{k} {v:.0f}' for k, v in cv_res['mean'].items())}. {mejor} es el mejor "
            f"pero con mucha varianza entre folds (±{cv_res.loc[mejor, 'std']:.0f}): el uplift individual es una señal débil; "
            f"en prueba, el top 30 % tiene efecto +{up_top30 * 100:.1f} pp (IC {lo30 * 100:.1f}–{hi30 * 100:.1f}); "
            f"Spearman predicho vs. observado por decil = {rho:.2f}.",
            f"Cambio de enfoque, de \"¿contacto o no?\" a \"¿qué correo a quién?\": la política personalizada logra "
            f"{politicas['Política del modelo']:.1f} % de visita (IPW) vs. {politicas['Todos: correo hombres']:.1f} % enviando "
            f"a todos el de hombres y {politicas['Nadie']:.1f} % sin correo. "
            + ("Conclusión: personalizar no supera a la mejor regla simple; recomendación = correo de hombres a todos."
               if politicas['Política del modelo'] - politicas['Todos: correo hombres'] < 0.3
               else "La personalización agrega valor sobre la mejor regla simple."),
        ],
        "metrics": [
            {"label": "Efecto promedio (IC 95 %)", "value": f"+{ate * 100:.1f} pp [{ate_lo * 100:.1f}–{ate_hi * 100:.1f}]"},
            {"label": "Efecto top 30 %", "value": f"+{up_top30 * 100:.1f} pp"},
            {"label": f"Qini CV {mejor}", "value": f"{cv_res.loc[mejor, 'mean']:.0f} ± {cv_res.loc[mejor, 'std']:.0f}"},
            {"label": "|SMD| máximo", "value": f"{max_smd:.3f}"},
        ],
        "charts": [
            line("qini", "Curva Qini en prueba: visitas incrementales acumuladas", curvas,
                 x_label="% de clientes contactados (ordenados por score)", y_label="Visitas incrementales"),
            bars("cvq", "Qini por learner (CV 5-fold, media)", cv_res.index, cv_res["mean"].round(1).values,
                 fmt="num", subtitle="Desviación por fold en la tabla de validación"),
            bars("dec", "Uplift observado por decil de uplift predicho", dec["decil"], dec["obs"].round(2),
                 fmt="pp", subtitle="Puntos porcentuales de visita: tratados − control"),
            bars("segm", "Efecto observado por segmento (pp de visita)", seg["Segmento"], seg["Efecto (pp)"].values,
                 fmt="pp", subtitle="Diferencia de medias dentro del segmento; IC en la tabla"),
            bars("pol", "Valor de la política: % de visita esperado (IPW)", list(politicas), [round(v, 2) for v in politicas.values()],
                 fmt="pct", subtitle="Política del modelo asigna: " + ", ".join(f"{k} {v:.0%}" for k, v in mix.items())),
            bars("smd", "Balance de la aleatorización (|SMD| < 0.1 = balanceado)", balance.index,
                 balance.round(4).values, fmt="dec"),
        ],
        "tables": [table("Efectos por brazo del experimento", brazo_tbl),
                   table("Validación cruzada del Qini por learner", cv_tbl),
                   table("Efecto por segmento con IC 95 %", seg),
                   table("Calibración del uplift por decil", dec_tbl)],
        "card": model_card(
            "Meta-learners de uplift (S, T y transformed outcome) con LightGBM poco profundo; política de "
            "tres brazos con S-learner multi-tratamiento.",
            ["Chequeo de aleatorización con diferencias de medias estandarizadas.",
             "Qini con validación cruzada estratificada (tratamiento × respuesta) de 5 folds.",
             "Calibración del uplift por decil y correlación de Spearman.",
             "Valor de políticas fuera de muestra con ponderación por probabilidad inversa (IPW)."],
            ["Aleatorización correcta (verificada) y sin interferencia entre clientes (SUTVA).",
             "La visita es un buen proxy de valor; el gasto es muy escaso (≈1 % convierte)."],
            ["El uplift individual es ruidoso: las decisiones se toman por segmentos o deciles, no por cliente.",
             "Experimento de 2008 de un minorista; el patrón se debe revalidar con un piloto propio."]),
    }
    sample = pd.DataFrame({"uplift": uplift.round(4), "visita": y_te, "tratamiento": t_te})
    export("p03-uplift", result, sample.sample(5000, random_state=1))
    return result


if __name__ == "__main__":
    run()