notebooks/simpsons-paradox-examples-exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9e6015f0",
   "metadata": {},
   "source": [
    "# Simpson's paradox\n",
    "\n",
    "## Example from Causality Primer\n",
    "\n",
    "From [Primer](http://bayes.cs.ucla.edu/PRIMER/primer-ch1.pdf), table 1.1. is\n",
    "\n",
    "Subpopulations | No Treatment | Treatment \n",
    "------------ | ------------- | ------\n",
    "Male |  234 of 270 recover (87%) | 81 of 87 recover (93%) \n",
    "Female | 55 of 80 recover (69%) | 192 of 263 recover (73%)\n",
    "Total | 289 of 350 (83%)) | 273 of 350 (78%)\n",
    "\n",
    "\n",
    "We exchange the rows due to keys having to be in alphabetical order for the fake-data-for-learning package (see [this issue](https://github.com/munichpavel/fake-data-for-learning/issues/8)):\n",
    "\n",
    "Subpopulations | No Treatment | Treatment \n",
    "------------ | ------------- | ------\n",
    "Female | 55 of 80 recover (69%) | 192 of 263 recover (73%)\n",
    "Male |  234 of 270 recover (87%) | 81 of 87 recover (93%) \n",
    "Total | 289 of 350 (83%)) | 273 of 350 (78%)\n",
    "\n",
    "\n",
    "\n",
    "## Equations for Simpson paradox\n",
    "\n",
    "Keeping with the above example, we consider the counts above as being derived from a space of counts along dimensions $(\\mathrm{RECOVERED}, \\mathrm{GENDER}, \\mathrm{TREATED})$ of $\\mathbb{N}^2 \\times \\mathbb{N}^2  \\times \\mathbb{N}^2$:\n",
    "\n",
    "$$\\begin{align*}\n",
    "u_{000} &= \\textrm{Count of non-recovered females who received no treatment} \\\\\n",
    "u_{100} &= \\textrm{Count of recovered females who received no treatment} \\\\\n",
    "u_{010} &= \\textrm{Count of non-recovered males who received no treatment} \\\\\n",
    "u_{110} &= \\textrm{Count of recovered males who received no treatment} \\\\\n",
    "u_{001} &= \\textrm{Count of non-recovered females who received treatment} \\\\\n",
    "u_{101} &= \\textrm{Count of recovered females who received treatment} \\\\\n",
    "u_{011} &= \\textrm{Count of non-recovered males who received treatment} \\\\\n",
    "u_{111} &= \\textrm{Count of recovered males who received treatment}\n",
    "\\end{align*}$$\n",
    "\n",
    "To denote sums along one or more of the count dimensions of RECOVERED, GENDER, TREATED, we substitute the dimension along which we sum by a $+$. So $u_{00+}=147$ is the count of recovered females, regardless of treatment. Hence we can re-write the above table as \n",
    "\n",
    "Subpopulations | No Treatment | Treatment \n",
    ":----------: | :-----------: | :----: \n",
    "Female |  $u_{100}$ of $u_{+00}$ recover (ratio $\\frac{u_{100}}{u_{+00}}$) | $u_{101}$ of $u_{+01}$ recover (ratio $\\frac{n_{101}}{n_{+01}})$ \n",
    "Male | $u_{110}$ of $u_{+10}$  recover (ratio $\\frac{u_{110}}{u_{+10}})$ | $u_{111}$ of $u_{+11}$ recover (ratio $\\frac{u_{111}}{u_{+11}})$\n",
    "Total |  $u_{1+0}$ of $u_{++0}$  recover $\\left(\\frac{u_{1+1}}{u_{++0}}\\right)$ |  $n_{1+1}$ of $u_{++1}$  recover (ratio $\\frac{u_{1+1}}{u_{++1}})$ \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a76eb1d",
   "metadata": {},
   "source": [
    "## Example: Contingency tables\n",
    "\n",
    "For computation, the $u_{ijk}$ lend themselves to representing experiment outcomes grouped by sub-populations as contigency tables; see [Lectures on Algebraic Statistics](http://math.berkeley.edu/~bernd/owl.pdf), Section 1.1]. \n",
    "\n",
    "In this example, the contingency table is a 2x2x2 matrices of event counts. Using [xarray](http://xarray.pydata.org/en/stable/), we can define this matrix in a convenient form for defining, accessing and performing operations (e.g. the marginal sums of the $+$ notation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1d03d303",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "\n",
    "from fake_data_for_learning import utils as ut\n",
    "from risk_learning.simpson import compute_margin, transform_data_array_component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "04b73840",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check against counts from contingency table above:\n",
      "\n",
      "Counts for event(s) {'recovered': 1, 'gender': 'female', 'treated': 0}: \n",
      "55\n",
      "Counts for event(s) {'recovered': 1, 'gender': 'female', 'treated': 1}: \n",
      "192\n",
      "Counts for event(s) {'recovered': 1, 'gender': 'male'}: \n",
      "[234  81]\n",
      "Counts for event(s) {'recovered': 1, 'gender': 'male'}: \n",
      "[234  81]\n"
     ]
    }
   ],
   "source": [
    "U = [\n",
    "    [\n",
    "        [25, 71],\n",
    "        [36, 6]\n",
    "        \n",
    "    ], \n",
    "    [\n",
    "         [55, 192],\n",
    "        [234, 81]\n",
    "    ]\n",
    "       \n",
    "]\n",
    "gender_vals = ['female', 'male']\n",
    "treated_vals = [0, 1]\n",
    "recovered_vals = [0, 1]\n",
    "contingency = xr.DataArray(\n",
    "    U, \n",
    "    dims=('recovered', 'gender', 'treated'),\n",
    "    coords=[recovered_vals, gender_vals, treated_vals]\n",
    ")\n",
    "\n",
    "print(\"Check against counts from contingency table above:\\n\")\n",
    "events = [\n",
    "    dict(recovered=1, gender='female', treated=0),\n",
    "    dict(recovered=1, gender='female', treated=1),\n",
    "    dict(recovered=1, gender='male'),\n",
    "    dict(recovered=1, gender='male'),\n",
    "\n",
    "]\n",
    "\n",
    "for event in events:\n",
    "    marginal_counts = contingency.sel(**event)\n",
    "    print(f'Counts for event(s) {event}: \\n{marginal_counts.values}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70ce8203",
   "metadata": {},
   "source": [
    "## Exercise: Recovery rates with `xarray`\n",
    "\n",
    "Difficulty: (*)\n",
    "\n",
    "Determine if the recovery rate $\\left(=\\frac{\\textrm{subpopulation-size-if-recovered}}{\\textrm{subpopulation-size}}\\right)$ is higher or lower for the following subpopulations\n",
    "\n",
    "1. females who do not recover vs males who do not recover\n",
    "1. females who did not receive treatment and recover vs males who did not receive treatment and recover"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b109bd",
   "metadata": {},
   "source": [
    "## Exercise: No-simpson\n",
    "\n",
    "Difficulty: (**)\n",
    "\n",
    "Show that if all sub-population counts $u_{ijk}$ in the above $(\\mathrm{RECOVERED}, \\mathrm{GENDER}, \\mathrm{TREATED})$ of $\\mathbb{N}^2 \\times \\mathbb{N}^2  \\times \\mathbb{N}^2$ are all equal to a constant $u_{i} > 0$ for all $j,k$, then Simpson's paradox cannot occur.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26dbc1b",
   "metadata": {},
   "source": [
    "## Exercise: Simple linear coordinate change about a Simpson point\n",
    "\n",
    "Difficulty: (**)\n",
    "\n",
    "Assume $p_{ijk}$ as before ($i,j,k \\in \\{0,1\\}$) such that the sub-population and total population (in)equalities are \n",
    "\n",
    "$$\\begin{align}\n",
    "p_{101}p_{+00} -p_{100}p_{+01} > 0 \\\\\n",
    "p_{111}p_{+10} -p_{110}p_{+11} = 0 \n",
    "\\end{align}$$\n",
    "\n",
    "$$\\begin{equation}\n",
    "p_{1+1}p_{++0} - p_{1+0}p_{++1}  = 0\n",
    "\\end{equation}$$\n",
    "\n",
    "For $\\epsilon > 0$, consider the coordinate change\n",
    "\n",
    "$$\\begin{align}\n",
    "p_{101} \\mapsto p_{101} + \\epsilon \\\\\n",
    "p_{100} \\mapsto p_{100} - \\epsilon\n",
    "\\end{align}$$\n",
    "\n",
    "where $\\epsilon$ is chosen small enough that $p_{101} + \\epsilon < 1$ and $p_{100} - \\epsilon > 0$.\n",
    "\n",
    "This change makes the first subpopulation inequality further from zero, potentially increasing the ''paradox,'' as the treated female subpopulation fares even better.\n",
    "\n",
    "What happens to the other two (in-)equalities? Give your answer in terms of the original $p_{ijk}$ and $\\epsilon$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98d0c55a",
   "metadata": {},
   "source": [
    "## Exercise: Simpson relative-risk in terms of counts\n",
    "\n",
    "Difficulty: (*)\n",
    "\n",
    "Recall that an the relative risk for two sub-population groups is the probability of an outcome in an exposed group (=treated) to the probability of the outcome in an unexposed group (=untreated).\n",
    "\n",
    "In terms of relative risk, Simpson's paradox roughly occurs when the relative risk of one or more subpopulations is less than (resp. greater than) 1, while the trend is reversed in the total population, or no trend exists, i.e. relative risk values are 1.\n",
    "\n",
    "Assuming that all terms we need to divide by are non-zero, show that, for the female, male and total subgroups, respectively, the relative risks with event $recover=1$ and exposure event $treated=1$,\n",
    "\n",
    "$$\\begin{align}\n",
    "rr_{gender=0} &= \\frac{ u_{101} u_{+00} }{ u_{100} u_{+01} } \\\\\n",
    "rr_{gender=1} &= \\frac{ u_{111} u_{+10} }{ u_{110} u_{+11} } \\\\\n",
    "rr &= \\frac{ u_{1+1} u_{++0} }{ u_{1+0} u_{++1} }\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3388e116",
   "metadata": {},
   "source": [
    "## Exercise: more extreme Simpson case of female subpopulation examples above\n",
    "\n",
    "Difficulty: (***)\n",
    "\n",
    "Consider the count transformation\n",
    "\n",
    "$$\\begin{align}\n",
    "\\phi_{101}: u_{101} &\\mapsto u_{101} + \\epsilon \\\\\n",
    "\\phi_{111}: u_{111} &\\mapsto u_{111} + \\nu\n",
    "\\end{align}$$\n",
    "\n",
    "And assume that the gender x treatment sub-population counts remain fixed, i.e. the $u_{+jk}$ do not change (by changing the non-recovery subpopulation counts, e.g. $u_{100} \\mapsto u_{100} - \\epsilon$.)\n",
    "\n",
    "Find a $\\phi_{110}: u_{110} \\mapsto \\phi_{110}(u_{110})$ in terms of $\\nu$ that preserves the relative-risk $rr_{gender=1}$. Hint: finding a solution with the concrete numbers of the example \"Better recovery with treatment for female subpopulation, all else the same recovery\" may help.\n",
    "\n",
    "What contraints are needed on $\\epsilon$ and $\\nu$ for the transformation to be valid? Recall our assumption that sub-population counts $u_{+jk}$ remain fixed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "933a86e0",
   "metadata": {},
   "source": [
    "## Exercise: more extreme Simpson case of female subpopulation examples above, continued\n",
    "\n",
    "Difficulty: (**)\n",
    "\n",
    "Consider the count transformation\n",
    "\n",
    "$$\\begin{align}\n",
    "\\phi_{101}: u_{101} &\\mapsto u_{101} + \\epsilon \\\\\n",
    "\\phi_{111}: u_{111} &\\mapsto u_{111} + \\nu \\\\\n",
    "\\phi_{110}: u_{110} &\\mapsto \\phi_{110}(u_{110}) \\\\\n",
    "\\phi_{100}: u_{100} &\\mapsto \\phi_{100} - \\delta\n",
    "\\end{align}$$\n",
    "\n",
    "where we assume $\\epsilon, \\delta > 0$, and $\\phi_{110}$ is from your solution to the previous exercise. We also keep the assumption that the $u_{+jk}$ are fixed under all transformations.\n",
    "\n",
    "Determine conditions on $\\delta$ in terms of $\\epsilon$ and $\\nu$ to ensure that the total population relative-risk $rr$ is preserved.\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "9c43a4003e844a50c4c7376b71826f3a9eac595bfd87f8f1f8e71e580913501f"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}