src/semi_synthetic.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\Large \\textbf{Deriving a lower bound for hidden confounding in obs studies (Part 2)}$ "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we demonstrate an example of our testing procedure using the semi-synthetic experiments, specifically applied to the Hillstrom dataset. This approach follows the methodology discussed in the tutorial for synthetic experiments, but with a focus on a real-world dataset. We subsample from the Hillstrom randomized trial to create an observational dataset that exhibits a known, true confounding strength.\n",
    "\n",
    "For a comprehensive understanding of these semi-synthetic experiments, please refer to the appendix in our paper. In this example, we use the observational study as the target population. We study two confounders within the Hillstrom dataset: \"history,\" which is a strong confounder highly correlated with the outcome, and \"mens,\" a weaker confounder with a lesser correlation to the outcome. This example highlights how our test can effectively discern between strong and weak confounders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Required imports\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import bootstrap\n",
    "from sklearn.compose import ColumnTransformer\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.preprocessing import MinMaxScaler, OneHotEncoder\n",
    "\n",
    "from test_confounding.ate_bounds.ate_bounds import BootstrapSensitivityAnalysis\n",
    "from test_confounding.datasets.semi_synthetic.fetch_hillstrom import (\n",
    "    CAT_COVAR_HILLSTROM,\n",
    "    NUM_COVAR_HILLSTROM,\n",
    "    load_fetch_hillstrom_data,\n",
    ")\n",
    "from test_confounding.datasets.semi_synthetic.sampling import (\n",
    "    resample_data_with_confound_func,\n",
    "    subsample_df,\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we apply the testing procedure to the strong confounder \"history\". We identify a lower bound $\\Gamma_{LB}=4.5$ for the true confounding strength $\\Gamma^{*}=5.0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==== DATASET CHECKS ====\n",
      "Num observations =  64000\n",
      "P(T=1) =  0.66709375\n",
      "F-statistic(C, Y) =  1.0763267714114924\n",
      "P-value: =  2.97898074146969e-11\n",
      "RCT_ATE = 0.596796066727898\n",
      "Mean(Y) =  1.05090828125\n",
      "\n",
      "Feature importance for each feature:\n",
      "C = 0.54428333\n",
      "recency = 0.18103991\n",
      "mens = 0.11426223\n",
      "womens = 0.032999776\n",
      "zip_code = 0.04475624\n",
      "newbie = 0.040160246\n",
      "channel = 0.042498257\n",
      "\n",
      "Original data num. examples =  56319\n",
      "Downsampled data num. examples =  20728\n",
      "RCT num. samples:  7681\n",
      "\n"
     ]
    }
   ],
   "source": [
    "seed = 50\n",
    "confound_func = {\"para_form\": \"adversarial_with_conf\", \"true_gamma\": 5.0} #true conf strenght 5.0\n",
    "\n",
    "obs_data_pre_conf, rct_data = load_fetch_hillstrom_data(\n",
    "    conf_var=\"history\",\n",
    "    support_var=\"zip_code\",\n",
    "    split_data=True,\n",
    "    support_feature_values=[0.0, 1.0],\n",
    "    proportion_full_support=0.8,\n",
    "    seed=52,\n",
    "    target_col=\"spend\",\n",
    ")\n",
    "\n",
    "obs_data = resample_data_with_confound_func(\n",
    "    obs_data_pre_conf,\n",
    "    confound_func_params=confound_func,\n",
    "    seed=seed,\n",
    ")\n",
    "\n",
    "numeric_covariates = NUM_COVAR_HILLSTROM\n",
    "categorical_covariates = CAT_COVAR_HILLSTROM\n",
    "\n",
    "# Define transformer for encoding and normalization\n",
    "transformer = ColumnTransformer(\n",
    "    [\n",
    "        (\n",
    "            \"encoder\",\n",
    "            OneHotEncoder(handle_unknown=\"ignore\", sparse_output=False),\n",
    "            [col for col in categorical_covariates if col != \"history\"],\n",
    "        ),\n",
    "        (\n",
    "            \"normalizer\",\n",
    "            MinMaxScaler(),\n",
    "            [col for col in numeric_covariates if col != \"history\"],\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "print(\"RCT num. samples: \", rct_data.shape[0])  # type: ignore\n",
    "print()\n",
    "\n",
    "# Transform the features into one-hot encoding\n",
    "x_obs_raw = obs_data.drop([\"Y\", \"T\", \"C\"], axis=1)\n",
    "x_rct_raw = rct_data.drop([\"Y\", \"T\", \"C\"], axis=1)\n",
    "x_obs_encoded = transformer.fit_transform(x_obs_raw)\n",
    "x_rct_encoded = transformer.transform(x_rct_raw)\n",
    "\n",
    "t_rct, y_rct, x_rct = (\n",
    "    np.array(rct_data[\"T\"].values),\n",
    "    np.array(rct_data[\"Y\"].values),\n",
    "    x_rct_encoded,\n",
    ")\n",
    "t_obs, y_obs, x_obs = (\n",
    "    np.array(obs_data[\"T\"].values),\n",
    "    np.array(obs_data[\"Y\"].values),\n",
    "    x_obs_encoded,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile functions are now trained for QB. Starting bootstrap.\n",
      "Elapsed time for 50 bootstrap samples: 92.77 seconds\n"
     ]
    }
   ],
   "source": [
    "n_bootstrap = 50\n",
    "user_conf = [3.5, 4.0, 4.5, 5.0, 5.5]\n",
    "\n",
    "n_rct = x_rct.shape[0]\n",
    "n_obs = x_obs.shape[0]\n",
    "alpha_trim = 0.001\n",
    "x_all = np.concatenate((x_rct, x_obs), axis=0)  \n",
    "s_all = np.concatenate((np.ones(n_rct), np.zeros(n_obs)))\n",
    "\n",
    "clf = RandomForestClassifier(random_state=seed)\n",
    "clf.fit(x_all, s_all)\n",
    "pi_s = clf.predict_proba(x_all)[:, 1]\n",
    "O_idx = pi_s > alpha_trim\n",
    "\n",
    "x_obs_trimm, t_obs_trimm, y_obs_trimm = (\n",
    "    x_obs[O_idx[n_rct:]],  \n",
    "    t_obs[O_idx[n_rct:]],\n",
    "    y_obs[O_idx[n_rct:]],\n",
    ")\n",
    "_, t_rct_trimm, y_rct_trimm = (\n",
    "    x_rct[O_idx[:n_rct]],  \n",
    "    t_rct[O_idx[:n_rct]],\n",
    "    y_rct[O_idx[:n_rct]],\n",
    ")\n",
    "\n",
    "mask = np.logical_and(O_idx, s_all)\n",
    "rct_to_obs_ratio = s_all[O_idx].sum() / (s_all[O_idx].size - s_all[O_idx].sum())\n",
    "ys = (\n",
    "    2\n",
    "    * (y_rct_trimm * t_rct_trimm - y_rct_trimm * (1 - t_rct_trimm))\n",
    "    * (1 - pi_s[mask])\n",
    "    / pi_s[mask]\n",
    ")\n",
    "bootstrap_rct = bootstrap(\n",
    "    (ys,), np.mean, n_resamples=n_bootstrap, axis=0, random_state=seed\n",
    ")\n",
    "std_rct = bootstrap_rct.standard_error\n",
    "var_rct = np.power(std_rct, 2) * (rct_to_obs_ratio**2)\n",
    "mean_rct = rct_to_obs_ratio * ys.mean()\n",
    "\n",
    "bootstrap_sa = BootstrapSensitivityAnalysis(\n",
    "    sa_name=\"QB\",\n",
    "    inputs=x_obs_trimm,\n",
    "    treatment=t_obs_trimm,\n",
    "    outcome=y_obs_trimm,\n",
    "    gammas=user_conf,\n",
    "    seed=seed,\n",
    ")\n",
    "bounds_dist = bootstrap_sa.bootstrap(num_samples=n_bootstrap, fast_quantile=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from test_confounding.test import run_multiple_ate_hypothesis_test\n",
    "\n",
    "results_dict_ate_strong_conf = run_multiple_ate_hypothesis_test(mean_rct = mean_rct, var_rct = var_rct, bounds_dist = bounds_dist, alpha = 5.0, user_conf = user_conf, verbose = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test Results:\n",
      "------------------------------\n",
      "Γ: 3.5\n",
      "  - Reject Null Hypothesis: Yes\n",
      "  - Test Statistic: -2.6669749399580023\n",
      "------------------------------\n",
      "Γ: 4.0\n",
      "  - Reject Null Hypothesis: Yes\n",
      "  - Test Statistic: -2.227974144909513\n",
      "------------------------------\n",
      "Γ: 4.5\n",
      "  - Reject Null Hypothesis: No\n",
      "  - Test Statistic: -1.8235700123565766\n",
      "------------------------------\n",
      "Γ: 5.0\n",
      "  - Reject Null Hypothesis: No\n",
      "  - Test Statistic: N/A\n",
      "------------------------------\n",
      "Γ: 5.5\n",
      "  - Reject Null Hypothesis: No\n",
      "  - Test Statistic: N/A\n",
      "------------------------------\n",
      "Γ_{LB} (Lower Bound): 4.5\n"
     ]
    }
   ],
   "source": [
    "def display_gamma_info(gamma_dict, user_conf): \n",
    "    print(\"Test Results:\")\n",
    "    print(\"-\" * 30)\n",
    "    \n",
    "    for gamma in user_conf:\n",
    "        if gamma in gamma_dict:\n",
    "            value = gamma_dict[gamma]\n",
    "            reject = \"Yes\" if value['reject'] else \"No\"\n",
    "            test_statistic = value['test_statistic']\n",
    "        else:\n",
    "            # If gamma is not in gamma_dict, assume hypothesis is accepted\n",
    "            reject = \"No\"\n",
    "            test_statistic = \"N/A\"\n",
    "\n",
    "        print(f\"Γ: {gamma}\")\n",
    "        print(f\"  - Reject Null Hypothesis: {reject}\")\n",
    "        print(f\"  - Test Statistic: {test_statistic}\")\n",
    "        print(\"-\" * 30)\n",
    "    \n",
    "    gamma_effective = gamma_dict.get('gamma_effective', 'Not available')\n",
    "    print(f\"Γ_{{LB}} (Lower Bound): {gamma_effective}\")\n",
    "\n",
    "\n",
    "display_gamma_info(results_dict_ate_strong_conf, user_conf)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we test for the weak confounder \"mens\". We determine a lower bound $\\Gamma_{LB}=2.5$ for the true confounding strength $\\Gamma^{*}=5.0$, demonstrating that a confounder-outcome correlation that is lower results in reduced power for the test. An intuition for this phenomenon is provided in the appendix of the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==== DATASET CHECKS ====\n",
      "Num observations =  64000\n",
      "P(C=1) =  0.55103125\n",
      "P(T=1) =  0.66709375\n",
      "Corr(C, Y) =  0.008599413005626647\n",
      "P-value =  0.02959294801554325\n",
      "RCT_ATE = 0.596796066727898\n",
      "Mean(Y) =  1.05090828125\n",
      "\n",
      "Feature importance for each feature:\n",
      "C = 0.18518682\n",
      "recency = 0.14364073\n",
      "womens = 0.027128274\n",
      "zip_code = 0.061120052\n",
      "newbie = 0.035805676\n",
      "channel = 0.038889784\n",
      "history = 0.50822866\n",
      "\n",
      "Original data num. examples =  56319\n",
      "Downsampled data num. examples =  20678\n",
      "RCT num. samples:  7681\n",
      "\n"
     ]
    }
   ],
   "source": [
    "obs_data_pre_conf, rct_data = load_fetch_hillstrom_data(\n",
    "    conf_var=\"mens\",\n",
    "    support_var=\"zip_code\",\n",
    "    split_data=True,\n",
    "    support_feature_values=[0.0, 1.0],\n",
    "    proportion_full_support=0.8,\n",
    "    seed=52,\n",
    "    target_col=\"spend\",\n",
    ")\n",
    "\n",
    "obs_data = resample_data_with_confound_func(\n",
    "    obs_data_pre_conf,\n",
    "    confound_func_params=confound_func,\n",
    "    seed=seed,\n",
    ")\n",
    "\n",
    "# Define transformer for encoding and normalization\n",
    "transformer = ColumnTransformer(\n",
    "    [\n",
    "        (\n",
    "            \"encoder\",\n",
    "            OneHotEncoder(handle_unknown=\"ignore\", sparse_output=False),\n",
    "            [col for col in categorical_covariates if col != \"mens\"],\n",
    "        ),\n",
    "        (\n",
    "            \"normalizer\",\n",
    "            MinMaxScaler(),\n",
    "            [col for col in numeric_covariates if col != \"mens\"],\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "print(\"RCT num. samples: \", rct_data.shape[0])  # type: ignore\n",
    "print()\n",
    "\n",
    "# Transform the features into one-hot encoding\n",
    "x_obs_raw = obs_data.drop([\"Y\", \"T\", \"C\"], axis=1)\n",
    "x_rct_raw = rct_data.drop([\"Y\", \"T\", \"C\"], axis=1)\n",
    "x_obs_encoded = transformer.fit_transform(x_obs_raw)\n",
    "x_rct_encoded = transformer.transform(x_rct_raw)\n",
    "\n",
    "t_rct, y_rct, x_rct = (\n",
    "    np.array(rct_data[\"T\"].values),\n",
    "    np.array(rct_data[\"Y\"].values),\n",
    "    x_rct_encoded,\n",
    ")\n",
    "t_obs, y_obs, x_obs = (\n",
    "    np.array(obs_data[\"T\"].values),\n",
    "    np.array(obs_data[\"Y\"].values),\n",
    "    x_obs_encoded,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile functions are now trained for QB. Starting bootstrap.\n",
      "Elapsed time for 50 bootstrap samples: 56.20 seconds\n"
     ]
    }
   ],
   "source": [
    "user_conf = [1.0, 1.5, 2.0, 2.5]\n",
    "\n",
    "n_rct = x_rct.shape[0]\n",
    "n_obs = x_obs.shape[0]\n",
    "alpha_trim = 0.001\n",
    "x_all = np.concatenate((x_rct, x_obs), axis=0)  \n",
    "s_all = np.concatenate((np.ones(n_rct), np.zeros(n_obs)))\n",
    "\n",
    "clf = RandomForestClassifier(random_state=seed)\n",
    "clf.fit(x_all, s_all)\n",
    "pi_s = clf.predict_proba(x_all)[:, 1]\n",
    "O_idx = pi_s > alpha_trim\n",
    "\n",
    "x_obs_trimm, t_obs_trimm, y_obs_trimm = (\n",
    "    x_obs[O_idx[n_rct:]],  \n",
    "    t_obs[O_idx[n_rct:]],\n",
    "    y_obs[O_idx[n_rct:]],\n",
    ")\n",
    "_, t_rct_trimm, y_rct_trimm = (\n",
    "    x_rct[O_idx[:n_rct]],  \n",
    "    t_rct[O_idx[:n_rct]],\n",
    "    y_rct[O_idx[:n_rct]],\n",
    ")\n",
    "\n",
    "mask = np.logical_and(O_idx, s_all)\n",
    "rct_to_obs_ratio = s_all[O_idx].sum() / (s_all[O_idx].size - s_all[O_idx].sum())\n",
    "ys = (\n",
    "    2\n",
    "    * (y_rct_trimm * t_rct_trimm - y_rct_trimm * (1 - t_rct_trimm))\n",
    "    * (1 - pi_s[mask])\n",
    "    / pi_s[mask]\n",
    ")\n",
    "bootstrap_rct = bootstrap(\n",
    "    (ys,), np.mean, n_resamples=n_bootstrap, axis=0, random_state=seed\n",
    ")\n",
    "std_rct = bootstrap_rct.standard_error\n",
    "var_rct = np.power(std_rct, 2) * (rct_to_obs_ratio**2)\n",
    "mean_rct = rct_to_obs_ratio * ys.mean()\n",
    "\n",
    "bootstrap_sa = BootstrapSensitivityAnalysis(\n",
    "    sa_name=\"QB\",\n",
    "    inputs=x_obs_trimm,\n",
    "    treatment=t_obs_trimm,\n",
    "    outcome=y_obs_trimm,\n",
    "    gammas=user_conf,\n",
    "    seed=seed,\n",
    ")\n",
    "bounds_dist = bootstrap_sa.bootstrap(num_samples=n_bootstrap, fast_quantile=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test Results:\n",
      "------------------------------\n",
      "Γ: 1.0\n",
      "  - Reject Null Hypothesis: Yes\n",
      "  - Test Statistic: -2.971132554321786\n",
      "------------------------------\n",
      "Γ: 1.5\n",
      "  - Reject Null Hypothesis: Yes\n",
      "  - Test Statistic: -2.6057077484500066\n",
      "------------------------------\n",
      "Γ: 2.0\n",
      "  - Reject Null Hypothesis: Yes\n",
      "  - Test Statistic: -2.0593676567374417\n",
      "------------------------------\n",
      "Γ: 2.5\n",
      "  - Reject Null Hypothesis: No\n",
      "  - Test Statistic: -1.436799711927924\n",
      "------------------------------\n",
      "Γ_{LB} (Lower Bound): 2.5\n"
     ]
    }
   ],
   "source": [
    "results_dict_ate_weak_conf = run_multiple_ate_hypothesis_test(mean_rct = mean_rct, var_rct = var_rct, bounds_dist = bounds_dist, alpha = 5.0, user_conf = user_conf, verbose = False)\n",
    "display_gamma_info(results_dict_ate_weak_conf, user_conf)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cate",
   "language": "python",
   "name": "python3"
  },
  "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.11.5 | packaged by conda-forge | (main, Aug 27 2023, 03:33:12) [Clang 15.0.7 ]"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "cf21b1a4de40f179e989ee5efafd4a121b190d1b06ef457fbf4e132f1ff3fa9c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}