causal-inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Neyman \"potential outcomes\" model for causal inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are $N$ subjects and $T$ possible treatments.\n",
    "Each subject is represented by a ticket. \n",
    "Ticket $j$ lists $T$ numbers, $(x_{j0}, \\ldots, x_{jT-1})$.\n",
    "The value $x_{jt}$ is the response subject $j$ will have if assigned to treatment $t$.\n",
    "(Treatment $0$ is commonly _control_ or _placebo_.)\n",
    "\n",
    "This mathematical set up embodies the _non-interference_ assumption: \n",
    "subject $j$'s response is assumed to depend only on which treatment subject $j$ receives, and not\n",
    "on the treatment any other subject receives.\n",
    "(That is not a good assumption in situations like vaccine trials, where whether one subject\n",
    "becomes infected may depend on which other subjects are vaccinated, if subjects\n",
    "may come in contact with each other.)\n",
    "\n",
    "This model is also called the _potential outcomes_ model, because it starts with the\n",
    "_potential_ outcomes each subject will have to each treatment. \n",
    "Assigning a subject to a \n",
    "treatment just reveals the potential outcome that corresponds to that treatment, for that subject. \n",
    "This model was introduced by Jerzy Neyman, the founder of the U.C. Berkeley Department of Statistics, in a 1923 paper in Polish [translated into English in 1990](https://projecteuclid.org/journals/statistical-science/volume-5/issue-4/On-the-Application-of-Probability-Theory-to-Agricultural-Experiments-Essay/10.1214/ss/1177012031.full).\n",
    "It was popularized by Donald Rubin in the 1970s and 1980s.\n",
    "\n",
    "There are generalizations of this model, including one in which the \"potential outcomes\" are random, rather than deterministic, but their distributions are fixed before assignment to treatment: if subject $j$ is assigned treatment $t$, a draw from the distribution $\\mathbb{P}_{jt}$ is observed. \n",
    "Draws for different subjects are independent.\n",
    "We shall see an example of this when we discuss [nuisance parameters](./nuisance.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Intention-to-treat (ITT) versus per-protocol\n",
    "\n",
    "An intention-to-treat analysis compares outcomes in the subjects assigned to treatment to the subjects assigned to control, whether or not the subjects assigned to a treatment received it.\n",
    "It can be though of as measuring the effect of _prescribing_ treatment, in contrast to the effect of _receiving_ treatment.\n",
    "It is common in experiments on humans that some subjects do not receive the treatment they were randomly assigned.\n",
    "Such _non-compliance_ may result from subjects breaking the blinding (e.g., noticing that their \"treatment\" is\n",
    "a sugar pill and seeking other treatment rather than receiving only the placebo) or seeking alternative treatments.\n",
    "For instance, this was an issue in the TOGETHER study of the effect of ivermectin on COVID-19, in which a substantial\n",
    "fraction of patients assigned to receive a placebo bought ivermectin over the counter and took it.\n",
    "Non-compliance may also result from subjects simply not doing as they are told.\n",
    "In general, ITT is the preferred analysis, but especially when compliance rates are low, _per-protocol_ (PP) analyisis, which compares outcomes for subjects who actually received the treatment to outcomes for subjects who actually received placebo (and no other treatment) may be more informative.\n",
    "\n",
    "For instance, in [a randomized study of the effect of colonoscopy screenings on colorectal cancer and cancer deaths published in 2022](https://www.nejm.org/doi/full/10.1056/NEJMoa2208375),\n",
    "only 42% of subjects who were \"invited\" to receive a colonoscopy actually underwent the procedure. \n",
    "The apparent effect of \"inviting\" subjects to receive a colonoscopy was small:\n",
    "the risk of death from colorectal cancer at 10 years was 0.28% in the invited group and 0.31% in the control group,\n",
    "an ITT relative risk of $0.28/0.31 = 0.9$.\n",
    "This led [major news outlets to report](https://www.statnews.com/2022/10/09/in-gold-standard-trial-colonoscopy-fails-to-reduce-rate-of-cancer-deaths/) that colonoscopies do not prevent\n",
    "cancer deaths. \n",
    "A better description might be that _recommending_ colonoscopies does not appear to have a large\n",
    "effect on deaths.\n",
    "In the same study, the 10-year risk of death among subjects who actually received colonoscopies was 0.15%, a \n",
    "PP relative risk of $0.15/0.3 = 0.5$, substantially stronger evidence that _receiving_ a colonoscopy has value.\n",
    "\n",
    "However, for PP it is crucial to know _why_ subjects did not receive their assigned treatment, since self-selection \n",
    "and other processes can lead to substantial biases in estimates of the effect of treatment.\n",
    "For example, imagine a trial of a treatment for a disease.\n",
    "Suppose that treatment is effective, and that subjects who were assigned a placebo and and get worse break\n",
    "protocol and take an active treatment instead.\n",
    "Then the PP analysis will be biased against finding that the treatment is effective.\n",
    "Conversely, suppose that treatment actually makes some people worse, and that those people stop taking\n",
    "the treatment.\n",
    "Then a PP analysis will tend to hide the fact that treatment hurts some subjects.\n",
    "\n",
    "### Censoring and missing data\n",
    "\n",
    "It is common in some kinds of experiments that the outcomes for some subjects cannot be observed, or can be observed\n",
    "only partially.\n",
    "For instance, in experiments on people that involve followup over time, some subjects might simply not show up\n",
    "for one or more followup examinations, resulting in _missing data_.\n",
    "\n",
    "In experiments involving how long something takes, for instance, time to \"failure\" or \"death,\" some subjects might\n",
    "not have experienced \"failure\" by the time the experiment ends.\n",
    "The experimenter knows that such subjects are \"still alive\" when the experiment ended, but they do not\n",
    "know how long it will be until those subjects \"die.\" \n",
    "That is, rather than observing the failure time itself,\n",
    "the experimenter observes that the failure time is greater than some value.\n",
    "This is an example of _(right) censoring_.\n",
    "There are special techniques for dealing with missing data (including subjects lost to followup) and censoring.\n",
    "\n",
    "Below, we assume that compliance is perfect, that no subject drops out of the experiment, and that\n",
    "all outcomes can be observed (no missing data and no censoring)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Null hypotheses for the Neyman model\n",
    "\n",
    "The _strong_ null hypothesis is that subject by subject, the effect of\n",
    "all $T$ treatments is the same.\n",
    "That is,\n",
    "\\begin{equation}\n",
    "x_{j0} = x_{j1} = \\cdots = x_{jT-1}, \\;\\; j=1, \\ldots, N.\n",
    "\\end{equation}\n",
    "Different subjects may have different responses ($x_{jt}$ might not equal $x_{kt}$ if $j \\ne k$), but each subject's response is the same regardless of the treatment assigned \n",
    "to that subject.\n",
    "This is the null hypothesis Fisher considered in _The Design of Experiments_ and which\n",
    "he generally considered the \"correct\" null in practice.\n",
    "\n",
    "The _weak_ null hypothesis is that on average across subjects, all treatments have the same effect. \n",
    "That is,\n",
    "\\begin{equation}\n",
    "\\frac{1}{N} \\sum_{j=1}^N x_{j0} = \\frac{1}{N} \\sum_{j=1}^N x_{j1} = \\ldots = \\frac{1}{N} \\sum_{j=1}^N x_{jT-1}.\n",
    "\\end{equation}\n",
    "Much of Neyman's work on experiments involves this null hypothesis.\n",
    "The statistical theory is more complex for the weak null hypothesis than for the strong null,\n",
    "and most results about the weak null are asymptotic, while the strong null is more amenable to exact or conservative\n",
    "inferences.\n",
    "\n",
    "The strong null is indeed a stronger hypothesis than the weak null, because if the strong null is true, the weak null must also be true: if $T$ lists are equal, element by element, then their means are equal. \n",
    "(However, see Ding, P., 2017. A Paradox from Randomization-Based Causal Inference, _Statist. Sci. 32_, 331-345. 10.1214/16-STS571.)\n",
    "\n",
    "The converse is not true: the weak null can be true even if the strong null is false.\n",
    "For example, for $T=2$ and $N=2$, we might have potential responses $(0, 1)$ for subject 1 and $(1,0)$ for subject 2. \n",
    "The effect of treatment is to increase subject 1's response from 0 to 1 and to decrease subject 2's response from 1 to 0.\n",
    "The treatment affects both subjects, but the average effect of treatment is the same: the average response across subjects is 1/2, with or without treatment.\n",
    "\n",
    "The _effect of treatment $t$ on subject $j$_ (compared to control) is \n",
    "\\begin{equation}\n",
    "\\tau_{jt} := x_{jt} - x_{j0}, \\;\\;t=1, \\ldots, T-1.\n",
    "\\end{equation}\n",
    "The _average effect of treatment $t$_ (compared to control) is\n",
    "\\begin{equation}\n",
    "\\bar{\\tau}_t := \\frac{1}{N} \\sum_{j=1}^N (x_{jt} - x_{j0}), \\;\\; t=1, \\ldots, T-1.\n",
    "\\end{equation}\n",
    "\n",
    "The strong null hypothesis for treatment $t > 0$ is that $\\tau_{jt} = 0$ for all $j = 1, \\ldots, N$.\n",
    "(Treatment $t$ has no effect on any subject in the study.)\n",
    "The \"overall\" or \"grand\" or \"omnibus\" strong null hypothesis is that \n",
    "$\\tau_{jt} = 0$ for all $j = 1, \\ldots, N$ and all $t = 1, \\ldots, T-1$.\n",
    "(No treatment makes any difference to any subject in the study.)\n",
    "Equivalently, the omnibus strong null hypothesis is the intersection of the strong null hypotheses for \n",
    "all the individual treatments.\n",
    "\n",
    "The weak null hypothesis for treatment $t > 0$ is that $\\bar{\\tau}_t = 0$.\n",
    "(On average, treatment $t$ has no effect on the subjects in the study.)\n",
    "The \"overall\" or \"grand\" or \"omnibus\" weak null hypothesis is that \n",
    "$\\bar{\\tau}_t = 0$ for all $t=1, \\ldots, T-1$.\n",
    "(No treatment affects the average response of the subjects in the study.)\n",
    "Equivalently, the omnibus weak null hypothesis is the intersection of the weak null hypotheses for \n",
    "all the individual treatments.\n",
    "\n",
    "\n",
    "In general we have to make assumptions about how treatment affects responses in order to make rigorous inferences about the average treatment effect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing the strong null hypothesis\n",
    "\n",
    "Under the strong null that the treatment makes no difference whatsoever--as if \n",
    "the response had been predetermined before assignment to treatment or control--the null distribution of any test statistic is completely determined once the data have\n",
    "been observed: we know what the data would have been for any other random assignment, namely, the same. \n",
    "And we know the chance that each of those possible datasets would have resulted from\n",
    "the experiment, since we know how subjects were assigned at random to treatments.\n",
    "\n",
    "For alternatives that allow us to find $x_{j0}$ from $x_{jt}$ and vice versa,\n",
    "the alternative also completely determines the \n",
    "probability distribution of any test statistic, once the data have been observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternative hypotheses in the Neyman model\n",
    "\n",
    "#### Constant shift\n",
    "\n",
    "For example, if we assume that the effect of treatment is to shift every subject's response by the same amount, then we can use a test of the strong null to make inferences about that constant effect.\n",
    "In symbols, this alternative states that there are some numbers $\\{\\delta_t\\}_{t=1}^T$\n",
    "such that $x_{jt} = x_{j0}+\\delta_t$ for all subjects $j$.\n",
    "\n",
    "Again, once the original data are observed, this hypothesis completely specifies the probability distribution of the data: we know what subject $j$'s response would have been had the subject been assigned any other treatment. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Other tractable alternative hypotheses\n",
    "\n",
    "Similarly, it is straightforward to test against any fixed alternative set of values\n",
    "of $\\{\\tau_{jt}\\}$ using methods for testing the strong null.\n",
    "It is also possible to test against the alternative $x_{jt} = f_t(x_{j0})$ for some strictly monotonic (and thus invertible) functions $\\{f_t\\}$. \n",
    "A simple example is that treatment multiplies the baseline response by a nonzero constant.\n",
    "\n",
    "In some contexts, it can be reasonable to assume that treatment can only help, that is that $x_{jt} \\ge x_{j0}$ for all $t$, without specifying a functional relationship between them. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary treatments\n",
    "\n",
    "Suppose $T=2$: we are comparing two treatments (typically, an \"active\" treatment and a placebo or other control). \n",
    "Let $\\boldsymbol{1}_N$ be the $N$-dimensional vector of 1s, and let $\\boldsymbol{0}_N$ be the $N$-dimensional zero vector.\n",
    "We will omit the subscript $N$ where the dimension is clear from context.\n",
    "\n",
    "Recall that the _individual treatment effect for subject $j$_ is \n",
    "\\begin{equation}\n",
    "\\tau_j := x_{1j}-x_{0j}.\n",
    "\\end{equation}\n",
    "Let $\\tau := (\\tau_1, \\ldots \\tau_N)$.\n",
    "The _average treatment effect_ for the subjects in the experiment is\n",
    "\\begin{equation}\n",
    "\\bar{\\tau} := \\frac{1}{N} \\sum_{j=1}^N \\tau_j = \\boldsymbol{1} \\cdot \\tau/N = \\frac{1}{N} \\sum_{j=1}^N (x_{1j}-x_{0j}) = \\bar{x}_1 - \\bar{x}_0.\n",
    "\\end{equation}\n",
    "(Whether $\\bar{\\tau}$ says anything about any larger population \n",
    "depends on how the subjects in the experiment were selected--i.e., on whether the experimental subjects are representative \n",
    "of that larger population.)\n",
    "The _strong_ null hypothesis is $\\tau = \\boldsymbol{0}$, i.e.,\n",
    "\\begin{equation}\n",
    "\\tau_1 = \\tau_2 = \\ldots = \\tau_N = 0.\n",
    "\\end{equation}\n",
    "This is a composite hypothesis, because it does not specify the \"baseline\" responses $\\{x_{j0}\\}$ of the subjects,\n",
    "while the distribution of the data depends on $\\{x_{j0}\\}$ in addition to the individual treatment\n",
    "effects $\\{\\tau_{jt}\\}$.\n",
    "The baseline responses are _nuisance parameters_ if the goal is to make inferences about treatment effects.\n",
    "\n",
    "For any vector $\\xi \\in \\mathbb{R}^N$, let $H_\\xi$ denote the hypothesis that $\\tau = \\xi$.\n",
    "For any set $\\Xi \\subset \\mathbb{R}^N$, let $H_\\Xi$ denote the hypothesis that $\\tau \\in \\Xi$.\n",
    "Equivalently, $H_\\Xi = \\cup_{\\xi \\in \\Xi} H_\\xi$.\n",
    "Then the _strong_ null is $H_\\boldsymbol{0}$ and the \n",
    "_weak_ null is $H_{\\Xi_0}$, \n",
    "where\n",
    "\\begin{equation}\n",
    "\\Xi_0 := \\left \\{ \\xi \\in \\mathbb{R}^N : \\boldsymbol{1}_N \\cdot \\xi/N = 0 \\right \\}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference about the average treatment effect\n",
    "Suppose we assign $n$ subjects at random\n",
    "to treatment 0 (_control_ or _placebo_) and the other $m = N-n$ to treatment 1 (_treatment_ or _active treatment_).\n",
    "\n",
    "Then testing the strong null hypothesis is identical to the _two-sample problem_:\n",
    "under the strong null, each subject's response would have been the same, regardless\n",
    "of treatment.\n",
    "If we know the mechanism for allocating subjects to treatments (e.g., by _simple random sampling_, by _Bernoulli sampling_,\n",
    "using a _blocked_ design, etc.) we can calculate or simulate the distribution of any test statistic\n",
    "whatsoever. \n",
    "The question is which test statistic to use.\n",
    "The answer depends on the alternatives against which we want to have power:\n",
    "if we think a particular set of alternatives is most plausible or most important, we can choose the test statistic to\n",
    "maximize some summary of power over that set of alternatives, such as the minimum or a weighted average\n",
    "(the weights can be thought of as a prior probability distribution on alternatives)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now consider the _weak_ null hypothesis and the average treatment effect.\n",
    "\n",
    "Let $W_j$ denote the treatment assigned to subject $j$: $W_j = 0$ if subject $j$ is assigned to control\n",
    "and $W_j = 1$ if subject $j$ is assigned to the active treatment.\n",
    "Assume for now that the assignment is by simple random sampling\n",
    "or by Bernoulli sampling. \n",
    "For Bernoulli sampling, conditioning on the number of subjects assigned to active treatment, $n$, \n",
    "makes the math the same as it is for simple random sampling, where $n$ is fixed in advance.\n",
    "\n",
    "The random variables $\\{W_j\\}$ are identically distributed but not independent (given $n$),\n",
    "because $\\sum_{j=1}^N W_j = n$.\n",
    "Because $\\{W_j\\}$ are identically distributed, they share the same expected value, and\n",
    "thus $\\mathbb{E} \\sum_{j=1}^N W_j = N \\mathbb{E} W_j$ for each $j$.\n",
    "Since $\\sum_{j=1}^N W_j = n$ is constant, $\\mathbb{E} \\sum_{j=1}^N W_j = n$.\n",
    "Hence, $\\mathbb{E} W_j = n/N$, $j = 1, \\ldots, N$.\n",
    "\n",
    "Here is a bit more notation:\n",
    "\\begin{equation}\n",
    "X_0^* := \\sum_{j=1}^N (1-W_j)x_{j0}\n",
    "\\end{equation}\n",
    "is the sum of the responses of the subjects assigned to treatment 0\n",
    "and\n",
    "\\begin{equation}\n",
    "\\bar{X}_0 := X_0^*/m\n",
    "\\end{equation}\n",
    "is their mean response.\n",
    "\\begin{equation}\n",
    "X_1^* := \\sum_{j=1}^N W_j x_{j1} \n",
    "\\end{equation}\n",
    "is the sum of the responses of the subjects assigned to treatment 1 and\n",
    "\\begin{equation}\n",
    "\\bar{X}_1 := X_1^*/n,\n",
    "\\end{equation}\n",
    "is their mean response.\n",
    "\n",
    "The statistics $\\bar{X}_0$ and $\\bar{X}_1$  are unbiased estimates of \n",
    "the corresponding population parameters,\n",
    "\\begin{equation}\n",
    "\\bar{x}_0 := \\frac{1}{N} \\sum_{j=1}^n x_{j0}\n",
    "\\end{equation}\n",
    "and\n",
    "\\begin{equation}\n",
    "\\bar{x}_1 := \\frac{1}{N} \\sum_{j=1}^n x_{j1}.\n",
    "\\end{equation}\n",
    "Therefore, $\\hat{\\tau} := \\bar{X}_1 - \\bar{X}_0$ is an unbiased estimate of the average treatment effect $\\bar{\\tau}$. \n",
    "One could base a permutation test of the strong null using as the test statistic any of $X_0^*$, $X_1^*$,\n",
    "$\\bar{X}_0$, $\\bar{X}_1$, or $\\hat{\\tau}$: all lead to equivalent tests, in the sense that\n",
    "they reject for the same data.\n",
    "\n",
    "See Wu and Ding, 2021. Randomization Tests for Weak Null Hypotheses in Randomized Experiments,\n",
    "_JASA_, _116_. https://www.tandfonline.com/doi/abs/10.1080/01621459.2020.1750415"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two treatments, binary responses.\n",
    " \n",
    "Imagine testing whether a vaccine prevents a disease.\n",
    "We assign a random sample of $n$ of the $N$ subjects to receive treatment 1;\n",
    "the other $m = N-n$ receive a placebo, treatment 0.\n",
    "Let $W_j$ denote the treatment assigned to subject $j$.\n",
    "We observe $W_j$ and\n",
    "\\begin{equation}\n",
    "X_j := (1-W_j) x_{j0} + W_j x_{j1}, \\;\\; j=1, \\ldots, N.\n",
    "\\end{equation}\n",
    "These are random variables, but (in this model)\n",
    "the only source of randomness is $\\{W_j\\}$, the \n",
    "treatment assignment variables.\n",
    "\n",
    "The total number of infections among the vaccinated is\n",
    "\\begin{equation}\n",
    "   X_1^* := \\sum_{j=1}^N W_j x_{j1}\n",
    "\\end{equation}\n",
    "and the total among the unvaccinated is\n",
    "\\begin{equation}\n",
    "   X_0^* := \\sum_{j=1}^N (1-W_j) x_{j0}.\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Under the strong null that the vaccine makes no difference whatsoever--as if whether a subject would become ill was predetermined before the assignment to treatment or control--the distribution of the number of infections among the vaccinated would have a hypergeometric distribution with parameters $N$, $G=X_0^*+X_1^*$, and $n=n$.\n",
    "Testing the strong null using this hypergeometric distribution yields [Fisher's Exact Test](./fisher-exact.ipynb).\n",
    "\n",
    "This approach can be generalized by considering the sample sizes $n$ and $m$ and/or the total number of infections $X_0^*+X_1^*$ to be random, then conditioning on their observed values to get a conditional test.\n",
    "\n",
    "If responses are binary, each $x_{jt}$ is $0$ or $1$, and each $\\tau_j$ is either $-1$, $0$, or $1$.\n",
    "Thus the only possible values of $\\bar{\\tau}$ are multiples of $1/N$.\n",
    "The largest and smallest possible values are $\\bar{\\tau}=-1$ (which\n",
    "occurs if $x_{j1}=0$ and $x_{j0}=1$ for all $j$) and\n",
    "$\\bar{\\tau}=1$ (which occurs if $x_{j1}=1$ and $x_{j0}=0$ for all $j$), a range of 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The core question\n",
    "\n",
    "What can we say about the average treatment effect from the data?\n",
    "\n",
    "The potential responses of the subjects in the experiment can be summarized by four numbers, \n",
    "the number of subjects with each possible combination of potential outcomes.\n",
    "Let $N_{ik}$ be the number of subjects $j$ for whom \n",
    "$x_{j0} = i$ and $x_{j1} = k$, for $i, k \\in \\{0, 1\\}$.\n",
    "That is, $N_{00}$ is the number of subjects whose response is \"0\" whether they are assigned to treatment or to control, while $N_{01}$ is the number of subjects whose response would be\n",
    "\"0\" if assigned to control and \"1\" if assigned to treatment, etc.\n",
    "Of course, $N = N_{00} + N_{01} + N_{10} + N_{11}$.\n",
    "(**Note**: this notation has the indices in the opposite order from that in [Li and Ding (2016)](https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6764).)\n",
    "These four numbers tell us everything relevant about the study population. \n",
    "Let $\\mathbf{N} := (N_{ij})_{i, j \\in \\{0, 1\\}}$ be the 2x2 _summary potential outcome table_.\n",
    "\n",
    "Let $N_{1+} := N_{10} + N_{11}$ be the number of subjects whose response would be \"1\" if assigned to control,\n",
    "and let $N_{+1} := N_{01} + N_{11}$ be the number of subjects whose response would be \"1\" if \n",
    "assigned to treatment.\n",
    "(**Note**: this notation also has the indices in the opposite order from that in \n",
    "[Li and Ding (2016)](https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6764).)\n",
    "\n",
    "Now $\\sum_{j=1}^N x_{j0} = N_{1+}$ and $\\sum_{j=1}^N x_{j1} = N_{+1}$, so\n",
    "the average treatment effect can be written\n",
    "\\begin{equation}\n",
    "\\bar{\\tau} = \\frac{1}{N} \\sum_{j=1}^N (x_{j1}-x_{j0}) = \\frac{1}{N} (N_{+1} - N_{1+}),\n",
    "\\end{equation}\n",
    "the average response among the treated minus the average among the controls."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the data have been observed, we know that $N_{1+}$ is at least $X_0^*$\n",
    "(it was that big for the control group alone, even if none of the $n$ treated subjects would have\n",
    "had a control response of 1)\n",
    "and at most $X_0^* + n$ (if all $n$ of the treated subjects would have had a response of 1 if they\n",
    "had been assigned to control instead of the active treatment).\n",
    "Similarly, we know that $N_{+1}$ is at least $X_1^*$\n",
    "(we saw that many 1s in the active treatment group)\n",
    "and at most $X_1^* + m$ (if every subject assigned to control would have responded 1 if they had received the\n",
    "active treatment).\n",
    "Their difference is thus between\n",
    "\\begin{equation}\n",
    "X_1^* - X_0^* - n \n",
    "\\end{equation}\n",
    "and\n",
    "\\begin{equation}\n",
    "X_1^* - X_0^* + m,\n",
    "\\end{equation}\n",
    "a range of $n+m = N$.\n",
    "Thus\n",
    "\\begin{equation}\n",
    "\\bar{\\tau} \\in \\{ (X_1^* - X_0^* - n)/N, (X_1^* - X_0^* - n + 1)/N, \\ldots, (X_1^* - X_0^* + m)/N\\},\n",
    "\\end{equation}\n",
    "which has range $1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First (conservative) approach to confidence bounds for $\\bar{\\tau}$: Bonferroni simultaneous confidence intervals\n",
    "\n",
    "A tuple of confidence set procedures $(\\mathcal{I}_i )_{i=1}^C$ for a corresponding tuple of parameters $(\\theta_i)_{i=1}^C \\in \\Theta$ has _simultaneous coverage probability_ $1-\\alpha$ if\n",
    "\\begin{equation}\n",
    "\\mathbb{P}_{\\theta_1, \\ldots, \\theta_C} \\bigcap_{i=1}^C ( \\mathcal{I}_i \\ni \\theta_i  ) \\ge 1-\\alpha,\n",
    "\\end{equation}\n",
    "whatever the true values $(\\theta_i) \\in \\Theta$ are.\n",
    "\n",
    "Bonferroni's inequality, also called _the union bound_, says that for any collection of events $\\{A_i\\}$,\n",
    "\\begin{equation}\n",
    "\\mathbb{P} \\left ( \\bigcup_i A_i \\right ) \\le \\sum_i \\mathbb{P}(A_i).\n",
    "\\end{equation}\n",
    "It follows that if $\\mathcal{I}_i$ is a confidence interval procedure for $\\theta_i$ with coverage probability $1-\\alpha_i$, $i=1, \\ldots, C$, then the tuple of $C$\n",
    "confidence set procedures $(\\mathcal{I}_i)_{i=1}^C$ has simultaneous coverage probability\n",
    "not smaller than $1-\\sum_{i=1}^C \\alpha_i$.\n",
    "\n",
    "In the current situation, if we had simultaneous confidence sets for both $N_{1+}$ and $N_{+1}$, we could find a confidence set for $\\bar{\\tau}$, because \n",
    "$\\bar{\\tau} = (N_{+1} - N_{1+})/N$.\n",
    "\n",
    "Now, $X_1^*$ is the number of \"1\"s in a random sample of size $n$ from a population of size $N$ of which $N_{+1}$ are labeled \"1\" and the rest are labeled \"0.\"\n",
    "Similarly, $X_0^*$ is the number of \"1\"s in a random sample of size $m$ from a population of size $N$ of which $N_{1+}$ are labeled \"1\" and the rest are labeled \"0.\"\n",
    "That is, the probability distribution of $X_1^*$ is \n",
    "hypergeometric with parameters $N$, $G=N_{+1}$, and $n$,\n",
    "and the probability distribution of $X_0^*$ is \n",
    "hypergeometric with parameters $N$, $G=N_{1+}$, and $m$.\n",
    "We can thus use [methods for finding confidence bounds for the hypergeometric $G$](./confidence-sets.ipynb)\n",
    "to find confidence bounds for $N_{+1}$ and $N_{1+}$ from $X_1^*$ and $X_0^*$.\n",
    "The variables $X_1^*$ and $X_0^*$ are dependent, but Bonferroni's inequality holds for dependent events.\n",
    "\n",
    "To construct an _upper_ 1-sided $1-\\alpha$ confidence bound for $\\bar{\\tau}$, we can find an _upper_ 1-sided $1-\\alpha/2$ confidence bound for $N_{+1}$ (using the hypergeometric distribution), subtract a _lower_ \n",
    "1-sided $1-\\alpha/2$ confidence bound for $N_{1+}$, and divide the result by $N$.\n",
    "\n",
    "To construct a _lower_ 1-sided $1-\\alpha$ confidence bound for $\\bar{\\tau}$, we can find a _lower_ 1-sided $1-\\alpha/2$ confidence bound for $N_{+1}$, subtract an _upper_ 1-sided $1-\\alpha/2$ confidence bound for $N_{1+}$, and divide the result by $N$.\n",
    "\n",
    "To construct a 2-sided confidence interval for $\\bar{\\tau}$, we can find a \n",
    "2-sided $1-\\alpha/2$ confidence bound for $N_{+1}$ and a 2-sided $1-\\alpha/2$ confidence bound for $N_{1+}$. \n",
    "The lower endpoint of the $1-\\alpha$ confidence interval for $\\bar{\\tau}$ \n",
    "is the lower endpoint of the 2-sided interval for $N_{+1}$  minus the upper\n",
    "endpoint of the 2-sided interval for $N_{1+}$, divided by $N$.\n",
    "The upper endpoint of the $1-\\alpha$ confidence interval for $\\bar{\\tau}$ \n",
    "is the upper endpoint of the 2-sided interval for $N_{+1}$  minus the lower\n",
    "endpoint of the 2-sided interval for $N_{1+}$, divided by $N$.\n",
    "\n",
    "This approach has a number of advantages: it is conceptually simple, conservative, and\n",
    "only requires the ability to compute confidence intervals for $G$ for\n",
    "hypergeometric distributions.\n",
    "But the intervals will in general be quite conservative, i.e., unnecessarily wide. \n",
    "(See Table 1 of [Li and Ding (2016)](https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6764).)\n",
    "We now examine a sharper approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Second approach: testing all tables of potential outcomes\n",
    "\n",
    "After the randomization, for each subject $j$, we observe either $x_{j0}$ or\n",
    "$x_{j1}$.\n",
    "At that point, we know $N$ of the $2N$ numbers $\\{x_{jt}\\}_{j=1}^N{}_{t=0}^1$.\n",
    "The other $N$ numbers--the responses that were not observed--can be any combination of 0s and 1s: there are $2^N$ possibilities in all.\n",
    "But not all of those yield distinguishable sets of responses: what matters is how many subjects\n",
    "have each possible pattern of potential outcomes: (0, 0), (0, 1), (1, 0), and (1, 1).\n",
    "\n",
    "The average treatment effect $\\bar{\\tau}$ is determined by\n",
    "the summary potential outcome table $\\mathbf{N}$, i.e., the\n",
    "four numbers, $N_{00}$, $N_{01}$, $N_{10}$, and $N_{11}$, which to sum to $N$. \n",
    "How many possible values are there for those four numbers?\n",
    "The total number of ways there are of partitioning $N$ items into 4 groups can be found by Feller's \"bars and stars\" argument (see the [notes on nuisance parameters](./nuisance.ipynb));\n",
    "the answer is $\\binom{N+3}{3} = (N+3)(N+2)(N+1)/6$. \n",
    "This is $O(N^3)$ tables\n",
    "(see [Rigdon and Hudgens (2015)](https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6384)).\n",
    "\n",
    "But many of those tables are incompatible with the observed data.\n",
    "For instance, we know that $ N_{1+} \\ge X_0^*$ and $N_{+1} \\ge X_1^*$.\n",
    "[Li and Ding (2016)](https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6764) show that taking into account \n",
    "the observed data reduces the number of tables that must be considered to $O(N^2)$, greatly speeding the computation.\n",
    "\n",
    "More recent work cuts the number of tables from $O(N^2)$ to $O(N \\ln N)$\n",
    "(Aronow, P.M., H. Chang, and P. Lopatto, 2022?. Fast computation of exact confidence intervals for randomized experiments with binary outcomes. https://lopat.to/permutation.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A permutation approach**\n",
    "\n",
    "Here is the Li and Ding (2016) approach.\n",
    "Together, $N_{00}$, $N_{01}$, $N_{10}$, and $N_{11}$ determine the sampling distribution of any statistic, \n",
    "through the random allocation of subjects to treatments.\n",
    "(This is true whether the allocation to treatment is by simple random sampling or Bernoulli sampling;\n",
    "for a blocked design, the distribution also depends on how many subjects of each kind are in each block.)\n",
    "To test the null hypothesis $H_0: \\mathbf{N} = \\mathbf{N}_0$, \n",
    "we can define some function $T$ of $X_0^*$ and $X_1^*$ and reject $H_0$ if the observed value of $T$ is in the \n",
    "tail of the probability distribution corresponding to $H_0$.\n",
    "\n",
    "What should we use for $T$? Since we are interested in $\\bar{\\tau}$, we might base the test on\n",
    "\\begin{equation}\n",
    "\\hat{\\tau} := \\bar{X}_1 - \\bar{X}_0 = \\frac{X_1^*}{n} - \\frac{X_0^*}{m}.\n",
    "\\end{equation}\n",
    "This is an unbiased estimator of $\\bar{\\tau}$: $X_1^*/n$ is the sample mean of a simple random sample of size $n$ from $\\{x_{j1}\\}_{j=1}^N$, so it is unbiased for\n",
    "$N_{+1}/N$ and \n",
    "$X_0^*/m$ is the sample mean of a simple random sample of size $m$ from $\\{x_{j0}\\}_{j=1}^N$, so it is unbiased for\n",
    "$N_{1+}/N$. \n",
    "Their difference is thus unbiased for $N_{+1}/N - N_{1+}/N = \\bar{\\tau}$.\n",
    "\n",
    "Rigdon and Hudgens (2015)  and Li and Ding (2016) take the test statistic to be $T=|\\hat{\\tau} - \\bar{\\tau}|$ (they look at other things, too)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that $N_{1+} := N_{10} + N_{11}$, $N_{+1} := N_{01} + N_{11}$,\n",
    "and \n",
    "\\begin{equation}\n",
    "\\bar{\\tau} = \\frac{N_{+1} - N_{1+}}{N} = \\frac{N_{01}-N_{10}}{N}.\n",
    "\\end{equation}\n",
    "For any 2x2 summary table $\\mathbf{M}$ of counts of potential outcomes, define\n",
    "\\begin{equation}\n",
    "\\bar{\\tau}(\\mathbf{M}) := \\frac{M_{01} - M_{10}}{\\sum_{j,k} M_{jk}}.\n",
    "\\end{equation}\n",
    "If we do not reject the simple hypothesis $\\mathbf{N} = \\mathbf{N}_0$, then we cannot reject the hypothesis $\\bar{\\tau} = \\bar{\\tau}(\\mathbf{N}_0) := \\bar{\\tau}_0$.\n",
    "We can reject the composite hypothesis $\\bar{\\tau} = \\bar{\\tau}_0$ if we can reject the simple hypothesis $\\mathbf{N} = \\mathbf{M}$\n",
    "for _every_ summary table $\\mathbf{M}$ for which $\\bar{\\tau}(\\mathbf{M}) = \\bar{\\tau}_0$.\n",
    "\n",
    "We can calibrate the test of the hypothesis $\\mathbf{N} = \\mathbf{N}_0$\n",
    "by finding out the probability distribution of $T$ under $H_0$:\n",
    "Given the summary potential outcome table $\\mathbf{N}_0$ (the counts of subjects with each possible combination of potential outcomes), \n",
    "we can construct a full table of potential outcomes that is consistent\n",
    "with that table.\n",
    "For each of the $\\binom{N}{n}$ equally likely possible treatment assignments, we can find\n",
    "the corresponding value of $T$ by allocating the subjects accordingly, finding $\\hat{\\tau}$, subtracting $\\bar{\\tau}(\\mathbf{N}_0)$, and\n",
    "taking the absolute value.\n",
    "\n",
    "When $N$ is large and $n$ is not close to $0$ or $N$, it is impractical to construct all $\\binom{N}{n}$ possible treatment assignments.\n",
    "Then, we might approximate the probability distribution by simulation: actually drawing\n",
    "pseudo-random samples of size $n$ from the subjects, considering them to be the treatment group, calculating $\\hat{\\tau}$, etc.\n",
    "As discussed in the chapter on [testing](./tests.ipynb), that can be viewed as an approximation to a theoretical $P$-value or as\n",
    "the exact $P$-value for a randomized test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Enumerating all 2x2 tables arithmetically consistent with the data**\n",
    "\n",
    "The table $\\mathbf{N}$ summarizing potential outcomes is constrained by the data. \n",
    "For instance, $N_{10} + N_{11} \\ge X_0^*$ and $N_{01} + N_{11} \\ge X_1^*$.\n",
    "There are other constraints, too, as we shall see.\n",
    "\n",
    "Define \n",
    "\\begin{equation}\n",
    "  n_{wk} := \\sum_{j=1}^N 1_{W_j=w, x_{jw}=k}\n",
    "\\end{equation}\n",
    "for $w \\in \\{0, 1\\}$ and $k \\in \\{0, 1\\}$.\n",
    "That is, $n_{00} = m-X_0^*$, $n_{01} = X_0^*$, $n_{10} = n-X_1^*$, and\n",
    "$n_{11} = X_1^*$.\n",
    "Clearly, $\\sum_{w,k} n_{wk} = N$: these numbers count the elements in \n",
    "a partition of the subjects.\n",
    "Using this notation, $\\hat{\\tau} = n_{11}/n - n_{01}/m$.\n",
    "Let $\\mathbf{n} := (n_{wk})_{w, k \\in \\{0, 1\\}}$ be these four numbers as a matrix, the _observed outcome\n",
    "matrix_, the number of subjects assigned to treatment $w$ whose response was $k$, for $w, k \\in \\{0, 1\\}$.\n",
    "\n",
    "[Rigdon and Hudgens (2014)](https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6384) show that it suffices to examine \n",
    "$n_{RH} := (n_{11} + 1)(n_{10} + 1)(n_{01} + 1)(n_{00} + 1)$ 2x2 tables to find confidence bounds\n",
    "for $\\bar{\\tau}$ using $\\hat{\\tau}$ as the test statistic.\n",
    "\n",
    "Their argument is as follows: \n",
    "Consider the $n_{11}$ subjects who were assigned to treatment $w=1$ and whose response was $x_{j1}=1$.\n",
    "Fix the unobserved values of the remaining $N-n_{11}$ subjects.\n",
    "As the unobserved responses of those $n_{11}$ vary, the value of $\\bar{\\tau}$ \n",
    "and the probability distribution of $T$ depend only on how many of them\n",
    "have $x_{j0}=0$ and how many have $x_{j0}=1$.\n",
    "The number $m$ for which $x_{j0}=0$ could be 0, 1, $\\ldots$, or $n_{11}$;\n",
    "the other $n_{11}-m$ have $x_{j0}=1$. \n",
    "There are thus only $n_{11}+1$ tables that need to be examined,\n",
    "given the unobserved values in the other three groups.\n",
    "By the same argument, as the unobserved values of the $n_{01}$\n",
    "subjects vary, holding constant the unobserved values for the rest, there are at\n",
    "most $n_{01}+1$ distinct values of $\\bar{\\tau}$ and distinct probability distributions for\n",
    "$T$. \n",
    "The same goes for $n_{10}$ and $n_{11}$. \n",
    "By the _fundamental rule of counting_,\n",
    "the total number of tables that give rise to distinguishable probability distributions\n",
    "for $T$ is thus at most $n_{RH}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Li and Ding (2016)](https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6764) improve on that result.\n",
    "They prove that\n",
    "a table $\\mathbf{N} = (N_{00}, N_{01}, N_{10}, N_{11})$ is consistent with the observed values $(n_{wk})$ iff\n",
    "\\begin{equation}\n",
    "\\max \\{0,n_{11}-N_{01}, N_{11}-n_{01}, N_{10}+N_{11} - n_{10} - n_{01} \\}\n",
    "  \\le\n",
    "\\min \\{N_{11}, n_{11}, N_{10}+N_{11}-n_{01}, N - N_{01} - n_{01} - n_{10} \\}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "**The algorithm**\n",
    "\n",
    "The overall approach of Li and Ding (2016) is as follows:\n",
    "\n",
    "The values $N$, $(n_{wk})$, and $\\alpha$ are given. \n",
    "Set\n",
    "\\begin{equation}\n",
    "\\hat{\\tau}^* := \\frac{n_{11}}{n} - \\frac{n_{01}}{m},\n",
    "\\end{equation}\n",
    "the observed value of the unbiased estimate of $\\bar{\\tau}$.\n",
    "\n",
    "+ for each table $\\mathbf{N}$ that is algebraically consistent with the observed values $(n_{wk})$:\n",
    "    - find $\\bar{\\tau}(\\mathbf{N}) = \\frac{N_{01}-N_{10}}{N}$.\n",
    "    - calculate $t = |\\hat{\\tau}^* - \\bar{\\tau}(\\mathbf{N})|$\n",
    "    - create a \"full\" list of potential outcomes $(x_{ik})$ consistent with the summary table $\\mathbf{N}$.\n",
    "    - find or simulate the sampling distribution of $|\\hat{\\tau} - \\bar{\\tau}|$ using $(x_{ik})$\n",
    "    - if $t$ is above the $1-\\alpha$ percentile of the sampling distribution, pass; otherwise include $\\bar{\\tau}(\\mathbf{N})$ in the confidence set\n",
    "+ report the smallest and largest values of $\\bar{\\tau}$ in the confidence set as the endpoints of the confidence interval\n",
    "    \n",
    "Below, the main step--generating tables that are algebraically consistent with the data--is implemented in python\n",
    "(`N_generator()`), as is the step of\n",
    "creating a table of potential outcomes consistent with $\\mathbf{N}$ (`potential_outcomes()`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from itertools import filterfalse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N_generator(N: int, n00: int, n01: int, n10: int, n11: int) -> tuple:\n",
    "    ''' \n",
    "    generate tables algebraically consistent with data from an experiment with binary outcomes\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    N : int\n",
    "        number of subjects\n",
    "    n00 : int\n",
    "        number of subjects assigned to treatment 0 who had outcome 0\n",
    "    n01 : int\n",
    "        number of subjects assigned to treatment 0 who had outcome 0\n",
    "    n10 : int\n",
    "        number of subjects assigned to treatment 1 who had outcome 0\n",
    "    n11 : int\n",
    "        number of subjects assigned to treatment 1 who had outcome 1\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Nt : list of 4 ints \n",
    "        N00, subjects with potential outcome 0 under treatments 0 and 1\n",
    "        N01, subjects with potential outcome 0 under treatment 0 and 1 under treatment 1\n",
    "        N10, subjects with potential outcome 1 under treatment 0 and 0 under treatment 1\n",
    "        N11, subjects with potential outcome 1 under treatments 0 and 1\n",
    "    '''\n",
    "    for i in range(min(N-n00, N-n10)+1):               # allocate space for the observed 0 outcomes, n00 and n10\n",
    "        N11 = i                                           \n",
    "        for j in range(max(0, n01-N11), N-n00-N11+1):    # N11+N10 >= n01; N11+N10+n00 <= N\n",
    "            N10 = j                                        \n",
    "            for k in range(max(0, n11-N11), min(N-n10-N11, N-N11-N10)+1): \n",
    "                                                       # N11+N01 >= n11; N11+N01+n10 <= N; no more than N subjects\n",
    "                N01 = k                                  \n",
    "                N00 = N-N11-N10-N01                  \n",
    "                if filter_table([N00, N01, N10, N11], n00, n01, n10, n11):\n",
    "                    yield [N00, N01, N10, N11]\n",
    "                else:\n",
    "                    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_table(Nt: list, n00: int, n01: int, n10: int, n11: int) -> bool:\n",
    "    '''\n",
    "    check whether summary table Nt of binary outcomes is consistent with observed counts\n",
    "    \n",
    "    implements the test in Theorem 1 of Li and Ding (2016):\n",
    "    \n",
    "       \\max \\{0,n_{11}-N_{01}, N_{11}-n_{01}, N_{10}+N_{11}-n_{10}-n_{01} \\}\n",
    "           \\le\n",
    "       \\min \\{N_{11}, n_{11}, N_{10}+N_{11}-n_{01}, N-N_{01}-n_{01}-n_{10} \\}\n",
    "    \n",
    "    Parameters:\n",
    "    ----------\n",
    "    Nt : list of four ints\n",
    "        the table of counts of subjects with each combination of potential outcomes, in the order\n",
    "        N_00, N_01, N_10, N_11\n",
    "       \n",
    "    n01 : int\n",
    "        number of subjects assigned to control whose observed response was 1\n",
    "\n",
    "    n11 : int\n",
    "        number of subjects assigned to treatment whose observed response was 1\n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    ok : bool\n",
    "        True if table is consistent with the data\n",
    "    '''\n",
    "    N = np.sum(Nt)   # total subjects\n",
    "    '''\n",
    "    '''\n",
    "    return max(0,n11-Nt[1], Nt[3]-n01, Nt[2]+Nt[3]-n10-n01) <= min(Nt[3], n11, Nt[2]+Nt[3]-n01, N-Nt[1]-n01-n10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def potential_outcomes(Nt: list) -> np.array:\n",
    "    '''\n",
    "    make a 2xN table of potential outcomes from the 2x2 summary table Nt\n",
    "    \n",
    "    Parameters\n",
    "    ----------   \n",
    "    Nt : list of 4 ints\n",
    "        N00, N01, N10, N11\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    po : Nx2 table of potential outcomes consistent with Nt\n",
    "    '''\n",
    "    return np.reshape(np.array([0,0]*Nt[0]+[0,1]*Nt[1]+[1,0]*Nt[2]+[1,1]*Nt[3]), [-1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0],\n",
       "       [0, 0],\n",
       "       [0, 0],\n",
       "       [0, 0],\n",
       "       [0, 0],\n",
       "       [0, 1],\n",
       "       [0, 1],\n",
       "       [0, 1],\n",
       "       [0, 1],\n",
       "       [1, 0],\n",
       "       [1, 0],\n",
       "       [1, 0],\n",
       "       [1, 1],\n",
       "       [1, 1]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test\n",
    "Nt = [5, 4, 3, 2]\n",
    "potential_outcomes(Nt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hypothetical experiment\n",
    "N = 10\n",
    "n = 5\n",
    "n00 = 3\n",
    "n01 = 2\n",
    "n10 = 1\n",
    "n11 = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[4, 4, 2, 0],\n",
       " [3, 5, 2, 0],\n",
       " [2, 6, 2, 0],\n",
       " [1, 7, 2, 0],\n",
       " [3, 4, 3, 0],\n",
       " [2, 5, 3, 0],\n",
       " [1, 6, 3, 0],\n",
       " [0, 7, 3, 0],\n",
       " [4, 4, 1, 1],\n",
       " [3, 5, 1, 1],\n",
       " [2, 6, 1, 1],\n",
       " [1, 7, 1, 1],\n",
       " [4, 3, 2, 1],\n",
       " [3, 4, 2, 1],\n",
       " [2, 5, 2, 1],\n",
       " [1, 6, 2, 1],\n",
       " [0, 7, 2, 1],\n",
       " [3, 3, 3, 1],\n",
       " [2, 4, 3, 1],\n",
       " [1, 5, 3, 1],\n",
       " [0, 6, 3, 1],\n",
       " [4, 4, 0, 2],\n",
       " [3, 5, 0, 2],\n",
       " [2, 6, 0, 2],\n",
       " [1, 7, 0, 2],\n",
       " [4, 3, 1, 2],\n",
       " [3, 4, 1, 2],\n",
       " [2, 5, 1, 2],\n",
       " [1, 6, 1, 2],\n",
       " [0, 7, 1, 2],\n",
       " [4, 2, 2, 2],\n",
       " [3, 3, 2, 2],\n",
       " [2, 4, 2, 2],\n",
       " [1, 5, 2, 2],\n",
       " [0, 6, 2, 2],\n",
       " [3, 2, 3, 2],\n",
       " [2, 3, 3, 2],\n",
       " [1, 4, 3, 2],\n",
       " [0, 5, 3, 2],\n",
       " [4, 3, 0, 3],\n",
       " [3, 4, 0, 3],\n",
       " [2, 5, 0, 3],\n",
       " [1, 6, 0, 3],\n",
       " [4, 2, 1, 3],\n",
       " [3, 3, 1, 3],\n",
       " [2, 4, 1, 3],\n",
       " [1, 5, 1, 3],\n",
       " [0, 6, 1, 3],\n",
       " [4, 1, 2, 3],\n",
       " [3, 2, 2, 3],\n",
       " [2, 3, 2, 3],\n",
       " [1, 4, 2, 3],\n",
       " [0, 5, 2, 3],\n",
       " [3, 1, 3, 3],\n",
       " [2, 2, 3, 3],\n",
       " [1, 3, 3, 3],\n",
       " [0, 4, 3, 3],\n",
       " [4, 2, 0, 4],\n",
       " [3, 3, 0, 4],\n",
       " [2, 4, 0, 4],\n",
       " [1, 5, 0, 4],\n",
       " [4, 1, 1, 4],\n",
       " [3, 2, 1, 4],\n",
       " [2, 3, 1, 4],\n",
       " [1, 4, 1, 4],\n",
       " [0, 5, 1, 4],\n",
       " [4, 0, 2, 4],\n",
       " [3, 1, 2, 4],\n",
       " [2, 2, 2, 4],\n",
       " [1, 3, 2, 4],\n",
       " [0, 4, 2, 4],\n",
       " [3, 0, 3, 4],\n",
       " [2, 1, 3, 4],\n",
       " [1, 2, 3, 4],\n",
       " [0, 3, 3, 4],\n",
       " [4, 1, 0, 5],\n",
       " [3, 2, 0, 5],\n",
       " [2, 3, 0, 5],\n",
       " [1, 4, 0, 5],\n",
       " [4, 0, 1, 5],\n",
       " [3, 1, 1, 5],\n",
       " [2, 2, 1, 5],\n",
       " [1, 3, 1, 5],\n",
       " [0, 4, 1, 5],\n",
       " [3, 0, 2, 5],\n",
       " [2, 1, 2, 5],\n",
       " [1, 2, 2, 5],\n",
       " [0, 3, 2, 5],\n",
       " [4, 0, 0, 6],\n",
       " [3, 1, 0, 6],\n",
       " [2, 2, 0, 6],\n",
       " [1, 3, 0, 6],\n",
       " [3, 0, 1, 6],\n",
       " [2, 1, 1, 6],\n",
       " [1, 2, 1, 6],\n",
       " [0, 3, 1, 6]]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[Nt for Nt in N_generator(N, n00, n01, n10, n11)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Regeneron data from \n",
    "# https://investor.regeneron.com/news-releases/news-release-details/phase-3-prevention-trial-showed-81-reduced-risk-symptomatic-sars\n",
    "n=753\n",
    "m=752\n",
    "N=n+m\n",
    "n01 = 59\n",
    "n11 = 11\n",
    "n00 = m-n01\n",
    "n10 = n-n11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other measures of the effect of treatment\n",
    "\n",
    "The average treatment effect is one of many ways of quantifying the effect of treatment. \n",
    "For vaccines, a common measure of effect size is _vaccine efficacy_, defined as\n",
    "\\begin{equation}\n",
    "\\mbox{VE} := \\frac{\\mbox{risk among unvaccinated} - \\mbox{risk among vaccinated}}{\\mbox{risk among unvaccinated}}\n",
    "= 1 - \\frac{\\mbox{risk among vaccinated}}{\\mbox{risk among unvaccinated}}.\n",
    "\\end{equation}\n",
    "This is the naive estimate of the fraction of infections that vaccinating everyone would prevent.\n",
    "\n",
    "Here, _risk_ is the fraction of subjects who become infected. \n",
    "In the notation we've been using,\n",
    "\\begin{equation}\n",
    "\\mbox{VE}(\\mathbf{N}) = \\frac{N_{1+}/N - N_{+1}/N}{N_{1+}/N} = \\frac{N_{1+} - N_{+1}}{N_{1+}} = 1 - \\frac{ N_{+1}}{N_{1+}}\n",
    "\\end{equation}\n",
    "if $N_{1+} > 0$.\n",
    "\n",
    "We can make confidence intervals for vaccine efficiacy in the same way as for the average treatment effect:\n",
    "\n",
    "+ Choose a reasonable test statistic $T$. Here, we might use the _plug-in estimator_ of VE, which is biased, but\n",
    "still a sensible choice: \n",
    "\\begin{equation}\n",
    "\\widehat{\\mbox{VE}} := \\left \\{ \\begin{array}{ll}\n",
    "    1 - \\bar{X}_1^*/\\bar{X}_0^*, & \\bar{X}_0^* > 0 \\\\\n",
    "    0, & \\bar{X}_0^* = 0.\n",
    "    \\end{array}\n",
    "    \\right .\n",
    "\\end{equation}\n",
    "+ For each summary table of potential outcomes, $\\mathbf{N}$, test the hypothesis that the observations came from that table using the randomization distribution of $T$: reject if the observed value of $T$ exceeds the $1-\\alpha$ percentile of that distribution.\n",
    "+ If the hypothesis is not rejected, include $\\mbox{VE}(\\mathbf{N})$ in the confidence set.\n",
    "\n",
    "A similar approach works with any measure of effect size in a randomized experiment with binary outcomes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What's special about binary outcomes?\n",
    "\n",
    "The approach above uses the assumption that the potential outcomes can be only 0 or 1.\n",
    "If the potential outcomes take values in a known finite set, a similar approach works, but if there are more than 2 possible values, constructing all tables algebraically consistent with the data is in general harder, and the number of such tables will be larger. \n",
    "\n",
    "If the potential outcomes are bounded with known bounds (but not necessarily discrete with known support), the approach can be generalized using nonparametric methods for estimating the mean of bounded, finite populations (see other chapters in the notes for such methods, including [Confidence bounds via the Chebychev and Hoeffding Inequalities](./hoeffding.ipynb), [Lower confidence bounds for the mean via Markov's Inequality and methods based on the Empirical Distribution](./markov.ipynb), [Penny sampling](./pennySampling.ipynb), [Wald's Sequential Probability Ratio Test (SPRT)](./sprt.ipynb), [Kaplan-Wald Confidence Bound for a Nonnegative Mean](./kaplanWald.ipynb), and [Method shootout](./shootout.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests for bounded treatment effects\n",
    "\n",
    "This section is based on \n",
    "Caughey, D., A. Dafoe, X. Li, and L. Miratrix, 2021. Randomization Inference beyond the Sharp Null: Bounded Null Hypotheses and Quantiles of Individual Treatment Effects https://arxiv.org/abs/2101.09195.\n",
    "\n",
    "That paper shows that for testing based on a difference in means, the permutation test based on the strong null in fact gives valid tests of \"bounded\" nulls.\n",
    "For instance, the usual one-sided test against the alternative that treatment increases response is valid for the null hypothesis that treatment has a non-positive effect for every subject.\n",
    "And a confidence interval for the shift alternative, based on the assumption that all treatment effects are equal, is in fact a confidence interval for the maximum treatment effect even when the effects are not\n",
    "necessarily equal.\n",
    "\n",
    "The basic idea is simple: recall that the treatment effect for the $j$th subject is $\\tau_j := x_{j1}-x_{j0}$.\n",
    "The key observation is that if the test statistic has a particular kind of monotonicity, then if we reject the hypothesis $H_\\delta: \\tau_j = \\delta_j \\;\\forall j$, the test would also reject the hypothesis $H_{\\ge \\delta}: \\tau_j \\ge \\delta_j \\;\\forall j$.\n",
    "A hypothesis of the form $H_{\\ge \\delta}: \\tau_j \\ge \\delta_j \\;\\forall j$ or\n",
    "$H_{\\le \\delta}: \\tau_j \\le \\delta_j \\;\\forall j$ is called a _bounded null hypothesis_.\n",
    "\n",
    "Suppose the ordered treatment\n",
    "effects are $\\tau_{(1)} \\le \\tau_{(2)} \\le \\cdots \\le \\tau_{(N)}$.\n",
    "Caughey et al. show how to test the hypotheses $\\tau_{(k)} \\le c$ for any $k \\in \\{1, \\ldots, N\\}$ and any $c \\in \\Re$, and to find simultaneous confidence intervals for $\\{\\tau_{(k)}\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As above, let $W_j$ be the treatment assigned to the $j$th subject. \n",
    "Let $\\mathbf{W}$ denote the $N$-vector of treatments.\n",
    "A method of assigning treatments is _exchangeable_ iff \n",
    "$\\mathbb{P} \\{\\mathbf{W} = \\mathbf{w}\\} = \\mathbb{P} \\{\\mathbf{W} = \\mathbf{w}'\\}$ for every permutation $\\mathbf{w}'$ of the elements of $\\mathbf{w}$. \n",
    "Assigning a simple random sample of subjects to treatment is exchangeable; so is Bernoulli assignment (each subject is assigned to active treatment or control by independent toss of a biased coin, with the\n",
    "same chance for each subject).\n",
    "\n",
    "We observe $\\mathbf{W} := (W_j)_{j=1}^N$ and $\\mathbf{X} := (X_j)_{j=1}^N$, where\n",
    "\\begin{equation}\n",
    "X_j := (1-W_j) x_{j0} + W_j x_{j1}, \\;\\; j=1, \\ldots, N.\n",
    "\\end{equation}\n",
    "Let $\\mathbf{x}_0 := (x_{j0})_{j=1}^N$ be the potential outcomes under control, and\n",
    "let $\\mathbf{x}_1 := (x_{j1})_{j=1}^N$ be the potential outcomes under active treatment.\n",
    "Let $\\mathbf{1}$ be the $N$-vector of ones.\n",
    "Let $\\circ$ denote elementwise multiplication of vectors, so (for example) we can write\n",
    "\\begin{equation}\n",
    "\\mathbf{X} = (\\mathbf{1} - \\mathbf{W}) \\circ \\mathbf{x}_0 + \\mathbf{W} \\circ \\mathbf{x}_1.\n",
    "\\end{equation}\n",
    "Let $\\mathbf{\\tau} := (\\tau_j)_{j=1}^N$ and $\\mathbf{\\delta} := (\\delta_j)_{j=1}^N$\n",
    "If $\\tau = \\delta$, then $\\mathbf{x}_1 = \\mathbf{x}_0 + \\delta$,\n",
    "so the implied set of potential outcomes under the control is\n",
    "\\begin{equation}\n",
    "\\mathbf{x}_0^{\\mathbf{W},\\delta} = \\mathbf{X} - \\mathbf{W} \\circ \\mathbf{\\delta}\n",
    "\\end{equation}\n",
    "and the implied set of potential outcomes under active treatment is \n",
    "\\begin{equation}\n",
    "\\mathbf{x}_1^{\\mathbf{W},\\delta} = \\mathbf{X} + (\\mathbf{1} - \\mathbf{W}) \\circ \\mathbf{\\delta}.\n",
    "\\end{equation}\n",
    "\n",
    "To test the null $H_\\delta$, we will use test statistics of the form\n",
    "$t(\\mathbf{W}, \\mathbf{x}_0^{\\mathbf{W},\\delta})$, and reject for large values of $t$.\n",
    "\n",
    "The randomization $P$-value for the hypothesis $\\tau = \\delta$ is\n",
    "\\begin{equation}\n",
    "P_{\\mathbf{w},\\delta} := \\mathbb{P}_\\delta \\{ t(\\mathbf{W}, \\mathbf{x}_0^{\\mathbf{w},\\delta}) \\ge\n",
    "t(\\mathbf{w}, \\mathbf{x}_0^{\\mathbf{w},\\delta}),\n",
    "\\end{equation}\n",
    "where $\\mathbf{w} \\in \\{0, 1\\}^N$ is the actual assignment and $W$ is generated at random according to\n",
    "the experimental design."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A test statistic $t$ is _effect-increasing_ if for any $\\mathbf{z} \\in \\{0, 1\\}^N$ and any $\\mathbf{y}$,\n",
    "$\\eta$, and $\\zeta$ in $\\Re^N$ with $\\eta \\ge 0 \\ge \\zeta$,\n",
    "\\begin{equation}\n",
    "t(\\mathbf{z}, \\mathbf{y} + \\mathbf{z}\\circ \\eta + (\\mathbf{1} - \\mathbf{z}) \\circ \\zeta) \\ge t(\\mathbf{z},\\mathbf{y}) .\n",
    "\\end{equation}\n",
    "That is, the test statistic increases if any control responses decrease or any treatment responses increase.\n",
    "\n",
    "A test statistic $t$ is _differential-increasing_ if for any \n",
    "$\\mathbf{z}, \\mathbf{a} \\in \\{0, 1\\}^N$ and any $\\mathbf{y}$,\n",
    "$\\eta$, in $\\Re^N$ with $\\eta \\ge 0$,\n",
    "\\begin{equation}\n",
    "t(\\mathbf{z}, \\mathbf{y} + \\mathbf{a}\\circ \\eta) − t(\\mathbf{z}, \\mathbf{y}) \\le\n",
    "   t(\\mathbf{a},\\mathbf{y}+\\mathbf{a}\\circ \\eta )−t(\\mathbf{a},\\mathbf{y}).\n",
    "\\end{equation}\n",
    "That is, adding something (positive) to a subset of the responses increases the test statistic the most\n",
    "when the increases are to responses of subjects assigned to treatment.\n",
    "\n",
    "A test statistic $t$ is _distribution-free_ if for any $\\mathbf{y}$, $\\mathbf{y}' \\in \\Re^N$,\n",
    "\\begin{equation}\n",
    "t(W, \\mathbf{y}) \\sim t(W, \\mathbf{y}'),\n",
    "\\end{equation}\n",
    "where $W$ is generated at random following the treatment assignment mechanism.\n",
    "That is, the probability distribution of the test statistic does not depend on the observed\n",
    "data. (Many rank statistics are distribution-free.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem** (Caughey et al. 2021 Theorem 1)\n",
    "\n",
    "(a) If the test statistic $t$ is either differential-increasing or effect-increasing, \n",
    "then for any constant $\\delta \\in \\Re^N$, the corresponding randomization p-value \n",
    "for the strong null $H_\\delta$ is also valid for testing the bounded null \n",
    "$H_{\\le \\delta}$. I.e., under $H_{\\le \\delta}$, $\\mathbb{P} \\{P_{\\mathbf{w},\\delta} \\le \\alpha \\} \\le \\alpha$ for\n",
    "all $\\alpha \\in [0, 1]$. \n",
    "\n",
    "(b) If the test statistic $t$ is differential-increasing, or it is both effect-increasing and \n",
    "distribution-free, then for any possible treatment assignment $w \\in \\{0, 1\\}^N$ \n",
    "the corresponding randomization $P$-value, $P_{\\mathbf{w},\\delta}$ \n",
    "is a monotone increasing function of $\\delta$: \n",
    "if $\\delta$, $\\bar{\\delta} \\in \\Re^N$ with $\\delta \\le \\bar{\\delta}$, \n",
    "then $P_{\\mathbf{w},\\delta} \\le P_{\\mathbf{w}, \\bar{\\delta}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sketch proof: If $\\tau \\le \\delta$, then\n",
    "\\begin{equation}\n",
    "\\mathbf{x}_0^{\\mathbf{W},\\delta} - \\mathbf{x}_0 = \\mathbf{W} \\circ (\\tau - \\delta) \\le 0\n",
    "\\end{equation}\n",
    "and \n",
    "\\begin{equation}\n",
    "\\mathbf{x}_1^{\\mathbf{W},\\delta} - \\mathbf{x}_1 = (\\mathbf{1} - \\mathbf{W}) \\circ (\\delta - \\tau) \\ge 0.\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}