permute-classical.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7ceae016-25bd-4ca2-b11d-9d257aedcdb9",
   "metadata": {},
   "source": [
    "# Permutation tests based on ranks, and calibrating parametric tests using permutations\n",
    "\n",
    "This notebook explores some common nonparametric tests based on ranks, and shows that (depending\n",
    "on the null hypothesis) parametric tests can be calibrated nonparametrically to have the correct significance level when the assumptions of the parametric tests fail."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d872a33e-a186-4f53-a194-8b093396066b",
   "metadata": {},
   "source": [
    "## Notation\n",
    "\n",
    "Given a set of real numbers $\\{x_j\\}_{j=1}^N$, $x_{(i)}$ denotes the $i$th smallest element of the set. \n",
    "That is, \n",
    "\\begin{equation}\n",
    "x_{(1)} \\le x_{(2)} \\le \\ldots \\le x_{(N)}.\n",
    "\\end{equation}\n",
    "\n",
    "The _(mid-)rank_ $r_j$ of the $j$th observation, $x_j$, is $\\#\\{ k: x_k < x_j \\} + (\\#\\{ k: x_k = x_j \\}+1)/2$.\n",
    "This assigns tied observations the average of the ranks they would have had if they had been distinct.\n",
    "\n",
    "For example, for the set $\\{1, 4, 2, 6\\}$ the corresponding (mid-)ranks are $1, 3, 2, 4$; for the set\n",
    "$\\{1, 4, 2, 6, 2, 2\\}$, the mid-ranks are  $1, 5, 3, 6, 3, 3$;\n",
    "and for the set $\\{1, 4, 2, 6, 2, 2, 2\\}$, the mid-ranks are  $1, 6, 3.5, 7, 3.5, 3.5, 3.5$.\n",
    "\n",
    "The sum of the (mid-)ranks of $N$ observations is $N(N+1)/2$.\n",
    "\n",
    "When the data are random variables $\\{X_j \\}_{j=1}^N$, the (random) ranks will be written using uppercase letters, $\\{R_j\\}_{j=1}^N$.\n",
    "When there are only two groups, control and treatment, $\\{S_j\\}_{j=1}^n$ will denote the ranks of the responses\n",
    "in the treatment group.\n",
    "\n",
    "Many of the methods described in this notebook can be performed using the original data or using\n",
    "the ranks of the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcea2117-596b-4d36-ac00-24f93bbf1494",
   "metadata": {},
   "source": [
    "### Why use ranks instead of the original data?\n",
    "\n",
    "+ In some problems, only ranks are observed. This is common when it is possible to make comparative judgements about outcomes, but not absolute judgements. For example, a doctor might be able to rank a collection of patients according to the severity with which they suffer from some disease they share, but might _not_ be able to rate the severity on an objective quantitative scale&mdash;and certainly not an equal-increment scale, e.g., a scale for which 5 is \"better than\" 4 by the same amount that 2 is \"better than\" 1. Similarly, a consumer might be able to order a collection of items by preference: she prefers A to B, B to C and C to D, for example. Yet she might not be able to quantify the strength of her preferences. Having techniques for dealing with ranks is then useful. The same issue occurs with measurements on Likert scales.\n",
    "+ If the measurements are \"contaminated\" by outliers, using ranks can decrease the sensitivity to the outliers, potentially increasing power.\n",
    "+ If the original data have no ties, the null distribution of the rank-based statistic is \"universal\": can tabulate it once and for all. (This is less important than it used to be because of increases in computing power.)\n",
    "+ Often, the normal approximation to the null distribution of the rank-based statistic is more accurate than the normal approximation to the null distribution of the test statistic for raw data. (This is less important than it used to be because of increases in computing power.)\n",
    "\n",
    "### Why use the original data instead of ranks?\n",
    "\n",
    "+ The original units are more meaningful for things like treatment effect, estimating shifts, etc.\n",
    "+ Tests using the original data may have more power.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58de0b71-6748-4726-b4e2-a09a09c1d93e",
   "metadata": {},
   "source": [
    "## Notation review\n",
    "\n",
    "As in the notebook on [causal inference](./causal-inference.ipynb), there are $T$ possible treatments (or $T$ groups being compared), numbered 0 to $T-1$, \n",
    "and $W_j$ is the treatment that subject $j$ receives (or the group that subject $j$ belongs to).\n",
    "Subject $j$'s response if assigned to treatment $t$ is $x_{jt}$, $j=1, \\ldots, N$, $t=0, \\ldots, T-1$.\n",
    "The number of subjects assigned to treatment $t$ is $\\sum_{j=1}^N 1(W_j = t)$.\n",
    "(When there are only two treatments, $t=0$ and $t=1$, $n :=\\sum_j W_j$ is the number of subjects \n",
    "assigned to treatment 1.)\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}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b629b22b-4bbe-4f0d-84eb-31abc9229ef0",
   "metadata": {},
   "source": [
    "## The two-sample problem\n",
    "\n",
    "We first specialize to the case of $T=2$ groups or treatments. \n",
    "There are $N$ items in all, $n$ in the \"treatment\" group and $m := N-n$ in the \"control\" group.\n",
    "Let $\\{X_j\\}_{j=1}^m$ denote the data for the control group and $\\{Y_j\\}_{j=1}^n$ denote the\n",
    "data for the treatment group.\n",
    "\n",
    "In the ***randomization*** version of the two-sample problem, the two groups are provided by nature, and the \n",
    "question is whether the groups are \"statistically distinguishable,\" in the sense that\n",
    "they don't look like a single group that was partitioned randomly. \n",
    "In the context of causal inference,\n",
    "the groups started as a single group of subjects, randomly partitioned into two groups\n",
    "that receive different treatments, and the question is whether their responses are statistically\n",
    "distinguishable.\n",
    "\n",
    "There is also a ***population*** version of the two-sample problem, which asks whether\n",
    "the two groups look like IID samples from the same parent population, i.e., whether there is some\n",
    "$F$ for which $\\{X_1, \\ldots, X_m, Y_1, \\ldots, Y_n \\}$ are IID $F$.\n",
    "If the two groups are IID samples from the same distribution, then, conditional on the set of $N$ observed values,\n",
    "every subset of $n$ of the \n",
    "$N$ values is equally likely to be the group of size $n$ and every subset of $m$\n",
    "of the $N$ values is equally likely to be the group of size $m$, so tests that work for\n",
    "the randomization version of the 2-sample problem also work for the population version\n",
    "of the 2-sample problem.\n",
    "\n",
    "\n",
    "Lehmann (1998, p.&nbsp;64&ndash;65) identifies five models for five problems, all of which lead to the same (conditional) null distribution of ranks:\n",
    "\n",
    "+ Randomization model for comparing two treatments. $N$ subjects are given and fixed; $n$ are assigned at random to treatment and $m= N-n$ to control.\n",
    "    \n",
    "\n",
    "+ Population model for comparing two treatments. $N$ subjects are a drawn as a\n",
    "simple random sample from a much larger population;\n",
    "$n$ are assigned at random to treatment and\n",
    "$m = N-n$ to control.\n",
    "Condition on the configuration of ties if there are ties among the data.\n",
    "\n",
    "\n",
    "+ Comparing two sub-populations using a sample from each.\n",
    "A simple random sample of $n$ subjects is drawn from one much larger population\n",
    "with many more than $N$ members, and\n",
    "a simple random sample of $m$ subjects is drawn from another population with\n",
    "many more than $m$ members.\n",
    "Condition on the configuration of ties if there are ties among the data.\n",
    "\n",
    "\n",
    "+ Comparing two sub-populations using a sample from the pooled population.\n",
    "A simple random sample of $N$ subjects is drawn from the pooled population,\n",
    "giving random samples from the two populations but with random sample sizes.\n",
    "Condition on the sample sizes and, if there are ties, on the configuration of the ties.\n",
    "\n",
    "\n",
    "+ Comparing two sets of measurements. Independent sets of $n$ and $m$\n",
    "measurements come from two sources.\n",
    "Condition on the configuration of ties if there are ties among the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2be9b65-9c08-4613-b994-48b854159a48",
   "metadata": {},
   "source": [
    "### Testing for a difference in location"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5e51408-2490-4dff-9ffc-2dc2b4f351da",
   "metadata": {},
   "source": [
    "#### The Wilcoxon rank-sum test.\n",
    "\n",
    "Define the sum of the ranks of the responses in the treatment group, group 1:\n",
    "\\begin{equation}\n",
    "W_Y := \\sum_j W_j R_j.\n",
    "\\end{equation}\n",
    "(Note that $W$ is being used in two very different ways in the notation.)\n",
    "This is the _Wilcoxon rank sum statistic_.\n",
    "Under the null, every subset of $n$ of the responses is equally likely to be the\n",
    "treatment group.\n",
    "The expected value of $W_Y$ under the null is $n(N+1)/2$.\n",
    "Moreover, under the null, the distribution of $W_Y$ is symmetric.\n",
    "If treatment tends to increase responses, $W_Y$ will tend to be larger than its median under the null; \n",
    "if treatment tends to decrease responses, $W_Y$ will tend to be smaller than its median under the null.\n",
    "\n",
    "To test against the alternative that treatment tends to increase responses,\n",
    "we would reject for large values of $W_Y$.\n",
    "To test against the alternative that treatment tends to decrease responses,\n",
    "we would reject for small values of $W_Y$.\n",
    "The critical value of the test is set using the probability distribution of\n",
    "$W_Y$ on the assumption that the strong null hypothesis is true.\n",
    "For a level-alpha test against the alternative that treatment increases\n",
    "responses, we would find the smallest $c$ such that, if the\n",
    "strong null is true, $\\mathbb{P}(W_Y \\ge c) \\le \\alpha$.\n",
    "We would then reject the strong null if the observed value of\n",
    "$W_Y$ is $c$ or greater.\n",
    "Recall that $1 + 2 + \\cdots + k = k(k+1)/2$.\n",
    "\n",
    "If the treated subjects have the smallest possible ranks, 1 to $n$, then\n",
    "$W_Y = 1 + 2 + \\cdots + n = n(n+1)/2$.\n",
    "If the treated subjects have the largest possible ranks, $N - n+1$ to $N$, then \n",
    "\\begin{eqnarray}\n",
    "W_Y & = & (N-n+1) + (N-n+2) + \\cdots + N \\\\\n",
    "    & = & (N-n) + 1 + (N-n) + 2 + \\cdots (N-n) + n \\\\\n",
    "    & = & n(N-n) + (1 + 2 + \\cdots + n) \\\\\n",
    "    & = & n(N-n) + n(n+1)/2.\n",
    "\\end{eqnarray}\n",
    "All the integers between $n(n+1)/2$ and $n(N-n) + n(n+1)/2$ are possible values of $W_Y$.\n",
    "The null distribution of $W_Y$ under the strong null hypothesis is symmetric about\n",
    "$n(N+1)/2$.\n",
    "\n",
    "(In the randomization model, each subset of $n$ of the $N$ ranks is equally likely to be the \n",
    "ranks of the treatment group.  \n",
    "The probability that the treatment ranks are $\\{1, 2, \\ldots , n\\}$ is equal to the \n",
    "probability that the treatment ranks are $\\{N, N-1, \\ldots, N-n+1\\}$.  \n",
    "More generally, consider re-labeling the $j$th-ranked observation to be the $N-j+1$st-ranked observation.\n",
    "The labels of the observations would still be $\\{1, 2, \\ldots , N\\}$, so the probability \n",
    "distribution of the sum of the labels on the treated subjects would be the same \n",
    "as the probability distribution of $W_Y$ under the strong null  hypothesis.  \n",
    "However, if the sum of the treatment ranks was $W_Y$,\n",
    "the sum of the new labels would be $n(N+1)-W_Y$. \n",
    "Thus $\\mathbb{P}\\{W_Y = k\\} = \\mathbb{P}\\{W_Y = n(N+1) - k \\}$. \n",
    "That is, the probability distribution of $W_Y$ is symmetric about $n(N+1)/2$.  \n",
    "This argument is essentially that in Lehmann, E.L., 1998. _Nonparametrics: Statistical \n",
    "Methods Based on Ranks_. Upper Saddle River, N.J.: Prentice Hall, pp 12-13.\n",
    "\n",
    "\n",
    "Thus the expected value of $W_Y$ under the null hypothesis is $n(N+1)/2$.\n",
    "(This also follows from the fact that the expected value of the sample sum of\n",
    "a simple random sample of size $n$ from a box of $N$ numbers is\n",
    "equal to $n$ times the mean of the $N$ numbers; the mean of the\n",
    "integers 1 to $N$ is $(N+1)/2$.)\n",
    "The variance of $W_Y$ under the strong null hypothesis\n",
    "is $mn(N+1)/12$.\n",
    "(This follows from the fact that the variance of the sample sum of a simple random sample of size $n$ from a list of $N$ \n",
    "numbers is $(N-n)n(\\mbox{variance of list})/(N-1) = mn(\\mbox{variance of list})/(N-1)$. \n",
    "The variance of the list of the integers 1 to $N$ is $(N^2-1)/12 = (N-1)(N+1)/12$, \n",
    "so the variance of $W_Y$ is $mn(N-1)(N+1)/(12(N-1)) = mn(N+1)/12$.)\n",
    "\n",
    "Define $W_X$ to be the sum of the control ranks.\n",
    "Because the treatment ranks and control ranks together comprise all the ranks,\n",
    "$ W_X + W_Y = 1 + 2 + \\cdots + N = N(N+1)/2$.\n",
    "Thus, $W_X  = N(N+1)/2 - W_Y$.\n",
    "\n",
    "It follows that the (null) expected value of $W_X$ is\n",
    "$N(N+1)/2 - n(N+1)/2 = m(N+1)/2$, and that the (null) variance of $W_X$ is also\n",
    "$mn(N+1)/12$.\n",
    "When $n$ and $m$ are both large, the normal approximations to the null distributions of \n",
    "$W_X$ and to $W_Y$ tend to be accurate.\n",
    "\n",
    "The Wilcoxon rank-sum test is exactly what you get if you replace each observation by its\n",
    "rank, then do a permutation test using the sum of the (ranks of the) responses in the treatment group.\n",
    "Thus, the code in the previous chapter can be used to simulate the null distribution,\n",
    "by replacing the raw data with their ranks.\n",
    "However, the normal approximation to the null distribution\n",
    "of $W_Y$ can be better than the normal approximation to the null\n",
    "distribution of the sample sum of the responses to treatment, because the ranks\n",
    "are evenly spread out, whereas the raw responses can be highly skewed or multimodal.\n",
    "\n",
    "### Mann-Whitney Statistics\n",
    "\n",
    "Define $W_{XY} := W_Y - n(n+1)/2$.\n",
    "This is $W_Y$ minus its minimum possible value.\n",
    "Let $W_{YX} := W_X - m(m+1)/2$, $W_X$ minus its minimum possible value.\n",
    "Under the strong null hypothesis, the probability distribution of $W_{XY}$\n",
    "is the same as the probability distribution of $W_{YX}$,\n",
    "a consequence of the symmetry of the probability distribution of $W_Y$.\n",
    "The statistics $W_{XY}$ and $W_{YX}$ are called the _Mann-Whitney_ statistics.\n",
    "\n",
    "Let $\\{X_1, \\ldots, X_m\\}$ denote the $m$ control responses and let\n",
    "$\\{Y_1, \\ldots, Y_n\\}$ denote the $n$ treatment responses.\n",
    "Consider $\\#\\{(i, j) : 1 \\le i \\le m, \\; 1 \\le j \\le n \\mbox{ and } X_i < Y_j \\}$.\n",
    "This is the number of (control, treatment) pairs such that the control\n",
    "response is less than the treatment response.\n",
    "\n",
    "Let $\\{S_j: j = 1, \\ldots , n\\}$ be the ranks of the treatment responses.\n",
    "Let $\\{S_{(j)}: j = 1, \\ldots , n\\}$ be those ranks in increasing order.\n",
    "Partition the set of (control, treatment)\n",
    "pairs for which the control response is less than the treatment response on the value of the treatment\n",
    "response: \n",
    "The total number of\n",
    "such pairs for which the control response is less than the treatment response\n",
    "is the number of pairs where the control response is less than the smallest treatment\n",
    "response, plus the number of pairs where the control response is less than the second-smallest\n",
    "treatment response, and so on. This will help us count the pairs.\n",
    "\n",
    "How many responses are less than $S_{(1)}$, the rank of the smallest\n",
    "treatment response?\n",
    "By definition, there are $S_{(1)}-1$ of them&mdash;all of which are control responses.\n",
    "The number of response values that are less than $S_{(2)}$  is\n",
    "$S_{(2)}-1$, one of which is $S_{(1)}$, so the\n",
    "total number of control responses that are less than $S_{(2)}$ is\n",
    "$S_{(2)}-2$.\n",
    "The total number of control responses that are less than $S_{(j)}$\n",
    "is $S_{(j)} - j$, so the total number of (control, treatment)\n",
    "pairs with the control response less than the treatment response is\n",
    "\\begin{equation}\n",
    "S_{(1)}-1 + S_{(2)} - 2 + \\cdots +\n",
    "    S_{(n)} - n = W_Y - (1 + 2  + \\cdots + n) =\n",
    "    W_Y - n(n+1)/2 = W_{XY}.\n",
    "\\end{equation}\n",
    "Thus \n",
    "\\begin{equation}\n",
    "W_{XY} = \\#\\{(i, j) : 1 \\le i  \\le m, \\; 1 \\le j \\le n \\mbox{ and } X_i < Y_j \\}.\n",
    "\\end{equation}\n",
    "That is, the Mann-Whitney statistic $W_{XY}$ (and the Wilcoxon rank sum\n",
    "$W_Y$, up to an additive constant)\n",
    "measures the number of (control, treatment) pairs for which the treatment response\n",
    "is at least as large as the control response.\n",
    "The larger the positive effect of treatment, the larger the Mann-Whitney and Wilcoxon\n",
    "rank sum statistics tend to be.\n",
    "\n",
    "**Equivalent tests.**\n",
    "Multiplying the test statistic by a (positive) constant produces an equivalent test (_equivalent_ means\n",
    "that the tests reject the null for exactly the same data).\n",
    "Hence using $\\bar{W}_Y := W_Y/n$ produces an equivalent test.\n",
    "\n",
    "Since $\\sum_j r_j = N(N+1)/2$, as $W$ increases, the sum (and the mean) of the control ranks $W_X := \\sum_j (1-W_j) r_j$ decreases monotonically, so the Wilcoxon rank-sum test is equivalent to a test based on the difference between the mean ranks for treatment and control, $W_Y/n - W_X/m$.\n",
    "\n",
    "The Mann-Whitney test and the Wilcoxon rank-sum test are equivalent."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c43e1aa2-1811-4d11-a291-241b56deea6e",
   "metadata": {},
   "source": [
    "#### Permutation $t$-test.\n",
    "\n",
    "Use the $t$ statistic as the test statistic, but calibrate it using the permutation distribution\n",
    "rather than Student's t distribution.\n",
    "See [introduction to permutation tests](./permute-intro.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e46cdc9-10ac-49a1-8c0f-c21e4619b7dc",
   "metadata": {},
   "source": [
    "## Confidence bounds for a shift in the two-sample problem\n",
    "\n",
    "The strong null is that each subject would have the same response regardless of the treatment assigned to the\n",
    "subject.\n",
    "One alternative hypothesis is that the effect of treatment is to shift each subject's response by the same amount,\n",
    "$\\Delta$. \n",
    "In applications, there is generally little reason to belive that treatment will affect every subject by exactly\n",
    "the same amount: this is a very artificial alternative.\n",
    "But it is a _strong_ alternative, in the sense that it lets us find the sampling distribution of any statistic from the observed data.\n",
    "\n",
    "The approach here is based on the duality between tests and confidence sets.\n",
    "Suppose we have a family of level $\\alpha$ tests of the null hypotheses\n",
    "$H_{d}$: $\\Delta=d$ for all real $d$.\n",
    "Let $S(Y)$ denote the set of values of $d$ for which\n",
    "the test of the corresponding hypothesis $\\Delta = d$ would not rejected\n",
    "the null hypothesis if the data were $y$ (here $y$ is generic; in the problems we have been considering,\n",
    "$y$ specifies all the treatment responses and all the control responses).\n",
    "Then $S(Y)$ is a $1-\\alpha$ confidence set for $\\Delta$.\n",
    "\n",
    "\n",
    "**Proof**.\n",
    "If $d$ is the true value of $\\Delta$, then the probability (under $d$)\n",
    "that the corresponding test rejects the hypothesis $\\Delta = d$ is\n",
    "at most $\\alpha$.\n",
    "The probability that the test does not reject is at least $1-\\alpha$, so with probability\n",
    "$1-\\alpha$ under $\\Delta$, the value $\\Delta$ is in $S(Y)$.\n",
    "That is,\n",
    "\\begin{equation}\n",
    "    \\mathbb{P}_d \\{S(Y) \\ni \\Delta \\} \\ge 1-\\alpha.\n",
    "\\end{equation}\n",
    "The shape and size of the confidence set are tied to the family of hypothesis tests\n",
    "that are being \"inverted.\"\n",
    "No matter what, the procedure gives a $1-\\alpha$ confidence set, but it\n",
    "can be large or bizarre (e.g., consist of disjoint pieces) according to the\n",
    "nature of the tests.\n",
    "The duality between tests and confidence sets can be exploited to produce\n",
    "confidence sets with desirable properties; see, for example,\n",
    "Benjamini and Stark (1996. Non-equivariant simultaneous confidence intervals\n",
    "less likely to contain zero, _J. Amer. Stat. Assoc._, _91_, 329&ndash;337);\n",
    "Benjamini, Hochberg and Stark (1998. Confidence Intervals with more Power\n",
    "to determine the Sign: Two Ends constrain the Means, _J. Amer. Stat. Assoc._,\n",
    "_93_, 309&ndash;317); or Evans, Hansen and Stark (2005.\n",
    "Minimax Expected Measure Confidence Sets for Restricted Location Parameters,\n",
    "_Bernoulli, 11,_ 571&ndash;590).\n",
    "\n",
    "\n",
    "\n",
    "To use this result to find a confidence set for $\\Delta$, we need a family of level $\\alpha$ tests of the hypotheses $\\Delta = d$.\n",
    "\n",
    "Consider the Mann-Whitney statistic $W_{XY}$.\n",
    "Under the strong null hypothesis, the expected value of $W_{XY}$\n",
    "is $m$$n$/2.\n",
    "Suppose that the critical values for the two-sided test of the hypothesis\n",
    "that $\\Delta=0$ are $mn/2-c$  and $mn/2+c$.\n",
    "Note that under the hypothesis that $\\Delta = d$, the distribution of\n",
    "the Mann-Whitney statistic for the data $\\{X_1, \\ldots, X_m, Y_1-d, \\ldots, Y_n-d \\}$,\n",
    "which we denote $W_{XY-d}$, is the same as the distribution of\n",
    "$W_{XY}$ under the original strong null hypothesis.\n",
    "A two-sided level-$\\alpha$ test of the hypothesis $\\Delta=d$ thus would reject\n",
    "if $W_{XY-d} \\le mn/2-c$ or $W_{XY-d} \\ge mn/2+c$.\n",
    "This is a family of tests on which we can base the confidence set.\n",
    "\n",
    "We can find a confidence interval for $\\Delta$ using these tests as follows:\n",
    "Let $d$ be the estimated treatment effect, the median of the $mn$ differences $d_{ij}$.\n",
    "Subtract $d$ from the treatment observations; compute the\n",
    "ranks of this new data set; perform a two-sided Wilcoxon rank sum test of\n",
    "the null hypothesis of no treatment effect using these new data.\n",
    "If the test does not reject, include $d$ in the confidence set.\n",
    "Increase $d$ systematically until the test rejects to get the\n",
    "upper endpoint of the confidence interval; decrease $d$ systematically\n",
    "until the test rejects to get the lower endpoint of the confidence interval.\n",
    "\n",
    "It is clear that the only values of $d$ that could be endpoints of the confidence\n",
    "interval are those that cause the ranks of the treatment observations to\n",
    "change&mdash;only those values\n",
    "change the value of $W_{XY}$.\n",
    "We can exploit this to find a confidence interval for $\\Delta$ with less trial-and-error\n",
    "searching.\n",
    "\n",
    "To change the value of $W_{XY-d}$, we need to change $d$ by enough to change the number of\n",
    "pairs $(i, j)$ such that $X_i < Y_j-d$.\n",
    "The values of $d$ that change the ranks are precisely the differences\n",
    "$d_{ij} = Y_j-X_i$.\n",
    "\n",
    "Let $d_{(1)}, \\ldots, d_{(mn)}$ \n",
    "denote the order statistics of the $mn$ differences $\\{d_{ij}\\}$,\n",
    "and define $d_{(0)} := -\\infty$ and $d_{(mn+1)} := \\infty$. \n",
    "Then $d_{(k)}  \\le d <  d_{(k+1)}$\n",
    "if and only if $W_{XY-d} = mn-k$.\n",
    "\n",
    "To see this, note that $W_{XY-d} \\le mn-k$ if $mn-k$ of the differences\n",
    "$\\{(Y_j - d - X_i)\\}$ are at least zero, that is, if $d_{(k)} \\le d$.\n",
    "The distribution of $W_{XY-d}$ when $d$ is the true shift\n",
    "does not depend on $d$: it  is the same as the distribution of $W_{XY}$ \n",
    "under the strong null hypothesis.\n",
    "Using the symmetry of the null distribution of $W_{XY}$, we see that for any $d$,\n",
    "\\begin{equation}\n",
    "\\mathbb{P}_d \\{ d_{(k)} \\le d < d_{(k+1)} \\} = \\mathbb{P}_0 \\{W_{XY} = mn - k \\} =\n",
    "    \\mathbb{P}_0 \\{W_{XY} = k \\}.\n",
    "\\end{equation}\n",
    "\n",
    "Thus the random differences $d_{ij}$\n",
    "carve the real line into random intervals with known probabilities of containing\n",
    "the true value of $d$ (if the shift alternative is true!).\n",
    "To find a confidence interval for $d$, we can start with (one of) the interval\n",
    "$[d_{(k)}, d_{(k+1)})$ \n",
    "with the largest chance of containing $d$, then append intervals in\n",
    "decreasing order of their probability of containing $d$ until the\n",
    "sum of their probabilities is at least $1-\\alpha$.\n",
    "Because of the unimodality of the distribution,\n",
    "these intervals will be contiguous; let\n",
    "$[d_{(\\ell)}, d_{(u)})$ \n",
    "denote the interval.\n",
    "Alternatively, we could choose the intervals to try to have symmetric\n",
    "non-coverage probabilities; that is, so that\n",
    "$\\mathbb{P}_d \\{ d_{(\\ell)} > d \\}$  is as close as possible to\n",
    "$\\mathbb{P}_d \\{ d_{(u)} < d \\}$.\n",
    "(I don't know whether these two approaches differ; I think the former corresponds to\n",
    "inverting the two-sided Wilcoxon test.)\n",
    "All the probability calculations are implicit in calculating the distribution of\n",
    "$W_{XY}$ under the strong null hypothesis.\n",
    "In fact, one need not even calculate the entire distribution; see Lehmann\n",
    "(1998, pp.&nbsp;92ff) for shortcuts."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ee2f501-ebc4-4c0b-923b-72f0a0e4357c",
   "metadata": {},
   "source": [
    "### Testing for a difference in dispersion: the Siegel-Tukey test\n",
    "\n",
    "Assign \"ranks\" differently: smallest observation has \"rank\" 1, largest has \"rank\" 2, second-smallest has \"rank\" 3,\n",
    "etc. Then apply Wilcoxon the rank-sum test to these \"ranks.\" The null distribution of the sum is the same, but this\n",
    "statistic has power against the alternative that treatment affects dispersion (and not location).\n",
    "In contrast, the Wilcoxon rank-sum test has power against the alternative that treatment affects location, not dispersion.\n",
    "\n",
    "\n",
    "### Testing for a difference in distribution: the Smirnov test\n",
    "\n",
    "A common measure of the \"distance\" between two probability distributions $F$ and $G$ is the _Smirnov_ or _Kolmogorov-Smirnov_\n",
    "distance,\n",
    "\\begin{equation}\n",
    "\\| F-G \\|_\\infty := \\sup_{x \\in \\mathbb{R}} |F(x)-G(x)|.\n",
    "\\end{equation}\n",
    "A two-sample test sensitive to general differences (rather than changes in location and/or spread)\n",
    "can be based on that distance.\n",
    "\n",
    "Let $\\hat{F}$ be the empirical CDF of the responses to treatment and $\\hat{G}$ be the empirical\n",
    "CDF of the responses to control:\n",
    "\\begin{equation}\n",
    "\\hat{F}(x) := \\frac{1}{n} \\sum_{j=1}^n 1_{x \\ge X_j}\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "\\hat{G}(x) := \\frac{1}{m} \\sum_{j=1}^n 1_{x \\ge Y_j}.\n",
    "\\end{equation}\n",
    "The _Smirnov_ test uses $\\|\\hat{F} - \\hat{G}\\|$ as the test statistic.\n",
    "\n",
    "Note that the value of $\\|\\hat{F} - \\hat{G}\\|$ depends only on the relative ordering\n",
    "of the values $\\{X_j\\}$ and $\\{Y_j\\}$, not on their numerical values.\n",
    "Thus, replacing the original data by their ranks does not change the test statistic\n",
    "(or its distribution under the null hypothesis): the Smirnov test is a rank-based test.\n",
    "Moreover, if there are no ties, the null distribution of the test statistic depends only on $n$ and $m$.\n",
    "\n",
    "Here is a simulation of the Smirnov test statistic under the null hypothesis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "db512aaf-a987-457a-80d2-be12353abb3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "from scipy.stats import binom\n",
    "import pandas as pd\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display, HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "c3f0e04d-02cf-4ea6-968e-980beba69abd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5bed697255b8494b8e8754641ee5e608",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=10, description='m', max=300, min=3), IntSlider(value=10, description='n…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plotSmirnov(m, n)>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def ecdf(x):\n",
    "    '''\n",
    "       calculates the empirical cdf of data x\n",
    "       returns the unique values of x in ascending order and the cumulative probabity at those values\n",
    "       NOTE: This is not an efficient algorithm: it is O(n^2), where n is the length of x. \n",
    "       A better algorithm would rely on the Collections package or something similar and could work\n",
    "       in O(n log n)\n",
    "    '''\n",
    "    theVals = sorted(np.unique(x))\n",
    "    theProbs = np.array([sum(x <= v) for v in theVals])/float(len(x))\n",
    "    if (theVals[0] > 0.0):\n",
    "        theVals = np.append(0., theVals)\n",
    "        theProbs = np.append(0., theProbs)\n",
    "    return theVals, theProbs\n",
    "\n",
    "\n",
    "def plotSmirnov(m, n):\n",
    "    '''\n",
    "    plot the ecdfs and Smirnov statistic for samples from a uniform distribution\n",
    "    \n",
    "    '''\n",
    "    sam = np.random.uniform(size=n+m)\n",
    "    y, pry = ecdf(sam[0:n])\n",
    "    x, prx = ecdf(sam[n:])\n",
    "    fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "    ax.step(y, pry, color='b', where='post', label=f'ecdf(y): {n} data')\n",
    "    ax.step(x, prx, color='r', where='post', label=f'ecdf(x): {m} data')\n",
    "    ax.legend(loc='best')\n",
    "    diff = [abs(np.sum(sam[n:] <= zi)/m - np.sum(sam[0:n] <= zi)/n) for zi in sam]\n",
    "    dLoc = np.argmax(diff)\n",
    "    zm = sam[dLoc]\n",
    "    xd = np.sum(sam[n:] <= zm)/m\n",
    "    yd = np.sum(sam[0:n] <= zm)/n\n",
    "    ymin, ymax = min(xd, yd), max(xd, yd)\n",
    "    ax.axvline(x=zm, ymin=ymin, ymax=ymax, color='g', linewidth='2')\n",
    "    ax.text(0.5, 0.1, f'Smirnov statistic: {round(ymax-ymin, 3)}', color='g', weight='heavy')\n",
    "    plt.show()\n",
    "\n",
    "interact(plotSmirnov, m=widgets.IntSlider(min=3, max=300, step=1, value=10), n=widgets.IntSlider(min=3, max=300, step=1, value=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff07b801-e839-4bf8-b7b0-798df673c436",
   "metadata": {},
   "source": [
    "#### Wilcoxon signed rank test\n",
    "\n",
    "The _signed (mid)-rank_ of an observation $x_k \\in \\{x_j\\}_{j=1}^N$ is $\\mathrm{sgn}(x_k) r_k$, where\n",
    "$r_k$ is the rank of $|x_k|$ in the set $\\{|x_j|\\}_{j=1}^N$, and \n",
    "\\begin{equation}\n",
    "\\mathrm{sgn}(x) := \\left \\{ \n",
    "    \\begin{array}{lr}  \n",
    "    -1, & x < 0 \\cr\n",
    "    0, & x=0 \\cr\n",
    "    1, & x > 0. \n",
    "    \\end{array}\n",
    "    \\right .\n",
    "\\end{equation}\n",
    "\n",
    "Suppose there are $N$ pairs of subjects; in each pair, one subject is assigned to treatment and one to control, at random,\n",
    "with equal probability.\n",
    "Look at the differences between the treated and control subject in each pair.\n",
    "Add the signed ranks of those differences.\n",
    "To test against the alternative that treatment increases (decreases) responses, reject if the sum is sufficiently large\n",
    "(small).\n",
    "\n",
    "\n",
    "[MORE TO COME]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5b9e92f-8867-469a-b775-6b42cdb1a4b5",
   "metadata": {},
   "source": [
    "## The $k$-sample problem\n",
    "\n",
    "(Calling this $k$, but we are using $T$ to denote the number of groups.)\n",
    "Size of group $t$ is $n_t$, $t=0, \\ldots, T-1$.\n",
    "\n",
    "### Testing for differences in location\n",
    "\n",
    "#### The Kruskal-Wallis test\n",
    "\n",
    "\\begin{equation}\n",
    "    K := \\frac{12}{N(N+1)} \\sum_{t=0}^{T-1} n_t \\left ( \\bar{R}_t - (N+1)/2 \\right )^2,\n",
    "\\end{equation}\n",
    "where $\\bar{R}_t$ is the mean of the ranks of the observations in group $t$.\n",
    "(The mean of all $N$ (mid-)ranks is $(N+1)/2$.)\n",
    "\n",
    "Recall that the $F$-statistic is \n",
    "\n",
    "\\begin{equation}\n",
    " F := \\frac{\\mbox{ \"explained variance\" }}{\\mbox{\"unexplained variance\"}} =\n",
    " \\frac{\\sum_{t=0}^{T-1} n_t({\\bar {Y}}_{t\\cdot }-{\\bar {Y}})^{2}/(T-1)}{\\sum_{t=0}^{T-1}\\sum _{j=1}^{n_t} \\left( Y_{tj}-{\\bar {Y}}_{t\\cdot } \\right)^{2}/(N-T)}.\n",
    "\\end{equation}\n",
    "\n",
    "Thus, $K$ is essentially the $F$-statistic of the ranks, since the overall sum of squared differences\n",
    "from the mean rank is constant.\n",
    "\n",
    "The asymptotic null distribution (as all the group sizes approach infinity)\n",
    "is chi-square with $T-1$ degrees of freedom.\n",
    "\n",
    "#### The permutation $F$ test\n",
    "\n",
    "Apply the usual $F$ test, but calibrate it using the permutation distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f72fae84-fafe-4940-89af-86f75f91e0d3",
   "metadata": {},
   "source": [
    "## Confidence bounds for quantiles\n",
    "\n",
    "Suppose $\\{X_j\\}_{j=1}^N$ are real random variables that are IID $F$. \n",
    "We want a confidence bound for $F_\\alpha$, the $\\alpha$ quantile of $F$, i.e., \n",
    "the smallest value of $x$ for which $F(x) := \\mathbb{P} \\{X_j \\le x \\} \\ge \\alpha$.\n",
    "\n",
    "The approach illustrated above for finding a confidence interval for the\n",
    "shift $d$ is very similar to a nonparametric approach to\n",
    "finding confidence intervals for a quantile of a probability distribution.\n",
    "Let $\\{X_j\\}_{j=1}^N$ be IID real-valued random variables with continuous cdf $F$.\n",
    "For $0 < q < 1$, let $x_q := F^{-1}(q) = \\inf \\{x: F(x) \\ge q \\}$ denote the $q$th quantile\n",
    "of the distribiution. \n",
    "We shall find a confidence interval for $x_q$ that does not depend on\n",
    "any other assumptions about $F$.\n",
    "\n",
    "If $x_q=x$, the probability that $X_j \\le x$ is $q$.\n",
    "The data are independent, so $N(x) := \\#\\{j : X_j \\le x\\}$ \n",
    "has a binomial distribution with parameters $n=N$ and $p=q$.\n",
    "Let $p_{N,q}(k) := {{N} \\choose {k}} q^k (1-q)^{N-k} = \\mathbb{P}_{x_q=x}\\{ N(x) = k \\}$,\n",
    "$k=0, \\ldots, N$.\n",
    "\n",
    "Note that $\\mathbb{P}_{x_q=x} \\{ N(x) = k \\}$ does not depend on $x_q$.\n",
    "The probability distribution of $N(x)$ is unimodal; we could build a level $\\alpha$\n",
    "test of the hypothesis that $x_q=x$ by assigning to the rejection region the values of $k$ for which\n",
    "$_{n,q}(k)$ is smallest, subject to the constraint that the sum of the probabilities is at most $\\alpha$.\n",
    "\n",
    "To find a confidence interval for $x_q$, we could start\n",
    "with the $k$th order statistic $X_{(k)}$ as a trial\n",
    "value for $x$, where $k$ is as close as possible to $qn$, then work\n",
    "\"outward\" as before, increasing and decreasing $x$ until the\n",
    "test of the hypothesis $x_q = x$  rejects.\n",
    "\n",
    "We can streamline the approach in the same way as we did to find a confidence\n",
    "interval for the treatment effect $\\Delta$ in the shift model.\n",
    "The values of $x$ at which $N(x)$  changes\n",
    "are the order statistics $X_{(j)}$, so it suffices to\n",
    "consider them when searching for the endpoints of the confidence interval.\n",
    "Now $N(x) = k$ iff $X_{(k)} \\le x < X_{(k+1)}$, $k = 0, 1, \\ldots, N$,\n",
    "where $X_{(0)} := -\\infty$ and $X_{(N+1)} := \\infty$.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathbb{P}_{x_q} \\{ X_{(k)} \\le x_q < X_{(k+1)} \\} =  {{N} \\choose {k}} q^k (1-q)^{N-k}.\n",
    "\\end{equation}\n",
    "\n",
    "The order statistics carve the real line into intervals with known chances of\n",
    "containing the true value of $x_q$.\n",
    "Those chances are determined by the binomial distribution with parameters\n",
    "$n$ and $q$.\n",
    "We can find a confidence interval for $x_q$\n",
    "by starting with the interval\n",
    "$[X_{(k)}, X_{(k+1)})$ \n",
    "with the largest chance of containing $x_q$, then appending\n",
    "intervals in decreasing order of their probability of containing $x_q$\n",
    "until the sum of their probabilities is at least $1-\\alpha$.\n",
    "\n",
    "**Exercise.** Extend this derivation to distributions that are not necessarily continuous\n",
    "to obtain a conservative confidence interval for $x_q$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b908f539-f4a9-473b-bdc2-c4b72999a0be",
   "metadata": {},
   "source": [
    "## One-sample tests\n",
    "\n",
    "Typically used to test whether data are a sample from a symmetric distribution or a distribution with\n",
    "a specified median."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f612220d-d4bf-40ca-8520-fb168d02b89a",
   "metadata": {},
   "source": [
    "## When there is no random assignment\n",
    "\n",
    "### Imaginary randomization\n",
    "\n",
    "In many situations in which hypothesis tests are used, there is no actual random assignment.\n",
    "This happens frequently in observational studies involving humans.\n",
    "For instance, the following characteristics are not really assigned at random\n",
    "\n",
    "+ gender\n",
    "+ ethnicity\n",
    "+ socio-economic status (SES)\n",
    "+ educational attainment\n",
    "+ whether people are housed, and if so, where they live\n",
    "\n",
    "In such situations, the corresponding null hypothesis is typically that the characteristic/label is assigned \"as if\" at random, i.e., that the characteristic is an arbitrary label unconnected with anything else.\n",
    "\n",
    "For more discussion, see \n",
    "\n",
    "+ Freedman, D.A., and R. Berk, 2012. Statistical Assumptions as Empirical Commitments. https://www.cambridge.org/core/books/abs/statistical-models-and-causal-inference/statistical-assumptions-as-empirical-commitments/B1879A1919929D5948F6BB539F91D786\n",
    "+ Freedman, D.A., and D. Lane, 1983. A Nonstochastic Interpretation of Reported Significance Levels, _J. Business and Economic Statistics_, _1_, 292-298. https://www.tandfonline.com/doi/abs/10.1080/07350015.1983.10509354"
   ]
  }
 ],
 "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": 5
}