tests-combo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8e92f105-53db-4b35-8573-c55cdd7bb45b",
   "metadata": {},
   "source": [
    "# Combining Hypothesis Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f33abd7-84b3-4422-92ee-38850f063944",
   "metadata": {},
   "source": [
    "## Intersection-Union Hypotheses\n",
    "\n",
    "In many situations, a null hypothesis of interest can be represented as the intersection of simpler hypotheses. \n",
    "For instance, the hypothesis that a university does not discriminate in its graduate admissions might be represented as \n",
    "\n",
    "(does not discriminate in arts and humanities) $\\cap$ (does not discriminate in sciences) $\\cap$ (does not discriminate in engineering) $\\cap$ (does not discriminate in professional schools).\n",
    "\n",
    "In this example, the alternative hypothesis is a _union_, viz.,\n",
    "\n",
    "(discriminates in arts and humanities) $\\cup$ (discriminates in sciences) $\\cup$ (discriminates in engineering) $\\cup$ (discriminates in professional schools).\n",
    "\n",
    "Framing a test this way leads to an _intersection-union test_.\n",
    "The null hypothesis is the intersection\n",
    "\n",
    "\\begin{equation*} \n",
    "   H_0 := \\cap_{j=1}^n H_{0j}\n",
    "\\end{equation*}\n",
    "\n",
    "and the alternative is the union\n",
    "\n",
    "\\begin{equation*} \n",
    "   H_1 := \\cup_{j=1}^n H_{0j}^c.\n",
    "\\end{equation*}\n",
    "\n",
    "There can be good reasons for representating a null hypothesis as such an intersection. \n",
    "In the example just mentioned, the applicant pool might be quite different across disciplines, making it hard to judge at the aggregate level whether there is discrimination, while testing within each discipline is more straightforward (that is, _Simpson's Paradox_ can be an issue).\n",
    "\n",
    "Some hypotheses about multivariate distributions can be expressed as the intersection of hypotheses about each dimension separately. For instance, the hypothesis that a $n$-dimensional distribution has zero mean could be represented as \n",
    "\n",
    "(1st component has zero mean) $\\cap$ (2nd component has zero mean) $\\cap$ $\\cdots$ $\\cap$ ($n$th component has zero mean)\n",
    "\n",
    "The alternative is again a union:\n",
    "\n",
    "(1st component has nonzero mean) $\\cup$ (2nd component has nonzero mean) $\\cup$ $\\cdots$ $\\cup$ ($n$th component has nonzero mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe10e7d5-b94a-4a40-9b72-9e6a0ca8c89c",
   "metadata": {},
   "source": [
    "## Combinations of experiments and stratified experiments\n",
    "\n",
    "The same kind of issue arises when combining information from different experiments.\n",
    "For instance, imagine testing whether a drug is effective. We might have several randomized, controlled trials in different places, or a large experiment involving a number of centers, each of which performs its own randomization (i.e., the randomization is stratified).\n",
    "\n",
    "How can we combine the information from the separate (independent) experiments to test the null hypothesis that the drug is ineffective? \n",
    "\n",
    "Again, the overall null hypothesis is \"the drug doesn't help,\" which can be written as an intersection of hypotheses\n",
    "\n",
    "(drug doesn't help in experiment 1) $\\cap$ (drug doesn't help in experiment 2) $\\cap$ $\\cdots$  $\\cap$ (drug doesn't help in experiment $n$),\n",
    "\n",
    "and the alternative can be written as\n",
    "\n",
    "(drug helps in experiment 1) $\\cup$ (drug helps in experiment 2) $\\cup$ $\\cdots$  $\\cup$ (drug helps in experiment $n$),\n",
    "\n",
    "a union."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cf91e2d-78c7-439a-a6c4-df23210d83c5",
   "metadata": {},
   "source": [
    "## Combining evidence\n",
    "\n",
    "Suppose we have a test of each \"partial\" null hypothesis $H_{0j}$. Clearly, if the $P$-value for one of those tests is sufficiently small, that's evidence that the overall null $H_0$ is false.\n",
    "(For instance, we could use Bonferroni's inequality, aka the union bound.)\n",
    "\n",
    "But we will not know in advance whether any individual $P$-value will turn out to be small.\n",
    "Perhaps none of the individual $P$-values turns out to be small, but many are \"not large.\" \n",
    "Is there a way to combine $P$-values across tests to get stronger evidence about $H_0$?\n",
    "(Our procedure for combining $P$-values should be specified _before_ we look at the data,\n",
    "or we are committing _selective inference_, which requires special methods.\n",
    "There are other approaches to combining evidence; \n",
    "in particular, see [E-values and betting scores](./E-values.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "210ecdb5-7a4f-4f22-85f9-d0185e85934f",
   "metadata": {},
   "source": [
    "## $P$-value Combining functions\n",
    "\n",
    "Let $\\lambda$ be a $n$-vector of statistics such that the distribution of $\\lambda_j$\n",
    "if hypothesis $H_{0j}$ is true is known. \n",
    "We assume that smaller values of $\\lambda_j$ are stronger evidence against $H_{0j}$.\n",
    "For instance, $\\lambda_j$ might be the $P$-value of $H_{0j}$ for some test.\n",
    "\n",
    "Consider a function\n",
    "\n",
    "\\begin{equation*}  \\phi: [0, 1]^n \\rightarrow \\Re; \\lambda = (\\lambda_1, \\ldots, \\lambda_n) \\mapsto \\phi(\\lambda)\n",
    "\\end{equation*} \n",
    "with the properties:\n",
    "\n",
    "+ $\\phi$ is non-increasing in every argument, i.e., $\\phi( \\ldots, \\lambda_j, \\ldots) \\ge \\phi(( \\ldots, \\lambda_j', \\ldots)$ if $\\lambda_j \\le \\lambda_j'$, $j = 1, \\ldots, n$.\n",
    "\n",
    "+ $\\phi$ attains its maximum if any of its arguments equals 0.\n",
    "\n",
    "+ $\\phi$ attains its minimum if all of its arguments equal 1.\n",
    "\n",
    "+ for all $\\alpha > 0$, there exist finite functions $\\phi_-(\\alpha)$, $\\phi_+(\\alpha)$ such that if every partial null hypothesis $\\{H_{0j}\\}$ is true, \n",
    "\\begin{equation*} \\mathbb{P} \\{\\phi_-(\\alpha) \\le \\phi(\\lambda) \\le \\phi_+(\\alpha) \\} \\ge 1-\\alpha\\end{equation*}\n",
    "and $[\\phi_-(\\alpha), \\phi_+(\\alpha)] \\subset [\\phi_-(\\alpha'), \\phi_+(\\alpha')]$ if $\\alpha \\ge \\alpha'$.\n",
    "\n",
    "Then we can use $\\phi(\\lambda)$ as the basis of a test of $H_0 = \\cap_{j=1}^n H_{0j}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfabccb7-dbde-4f6b-ba82-e482f352bb54",
   "metadata": {},
   "source": [
    "### Fisher's combining function\n",
    "\n",
    "\\begin{equation*}  \\phi_F(\\lambda) := -2 \\sum_{j=1}^n \\ln(\\lambda_j).\\end{equation*}\n",
    "\n",
    "### Liptak's combining function\n",
    "\n",
    "\\begin{equation*}  \\phi_L(\\lambda) := \\sum_{j=1}^n \\Phi^{-1}(1-\\lambda_j),\\end{equation*}\n",
    "\n",
    "where $\\Phi^{-1}$ is the inverse standard normal CDF.\n",
    "\n",
    "### Tippet's combining function\n",
    "\n",
    "\\begin{equation*}  \\phi_T(\\lambda) := \\max_{j=1}^n (1-\\lambda_j).\\end{equation*}\n",
    "\n",
    "### Direct combination of test statistics\n",
    "\n",
    "\\begin{equation*}  \\phi_D := \\sum_{j=1}^n f_j(\\lambda_j), \\end{equation*}\n",
    "\n",
    "where $\\{ f_j \\}$ are suitable decreasing functions. For instance, if $\\lambda_j$ is the $P$-value for $H_{0j}$ corresponding to some test statistic $T_j$ for which larger values are stronger evidence against $H_{0j}$, we could use $\\phi_D = \\sum_j T_j$.\n",
    "\n",
    "### Combining $E$-values\n",
    "\n",
    "$E$-values are generally much easier to combine than $P$-values, but a combined $E$-value can be turned into a $P$-value.\n",
    "See [E-values and betting scores](./E-values.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98bd090d-f04d-4ccf-9568-515076ed93d1",
   "metadata": {},
   "source": [
    "## Fisher's combining function for independent $P$-values\n",
    "\n",
    "Suppose $H_0$ is true, that $\\lambda_j$ is the $P$-value of $H_{0j}$ for some pre-specified test, that the distribution of $\\lambda_j$ is continuous under $H_{0j}$, and that $\\{ \\lambda_j \\}$ are independent if $H_0$ is true.\n",
    "\n",
    "Then, if $H_0$ is true, $\\{ \\lambda_j \\}$ are IID $U[0,1]$.\n",
    "\n",
    "Under $H_{0j}$, the distribution of $-\\ln \\lambda_j$ is exponential(1):\n",
    "\n",
    "\\begin{equation*} \n",
    "   \\mathbb{P} \\{ -\\ln \\lambda_j \\le x \\} = \\mathbb{P} \\{ \\ln \\lambda_j \\ge -x \\} = \\mathbb{P} \\{ \\lambda_j \\ge e^{-x} \\} = 1 - e^{-x}.\n",
    "\\end{equation*}\n",
    "\n",
    "The distribution of 2 times an exponential is $\\chi_2^2$:\n",
    "the pdf of a chi-square with $k$ degrees of freedom is\n",
    "\n",
    "\\begin{equation*} \n",
    "   \\frac{1}{2^{k/2}\\Gamma(k/2)} x^{k/2-1} e^{-x/2}.\n",
    "\\end{equation*}\n",
    "\n",
    "For $k=2$, this simplifies to $e^{-x/2}/2$, the exponential density scaled by a factor of 2.\n",
    "\n",
    "Thus, under $H_0$, $\\phi_F(\\lambda)$ is the sum of $n$ independent $\\chi_2^2$ random variables. The distribution of a sum of independent chi-square random variables is a chi-square random variable with degrees of freedom equal to the sum of the degrees of freedom of the variables that were added.\n",
    "\n",
    "Hence, under $H_0$,\n",
    "\n",
    "\\begin{equation*} \n",
    "  \\phi_F(\\lambda) \\sim \\chi_{2n}^2,\n",
    "\\end{equation*}\n",
    "\n",
    "the chi-square distribution with $2n$ degrees of freedom.\n",
    "\n",
    "Let $\\chi_{k}^2(\\alpha)$ denote the $1-\\alpha$ quantile of the chi-square distribution\n",
    "with $k$ degrees of freedom.\n",
    "If we reject $H_0$ when\n",
    "\n",
    "\\begin{equation*} \n",
    "   \\phi_F(\\lambda) \\ge \\chi_{2n}^2(\\alpha),\n",
    "\\end{equation*}\n",
    "\n",
    "that yields a significance level $\\alpha$ test of $H_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "240d0fa5-6941-49b0-ba5a-ba3db923a529",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Simulate distribution of Fisher's combining function when all nulls are true\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "import numpy as np\n",
    "from numpy.polynomial import polynomial as P\n",
    "\n",
    "import scipy as sp\n",
    "from scipy.stats import chi2, binom\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display, HTML\n",
    "\n",
    "def plot_fisher_null(n=5, reps=10000):\n",
    "    U = sp.stats.uniform.rvs(size=[reps,n])\n",
    "    vals = np.apply_along_axis(lambda x: -2*np.sum(np.log(x)), 1, U)\n",
    "    fig, ax = plt.subplots(1, 1)\n",
    "    ax.hist(vals, bins=max(int(reps/40), 5), density=True, label=\"simulation\")\n",
    "    mxv = max(vals)\n",
    "    grid = np.linspace(0, mxv, 200)\n",
    "    ax.plot(grid, chi2.pdf(grid, df=2*n), 'r-', lw=3, label='chi-square pdf, df='+str(2*n))\n",
    "    ax.legend(loc='best')\n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9d93a84d-6456-4f66-85a0-f0d6cfa2e13c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2da85cc520864c9181efad2599304a13",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=5, description='n', max=200, min=1), IntSlider(value=10000, description=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_fisher_null(n=5, reps=10000)>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interact(plot_fisher_null, n=widgets.IntSlider(min=1, max=200, step=1, value=5) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1449a9c6-5d4f-4bac-bce7-7282053f0c72",
   "metadata": {},
   "source": [
    "## When $P$-values have atoms\n",
    "\n",
    "A real random variable $X$ is first-order stochastically larger than a real random variable $Y$ if for all $x \\in \\Re$,\n",
    "\n",
    "\\begin{equation*} \n",
    "   \\mathbb{P} \\{ X \\ge x \\} \\ge \\mathbb{P} \\{ Y \\ge x \\},\n",
    "\\end{equation*}\n",
    "with strict inequality for some $x \\in \\Re$.\n",
    "\n",
    "Suppose $\\{\\lambda_j \\}$ for $\\{ H_{0j}\\}$ satisfy\n",
    "\n",
    "\\begin{equation*} \n",
    "   \\mathbb{P} \\{ \\lambda_j \\le p  || H_{0j} \\} \\le p.\n",
    "\\end{equation*}\n",
    "\n",
    "This takes into account the possibility that $\\lambda_j$ does not have a continuous \n",
    "distribution under $H_{0j}$, ensuring that $\\lambda_j$ is still a _conservative_ $P$-value.\n",
    "\n",
    "Since $\\ln$ is monotone, it follows that for all $x \\in \\Re$\n",
    "\n",
    "\\begin{equation*} \n",
    "   \\mathbb{P} \\{ -2 \\ln \\lambda_j \\ge x \\} \\le \\mathbb{P} \\{ -2 \\ln U \\ge x \\}.\n",
    "\\end{equation*}\n",
    "\n",
    "That is, if $\\lambda_j$ does not have a continuous distribution, \n",
    "the a $\\chi_2^2$ variable is stochastically larger than the distribution of $-2\\ln \\lambda_j$.\n",
    "\n",
    "It turns out that $X$ is stochastically larger than $Y$ if and only if\n",
    "there is some probability space on which there exist \n",
    "two random variables, $\\tilde{X}$ and $\\tilde{Y}$ such that $\\tilde{X} \\sim X$,\n",
    "$\\tilde{Y} \\sim Y$, and $\\mathbb{P} \\{\\tilde{X} \\ge \\tilde{Y} \\} = 1$. \n",
    "(See, e.g., Grimmett and Stirzaker,_Probability and Random Processes_, 3rd edition,\n",
    "Theorem 4.12.3.)\n",
    "\n",
    "Let $\\{X_j\\}_{j=1}^n$ be IID $\\chi_2^2$ random variables,\n",
    "and let $Y_j := - 2 \\ln \\lambda_j$, $j=1, \\ldots, n$.\n",
    "\n",
    "Then there is some probability space \n",
    "for which we can define $\\{\\tilde{Y_j}\\}$ and $\\{\\tilde{X_j}\\}$ such that\n",
    "\n",
    "+ $(\\tilde{Y_j})$ has the same joint distribution as $(Y_j)$\n",
    "\n",
    "+ $(\\tilde{X_j})$ has the same joint distribution as $(X_j)$\n",
    "\n",
    "+ $\\tilde{X_j} \\ge \\tilde{Y_j}$ for all $j$ with probability one.\n",
    "\n",
    "Then\n",
    "\n",
    "+ $\\sum_j \\tilde{Y_j}$ has the same distribution as $\\sum_j Y_j = -2 \\sum_j \\ln \\lambda_j$,\n",
    "\n",
    "+ $\\sum_j \\tilde{X_j}$ has the same distribution as $\\sum_j X_j$ (namely, chi-square with $2n$ degrees of freedom),\n",
    "\n",
    "+ $\\sum_j \\tilde X_j  \\ge \\sum_j \\tilde{Y_j}$.\n",
    "\n",
    "That is, \n",
    "\n",
    "\\begin{equation*} \n",
    "  \\mathbb{P} \\left \\{-2 \\sum_j \\ln \\lambda_j \\ge \\chi_{2n}^2(\\alpha) \\right \\} \\le \\alpha.\n",
    "\\end{equation*}\n",
    "\n",
    "Thus, we still get a conservative hypothesis test if one or more of the $p$-values for the\n",
    "partial tests have atoms under their respective null hypotheses $\\{H_{0j}\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f44f371-988d-4b70-857a-7f5bc02b18a2",
   "metadata": {},
   "source": [
    "## Accounting for simulation error in combined $P$-values\n",
    "\n",
    "Suppose that the \"every-permutation\" $P$-value $\\lambda_j$ for $H_{0j}$ is estimated by $b_j$ simulations.\n",
    "How can we take the uncertainty of the simulation estimate into account?\n",
    "\n",
    "Here, we will pretend that the simulation itself is perfect: that the PRNG generates true IID $U[0,1]$ variables, that pseudo-random integers on $\\{0, 1, \\ldots, N\\}$ really are equally likely, and that pseudo-random samples or permutations really are equally likely, etc. \n",
    "\n",
    "The error we are accounting for is not the imperfection of the PRNG or other algorithms, just the uncertainty due to approximating a theoretical probability $\\lambda_j$ by an estimate via (perfect) simulation. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf8d1773-5afd-4464-a48f-7f8713b4bcdc",
   "metadata": {},
   "source": [
    "### A crude approach: simultaneous one-sided upper confidence bounds for every $\\lambda_j$\n",
    "\n",
    "Suppose we find, for each $j$, an upper confidence bound for $\\lambda_j$ (the \"every-permutation\" $P$-value in stratum $j$),\n",
    "for instance, by inverting binomial tests based on $\\# \\{k > 0: T(\\pi_k(x_0)) \\ge T(x_0) \\}$.\n",
    "\n",
    "Since $\\phi$ is monotonic in every coordinate, the upper confidence confidence bounds \n",
    "for $\\{ \\lambda_j \\}$ imply a lower confidence bound for $\\phi(\\lambda)$, which translates to an upper confidence bound for the combined $P$-value.\n",
    "\n",
    "What is the confidence level of the bound on the combined $P$-value? \n",
    "If the $P$-value estimates are independent, the joint coverage probability of a set of $n$ independent confidence bounds with confidence level $\\alpha$ is $1-(1-\\alpha)^n$, as we shall show.\n",
    "\n",
    "Let $A_j$ denote the event that the upper confidence bound for $\\lambda_j$ is greater than or equal to $\\lambda_j$, and suppose $\\mathbb{P} \\{A_j\\} = 1-\\alpha_j$.\n",
    "\n",
    "Regardless of the dependence among the events $\\{A_j \\}$, the chance that all of the confidence bounds cover their corresponding parameters can be bounded using Bonferroni's inequality:\n",
    "\n",
    "\\begin{equation*} \n",
    "\\mathbb{P} \\{ \\cap_j A_j \\} = 1 - \\mathbb{P} \\{ \\cup_j A_j^c \\} \\ge 1 - \\sum_j \\mathbb{P} \\{A_j^c \\} \n",
    "   = 1 - \\sum_j (1- \\mathbb{P} \\{A_j \\}) = 1 - \\sum_j \\alpha_j.\n",
    "\\end{equation*}\n",
    "\n",
    "If $\\{A_j\\}$ are independent, \n",
    "\n",
    "\\begin{equation*} \n",
    "\\mathbb{P} \\{ \\cap_j A_j \\} = \\prod_j \\mathbb{P} \\{ A_j \\} = \\prod_j (1-\\alpha_j).\n",
    "\\end{equation*}\n",
    "\n",
    "Both of those expressions tend to get small quickly as $n$ gets large;\n",
    "bounding $\\phi(\\lambda)$ by bounding the components of $\\lambda$ is inefficient.\n",
    "\n",
    "A far more efficient approach is to use the exact $P$-value based on the related randomized test.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66cf5a6f-560f-4353-910a-bfcd25704a80",
   "metadata": {},
   "source": [
    "### A sharper approach: use a related randomized test\n",
    "\n",
    "See [Simulation and permutation test $P$-values](./permute-simulate-p.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "975c0ee9-5de7-4d9b-b520-57816161fa92",
   "metadata": {},
   "source": [
    "## Dependent tests\n",
    "\n",
    "If $\\{ \\lambda_j \\}_{j=1}^n$ are dependent, the distribution of $\\phi_F(\\lambda)$ is no longer chi-square when the null hypotheses are true.\n",
    "Nonetheless, one can calibrate a test based on Fisher's combining function (or any other combining function) by simulation.\n",
    "This is commonly used in multivariate permutation tests involving dependent partial tests\n",
    "using \"lockstep\" permutations.\n",
    "\n",
    "See, e.g., Pesarin, F. and L. Salmaso, 2010. _Permutation Tests for Complex Data: Theory, Applications and Software_, Wiley, 978-0-470-51641-6."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e95a1366-8810-4ef8-98d0-68c7049738fe",
   "metadata": {},
   "source": [
    "We shall illustrate how the approach can be used to construct nonparametric multivariate tests from univariate tests to address for the two-sample problem (i.e., is there evidence that two samples come from different populations, or is it plausible\n",
    "that they are a single population randomly divided into two groups?).\n",
    "This is equivalent to testing whether treatment has an effect in a controlled, randomized study\n",
    "in which the subjects who receive treatment are a simple random sample of the study group,\n",
    "using the Neyman model for causal inference.\n",
    "(The null hypothesis is that treatment makes no difference whatsoever: each subject's response would\n",
    "be the same whether that subject was assigned to treatment or to control, and without regard for the assignment of\n",
    "other subjects to treatment or control.)\n",
    "\n",
    "We have $N$ subjects of whom $N_t$ are treated and $N_c$ are controls. \n",
    "Each subject has a vector of $n$ measurements. \n",
    "For each of these $n$ \"dimensions\" we have a test statistic $T_j(\\cdot, \\cdot)$ \n",
    "(for instance, the difference between the mean of the treated and the mean of the controls on that dimension or the Smirnov distance between the two groups--but we could use something else, and we don't have to use the same test statistic for different \"dimensions\"). \n",
    "\n",
    "Each $T_j$ takes the responses of the treated and the controls on dimension $j$ and yields a number. \n",
    "We will assume that larger values of $T_j$ are stronger evidence against the null hypothesis that for dimension $j$ treatment doesn't make any difference. \n",
    "\n",
    "Let $T$ denote the whole $n$-vector of test statistics. \n",
    "Let $t(0)$ denote the observed value of the \n",
    "$n$-vector of test statistics for the original data.\n",
    "\n",
    "Now consider randomly re-labelling $N_t$ of the $N$ subjects as treated and the remaining $N_c$ as controls\n",
    "by simple random sampling, so that all subsets of size $N_t$ are equally likely to be labeled \"treated.\"\n",
    "Each re-labelling carries the subject's entire $n$-vector of responses with it: \n",
    "the dimensions are randomized \"in lockstep.\"\n",
    "\n",
    "Let $t(k)$ denote the observed $n$-vector of test statistics for the $k$th random allocation. \n",
    "We randomize the data $K$ times in all, each yielding a $n$-vector $t(k)$ of observed values of the test statistics. \n",
    "This gives $K+1$ permutations in all, including the original data (the zeroth random allocation), for which the \n",
    "vector is $t(0)$. \n",
    "Let $t_j(k)$ denote the test statistic for dimension $j$ for the $k$th permutation.\n",
    "\n",
    "We now transform the $n$ by $(K+1)$ matrix $[t_j(k)]_{j=1}^n{}_{k=0}^K$ to a corresponding matrix of $P$-values for the univariate permutation tests:\n",
    "For $j=1, \\ldots , n$ and $k = 0, \\ldots , K$, define\n",
    "\\begin{equation*} \n",
    "   P_j(k) := \\frac{\\#\\{ \\ell \\in \\{0, \\ldots, K \\} : t_j(\\ell) \\ge t_j(k) \\}}{K+1}.  \n",
    "\\end{equation*}\n",
    "This is the simulated upper tail probability of the $k$th observed value of the $j$th test statistic \n",
    "under the null hypothesis. \n",
    "\n",
    "Think of the values of $P_j(k)$ as a matrix. \n",
    "Each column corresponds to a random permutation of the original data (the 0th column corresponds to the original data); each row corresponds to a dimension of measurement; each entry is a number between $1/(K+1)$ and 1.\n",
    "\n",
    "Now apply Fisher's combining function $\\phi_F$ (or Tippett's, or Stouffer's, or anything else) to each column of \n",
    "$n$ numbers. \n",
    "That gives $K+1$ values $\\{f(k): k=0, \\ldots , K\\}$, one for each permutation of the data. \n",
    "The overall \"Non-Parametric Combination of tests\" (NPC) $P$-value is\n",
    "\\begin{equation*} \n",
    " P_{\\mbox{NPC}} := \\frac{\\#\\{ k \\in \\{0, \\ldots, K\\} : f(k) \\ge f(0) \\}}{K+1}. \n",
    "\\end{equation*}\n",
    "This is the simulated lower tail probability of the observed value of the combining function\n",
    "applied to the partial tests under the randomization.\n",
    "\n",
    "Ultimately, this whole thing is just a univariate permutation test that uses a complicated test statistic \n",
    "that assigns one number to each permutation of the multivariate data.\n",
    "\n",
    "Also see [the permute Python package](http://statlab.github.io/permute/), and in particular,\n",
    "[the npc function](https://github.com/statlab/permute/blob/main/permute/npc.py)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "372b6a26-286d-48f4-aced-7ee36ffd6ea2",
   "metadata": {},
   "source": [
    "Pesarin and Salmaso (2010) give conditions under which nonparametric combinations of tests (NPC) yield consistent tests (tests that reject with probability approaching 1 as the sample size increases or as the number of dimensions increases, if the null is false).\n",
    "They show how to adjust NPC for multiple comparisons and how to use NPC with categorical data, repeated measurements,\n",
    "survival data (with censoring), shape data, propensity scores, etc., and how to use NPC to to test against alternatives like stochastic ordering.\n",
    "\n",
    "Here are some ways NPC can be used:\n",
    "\n",
    "+ Construct two-sided tests by combining one-sided tests\n",
    "+ Construct multivariate tests by combining univariate tests\n",
    "+ Construct union-intersection tests\n",
    "+ Construct multi-aspect tests from univariate tests, e.g., combining a test for a difference in means with a test for a difference in variances\n",
    "+ Testing for non-inferiority (\"bioequivalence\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5003c4aa-08c7-4a95-a930-68c73dfe975c",
   "metadata": {},
   "source": [
    "## Stratified Permutation Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60caed34-8041-4289-a2ff-210f157a9a31",
   "metadata": {},
   "source": [
    "Two examples: \n",
    "\n",
    "+ Boring, A., K. Ottoboni, and P.B. Stark, 2016. Student Evaluations of Teaching (Mostly) Do Not Measure Teaching Effectiveness, _ScienceOpen_, doi 10.14293/S2199-1006.1.SOR-EDU.AETBZC.v1\n",
    "\n",
    "+ Hessler, M.,  D.M. Pöpping, H. Hollstein, H. Ohlenburg, P.H. Arnemann, C. Massoth, L.M. Seidel, A. Zarbock & M. Wenk, 2018. Availability of cookies during an academic course session affects evaluation of teaching, _Medical Education, 52_, 1064–1072. doi 10.1111/medu.13627\n"
   ]
  }
 ],
 "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
}