Chapter 11 - Hypothesis Testing and p-values.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Hypothesis Testing and p-values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we partition the parameters space $\\Theta$ into two disjoint sets $\\Theta_0$ and $\\Theta_1$ and we wish to test\n",
    "\n",
    "$$\n",
    "H_0: \\theta \\in \\Theta_0\n",
    "\\quad \\text{versus} \\quad\n",
    "H_1: \\theta \\in \\Theta_1\n",
    "$$\n",
    "\n",
    "We call $H_0$ the **null hypothesis** and $H_1$ the **alternative hypothesis**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $X$ be a random variable and let $\\mathcal{X}$ be the range of $X$.  We test a hypothesis by finding an appropriate subset of outcomes $R \\subset \\mathcal{X}$ called the **rejection region**.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "X \\in R & \\Longrightarrow  \\text{reject } H_0 \\\\\n",
    "X \\notin R & \\Longrightarrow  \\text{retain (do not reject) } H_0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Usually the rejection region is of form\n",
    "\n",
    "$$ R = \\bigg\\{x: T(x) > c \\bigg\\}$$\n",
    "\n",
    "where $T$ is a **test statistic** and $c$ is a **critical value**.  The problem in hypothesis testing is to find an appropriate test statistic $T$ and and appropriate cutoff value $c$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|            | Retain Null   | Reject Null  |\n",
    "|------------|---------------|--------------|    \n",
    "|$H_0$ true  | $\\checkmark$  | Type I error |\n",
    "|$H_1$ true  | Type II error | $\\checkmark$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **power function** of a test with rejection region $R$ is defined by\n",
    "\n",
    "$$\\beta(\\theta) = \\mathbb{P}_\\theta(X \\in R)$$\n",
    "\n",
    "The **size** of a test is defined to be\n",
    "\n",
    "$$\\alpha = \\sup_{\\theta \\in \\Theta_0} \\beta(\\theta)$$\n",
    "\n",
    "A test is said to have **level $\\alpha$** if its size is less than or equal to $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- A hypothesis of the form $\\theta = \\theta_0$ is called a **simple hypothesis**.\n",
    "- A hypothesis of the form $\\theta > \\theta_0$ or $\\theta < \\theta_0$ is called a **composite hypothesis**.\n",
    "- A test of the form $H_0 : \\theta = \\theta_0$ versus $H_1 : \\theta \\neq \\theta_0$ is called a **two-sided test**.\n",
    "- A test of the form $H_0 : \\theta \\leq \\theta_0$ versus $H_1: \\theta > \\theta_0$ or $H_0: \\theta \\geq \\theta_0$ versus $H_1: \\theta < \\theta_0$ is called a **one-sided test**.\n",
    "\n",
    "The most common tests are two-sided."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finding most powerful tests is hard and, in many cases, most powerful tests don't even exist.  We will just consider three widely used tests: the Wald test, the $\\chi^2$ test, and the permutation test.  A fourth test, the likelihood ratio test, is discussed in the appendix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.1 The Wald Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider testing\n",
    "\n",
    "$$ H_0: \\theta = \\theta_0\n",
    "\\quad \\text{versus} \\quad\n",
    "H_1: \\theta \\neq \\theta_0$$\n",
    "\n",
    "Assume that $\\hat{\\theta}$ is asymptotically Normal:\n",
    "\n",
    "$$ \\frac{\\sqrt{n}(\\hat{\\theta} - \\theta_0)}{\\hat{\\text{se}}} \\leadsto N(0, 1) $$\n",
    "\n",
    "The size $\\alpha$ **Wald test** is: reject $H_0$ when $|W| > z_{\\alpha/2}$, where\n",
    "\n",
    "$$ W = \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.4**.  Asymptotically, the Wald test has size $\\alpha$, that is,\n",
    "\n",
    "$$ \\mathbb{P}_{\\theta_0} \\left(|Z| > z_{\\alpha/2} \\right) \\rightarrow \\alpha$$\n",
    "\n",
    "as $n \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  Under $\\theta = \\theta_0$, $(\\hat{\\theta} - \\theta) / \\text{se} \\leadsto N(0, 1)$, so the probability of rejecting when the null hypothesis $\\theta = \\theta_0$ is true is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}_{\\theta_0}(|W| > z_{\\alpha / 2}) &= \\mathbb{P}_{\\theta_0} \\left(\\frac{|\\hat{\\theta} - \\theta_0|}{\\text{se}} > z_{\\alpha/2} \\right) \\\\\n",
    "& \\rightarrow \\mathbb{P}_{\\theta_0}(| N(0, 1) | > z_{\\alpha/2}) \\\\\n",
    "& = \\alpha\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most texts define the Wald test with the standard error computed at $\\theta = \\theta_0$ rather than at the estimated value $\\hat{\\theta}$.  Both versions are valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.6**.  Suppose that the true value of $\\theta$ is $\\theta_* \\neq \\theta_0$.  The power $\\beta(\\theta_*)$, the probability of correctly rejecting the null hypothesis, is given (approximately) by\n",
    "\n",
    "$$ 1 - \\Phi \\left(\\frac{\\theta_0 - \\theta_*}{\\hat{\\text{se}}} + z_{\\alpha/2} \\right) + \\Phi \\left(\\frac{\\theta_0 - \\theta_*}{\\hat{\\text{se}}} - z_{\\alpha/2} \\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.10**.  The size $\\alpha$ Wald test rejects $H_0: \\theta = \\theta_0$ versus $H_1: \\theta \\neq \\theta_0$ if and only if $\\theta_0 \\notin C$ where\n",
    "\n",
    "$$ C = \\left(\\hat{\\theta} - \\hat{\\text{se}} z_{\\alpha/2}, \\; \\hat{\\theta} + \\hat{\\text{se}} z_{\\alpha / 2} \\right) $$\n",
    "\n",
    "Thus, testing the hypothesis is equivalent to checking whether the null value is in the confidence interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.2 p-values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that for every $\\alpha \\in (0, 1)$ we have a size $\\alpha$ test with rejection region $R_\\alpha$.  Then,\n",
    "\n",
    "$$ \\text{p-value} = \\inf \\Big\\{ \\alpha : T(X^n) \\in R_\\alpha \\Big\\} $$\n",
    "\n",
    "That is, the p-value is the smallest level at which we can reject $H_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Informally, the p-value is a measure of the evidence against $H_0$: the smaller the p-value, the stronger the evidence against $H_0$.  Typically, researchers use the following evidence scale:\n",
    "\n",
    "| p-value    | evidence                            |\n",
    "|------------|-------------------------------------|\n",
    "| under 1%   | very strong evidence against $H_0$  |\n",
    "| 1% to 5%   | strong evidence against $H_0$       |\n",
    "| 5% to 10%  | weak evidence against $H_0$         |\n",
    "| over 10%   | little or no evidence against $H_0$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Warning**: a large p-value is not strong evidence in favor of $H_0$. A large p-value can occur for two reasons:\n",
    "- $H_0$ is true, or\n",
    "- $H_0$ is false but the test has low power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The p-value is not the probability that the null hypothesis is true**.  We discuss quantities like $\\mathbb{P}(H_0 | \\text{Data})$ in the chapter on Bayesian inference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.12**.  Suppose that the size $\\alpha$ test is of the form\n",
    "\n",
    "$$ \\text{reject } H_0 \\text{ if and only if } T(X^n) \\geq c_\\alpha$$\n",
    "\n",
    "Then,\n",
    "\n",
    "$$ \\text{p-value} = \\sup_{\\theta \\in \\Theta_0} \\mathbb{P}_\\theta \\left(T(X^n) \\geq T(x^n) \\right)$$\n",
    "\n",
    "In words, the p-value is the probability (under $H_0$) of observing a value of the test statistic as or more extreme than what was actually observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a Wald test, $W$ has an approximate $N(0, 1)$ distribution under $H_0$.  Hence, the p-value is\n",
    "\n",
    "$$ \\text{p-value} \\approx \\mathbb{P}(|Z| > |w|) = 2\\mathbb{P}(Z < -|w|) = 2 \\Phi(Z < -|w|) $$\n",
    "\n",
    "where $Z \\sim N(0, 1)$ and $w = (\\hat{\\theta} - \\theta_0) / \\hat{\\text{se}}$ is the observed value of the test statistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.13**. If the test statistic has a continuous distribution, then under $H_0: \\theta = \\theta_0$ the p-value has an $\\text{Uniform}(0, 1)$ distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.3 The $\\chi^2$ distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $Z_1, \\dots, Z_k$ be independent, standard normals. Let $V = \\sum_{i=1}^k Z_i^2$.  Then we say that $V$ has a $\\chi^2$ distribution with $k$ degrees of freedom, written $V \\sim \\chi^2_k$.  The probability density of $V$ is\n",
    "\n",
    "$$ f(v) = \\frac{v^{(k/2) - 1}e^{-v/2}}{2^{k/2} \\Gamma(k / 2)} $$\n",
    "\n",
    "for $v > 0$.  It can be shown that $\\mathbb{E}(V) = k$ and $\\mathbb{V}(k) = 2k$.  We define the upper $\\alpha$ quantile $\\chi^2_{k, \\alpha} = F^{-1}(1 - alpha)$ where $F$ is the CDF.  That is, $\\mathbb{P}(\\chi^2_k > \\chi^2_{k, \\alpha}) = \\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.4 Pearson's $\\chi^2$ Test for Multinomial Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that $X = (X_1, \\dots, X_n)$ has a multinomial distribution if\n",
    "\n",
    "$$ f(x_1, \\dots, x_k; p) = \\begin{pmatrix}\n",
    "n \\\\\n",
    "x_1 \\dots x_k\n",
    "\\end{pmatrix} p_1^{x_1} \\dots p_k^{x_k}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "n \\\\\n",
    "x_1 \\dots x_k\n",
    "\\end{pmatrix} = \n",
    "\\frac{n!}{x_1! \\dots x_k!}\n",
    "$$\n",
    "\n",
    "The MLE of $p$ is $\\hat{p} = (\\hat{p_1}, \\dots, \\hat{p_k}) = \\left(X_1 / n, \\dots, X_k / n\\right)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $(p_{01}, \\dots, p_{0k})$ be some fixed test of probabilities and suppose we want to test\n",
    "\n",
    "$$ H_0: (p_1, \\dots, p_k) = (p_{01}, \\dots, p_{0k})\n",
    "\\quad \\text{versus} \\quad\n",
    "H_1: (p_1, \\dots, p_k) \\neq (p_{01}, \\dots, p_{0k})$$\n",
    "\n",
    "**Pearson's $\\chi^2$ statistic** is\n",
    "\n",
    "$$ T = \\sum_{j=1}^k \\frac{(X_j - np_{0j})^2}{np_{0j}} = \\sum_{j=1}^k \\frac{(O_j - E_j)^2}{E_j}$$\n",
    "\n",
    "where $O_j = X_j$ is the observed data and $E_j = \\mathbb{E}(X_j) = np_{0j}$ is the expected value of $X_j$ under $H_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.15**.  Under $H_0$, $T \\leadsto \\chi^2_{k - 1}$.  Hence the test:  reject $H_0$ if $T > \\chi^2_{k - 1, \\alpha}$ has asymptotic level $\\alpha$.  The p-value is $\\mathbb{P}(\\chi^2_k > t)$ where $t$ is the observed value of the test statistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.5 The Permutation Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $X_1, \\dots, X_m \\sim F_X$ and $Y_1, \\dots, Y_n \\sim F_Y$ are two independent samples and $H_0$ is the hypothesis that both samples are identically distributed.  More precisely we are testing\n",
    "\n",
    "$$ H_0: F_X = F_Y \\quad \\text{versus} \\quad F_1: F_X \\neq F_Y$$\n",
    "\n",
    "Let $T(x_1, \\dots, x_m, y_1, \\dots, y_n)$ be some test statistic, for example\n",
    "\n",
    "$$ T(X_1, \\dots, X_m, Y_1, \\dots, Y_n) = | \\overline{X}_m - \\overline{Y}_n |$$\n",
    "\n",
    "Let $N = m + n$ and consider forming all $N!$ permutations of the data $X_1, \\dots, X_m, Y_1, \\dots, Y_n$.  For exach permutation, compute the test statistic $T$.  Denote these values by $T_1, \\dots, T_{N!}$.  Under the null hypothesis, each of these values is equally likely.  The distribution $P_0$ that puts mass $1 / N!$ on each $T_j$ is called the **permutation distribution** of $T$.  Let $t_\\text{obs}$ be the observed value of the test statistic.  Assuming we reject when $T$ is large, the p-value is\n",
    "\n",
    "$$ \\text{p-value} = \\mathbb{P}_0(T > t_\\text{obs}) = \\frac{1}{N!} \\sum_{i=1}^{N!} I(T_j > t_\\text{obs}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually it is not practical to evaluate all $N!$ permutations.  We can approximate the p-value by sampling randomly from the set of permutations.  The fraction of times $T > t_\\text{obs}$ among these samples approximates the p-value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Algorithm for the Permutation Test**\n",
    "\n",
    "1. Compute the observed value of the test statistic $t_\\text{obs} = T(X_1, \\dots, X_m, Y_1, \\dots, Y_n)$.\n",
    "\n",
    "2. Randomly permute the data. Compute the statistic again using the permuted data.\n",
    "\n",
    "3. Repeat the previous step $B$ times and let $T_1, \\dots, T_B$ denote the resulting values.\n",
    "\n",
    "4. The approximate p-value is\n",
    "\n",
    "$$ \\frac{1}{B} \\sum_{i=1}^B I(T_j > t_\\text{obs}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In large samples, the permutation test usually gives similar results to a test that is based on large sample theory.  The permutation test is thus most useful for small samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.6 Multiple Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider $m$ hypothesis tests:\n",
    "\n",
    "$$ H_{0i} \\quad \\text{versus} \\quad H_{1i}, \\quad i = 1, \\dots, m$$\n",
    "\n",
    "and let $P_1, \\dots, P_m$ denote $m$ p-values for these tests."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The Bonferroni Method**:  Given p-values $P_1, \\dots, P_m$, reject null hypothesis $H_{0i}$ if $P_i < \\alpha / m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.19**.  Using the Bonferroni method, the probability of falsely rejecting any null hypothesis is less than or equal to $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  Let $R$ be the event that at least one null hypothesis is falsely rejected.  Let $R_i$ be the event that the $i$-th null hypothesis is falsely rejected.  Recall that if $A_1, \\dots, A_k$ are events then $\\mathbb{P} \\left( \\bigcup_{i=1}^k A_i \\right) \\leq \\sum_{i=1}^k \\mathbb{P}(A_i)$.  Hence,\n",
    "\n",
    "$$\\mathbb{P}(R) = \\mathbb{P}\\left( \\bigcup_{i=1}^k R_i \\right) \\leq \\sum_{i=1}^k \\mathbb{P}(R_i) = \\sum_{i=1}^k \\frac{\\alpha}{m} = \\alpha$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bonferroni Method is very conservative because it is trying to make it unlikely that you would make even one false rejection.  Sometimes, a more reasonable idea is to control the **false discovery rate** (FDR) which is defined as the mean number of false rejections divided by the number of rejections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|           |  $H_0$ not rejected | $H_0$ rejected  | Total |\n",
    "|-----------|---------------------|-----------------|-------|\n",
    "|$H_0$ true | $U$                 | $V$             | $m_0$ |\n",
    "|$H_0$ false| $T$                 | $S$             | $m_1$ |\n",
    "|Total      | $m - R$             | $R$             | $m$   |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the **false discovery proportion** (FDP):\n",
    "\n",
    "$$ \\text{FDP} = \\begin{cases}\n",
    "V / R & \\text{if } R > 0\\\\\n",
    "0     & \\text{if}  R = 0\n",
    "\\end{cases}$$\n",
    "\n",
    "The FDP is the proportion of rejections that are incorrect.  Next define $\\text{FDR} = \\mathbb{E}(\\text{FDP})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The Benjamini-Hochberg (BH) Method**\n",
    "\n",
    "1. Let $P_{(1)} < \\cdots < P_{(m)}$ denote the ordered p-values.\n",
    "\n",
    "2. Define\n",
    "\n",
    "$$\\ell_i = \\frac{i \\alpha}{C_m m},\n",
    "\\quad \\text{and} \\quad\n",
    "R = \\max \\bigg\\{ i: P_{(i)} < \\ell_i \\bigg\\}$$\n",
    "\n",
    "where $C_m$ is defined to be 1 if the p-values are independent and $C_m = \\sum_{i=1}^m (1/i)$ otherwise.\n",
    "\n",
    "3. Let $t = P_{(R)}$; we call $t$ the **BH rejection threshold**.\n",
    "\n",
    "4. Reject all null hypothesis $H_{0i}$ for which $P_i \\leq t$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.21 (Benjamini and Hocheberg)**.  If the procedure above is applied, then regardless of how many nulls are true and regardless of the distribution of p-values when the null hypothesis is false,\n",
    "\n",
    "$$\\text{FDR} = \\mathbb{E}(\\text{FDP}) \\leq \\frac{m_0}{m} \\alpha \\leq \\alpha $$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.7 Technical Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 11.7.1 The Neyman-Pearson Lemma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.23 (Neyman-Pearson)**.  Suppose we test $H_0: \\theta = \\theta_0$ versus $H_1: \\theta = \\theta_1$. Let\n",
    "\n",
    "$$ T = \\frac{\\mathcal{L}(\\theta_1)}{\\mathcal{L}(\\theta_0)} = \\frac{\\prod_{i=1}^n f(x_i; \\theta_1)}{\\prod_{i=1}^n f(x_i; \\theta_0)}$$\n",
    "\n",
    "Suppose we reject $H_0$ when $T > k$.  If we choose $k$ so that $\\mathbb{P}_{\\theta_0}(T > k) = \\alpha$ then this test is the most powerful, size $\\alpha$ test.  That is, among all tests with size $\\alpha$, this test maximizes power $\\beta(\\theta_1)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 11.7.2 Power of the Wald Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof of Theorem 11.6**.  \n",
    "\n",
    "Let $Z \\sim N(0, 1)$.  Then,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{Power} &= \\beta(\\theta_*) \\\\\n",
    "&= \\mathbb{P}_{\\theta_*}(\\text{Reject } H_0) \\\\\n",
    "&= \\mathbb{P}_{\\theta_*}\\left( \\frac{|\\hat{\\theta} - \\theta_0|}{\\hat{\\text{se}}} > z_{\\alpha/2} \\right) \\\\\n",
    "&= \\mathbb{P}_{\\theta_*}\\left( \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} > z_{\\alpha/2} \\right) \n",
    "+ \\mathbb{P}_{\\theta_*}\\left( \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} < -z_{\\alpha/2} \\right) \\\\\n",
    "&= \\mathbb{P}_{\\theta_*}(\\hat{\\theta} > \\theta_0 + \\hat{\\text{se}} z_{\\alpha/2})\n",
    "+ \\mathbb{P}_{\\theta_*}(\\hat{\\theta} < \\theta_0 - \\hat{\\text{se}} z_{\\alpha/2}) \\\\\n",
    "&= \\mathbb{P}_{\\theta_*}\\left( \\frac{\\hat{\\theta} - \\theta_*}{\\hat{\\text{se}}} > \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} + z_{\\alpha/2} \\right) \n",
    "+ \\mathbb{P}_{\\theta_*}\\left( \\frac{\\hat{\\theta} - \\theta_*}{\\hat{\\text{se}}} < \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} - z_{\\alpha/2} \\right) \\\\\n",
    "& \\approx \\mathbb{P}\\left(Z > \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} + z_{\\alpha/2} \\right) \n",
    "+ \\mathbb{P}\\left(Z < \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} - z_{\\alpha/2} \\right) \\\\\n",
    "&= 1 - \\Phi\\left( \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} + z_{\\alpha/2} \\right) \n",
    "+ \\Phi\\left( \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} - z_{\\alpha/2} \\right)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 11.7.3 The t-test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test $H_0: \\mu = \\mu_0$ where $\\mu$ is the mean, we can use the Wald test. When the data is assumed to be Normal and the sample size is small, it is common instead to use the **t-test**. A random variable $T$ as a *t-distribution with $k$ degrees of freedom* if it has density\n",
    "\n",
    "$$ f(t) = \\frac{\\Gamma\\left(\\frac{k+1}{2}\\right)}{\\sqrt{k \\pi} \\Gamma\\left(\\frac{k}{2}\\right) \\left(1 + \\frac{t^2}{k}\\right)^{(k+1)/2}}$$\n",
    "\n",
    "When the degrees of freedom $k \\rightarrow \\infty$, this tends to a Normal distribution.  When $k = 1$ it reduces to a Cauchy distribution.\n",
    "\n",
    "Let $X_1, \\dots, X_n \\sim N(\\mu, \\sigma^2)$ where $\\theta = (\\mu, \\sigma^2)$ are both unknown.  Suppose we want to test $\\mu = \\mu_0$ versus $\\mu \\neq \\mu_0$.  Let\n",
    "\n",
    "$$T = \\frac{\\sqrt{n}(\\overline{X}_n - \\mu_0)}{S_n}$$\n",
    "\n",
    "where $S_n^2$ is the sample variance. For large samples $T \\approx N(0, 1)$ under $H_0$.  The exact distribution of $T$ under $H_0$ is $t_{n-1}$.  Hence if we reject when $|T| > t_{n-1, \\alpha/2}$ then we get a size $\\alpha$ test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 11.7.4 The Likelihood Ratio Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\theta = (\\theta_1, \\dots, \\theta_q, \\theta_{q+1}, \\dots, \\theta_r)$ and suppose that $\\Theta_0$ consists of all parameter values $\\theta$ such that $(\\theta_{q+1}, \\dots, \\theta_r) = (\\theta_{0, q+1}, \\dots, \\theta_{0, r})$.\n",
    "\n",
    "Define the **likelihood ratio statistic** by\n",
    "\n",
    "$$ \\lambda \n",
    "= 2 \\log \\left(  \\frac{\\sup_{\\theta \\in \\Theta} \\mathcal{L}(\\theta)}{\\sup_{\\theta \\in \\Theta_0} \\mathcal{L}(\\theta)} \\right) \n",
    "= 2 \\log \\left(  \\frac{\\mathcal{L}(\\hat{\\theta})}{\\mathcal{L}(\\hat{\\theta_0})} \\right) $$\n",
    "\n",
    "where $\\hat{\\theta}$ is the MLE and $\\hat{\\theta_0}$ is the MLE when $\\theta$ is restricted to lie in $\\Theta_0$.  \n",
    "\n",
    "The **likelihood ratio test** is: reject $H_0$ when $\\lambda(x^n) > \\chi^2_{r-q, \\alpha}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 11.25**.  Under $H_0$,\n",
    "\n",
    "$$ 2 \\log \\lambda (x^n) \\leadsto \\chi^2_{r - q}$$\n",
    "\n",
    "Hence, asymptotically, the LR test level is $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.9 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.1**.  Prove Theorem 11.13."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the test statistic has a continuous distribution, then under $H_0: \\theta = \\theta_0$ the p-value has an $\\text{Uniform}(0, 1)$ distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  Let $T$ be the test statistic, and $F_T$ be its CDF under $H_0$.  Then\n",
    "\n",
    "$$ F_\\text{p-value}(p) = \\mathbb{P}(\\text{p-value} < p) = \\mathbb{P}(F(T) < p) = \\mathbb{P}(T < F^{-1}(p)) = F(F^{-1}(p)) = p$$\n",
    "\n",
    "so the CDF for the p-value is the same as the CDF for the $\\text{Uniform}(0, 1)$ distribution, therefore the p-value follows a $\\text{Uniform}(0, 1)$ distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.2**.  Prove Theorem 11.10."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The size $\\alpha$ Wald test rejects $H_0: \\theta = \\theta_0$ versus $H_1: \\theta \\neq \\theta_0$ if and only if $\\theta_0 \\notin C$ where\n",
    "\n",
    "$$ C = \\left(\\hat{\\theta} - \\hat{\\text{se}} z_{\\alpha/2}, \\; \\hat{\\theta} + \\hat{\\text{se}} z_{\\alpha / 2} \\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "$\\theta_0 \\in C$ if and only if\n",
    "\n",
    "$$\n",
    "\\hat{\\theta} - \\hat{\\text{se}} z_{\\alpha/2} < \\theta_0 < \\hat{\\theta} + \\hat{\\text{se}} z_{\\alpha/2} \\\\\n",
    "\\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} - z_{\\alpha/2} < 0 < \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} +  z_{\\alpha/2} \\\\\n",
    "- z_{\\alpha/2} < \\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} < z_{\\alpha/2} \\\\\n",
    "|Z| > z_{\\alpha/2}\n",
    "$$\n",
    "\n",
    "which is the criteria for the Wald test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.3**.  Let $X_1, \\dots, X_n \\sim \\text{Uniform}(0, \\theta)$ and let $Y = \\max \\{ X_1, \\dots, X_n \\}$.  We want to test\n",
    "\n",
    "$$ H_0: \\theta = 1/2 \\quad \\text{versus} \\quad H_1: \\theta > 1/2 $$\n",
    "\n",
    "The Wald test is not appropriate since $Y$ does not converge to a Normal.  Suppose we decide to test this hypothesis by rejecting $H_0$ when $Y > c$.\n",
    "\n",
    "**(a)** Find the power function.\n",
    "\n",
    "**(b)** What choice of $c$ will make the size of the test 0.05?\n",
    "\n",
    "**(c)** In a sample of size $n = 20$ with $Y = 0.48$ what is the p-value? What conclusion about $H_0$ would you make?\n",
    "\n",
    "**(d)** In a sample of size $n = 20$ with $Y = 0.52$ what is the p-value? What conclusion about $H_0$ would you make?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**\n",
    "\n",
    "**(a)** The power function for this test is\n",
    "\n",
    "$$\\beta(\\theta) = \\mathbb{P}_\\theta(X \\in R) = \\mathbb{P}_\\theta(Y > c) = 1 - \\mathbb{P}_\\theta(Y \\leq c) = 1 - \\prod_{i=1}^n \\mathbb{P}_\\theta(X_i \\leq c) = 1 - \\left(\\text{clip}\\left(\\frac{c}{\\theta}\\right)\\right)^n$$\n",
    "\n",
    "where we limit $c/\\theta$ to the $[0, 1]$ interval, $$\\text{clip}(t) = \\begin{cases} 1 & \\text{if } t > 1 \\\\ t & \\text{if } 0 \\leq t \\leq 1 \\\\ 0 & \\text{otherwise} \\end{cases}$$\n",
    "\n",
    "**(b)** There is only one value in $\\Theta_0$, $\\theta = 1/2$.  The size of the test then is $\\beta(1/2) = 1 - \\left(\\text{clip}\\left(2c\\right)\\right)^n$.  Equating it to 0.05 and solving, we get $c = 0.5 \\cdot 0.95^{1/n}$.\n",
    "\n",
    "**(c)** The test has size $\\alpha$ when $ 1 - \\left(\\text{clip}\\left(2c\\right)\\right)^n = \\alpha$, or $c = 0.5 \\cdot (1 - \\alpha)^{1/n}$.  The p-value occurs on the threshold of the rejection region, that is, the infimal $\\alpha$ such that $Y > c$, so $Y > 0.5 \\cdot (1 - \\alpha)^{1/n}$, or $\\alpha > 1 - (2Y)^n$.\n",
    "\n",
    "For $n = 20$, $Y = 0.48$, that reduces to $\\alpha > 1 - (2\\cdot 0.48)^{20} \\approx 0.557$.  This means that the test provides little or no evidence against $H_0$.\n",
    "\n",
    "**(d)** For $n = 20$, $Y = 0.52$, the clipping function gets stuck at its maximum -- even for $\\alpha = 0$ we get $Y > c = 0.05$, so the p-value is 0, and the test almost surely rejects the null hypothesis.  This is in agreement with the observation that, since $Y > \\theta_0$, the maximum sample was generated with value greater than $\\theta_0$, so it is certain that $\\theta_* > \\theta_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.4**. There is a theory that people can postpone their death until after an important event.  To test the theory, Phillips and King (1988) collected data on deaths around the Jewish holiday Passover.  Of 1919 deaths, 922 died the week before the holiday and 997 died the week after.  Think of this as a binomial and test the null hypothesis that $\\theta = 1/2$.  Report and interpret the p-value.  Also construct a confidence interval for $\\theta$.\n",
    "\n",
    "Reference:\n",
    "Philips, D.P and King, E.W. (1988).\n",
    "*Death takes a holiday: Mortality surrounding major social occasions.*\n",
    "The Lancet, 2, 728-732."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  Let the number of deaths be $X \\sim \\text{Binomial}(n, \\theta)$, with $n = 1919$.  We have a measurement for $X = 922$, and we want to test the null hypothesis:\n",
    "\n",
    "$$H_0:  \\theta = 1/2 \\quad \\text{versus} \\quad H_1: \\theta \\neq 1/2$$\n",
    "\n",
    "The MLE is $\\hat{\\theta} = X / n = 922 / 1919$; the estimated standard error is $\\sqrt{\\hat{\\theta}(1 - \\hat{\\theta})/n}$.  The Wald test statistic is $W = (\\hat{\\theta} - \\theta_0) / \\hat{\\text{se}}$, and the p-value is $\\mathbb{P}(|Z| > |W|) = 2(1 - \\Phi(|W|))$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated theta: \t 0.480\n",
      "Estimated SE: \t\t 0.011\n",
      "95% confidence interval: (0.458, 0.503)\n",
      "Wald statistic: \t -1.713\n",
      "p-value: \t\t 0.087\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from scipy.stats import norm\n",
    "\n",
    "n = 1919\n",
    "X = 922\n",
    "theta_hat = 922 / 1919\n",
    "se_hat = math.sqrt(theta_hat * (1 - theta_hat) / n)\n",
    "\n",
    "z_05 = norm.ppf(0.975)\n",
    "confidence_interval = (theta_hat - z_05 * se_hat, theta_hat + z_05 * se_hat)\n",
    "\n",
    "w = (theta_hat - 0.5) / se_hat\n",
    "p_value = 2 * (1 - norm.cdf(abs(w)))\n",
    "\n",
    "print('Estimated theta: \\t %.3f' % theta_hat)\n",
    "print('Estimated SE: \\t\\t %.3f' % se_hat)\n",
    "print('95%% confidence interval: (%.3f, %.3f)' % confidence_interval)\n",
    "print('Wald statistic: \\t %.3f' % w)\n",
    "print('p-value: \\t\\t %.3f' % p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A 95% confidence interval can be built based on the estimated value for $\\theta$ and its estimated standard error, from 45.8% to 50.3%.\n",
    "\n",
    "The Wald test produces a p-value of 8.7%, which represents only weak evidence against the null hypothesis -- inconclusive results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.5**. In 1861, 10 essays appeared in the New Orleans Daily Crescent.  They were signed \"Quintus Curtuis Snodgrass\" and some people suspected they were actually written by Mark Twain.  To investigate this, we will investigate the proportion of three letter words found in an author's work.\n",
    "\n",
    "From eight Twain essays we have:\n",
    "\n",
    "$$ \n",
    "0.225 \\quad 0.262 \\quad 0.217 \\quad 0.240 \\quad 0.230 \\quad 0.229 \\quad 0.235 \\quad 0.217\n",
    "$$\n",
    "\n",
    "From 10 Snodgrass essays we have:\n",
    "\n",
    "$$\n",
    "0.209 \\quad 0.205 \\quad 0.196 \\quad 0.210 \\quad 0.202 \\quad 0.207 \\quad 0.224 \\quad 0.223 \\quad 0.220 \\quad 0.201\n",
    "$$\n",
    "\n",
    "(source: Rice xxxx)\n",
    "\n",
    "**(a)** Perform a Wald test for equality of the means.  Use the nonparametric plug-in estimator. Report the p-value and a 95% confidence interval for the difference of means.  What do you conclude?\n",
    "\n",
    "**(b)** Now use a permutation test to avoid the use of large sample methods.  What is your conclusion?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [0.225, 0.262, 0.217, 0.240, 0.230, 0.229, 0.235, 0.217]\n",
    "Y = [0.209, 0.205, 0.196, 0.210, 0.202, 0.207, 0.224, 0.223, 0.220, 0.201]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated difference of means:\t 0.022\n",
      "Estimated SE: \t\t\t 0.006\n",
      "95% confidence interval:\t (0.010, 0.034)\n",
      "Wald statistic: \t\t 3.704\n",
      "Wald test p-value: \t\t 0.0002\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "X = np.array(X)\n",
    "Y = np.array(Y)\n",
    "\n",
    "x_hat = X.mean()\n",
    "y_hat = Y.mean()\n",
    "\n",
    "diff_hat = x_hat - y_hat\n",
    "se_hat = np.sqrt(X.var(ddof=1)/len(X) + Y.var(ddof=1)/len(Y))\n",
    "\n",
    "z_05 = norm.ppf(0.975)\n",
    "confidence_interval = (diff_hat - z_05 * se_hat, diff_hat + z_05 * se_hat)\n",
    "\n",
    "w = diff_hat / se_hat\n",
    "p_value = 2 * (1 - norm.cdf(abs(w)))\n",
    "\n",
    "print('Estimated difference of means:\\t %.3f' % diff_hat)\n",
    "print('Estimated SE: \\t\\t\\t %.3f' % se_hat)\n",
    "print('95%% confidence interval:\\t (%.3f, %.3f)' % confidence_interval)\n",
    "print('Wald statistic: \\t\\t %.3f' % w)\n",
    "print('Wald test p-value: \\t\\t %.4f' % p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Such a small p-value (0.02%) offers very strong evidence to reject the null hypothesis -- that is, that the series follow distributions with different means.  The writing styles are different according to this metric, even if the same writer is responsible for both."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f1152ef31b3f42ecacd6cf539829f41e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Permutation test p-value: \t\t 0.0005\n"
     ]
    }
   ],
   "source": [
    "# Permutation test using random shuffling\n",
    "\n",
    "B = 1000000\n",
    "full_series = np.concatenate([X, Y])\n",
    "nx = len(X)\n",
    "diff_boot_count = 0\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    np.random.shuffle(full_series)\n",
    "    xx, yy = full_series[:nx], full_series[nx:]\n",
    "    diff_boot = xx.mean() - yy.mean()\n",
    "    if diff_boot > diff_hat:\n",
    "        diff_boot_count += 1\n",
    "        \n",
    "p_value_perm = diff_boot_count / B\n",
    "\n",
    "print('Permutation test p-value: \\t\\t %.4f' % p_value_perm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This small p-value (0.05%) also offers very strong evidence that the two samples come from series with different means."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.6**.  Let $X_1, \\dots, X_n \\sim N(\\theta, 1)$.  Consider testing\n",
    "\n",
    "$$ H_0: \\theta = 0 \\quad \\text{versus} \\quad H_1: \\theta = 1$$\n",
    "\n",
    "Let the rejection region be $R = \\{ x^n : T(x^n) > c \\}$ where $T(x^n) = n^{-1} \\sum_{i=1}^n X_i$.\n",
    "\n",
    "**(a)** Find $c$ so that the test has size $\\alpha$.\n",
    "\n",
    "**(b)** Find the power under $H_1$, i.e., find $\\beta(1)$.\n",
    "\n",
    "**(c)** Show that $\\beta(1) \\rightarrow 1$ as $n \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)** The size of the rejection region is\n",
    "\n",
    "$$\\mathbb{P}_{\\theta = 0}(T > c) = \\mathbb{P}_{\\theta = 0}(\\sqrt{n}T / 1 > \\sqrt{n}c) = \\mathbb{P}(Z > \\sqrt{n} c) = 1 - \\Phi(\\sqrt{n}c)$$\n",
    "\n",
    "If this value is $\\alpha$, then \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\alpha &= 1 - \\Phi(\\sqrt{n}c) \\\\\n",
    "\\Phi(\\sqrt{n}c) &= 1 - \\alpha \\\\\n",
    "\\sqrt{n}c &= z_{\\alpha} \\\\\n",
    "c &= z_{\\alpha} / \\sqrt{n}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** We have $\\beta(1) = \\mathbb{P}_{\\theta = 1}(T > c) = \\mathbb{\\theta = 1}(\\sqrt{n}(T - 1)/1 > \\sqrt{n}(c - 1)) = \\mathbb{P}(Z > \\sqrt{n}(c - 1)) = 1 - \\Phi(\\sqrt{n}(c - 1)) = 1 - \\Phi(z_\\alpha - \\sqrt{n})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)** As $n \\rightarrow \\infty$, $\\beta(1) \\rightarrow 1 - \\Phi(-\\infty) = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.7**.  Let $\\hat{\\theta}$ be the MLE of a parameter theta and let $\\hat{\\text{se}} = \\{ n I(\\hat{\\theta})\\}^{-1/2}$ where $I(\\hat{\\theta})$ is the Fisher information.  Consider testing\n",
    "\n",
    "$$H_0: \\theta = \\theta_0 \\quad \\text{versus} \\quad H_1: \\theta \\neq \\theta_0$$\n",
    "\n",
    "Consider the Wald test with rejection region $R = \\{ x^n: |Z| > z_{\\alpha/2} \\}$ where $Z = (\\hat{\\theta} - \\theta_0) / \\hat{\\text{se}}$.  Let $\\theta_1 > \\theta_0$ be some alternative.  Show that $\\beta(\\theta_1) \\rightarrow 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "We have \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\beta(\\theta_1) &= \\mathbb{P}_{\\theta_1}(|Z| > z_{\\alpha/2}) \\\\\n",
    "&= \\mathbb{P}_{\\theta_1}(Z > z_{\\alpha/2}) + \\mathbb{P}_{\\theta_1}(Z < - z_{\\alpha/2})\\\\\n",
    "&= \\mathbb{P}_{\\theta_1}\\left(\\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} > z_{\\alpha/2}\\right)\n",
    "+ \\mathbb{P}_{\\theta_1}\\left(\\frac{\\hat{\\theta} - \\theta_0}{\\hat{\\text{se}}} < -z_{\\alpha/2}\\right) \\\\\n",
    "&= \\mathbb{P}_{\\theta_1}\\left(\\hat{\\theta} > \\hat{\\text{se}} z_{\\alpha/2} + \\theta_0 \\right)\n",
    "+ \\mathbb{P}_{\\theta_1}\\left(\\hat{\\theta} < -\\hat{\\text{se}} z_{\\alpha/2} + \\theta_0 \\right) \\\\\n",
    "&= \\mathbb{P}_{\\theta_1}\\left(\\frac{\\hat{\\theta} - \\theta_1}{\\hat{\\text{se}}} >  z_{\\alpha/2} + \\frac{\\theta_0 - \\theta_1}{\\hat{\\text{se}}} \\right)\n",
    "+ \\mathbb{P}_{\\theta_1}\\left(\\frac{\\hat{\\theta} - \\theta_1}{\\hat{\\text{se}}} <  -z_{\\alpha/2} + \\frac{\\theta_0 - \\theta_1}{\\hat{\\text{se}}} \\right) \\\\\n",
    "&= \\mathbb{P}_{\\theta_1}\\left(W >  z_{\\alpha/2} + \\frac{\\theta_0 - \\theta_1}{\\hat{\\text{se}}} \\right)\n",
    "+ \\mathbb{P}_{\\theta_1}\\left(W <  -z_{\\alpha/2} + \\frac{\\theta_0 - \\theta_1}{\\hat{\\text{se}}} \\right) \\\\\n",
    "& \\geq \\mathbb{P}_{\\theta_1}\\left(W >  z_{\\alpha/2} + \\frac{\\theta_0 - \\theta_1}{\\hat{\\text{se}}} \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "When $n \\rightarrow \\infty$, $\\text{\\hat{se}} \\rightarrow 0$.  Since $\\theta_1 > \\theta_0$, the right hand side of the inequality goes to $-\\infty$, and so the lower bound probability goes to 1.  Therefore $\\beta(\\theta_1) \\rightarrow 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.8**.  Here are the number of elderly Jewish and Chinese women who died just before and after the Chinese Harvest Moon Festival.\n",
    "\n",
    "| Week | Chinese | Jewish |\n",
    "|------|---------|--------|\n",
    "| -2   | 55      | 141    |\n",
    "| -1   | 33      | 145    |\n",
    "| 1    | 70      | 139    |\n",
    "| 2    | 49      | 161    |\n",
    "\n",
    "Compare the two mortality patterns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame({\n",
    "    'Chinese': [55, 33, 70, 49],\n",
    "    'Jewish':  [141, 145, 139, 161]\n",
    "}, index=[-2, -1, 1, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use a few different approaches.\n",
    "\n",
    "First:  $\\chi^2$ tests comparing, separately, Chinese and Jewish mortalities before and after the Chinese Harvest Moon Festival."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test statistic:\t\t\t\t4.643 \n",
      "95% percentile chi squared with 1 df:\t3.841\n",
      "p-value for different distributions:\t0.031\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import chi2\n",
    "\n",
    "chinese_before = df['Chinese'].loc[-2:-1].sum()\n",
    "chinese_after = df['Chinese'].loc[1:2].sum()\n",
    "\n",
    "mu_chinese_hat = (chinese_before + chinese_after) / 2\n",
    "t_pearson_chinese = (chinese_before - mu_chinese_hat)**2/mu_chinese_hat \\\n",
    "    + (chinese_after - mu_chinese_hat)**2/mu_chinese_hat\n",
    "\n",
    "chi2_95 = chi2.ppf(0.95, 1)\n",
    "p_value_chinese = 1 - chi2.cdf(t_pearson_chinese, 1)\n",
    "\n",
    "print('Test statistic:\\t\\t\\t\\t%.3f ' % t_pearson_chinese)\n",
    "print('95%% percentile chi squared with 1 df:\\t%.3f' % chi2_95)\n",
    "print('p-value for different distributions:\\t%.3f' % p_value_chinese)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test statistic:\t\t\t\t0.334 \n",
      "95% percentile chi squared with 1 df:\t3.841\n",
      "p-value for different distributions:\t0.563\n"
     ]
    }
   ],
   "source": [
    "jewish_before = df['Jewish'].loc[-2:-1].sum()\n",
    "jewish_after = df['Jewish'].loc[1:2].sum()\n",
    "\n",
    "mu_jewish_hat = (jewish_before + jewish_after) / 2\n",
    "t_pearson_jewish = (jewish_before - mu_jewish_hat)**2/mu_jewish_hat \\\n",
    "    + (jewish_after - mu_jewish_hat)**2/mu_jewish_hat\n",
    "\n",
    "chi2_95 = chi2.ppf(0.95, 1)\n",
    "p_value_jewish = 1 - chi2.cdf(t_pearson_jewish, 1)\n",
    "\n",
    "print('Test statistic:\\t\\t\\t\\t%.3f ' % t_pearson_jewish)\n",
    "print('95%% percentile chi squared with 1 df:\\t%.3f' % chi2_95)\n",
    "print('p-value for different distributions:\\t%.3f' % p_value_jewish)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\chi^2$ test statistics for the Chinese population *do* suggest, with a p-value of 3.1%, that the mortality rates are distinct before and after the event, while the statistics for the Jewish population do not suggest so (with a no-evidence p-value of 56.3%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the likelihood test to compare mortality in each group before and after the event:\n",
    "\n",
    "| Group   | Chinese | Jewish |\n",
    "|---------|---------|--------|\n",
    "| Before  | 88      | 286    |\n",
    "| After   | 119     | 300    |\n",
    "\n",
    "MLE under alternate hypothesis: $$\\hat{\\theta} = \\left( \\frac{88}{793}, \\frac{119}{793}, \\frac{286}{793}, \\frac{300}{793} \\right)$$\n",
    "\n",
    "MLE under null hypothesis:  $$\\hat{\\theta_0} =\\left( \\frac{103.5}{793}, \\frac{103.5}{793}, \\frac{293}{793}, \\frac{293}{793} \\right)$$\n",
    "\n",
    "Log-likelihood ratio statistic: \n",
    "\n",
    "$$ \\lambda = 2 \\log \\left( \\frac{\\mathcal{L}(\\hat{\\theta})}{\\mathcal{L}(\\hat{\\theta_0})} \\right)\n",
    "= 2 \\left( 88 \\log \\frac{88}{103.5} + 119 \\log \\frac{119}{103.5} + 286 \\log \\frac{286}{293} + 300 \\log \\frac{300}{293} \\right) \\approx 2.497\n",
    "$$\n",
    "\n",
    "The 95% percentile value for the $\\chi^2$ distribution with 1 degree of freedom is $\\approx 3.841$, so the likelihood test does not reject the null hypothesis at this confidence level.  The p-value for this test is 0.114, suggesting little or no evidence that, for both groups, the mortality rates are distinct."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's do 4 distinct sample binomial tests, one for each week:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>p-value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>-2</th>\n",
       "      <td>0.516380</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-1</th>\n",
       "      <td>0.021061</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.017937</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.388072</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     p-value\n",
       "-2  0.516380\n",
       "-1  0.021061\n",
       " 1  0.017937\n",
       " 2  0.388072"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from scipy.stats import binom_test\n",
    "\n",
    "theta_week1 = 55 / (55 + 141)\n",
    "theta_week2 = 33 / (33 + 145)\n",
    "theta_week3 = 70 / (70 + 139)\n",
    "theta_week4 = 49 / (49 + 161)\n",
    "\n",
    "theta_null = (55 + 33 + 70 + 49) / (55 + 33 + 70 + 49 + 141 + 145 + 139 + 161)\n",
    "\n",
    "p1 = binom_test(55, 55 + 141, theta_null, alternative=\"two-sided\")\n",
    "p2 = binom_test(33, 33 + 145, theta_null, alternative=\"two-sided\")\n",
    "p3 = binom_test(70, 70 + 139, theta_null, alternative=\"two-sided\")\n",
    "p4 = binom_test(49, 49 + 161, theta_null, alternative=\"two-sided\")\n",
    "\n",
    "results = pd.DataFrame({'p-value': [p1, p2, p3, p4]}, index=[-2, -1, 1, 2])\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tests suggest that there *is* strong evidence for distinct mortality rates between the Chinese and Jewish populations on the weeks immediately preceding and following the event (p-values of 2.1% and 1.8%), but not in the 2 weeks before or after the event (p-values of 51.6% and 38.8%)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.9**. A randomized, double-blind experiment was conducted to assess the effectiveness of several drugs for reducing post-operative nausea.  The data are as follows.\n",
    "\n",
    "|                        | Number of patients | Incidents of Nausea |\n",
    "|------------------------|--------------------|---------------------|\n",
    "| Placebo                | 80                 | 45                  |\n",
    "| Chlorpromazine         | 75                 | 26                  |\n",
    "| Dimenhydrate           | 85                 | 52                  |\n",
    "| Pentobarbital (100 mg) | 67                 | 35                  |\n",
    "| Pentobarbital (150 mg) | 85                 | 37                  |\n",
    "\n",
    "**(a)** Test each drug versus the placebo at the 5% level.  Also, report the estimated odds-ratios.  Summarize your findings.\n",
    "\n",
    "**(b)** Use the Bonferroni and the FDR method to adjust for multiple testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm\n",
    "\n",
    "df = pd.DataFrame({\n",
    "    'Treatment': ['Placebo', 'Chlorpromazine', 'Dimenhydrate', 'Pentobarbital (100 mg)', 'Pentobarbital (150 mg)'],\n",
    "    'Number of patients': [80, 75, 85, 67, 85],\n",
    "    'Incidents of Nausea': [45, 26, 52, 35, 37]\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['p_hat'] = df['Incidents of Nausea'] / df['Number of patients']\n",
    "df['variance'] = df['p_hat'] * (1 - df['p_hat']) / df['Number of patients']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_hat_placebo = df[df['Treatment'] == 'Placebo']['p_hat'][0]\n",
    "variance_placebo = df[df['Treatment'] == 'Placebo']['variance'][0]\n",
    "\n",
    "# odds = affected / normal = p / (1 - p)\n",
    "df['Odds'] = df['p_hat'] / (1 - df['p_hat'])\n",
    "\n",
    "# odds ratio = odds_placebo / odds_treatment\n",
    "df['Odds ratio'] = df[df['Treatment'] == 'Placebo']['Odds'][0] / df['Odds']\n",
    "\n",
    "df['Wald statistic'] = (df['p_hat'] - p_hat_placebo) / np.sqrt(df['variance'] + variance_placebo)\n",
    "df['p-value'] = 2 * (1 - norm.cdf(abs(df['Wald statistic'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Treatment</th>\n",
       "      <th>Odds ratio</th>\n",
       "      <th>Wald statistic</th>\n",
       "      <th>p-value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Chlorpromazine</td>\n",
       "      <td>2.423077</td>\n",
       "      <td>-2.764364</td>\n",
       "      <td>0.005703</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Dimenhydrate</td>\n",
       "      <td>0.815934</td>\n",
       "      <td>0.642987</td>\n",
       "      <td>0.520232</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Pentobarbital (100 mg)</td>\n",
       "      <td>1.175510</td>\n",
       "      <td>-0.486428</td>\n",
       "      <td>0.626664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Pentobarbital (150 mg)</td>\n",
       "      <td>1.667954</td>\n",
       "      <td>-1.646605</td>\n",
       "      <td>0.099639</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Treatment  Odds ratio  Wald statistic   p-value\n",
       "1          Chlorpromazine    2.423077       -2.764364  0.005703\n",
       "2            Dimenhydrate    0.815934        0.642987  0.520232\n",
       "3  Pentobarbital (100 mg)    1.175510       -0.486428  0.626664\n",
       "4  Pentobarbital (150 mg)    1.667954       -1.646605  0.099639"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[['Treatment', 'Odds ratio', 'Wald statistic', 'p-value']][df['Treatment'] != 'Placebo']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the 5% significance level, only Chlorpromazine changes the Nausea incidence, with an odds ratio of 2.42."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** There are 4 tests, so the significance ratio for the Bonferroni method will be $0.05 / 4 = 0.0125$.  Again, only Chlorpromazine rejects the null hypothesis.\n",
    "\n",
    "For the Benjamini-Hochfield method, the ordered p-values are $0.0057 < 0.0996 < 0.5202 < 0.6266$; $\\ell_i = i \\alpha / m = 0.0125 \\cdot i$, and $R = \\max \\{ i : P_{(i)} < \\ell_i \\} = 1$, so the rejection threshold again only rejects the null hypothesis for the first treatment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 11.9.10**.  Let $X_1, \\dots, X_n \\sim \\text{Poisson}(\\lambda)$.\n",
    "\n",
    "**(a)** Let $\\lambda_0 > 0$. Find the size $\\alpha$ Wald test for \n",
    "\n",
    "$$ H_0: \\lambda = \\lambda_0 \\quad \\text{versus} \\quad H_1: \\lambda \\neq \\lambda_0$$\n",
    "\n",
    "**(b)** (Computer Experiment)  Let $\\lambda_0 = 1$, $n = 20$ and $\\alpha = 0.05$.  Simulate $X_1, \\dots, X_n \\sim \\text{Poisson}(\\lambda_0)$ and perform the Wald test.  Repeat many times and count how often you reject the null.  How close is the type I error rate to 0.05?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)**\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\beta(\\lambda_0) &= \\mathbb{P}_{\\lambda_0}\\left( |W| > z_{\\alpha/2} \\right)  \\\\\n",
    "&= \\mathbb{P}_{\\lambda_0}\\left( \\Bigg| \\frac{\\hat{\\lambda} - \\lambda_0}{\\hat{\\text{se}}} \\Bigg| > z_{\\alpha/2} \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The MLE for the parameter is the mean of the data, $\\hat{\\lambda} = \\overline{X}$, while the estimate of the SE is $\\sqrt{1 / I_n(\\hat{\\lambda})} = 1 / \\sqrt{n\\overline{X}}$.  Replacing these, we get\n",
    "\n",
    "$$\n",
    "\\mathbb{P}_{\\lambda_0}\\left( \\big| (\\overline{X} - \\lambda_0)\\sqrt{n \\overline{X}} \\big| > z_{\\alpha/2} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm, poisson\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "lambda_0 = 1\n",
    "n = 20\n",
    "alpha = 0.05\n",
    "\n",
    "z = norm.ppf(1 - alpha / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rejection:  False\n"
     ]
    }
   ],
   "source": [
    "# Perform process once\n",
    "X = poisson.rvs(lambda_0, size=n)\n",
    "lambda_hat = X.mean()\n",
    "rejection = ((lambda_hat - lambda_0)**2) * (n * lambda_hat) > z**2\n",
    "print('Rejection: ', rejection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0dc88ab9f88c49dda002cbb42e258c66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Fraction of rejections: 0.052\n"
     ]
    }
   ],
   "source": [
    "# Perform process many times\n",
    "B = 1000000\n",
    "num_rejections = 0\n",
    "for i in tqdm_notebook(range(B)):\n",
    "    X = poisson.rvs(lambda_0, size=n)\n",
    "    lambda_hat = X.mean()\n",
    "    if ((lambda_hat - lambda_0)**2) * (n * lambda_hat) > z**2:\n",
    "        num_rejections += 1\n",
    "        \n",
    "fraction_rejections = num_rejections / B\n",
    "\n",
    "print('Fraction of rejections: %.3f' % fraction_rejections)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The measured fraction of rejections of the null hypothesis is $0.052$, which is very close to $\\alpha = 0.05$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}