criteo_benchmark.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "11f49bb1-f229-471b-b23f-6d954225d1c6",
   "metadata": {
    "tags": []
   },
   "source": [
    "**Notebook for the paper**<br/>\n",
    "**_A churn prediction dataset from the telecom sector: a new benchmark for uplift modeling_**<br/>\n",
    "_Anonymous authors_\n",
    "# Benchmark for uplift models on the Criteo dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "826c5bfc-7111-482b-a905-e0fff504fa01",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from math import sqrt\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\n",
    "from sklearn.preprocessing import normalize\n",
    "from os.path import join\n",
    "from scipy.stats import entropy\n",
    "\n",
    "from functions.wrappers import URFCWrapper, RandomForestWrapper, SLearnerWrapper\n",
    "from functions.benchmark import benchmark\n",
    "from functions.dataset import Dataset\n",
    "from functions.easy_ensemble import EasyEnsemble\n",
    "from functions.eval_measures import uplift_curve, calibrate_score"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e32c142b-fed0-41d4-80f6-52e1af263d02",
   "metadata": {},
   "source": [
    "### Data preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faa9e6e0-33a6-4df4-a253-a4c1b9fad657",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"../data/criteo-uplift-v2.1.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a97fec3c-53b8-43a4-aec9-b65099207326",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.iloc[np.random.choice(data.shape[0], replace=False, size=50000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a08847f3-8bea-4c3a-8212-ab4601616bc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "exclude_from_X = ['treatment', 'conversion', 'visit', 'exposure']\n",
    "columns_X = [c for c in data.columns if c not in exclude_from_X]\n",
    "dataset = Dataset(\n",
    "    X = data[columns_X].reset_index(drop=True).to_numpy(),\n",
    "    y = (data.visit == 1).to_numpy(),\n",
    "    t = (data.treatment == 1).to_numpy()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03b286cb-5489-4d89-ac10-d5bf81865bfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data normalization\n",
    "normalization = \"minmax\"\n",
    "if normalization == \"gaussian\":\n",
    "    m = dataset.X.mean(axis=0)\n",
    "    s = dataset.X.std(axis=0)\n",
    "    dataset.X = (dataset.X - m) / s\n",
    "elif normalization == \"minmax\":\n",
    "    M = dataset.X.max(axis=0)\n",
    "    m = dataset.X.min(axis=0)\n",
    "    dataset.X = (dataset.X - m) / (M - m)\n",
    "\n",
    "# Target variables\n",
    "dataset.t = np.array(dataset.t)\n",
    "dataset.y = np.array(dataset.y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c870c7aa-bd89-4905-833b-d791fe093774",
   "metadata": {},
   "outputs": [],
   "source": [
    "U = dataset.y[~dataset.t].mean() - dataset.y[dataset.t].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f99f972d-c561-4742-8ea8-9c37aea27e0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Uplift in this dataset: {:.2%}\".format(U))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d874d2c1-c183-4c28-8f2b-7ed2d9724198",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"../data/prepared_criteo.pickle\", 'wb') as f:\n",
    "    pickle.dump(dataset, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf51896f-a317-49ad-8c1d-832a0657a0c7",
   "metadata": {},
   "source": [
    "### Benchmark setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5549dbe4-dfed-417d-99f7-cd0874581632",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"../data/prepared_criteo.pickle\", \"rb\") as f:\n",
    "    dataset = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e15cb2c-8960-4a30-9b0f-5238540a2c08",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e2b6631-35d2-4160-8386-7bbaaa46fba1",
   "metadata": {},
   "outputs": [],
   "source": [
    "models = {}\n",
    "max_features = int(sqrt(dataset.X.shape[1]))\n",
    "n_estimators = 100\n",
    "max_depth = 20\n",
    "min_samples_leaf = 10\n",
    "n_folds = 8\n",
    "\n",
    "# RFeature with uplift random forest and EasyEnsemble\n",
    "models[\"urf_rfeature\"] = RFeature(\n",
    "    RandomForestClassifier(n_estimators=n_estimators, n_jobs=-1, random_state=seed),\n",
    "    EasyEnsemble(\n",
    "        URFCWrapper(\n",
    "            n_estimators=n_estimators,\n",
    "            max_depth=max_depth,\n",
    "            min_samples_leaf=min_samples_leaf,\n",
    "            min_samples_treatment=min_samples_leaf,\n",
    "            max_features=max_features,\n",
    "            n_jobs=1,\n",
    "            random_state=seed\n",
    "        ),\n",
    "        n_folds=n_folds,\n",
    "        n_jobs=-1,\n",
    "        random_state=seed,\n",
    "        verbose=False\n",
    "    ),\n",
    "    verbose=False\n",
    ")\n",
    "\n",
    "# Uplift random forest and EasyEnsemble without rfeature\n",
    "models[\"urf\"] = EasyEnsemble(\n",
    "    URFCWrapper(\n",
    "        n_estimators=n_estimators,\n",
    "        max_depth=max_depth,\n",
    "        min_samples_leaf=min_samples_leaf,\n",
    "        min_samples_treatment=min_samples_leaf,\n",
    "        max_features=max_features,\n",
    "        n_jobs=1,\n",
    "        random_state=seed\n",
    "    ),\n",
    "    n_folds=n_folds,\n",
    "    n_jobs=-1,\n",
    "    random_state=seed,\n",
    "    verbose=False\n",
    ")\n",
    "\n",
    "# RFeature with X-learner and EasyEnsemble\n",
    "models[\"xlearner_rfeature\"] = RFeature(\n",
    "    RandomForestClassifier(n_estimators=100, n_jobs=-1),\n",
    "    EasyEnsemble(\n",
    "        XClassifierWrapper(\n",
    "            outcome_learner=RandomForestClassifier(n_estimators=n_estimators, n_jobs=1),\n",
    "            effect_learner=RandomForestRegressor(n_estimators=n_estimators, n_jobs=1)\n",
    "        ),\n",
    "        n_folds=n_folds,\n",
    "        n_jobs=-1,\n",
    "        random_state=seed,\n",
    "        verbose=False\n",
    "    ),\n",
    "    verbose=True\n",
    ")\n",
    "\n",
    "# RFeature with X-learner and EasyEnsemble\n",
    "models[\"xlearner\"] = EasyEnsemble(\n",
    "    XClassifierWrapper(\n",
    "        outcome_learner=RandomForestClassifier(n_estimators=100, n_jobs=1),\n",
    "        effect_learner=RandomForestRegressor(n_estimators=100, n_jobs=1)\n",
    "    ),\n",
    "    n_folds=n_folds,\n",
    "    n_jobs=-1,\n",
    "    random_state=seed,\n",
    "    verbose=False\n",
    ")\n",
    "\n",
    "# Causal effect variational autoencoder\n",
    "models[\"cevae_rfeature\"] = RFeature(\n",
    "    RandomForestClassifier(n_estimators=100, n_jobs=-1, random_state=seed),\n",
    "    CEVAEWrapper(),\n",
    "    verbose=False\n",
    ")\n",
    "\n",
    "# Churn risk\n",
    "models[\"rf_rfeature\"] = RFeature(\n",
    "    RandomForestClassifier(n_estimators=100, n_jobs=-1, random_state=seed),\n",
    "    EasyEnsemble(\n",
    "        RandomForestWrapper(\n",
    "            n_jobs=1,\n",
    "            n_estimators=n_estimators,\n",
    "            max_depth=max_depth,\n",
    "            min_samples_leaf=min_samples_leaf,\n",
    "            max_features=max_features,\n",
    "            random_state=seed\n",
    "        ),\n",
    "        n_folds=n_folds,\n",
    "        n_jobs=-1,\n",
    "        random_state=seed,\n",
    "        verbose=False\n",
    "    ),\n",
    "    verbose=False\n",
    ")\n",
    "\n",
    "# Churn risk no reach\n",
    "models[\"rf\"] = EasyEnsemble(\n",
    "    RandomForestWrapper(\n",
    "        n_jobs=1,\n",
    "        n_estimators=n_estimators,\n",
    "        max_depth=max_depth,\n",
    "        min_samples_leaf=min_samples_leaf,\n",
    "        max_features=max_features,\n",
    "        random_state=seed\n",
    "    ),\n",
    "    n_folds=n_folds,\n",
    "    n_jobs=-1,\n",
    "    random_state=seed,\n",
    "    verbose=False\n",
    ")\n",
    "\n",
    "# RF S learner\n",
    "models[\"rf_slearner\"] = EasyEnsemble(\n",
    "    SLearnerWrapper(\n",
    "        n_jobs=1,\n",
    "        n_estimators=n_estimators,\n",
    "        max_depth=max_depth,\n",
    "        min_samples_leaf=min_samples_leaf,\n",
    "        max_features=max_features,\n",
    "        random_state=seed\n",
    "    ),\n",
    "    n_folds=n_folds,\n",
    "    n_jobs=-1,\n",
    "    random_state=seed,\n",
    "    verbose=False\n",
    ")\n",
    "\n",
    "models_to_run = [\"rf\", \"rf_slearner\", \"urf\"]\n",
    "models = {name: models[name] for name in models_to_run}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "616c41e4-3456-4a37-9631-a0864cb8b1a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_file = \"results/benchmark_criteo.pickle\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b493dce-acfc-4946-a12a-030fc10eee89",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_repeats = 10\n",
    "k_folds = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0114336-6ae8-48a4-b960-9ab84ee615d7",
   "metadata": {},
   "source": [
    "### Benchmark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77a99b73-7e5e-4c21-b9c4-728a0d79f37d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = benchmark(\n",
    "    dataset,\n",
    "    models=models,\n",
    "    k_folds=k_folds,\n",
    "    n_repeats=n_repeats,\n",
    "    seed=None,\n",
    "    verbose=False,\n",
    "    predict_params={\n",
    "        \"urf_rfeature\": {\"full_output\": True},\n",
    "        \"urf\": {\"full_output\": True}\n",
    "    }\n",
    ")\n",
    "if True:\n",
    "    with open(results_file, 'wb') as f:\n",
    "        pickle.dump(results, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "317aeaf1-10aa-47d9-bbde-bb2e68e71676",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(results_file, \"rb\") as f:\n",
    "    results = pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b2d453e-5b99-4391-bd04-dbf9a2864d27",
   "metadata": {},
   "source": [
    "### Estimating the mutual information and the estimator variance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7559181b-e9c6-42f5-88fd-5b3b7eeb6280",
   "metadata": {},
   "source": [
    "First, build the array of predictions of each model, to compute the variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9535bd7-d9aa-4858-8921-e5ea425eceac",
   "metadata": {},
   "outputs": [],
   "source": [
    "preds = {}\n",
    "for model_name in models:\n",
    "    all_pred = np.empty((dataset.X.shape[0], n_repeats))\n",
    "    for i, split in enumerate(results):\n",
    "        i_repeat = i // k_folds\n",
    "        pred = split[\"results\"][model_name][\"pred\"]\n",
    "        if model_name == \"urf\" or model_name == \"rf_slearner\":\n",
    "            pred = pred[\"control\"] - pred[\"target\"]\n",
    "        all_pred[split[\"test_indices\"], i_repeat] = pred\n",
    "    all_pred[np.isnan(all_pred)] = 0\n",
    "    all_pred[all_pred > 1e10] = 0 # Sometimes some models give huge scores\n",
    "    preds[model_name] = all_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c53e66c0-2a77-43ca-997d-d2557e61dd82",
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimator_variance(preds):\n",
    "    preds[preds <= -1] = -1\n",
    "    preds[preds >= 1] = 1\n",
    "    preds[np.isnan(preds)] = 0\n",
    "    return np.mean(np.var(preds, axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4adf0f86-e373-4f81-992c-c5a7895d1df9",
   "metadata": {},
   "outputs": [],
   "source": [
    "variances = {model_name: estimator_variance(preds[model_name]) for model_name in models}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5833b30-ce38-496f-924c-929334561759",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(\"Var_u = {:.2%}\".format(variances[\"urf\"]))\n",
    "#print(\"Var_p = {:.2%}\".format(variances[\"rf\"]))\n",
    "print(\"Var_s = {:.2e}\".format(variances[\"rf_slearner\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac9684fa-06ba-451f-81eb-087edd121b74",
   "metadata": {},
   "source": [
    "Then, estimate the individual counterfactual probabilities to estimate the mutual information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5ce7078-d382-4d1a-a596-bb3f490e39b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_0_prior = np.mean(dataset.y[~dataset.t])\n",
    "S_1_prior = np.mean(dataset.y[dataset.t])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e1302b6-30ae-4e17-9527-d386a966904f",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_0 = np.empty((dataset.X.shape[0], n_repeats))\n",
    "S_1 = np.empty((dataset.X.shape[0], n_repeats))\n",
    "\n",
    "for i, split in enumerate(results):\n",
    "    i_repeat = i // k_folds\n",
    "    S_0_split = split[\"results\"][\"rf_slearner\"][\"pred\"][\"control\"]\n",
    "    S_1_split = split[\"results\"][\"rf_slearner\"][\"pred\"][\"target\"]\n",
    "    #S_0_split = calibrate_score(S_0_split, S_0_prior)\n",
    "    #S_1_split = calibrate_score(S_1_split, S_1_prior)\n",
    "    S_0[split[\"test_indices\"], i_repeat] = S_0_split\n",
    "    S_1[split[\"test_indices\"], i_repeat] = S_1_split\n",
    "\n",
    "S_0[np.isnan(S_0)] = 0\n",
    "S_1[np.isnan(S_1)] = 0\n",
    "S_0[S_0 >= 1] = 1\n",
    "S_1[S_1 >= 1] = 1\n",
    "S_0[S_0 <= 0] = 0\n",
    "S_1[S_1 <= 0] = 0\n",
    "S_0 = np.mean(S_0, axis=1)\n",
    "S_1 = np.mean(S_1, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef4dc3e0-f634-40d8-8607-40a6345c2e6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calibrate two times for more precision\n",
    "S_0 = calibrate_score(calibrate_score(S_0, S_0_prior), S_0_prior)\n",
    "S_1 = calibrate_score(calibrate_score(S_1, S_1_prior), S_1_prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a95cd4ba-54eb-4519-8b11-f04c857e834c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def counterfactuals(S_0_x, S_1_x):\n",
    "    alpha_x = (1 - S_0_x) * (1 - S_1_x)\n",
    "    beta_x = S_0_x * (1 - S_1_x)\n",
    "    gamma_x = (1 - S_0_x) * S_1_x\n",
    "    delta_x = S_0_x * S_1_x\n",
    "    alpha = np.mean(alpha_x)\n",
    "    beta = np.mean(beta_x)\n",
    "    gamma = np.mean(gamma_x)\n",
    "    delta = np.mean(delta_x)\n",
    "    return (alpha, beta, gamma, delta)\n",
    "\n",
    "def mutual_information_marginal(S_x, S):\n",
    "    return entropy([S, 1 - S]), np.mean(entropy(np.vstack((S_x, 1 - S_x))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81c8682-f7a8-440f-b274-4da7e9ba9b46",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha, beta, gamma, delta = counterfactuals(S_0, S_1)\n",
    "print(\"alpha = {:.1%}\".format(alpha))\n",
    "print(\"beta  = {:.1%}\".format(beta))\n",
    "print(\"gamma = {:.1%}\".format(gamma))\n",
    "print(\"delta = {:.1%}\".format(delta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "967d0486-0eee-48e3-aa2b-0627ca071c44",
   "metadata": {},
   "outputs": [],
   "source": [
    "H_Y_0, H_Y_0_X = mutual_information_marginal(S_0, S_0_prior)\n",
    "H_Y_1, H_Y_1_X = mutual_information_marginal(S_1, S_1_prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c77abbfa-3367-49ae-90ac-d1457cad79d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"I_0 = {:.2%}\".format(1 - H_Y_0_X / H_Y_0))\n",
    "print(\"I_1 = {:.2%}\".format(1 - H_Y_1_X / H_Y_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "394db083-1541-4aa1-964f-731b04d4b920",
   "metadata": {},
   "source": [
    "### Plotting uplift curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4485c7af-3d57-4aae-aa87-9723ddb79a4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "curves = {}\n",
    "auucs = {}\n",
    "\n",
    "auucs[\"random\"] = (S_0_prior - S_1_prior) / 2\n",
    "N = dataset.X.shape[0]\n",
    "\n",
    "for model_name in models:\n",
    "    curves[model_name] = []\n",
    "    auucs[model_name] = []\n",
    "    for i, split in enumerate(results):\n",
    "        i_repeat = i // k_folds\n",
    "        pred = split[\"results\"][model_name][\"pred\"]\n",
    "        if model_name.startswith(\"urf\") or model_name == \"rf_slearner\":\n",
    "            pred = pred[\"control\"] - pred[\"target\"]\n",
    "        if model_name in (\"urf_rfeature\", \"urf\", \"cevae\", \"xlearner\", \"xlearner_rfeature\", \"rf_slearner\"):\n",
    "            pred = -pred\n",
    "        pred[np.isnan(pred)] = 0\n",
    "        pred[pred > 1e10] = 0 # Sometimes some models give huge scores\n",
    "        curve = uplift_curve(dataset.y[split[\"test_indices\"]], dataset.t[split[\"test_indices\"]], pred)\n",
    "        curves[model_name].append(curve)\n",
    "        auucs[model_name].append(-curve.profit.mean() / len(split[\"test_indices\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8fdd305-e0ec-4b55-a11e-0a38ec50455d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(\"AUUC_u = {:.2%}\".format(np.mean(auucs[\"urf\"])))\n",
    "#print(\"AUUC_p = {:.2%}\".format(np.mean(auucs[\"rf\"])))\n",
    "print(\"AUUC_p = {:.2%}\".format(np.mean(auucs[\"rf_slearner\"])))"
   ]
  }
 ],
 "metadata": {
  "environment": {
   "kernel": "conda-root-py",
   "name": "r-cpu.4-1.m89",
   "type": "gcloud",
   "uri": "gcr.io/deeplearning-platform-release/r-cpu.4-1:m89"
  },
  "kernelspec": {
   "display_name": "Python [conda env:root] *",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}