docs/Examples/Morris-Trials-2022.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9a493449",
   "metadata": {},
   "source": [
    "# Morris et al. (2022): Precision in Randomized Trials\n",
    "\n",
    "The following is a replication of the analyses described in Morris et al. (2022)., which uses data from Wilson E et al. (2017). Morris et al. examine three different approaches to adjust for prognostic factors in randomized trials. The three ways are direct adjustment, standardization, and inverse probability weighting. This example highlights how causal inference methods (like those in the Hernan & Robins (2023)) can be used in randomized trials.\n",
    "\n",
    "Here, we are not going to consider direct adjustment. Attention is further restricted to the any-test outcome. For finer details and comparisons between the covariate adjustment approaches, see the original publication.\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1a0e8811",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Versions\n",
      "NumPy:         1.25.2\n",
      "pandas:        1.4.1\n",
      "Delicatessen:  2.1\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import patsy\n",
    "\n",
    "import delicatessen\n",
    "from delicatessen import MEstimator\n",
    "from delicatessen.estimating_equations import ee_regression\n",
    "from delicatessen.utilities import inverse_logit\n",
    "\n",
    "print(\"Versions\")\n",
    "print(\"NumPy:        \", np.__version__)\n",
    "print(\"pandas:       \", pd.__version__)\n",
    "print(\"Delicatessen: \", delicatessen.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be76bfcf",
   "metadata": {},
   "source": [
    "Data is available from Supplement S1 of Wilson et al. (2017) (referenced below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "63c46420",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "d = pd.read_excel(\"data/GetTestedData.xls\", sheet_name=1)\n",
    "\n",
    "# Subsetting the desired columns\n",
    "cols = ['group',      # Randomized arm\n",
    "        'anytest',    # Outcome 1 (any test)\n",
    "        'anydiag',    # Outcome 2 (any diagnosis)\n",
    "        'gender',     # gender (male, female, transgender)\n",
    "        'msm',        # MSM, other\n",
    "        'age',        # age (continuous)\n",
    "        'partners',   # Number of partners in <12 months\n",
    "        'ethnicgrp']  # Ethnicity (5 categories)\n",
    "d = d[cols].copy()\n",
    "\n",
    "# Re-coding columns as numbers\n",
    "d['group_n'] = np.where(d['group'] == 'SH:24', 1, np.nan)\n",
    "d['group_n'] = np.where(d['group'] == 'Control', 0, d['group_n'])\n",
    "d['gender_n'] = np.where(d['gender'] == 'Female', 0, np.nan)\n",
    "d['gender_n'] = np.where(d['gender'] == 'Male', 1, d['gender_n'])\n",
    "d['gender_n'] = np.where(d['gender'] == 'Transgender', 2, d['gender_n'])\n",
    "d['msm_n'] = np.where(d['msm'] == 'other', 0, np.nan)\n",
    "d['msm_n'] = np.where(d['msm'] == 'msm', 1, d['msm_n'])\n",
    "d['msm_n'] = np.where(d['msm'] == '99', 2, d['msm_n'])\n",
    "d['partners_n'] = np.where(d['partners'] == '1', 0, np.nan)\n",
    "categories = ['2', '3', '4', '5', '6', '7', '8', '9', '10+']\n",
    "for index in range(len(categories)):\n",
    "    d['partners_n'] = np.where(d['partners'] == categories[index], \n",
    "                               index, d['partners_n'])\n",
    "\n",
    "d['ethnicgrp_n'] = np.where(d['ethnicgrp'] == 'White/ White British', 0, np.nan)\n",
    "d['ethnicgrp_n'] = np.where(d['ethnicgrp'] == 'Black/ Black British', 1, d['ethnicgrp_n'])\n",
    "d['ethnicgrp_n'] = np.where(d['ethnicgrp'] == 'Mixed/ Multiple ethnicity', 2, d['ethnicgrp_n'])\n",
    "d['ethnicgrp_n'] = np.where(d['ethnicgrp'] == 'Asian/ Asian British', 3, d['ethnicgrp_n'])\n",
    "d['ethnicgrp_n'] = np.where(d['ethnicgrp'] == 'Other', 4, d['ethnicgrp_n'])\n",
    "\n",
    "# Dropping old columns and renaming new ones\n",
    "d = d.drop(columns=['group', 'gender', 'msm', 'partners', 'ethnicgrp'])\n",
    "relabs = dict()\n",
    "for c in cols:\n",
    "    relabs[c + \"_n\"] = c\n",
    "\n",
    "d = d.rename(columns=relabs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05fc505d",
   "metadata": {},
   "source": [
    "As done in the main paper, we conduct a complete-case analysis (i.e., discard all the observations with missing outcomes).\n",
    "\n",
    "As `delicatessen` interacts with NumPy arrays, we will further extract the covariates from the complete-case data set. Here, we use `patsy`, which allows for R-like formulas, to make the covariate manipulations easier. Specifically, the `C(...)` functionality generates an array of indicator variables (saving us the trouble of hand-coding disjoint indicators).\n",
    "\n",
    "Finally, we generate some copies of the data set and set the corresponding `group` values to all-1 and all-0. These sets of covariates will be used to generate the predicted outcome under each treatment for the standardization approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bbb82143",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = d.dropna().copy()\n",
    "\n",
    "# Outcome\n",
    "y = np.asarray(ds['anytest'])\n",
    "\n",
    "# Treatment\n",
    "a = np.asarray(ds['group'])\n",
    "\n",
    "# Observed covariates\n",
    "prop_score_covs = patsy.dmatrix(\"age + C(gender) + C(msm) + C(partners) + C(ethnicgrp)\", \n",
    "                                ds)\n",
    "outcome_covs = patsy.dmatrix(\"group + age + C(gender) + C(msm) + C(partners) + C(ethnicgrp)\", \n",
    "                                ds)\n",
    "\n",
    "# Setting group=a and getting matrices\n",
    "dsa = ds.copy()\n",
    "dsa['group'] = 1\n",
    "outcome_a1_covs = patsy.dmatrix(\"group + age + C(gender) + C(msm) + C(partners) + C(ethnicgrp)\", \n",
    "                                dsa)\n",
    "dsa['group'] = 0\n",
    "outcome_a0_covs = patsy.dmatrix(\"group + age + C(gender) + C(msm) + C(partners) + C(ethnicgrp)\", \n",
    "                                dsa)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bea17d43",
   "metadata": {},
   "source": [
    "## Standardization\n",
    "\n",
    "For standardization (or g-computation), we stack several estimating equations together. Let $Y$ indicate the outcome, $A$ indicate the treatment, and $W$ are the prognostic factors included in the model (age, gender, sexual orientation, number of partners, ethnicity). The parameters are $\\theta$, which we further divide into the interest parameters ($\\mu_0, \\mu_1, \\mu_2, \\mu_3$) and the nuisance parameters ($\\beta$). As their name implies, nuisance parameters are not of interest but are necessary for estimation for our interest parameters. Importantly, by stacking the estimating equations, we can account for the uncertainty in the estimation of the nuisance parameters into the uncertainty in our interest parameters automatically.\n",
    "\n",
    "The estimating equations are\n",
    "$$\\psi(Y_i, A_i, W_i; \\theta) =\n",
    "\\begin{bmatrix} \n",
    "    (\\mu_1 - \\mu_0) - \\mu_2 \\\\ \n",
    "    \\log(\\mu_1 / \\mu_0) - \\mu_3 \\\\\n",
    "    \\hat{Y}_i^{a=1} - \\mu_1 \\\\\n",
    "    \\hat{Y}_i^{a=0} - \\mu_0 \\\\\n",
    "    (Y_i - \\text{expit}(g(A_i,W_i)^T \\beta)) g(A_i,W_i)_i\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $g(A_i,W_i)$ is a specific function to generate the design matrix, and $\\hat{Y}_i^{a} = \\text{expit}(g(a,W_i)^T \\beta)$. The first equation in the stack is the risk difference (RD) and the second is the log-transformed risk ratio (RR). Notice that these are defined in terms of other parameters. Specifically, parameters from the third and fourth equations, which correspond to the risk under all-1 and the risk under all-0. The predicted values of the outcome are generated using the parameters from the final equation, which is the logistic model used to estimate the probability of $Y_i$ conditional on $A_i,W_i$.\n",
    "\n",
    "An astute reader may notice that we could actually shorten this stack of equations. Rather than estimate $\\mu_2,\\mu_3$, we could have plugged the corresponding $\\hat{Y}_i^{a}$ into the first and second equations (thereby removing the need for the third and fourth). While this is correct, we prefer the process above for three reasons:\n",
    "1. The above stacked equations give us the correspond risk under each arm, which also may be of interest (what was the risk under a=0). Therefore, no extra steps are necessary to get this estimated parameter.\n",
    "2. This second argument is more technical. Here, we use a root-finding procedure to get $\\hat{\\theta}$. During this process we plug in different values for $\\hat{\\theta}$ until we find the (approximate) zeroes. For the RR, some values during the exploration of values may result in $\\hat{Y}_i^{a=0} = 0$. This causes a division by zero and may break the root-finding procedure. By instead calculating the average of $\\hat{Y}_i^{a=0}$ across all $i$'s, we reduce the opporuntity for a rogue zero.\n",
    "3. Related to the previous reason, the optimization procedure will run much faster since we only need to take the log of a single number as opposed to as an array. While this won't make much of a difference in run-time here, it is good practice (and may be important in problems with more complex estimating equations).\n",
    "\n",
    "Below we write out the estimating equation by-hand (except for the regression component, which we use a built-in functionality), estimate with our M-estimator, and then provide the results (similar to Table 3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f8b6e272",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    # Extracting parameters from theta for ease\n",
    "    risk_diff = theta[0]\n",
    "    log_risk_ratio = theta[1]\n",
    "    risk_a1 = theta[2]\n",
    "    risk_a0 = theta[3]\n",
    "    beta = theta[4:]\n",
    "    \n",
    "    # Estimating nuisance model (outcome regression)\n",
    "    ee_log = ee_regression(theta=beta,        # beta coefficients\n",
    "                           X=outcome_covs,    # X covariates\n",
    "                           y=y,               # y\n",
    "                           model='logistic')  # logit\n",
    "\n",
    "    # Generating predicted outcome values\n",
    "    ee_ya1 = inverse_logit(np.dot(outcome_a1_covs, beta)) - risk_a1\n",
    "    ee_ya0 = inverse_logit(np.dot(outcome_a0_covs, beta)) - risk_a0\n",
    "    \n",
    "    # Estimating interest parameters\n",
    "    ee_risk_diff = np.ones(y.shape[0])*(risk_a1 - risk_a0) - risk_diff\n",
    "    ee_risk_ratio = np.ones(y.shape[0])*np.log(risk_a1 / risk_a0) - log_risk_ratio    \n",
    "    \n",
    "    # Returning stacked estimating equations (order matters)\n",
    "    return np.vstack((ee_risk_diff,  # risk difference\n",
    "                      ee_risk_ratio, # risk ratio\n",
    "                      ee_ya1,        # risk a=1\n",
    "                      ee_ya0,        # risk a=0\n",
    "                      ee_log,))      # logistic model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5f148316",
   "metadata": {},
   "outputs": [],
   "source": [
    "estr = MEstimator(psi, \n",
    "                        init=[0, 0, 0.5, 0.5, ] + [0, ]*outcome_covs.shape[1])\n",
    "estr.estimate(solver='lm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c020844a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Risk Difference: 0.260 (0.021)\n",
      "Risk Ratio:      0.795 (0.075)\n"
     ]
    }
   ],
   "source": [
    "std = np.sqrt(np.diag(estr.variance))\n",
    "fmt = '{:.3f}'\n",
    "\n",
    "print(\"Risk Difference:\", fmt.format(estr.theta[0]), \"(\"+fmt.format(std[0])+\")\")\n",
    "print(\"Risk Ratio:     \", fmt.format(estr.theta[1]), \"(\"+fmt.format(std[1])+\")\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bbdde5f",
   "metadata": {},
   "source": [
    "At this point, you may look at the results for the risk ratio and notice some differences in the third decimal place. We reserve conversation regarding this point till the end.\n",
    "\n",
    "## Inverse Probability Weighting\n",
    "\n",
    "Now, we will consider the case of inverse probability weighting. Inverse probability weighting relies on a separate nuisance model. To clarify this, our nuisance parameters will be indicated by $\\alpha$ here.\n",
    "\n",
    "The estimating equations are\n",
    "$$\\psi(Y_i, A_i, W_i; \\theta) =\n",
    "\\begin{bmatrix} \n",
    "    (\\mu_3 - \\mu_4) - \\mu_1 \\\\ \n",
    "    \\log(\\mu_3 / \\mu_4) - \\mu_2 \\\\\n",
    "    \\frac{Y_i A_i}{\\text{expit}(g(W_i)^T \\beta)} - \\mu_3 \\\\\n",
    "    \\frac{Y_i (1-A_i)}{1 - \\text{expit}(g(W_i)^T \\beta)} - \\mu_4 \\\\\n",
    "    (A_i - \\text{expit}(g(W_i)^T \\beta)) g(W_i)\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $g(W_i)$ is a specific function to generate the design matrix. As before, the first and second equations are the risk difference and risk ratio, respectively. Third and fourth equations correspond to the risk under all-1 and the risk under all-0. The final estimating equation is the propensity score model and is used to generate the propensity scores given $W_i$.\n",
    "\n",
    "Again, we write out the estimating equation by-hand (except for the regression component, which we use a built-in functionality), estimate with our M-estimator, and then provide the results (similar to Table 3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7386a203",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    # Extracting parameters from theta for ease\n",
    "    risk_diff = theta[0]\n",
    "    log_risk_ratio = theta[1]\n",
    "    risk_a1 = theta[2]\n",
    "    risk_a0 = theta[3]\n",
    "    beta = theta[4:]\n",
    "    \n",
    "    # Estimating nuisance model (outcome regression)\n",
    "    ee_log = ee_regression(theta=beta,         # beta coefficients\n",
    "                           X=prop_score_covs,  # W covariates\n",
    "                           y=a,                # a\n",
    "                           model='logistic')      # logit\n",
    "    \n",
    "    # Generating predicted propensity score\n",
    "    prop_score = inverse_logit(np.dot(prop_score_covs, beta))\n",
    "   \n",
    "    # Calculating weighted pieces\n",
    "    ee_ya1 = (a*y) / prop_score - risk_a1\n",
    "    ee_ya0 = ((1-a)*y) / (1-prop_score) - risk_a0\n",
    "    \n",
    "    # Estimating interest parameters\n",
    "    ee_risk_diff = np.ones(a.shape[0])*(risk_a1 - risk_a0) - risk_diff\n",
    "    ee_risk_ratio = np.ones(a.shape[0])*np.log(risk_a1 / risk_a0) - log_risk_ratio    \n",
    "\n",
    "    # Returning stacked estimating equations (order matters)\n",
    "    return np.vstack((ee_risk_diff,  # risk difference\n",
    "                      ee_risk_ratio, # risk ratio\n",
    "                      ee_ya1,        # risk a=1\n",
    "                      ee_ya0,        # risk a=0\n",
    "                      ee_log,))      # logistic model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bff602f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "estr = MEstimator(psi, init=[0., 0., 0.5, 0.5] + [0, ]*prop_score_covs.shape[1])\n",
    "estr.estimate(solver='lm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1ddb5cdc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Risk Difference: 0.261 (0.021)\n",
      "Risk Ratio:      0.805 (0.075)\n"
     ]
    }
   ],
   "source": [
    "std = np.sqrt(np.diag(estr.variance))\n",
    "fmt = '{:.3f}'\n",
    "\n",
    "print(\"Risk Difference:\", fmt.format(estr.theta[0]), \"(\"+fmt.format(std[0])+\")\")\n",
    "print(\"Risk Ratio:     \", fmt.format(estr.theta[1]), \"(\"+fmt.format(std[1])+\")\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c6de2d1",
   "metadata": {},
   "source": [
    "Results are the same as Morris et al. up to the third decimal place.\n",
    "\n",
    "### Differences in Results\n",
    "As promised, let's return to the differences in the answers. Why does this occur? The most likely culprit is the differences in the procedures for finding $\\hat{\\theta}$ across software.\n",
    "\n",
    "`delicatessen` operates by a root-finding procedure. Essentially, we are looking for the values of $\\hat{\\theta}$ that lead to the estimating equations all being approximately zero. Specifically, I am using the Levenberg-Marquardt algorithm in these examples. `delicatessen` uses a root-finding procedure because it simplifies building stacked estimating equations. However, this comes at the cost of speed and robustness of other approaches.\n",
    "\n",
    "Other regression software tends to proceed under a different approach. Generally, most software aims to maximize the log-likelihood. As before, there are a variety of algorithms that can be used to find the maximize the likelihood. However, most algorithms take advantage of both the log-likelihood and the score function (the derivative of the log-likelihood) to find the maximum. This means they are quite robust to sparse data while also being fast. As there is sparse data (there are many categorical variables for number of partners), this is the most likely explanation for the small differences.\n",
    "\n",
    "### Conclusion \n",
    "\n",
    "We replicated selected analyses in Morris et al. (2022) to showcase the basics of stacking estimating equations and how they are implemented in `delicatessen`. While `delicatessen` has built-in functionalities for standardization and inverse probability weighting, we opted to present them by-hand here. By-hand takes some extra work but it provides a better idea of how `delicatessen` can be used to build estimating equations. For example, each of the above approaches could further be extended to account for missing outcomes under a weaker assumption.\n",
    "\n",
    "\n",
    "### References\n",
    "\n",
    "Morris TP, Walker AS, Williamson EJ, & White IR. (2022). \"Planning a method for covariate adjustment in individually-randomised trials: a practical guide\". *Trials*, 23:328\n",
    "\n",
    "Wilson E, Free C, Morris TP, Syred J, Ahamed I, Menon-Johansson AS et al.. (2017). \"Internet-accessed sexually transmitted infection (e-STI) testing and results service: a randomised, single-blind, controlled trial\". *PLoS Medicine*, 14(12), e1002479."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}