第 12 章假设检验理论

CHAPTER STATUS编写进度1%查看全书进度

Definition 12.1. 设\((X,\mathscr{A},\mathscr{P})\)是统计结构，则称\(\mathscr{P}\)的非空子集为假设（hypothesis）。若\((X,\mathscr{A},\mathscr{P})\)是参数结构，\(\Theta\)是参数空间，则\(\Theta\)的非空子集也被称之为假设。设\(f_1, f_2, \dots, f_{n}\)是从总体\(F\)中抽取的简单样本，判断\(F\)是否服从某一假设的问题是假设检验（hypothesis testing）问题。称要检验的假设\(\mathscr{P}_0(\Theta_0)\)为原假设（null hypothesis），记为\(H_0\)，\(\mathscr{P}_1\subseteq\mathscr{P}\setminus\mathscr{P}_0(\Theta_1\subseteq\Theta\setminus\Theta_0)\)被称之为备择假设（alternative hypothesis），记为\(H_1\)。检验指的是：

将\(X\)划分为互不相交的两个\(\mathscr{A}\)中的元素，即\(X=W\cup\overline{W},\;W\cap\overline{W}=\varnothing,\;W,\overline{W}\in\mathscr{A}\)；
作如下规定，若一次试验中\(\underset{i=1}{\overset{n}{\cap}}\{f_i\in W\}\)发生了，则拒绝原假设\(H_0\)，接受备择假设\(H_1\)；反之则接受原假设\(H_0\)，拒绝备择假设\(H_1\)。

称\(W\)为拒绝域（reject region）。

Definition 12.2. 在检验问题\((H_0,H_1)\)中，若\(H_0\)为真但由于样本的随机性拒绝了\(H_0\)，称此时犯的错误为第一类错误（type I error）或拒真错误；若\(H_0\)为伪但由于样本的随机性接受了\(H_0\)，称此时犯的错误为第二类错误（type II error）或取伪错误。

Definition 12.3. 设\((X,\mathscr{A},\mathscr{P})\)是统计结构，\(f_1, f_2, \dots, f_{n}\)为从总体\(F\)中抽取的简单样本，\((H_0,H_1)\)是检验问题，\(W\)为拒绝域。称满足条件：

\[\begin{equation*} \forall\;P\in\mathscr{P},\;g(P)=P(W) \end{equation*}\]

的函数\(g\)为检验问题\((H_0,H_1)\)的势函数（power function）。当\(F\in H_0\)时，\(g(F)\leqslant\alpha\)的检验被称之为显著性水平（significance level）为\(\alpha\)的检验。

note 12.1. 当\(F\in H_0\)时，\(g(F)\)为犯第一类错误的概率；当\(F\in H_1\)时，\(1-g(F)\)为犯第二类错误的概率。

Method 12.1. 假设检验基本步骤为：

建立假设检验问题；
选取统计量，由直观得出拒绝域形式；
选择检验的显著性水平；
确定具体的拒绝域；
做出判断。

Definition 12.4. 设\((H_0,H_1)\)是统计结构\((X,\mathscr{A},\mathscr{P})\)上的一个假设检验问题，\(f_1, f_2, \dots, f_{n}\)是从总体\(F\)中抽取的简单样本，利用\(f_1, f_2, \dots, f_{n}\)能够作出拒绝原假设的最小显著性水平被称为检验的p值（pvalue）。

note 12.2. 为什么需要\(p\)值的概念？由于不同的人对原假设的保守性程度不一样，他们选择的显著性水平就会不一样，对原假设越保守那所选取的\(\alpha\)就会越小，只关注是否拒绝原假设这一结论就会引起争议。如果能够给出p值让大家有一个大致的参照，虽然此时还是会有争议，但至少我们有了一个参考值。

如何求\(p\)值？为了得到拒绝原假设的最小显著性水平，那必然需要使得当前样本下的统计量值在拒绝域内且在拒绝域的边界上，否则拒绝域还能进一步缩小，即该显著性水平不是最小的。因为此时拒绝域的形式已经确定了，边界值也给出了，那么实际上我们就得到了一个具体的拒绝域，所以有：

\[\begin{equation*} p=\sup_{P\in H_0}g(P) \end{equation*}\]

如果\(p\)更小，那么就无法保证当\(F\in H_0\)的时候\(g(F)\leqslant\alpha\)；如果\(p\)更大则不最小。

note 12.3. 实际中我们根本不知道\(F\)到底是在\(H_0\)中还是在\(H_1\)中，假设检验本就是为了判断这个问题才出现的，为了得到显著性为\(\alpha\)的检验，实际上我们做的是找到\(X\)的一个划分使得：

\[\begin{equation*} \sup_{P\in H_0}g(P)\leqslant\alpha \end{equation*}\]

这就会导致当\(F\)真的在\(H_0\)中时，犯第一类错误的概率实际上有可能小于\(\alpha\)。这代表着检验还能进一步优化，由此引入下述随机化检验的概念。

Definition 12.5. 设\(\varphi\)是可测空间\((X,\mathscr{A})\)上的可测函数。若\(0\leqslant\varphi\leqslant1\)，则称\(\varphi\)是检验函数（test function）。当\(\varphi\)只取\(0,1\)两个值时，称对应的检验为非随机化检验（non-randomized testing），否则称为随机化检验（randomized testing）。\(\varphi\)的势函数为\(\operatorname{E}(\varphi)\)。

12.1 多重假设检验

Property 12.1.1. 若控制每一次假设检验的显著性水平为\(\alpha\)，\(n\)次独立假设检验中至少犯一次第一次错误的概率为\(1-(1-\alpha)^n\)。

证明. 此时一次第一类错误都没犯的概率为\((1-\alpha)^n\)。 ◻

设\(n\)个原假设中有\(n_0\)个为真，\(n_0\)个为假。记试验中被错误拒绝的假设个数为\(V\)，被正确拒绝的假设个数为\(S\)，总拒绝数为\(R=V+S\)。

Definition 12.6. 称\(P(V\geqslant1)\)为族错误率（family-wise error rate）。

Definition 12.7. 称：

\[\begin{equation*} \operatorname{E}\left(\frac{V}{\max\{R,1\}}\right) \end{equation*}\]

为错误检出率（false discovery rate）。

Method 12.2. （BonFerroni）

由性质 5.2.1(3)（次有限可加性）可知，将\(n\)个假设检验问题中单个假设检验的显著性水平设置为\(\dfrac{\alpha}{n}\)，即可保证FWER小于\(\alpha\)。

Method 12.3. （Benjamini–Hochberg）

设\(n\)个假设检验对应的\(p\)值为\(p_1,p_2\dots,p_n\)，将其从小到大排序为\(p_{(1)}\leqslant p_{(2)}\leqslant\cdots\leqslant p_{(n)}\)。给定目标水平\(\alpha\)，取最大的\(k\)使得：

\[\begin{equation*} p_{(k)}\leqslant\frac{k}{n}\alpha \end{equation*}\]

则拒绝\(H_{(1)},H_{(2)}\dots,H_{(k)}\)。若各真原假设对应的\(p\)值独立且服从\(U(0,1)\)，并且与其余\(p\)值独立，则此时的FDR不超过\(\alpha\)。

证明. 记\(I_0\)为真原假设的指标集合，\(m_0=|I_0|\)。对每个\(i\in I_0\)，记\(R^{(-i)}\)为删去\(p_i\)之后，其余\(n-1\)个\(p\)值按同一Benjamini–Hochberg规则所产生的拒绝数。由逐步上升法的结构可知，若第\(i\)个真原假设被错误拒绝且\(R=r\)，则必有

\[\begin{equation*} p_i\leqslant \frac{r}{n}\alpha \end{equation*}\]

并且

\[\begin{equation*} R^{(-i)}=r-1 \end{equation*}\]

于是

\[\begin{align*} \mathrm{FDR} &=E\left(\frac{V}{R\vee1}\right)\\ &=\sum_{i\in I_0}E\left(\frac{\mathbf{1}\{p_i\text{被拒绝}\}}{R\vee1}\right)\\ &=\sum_{i\in I_0}\sum_{r=1}^n\frac1rP\left(p_i\leqslant \frac{r}{n}\alpha,R=r\right)\\ &\leqslant \sum_{i\in I_0}\sum_{r=1}^n\frac1rP\left(p_i\leqslant \frac{r}{n}\alpha,R^{(-i)}=r-1\right) \end{align*}\]

若各真原假设对应的\(p\)值独立且服从\(U(0,1)\)，并且与其余\(p\)值独立，则

\[\begin{align*} \mathrm{FDR} &\leqslant \sum_{i\in I_0}\sum_{r=1}^n\frac1rP\left(p_i\leqslant \frac{r}{n}\alpha\right)P\left(R^{(-i)}=r-1\right)\\ &=\sum_{i\in I_0}\sum_{r=1}^n\frac1r\frac{r\alpha}{n}P\left(R^{(-i)}=r-1\right)\\ &=\sum_{i\in I_0}\frac{\alpha}{n}\sum_{r=1}^nP\left(R^{(-i)}=r-1\right)\\ &=\frac{m_0}{n}\alpha\\ &\leqslant \alpha \end{align*}\]

故Benjamini–Hochberg方法在独立情形下控制FDR ◻

Method 12.4. （Benjamini–Yekutieli）

设\(n\)个假设检验对应的\(p\)值为\(p_1,\dots,p_n\)，将其从小到大排序为

\[\begin{equation*} p_{(1)}\leqslant p_{(2)}\leqslant\cdots\leqslant p_{(n)} \end{equation*}\]

记

\[\begin{equation*} c_n=\sum_{j=1}^n\frac1j \end{equation*}\]

给定目标水平\(\alpha\)，取最大的\(k\)使得

\[\begin{equation*} p_{(k)}\leqslant \frac{k}{n\,c_n}\alpha \end{equation*}\]

则拒绝

\[\begin{equation*} H_{(1)},\dots,H_{(k)} \end{equation*}\]

证明：Benjamini–Yekutieli证明了如下上界：对任意依赖结构，若采用Benjamini–Hochberg型逐步上升法，并以

\[\begin{equation*} \alpha_i=\frac{i}{n}\gamma,\qquad i=1,\dots,n \end{equation*}\]

作为临界值，则有

\[\begin{equation*} \mathrm{FDR}\leqslant \frac{m_0}{n}\gamma c_n \end{equation*}\]

现在取

\[\begin{equation*} \gamma=\frac{\alpha}{c_n} \end{equation*}\]

则对应临界值恰为

\[\begin{equation*} \alpha_i=\frac{i}{n\,c_n}\alpha \end{equation*}\]

于是

\[\begin{align*} \mathrm{FDR} &\leqslant \frac{m_0}{n}\frac{\alpha}{c_n}c_n\\ &=\frac{m_0}{n}\alpha\\ &\leqslant \alpha \end{align*}\]

故Benjamini–Yekutieli方法在任意依赖结构下控制FDR

12.2 非参数假设检验

12.2.1 位置检验

位置检验分为两类：

单样本位置检验：给定一个数\(q_0\)，检验它与总体\(\pi\)分位点\(Q_{\pi}\)之间的大小关系；
多样本位置检验：检验多个总体位置参数的大小关系。

12.2.1.1 单样本位置检验

Method 12.5. （Sign Test）

设\((X,\mathscr{A},\mathscr{P})\)是统计结构，\(\mathbf{X}=(X_1, X_2, \dots, X_{n})\)为从总体\(F\)中抽取的简单样本，\(\mathscr{P}\mathbf{X}^{-1}\)为连续型随机变量的概率分布族，检验问题为：

\[\begin{equation*} \begin{cases} H_0:Q_\pi\leqslant q_0,\;H_1:Q_\pi>q_0 \\ H_0:Q_\pi\geqslant q_0,\;H_1:Q_\pi<q_0 \\ H_0:Q_\pi=q_0,\;H_1:Q_\pi\ne q_0 \end{cases} \end{equation*}\]

记\(\mathbf{X}\)中大于\(q_0\)的样本单元的个数为\(S^+\)，小于\(q_0\)的个数为\(S^-\)。令\(K\sim\operatorname{Binom}(S^++S^-,\pi)\)，则上述检验问题的\(p\)值分别为：

\[\begin{equation*} P_{H_0}(K\leqslant S^-),\quad P_{H_0}(K\geqslant S^-),\quad 2\min\{P_{H_0}(K\leqslant S^-),\;P_{H_0}(K\geqslant S^-)\} \end{equation*}\]

当\(n\to+\infty\)时，令\(Z=\dfrac{K-(S^++S^-)\pi}{\sqrt{(S^++S^-)\pi(1-\pi)}}\)，\(\Phi\)为标准正态分布的分布函数，则有如下对应的近似\(p\)值：

\[\begin{equation*} \Phi(Z),\quad1-\Phi(Z),\quad2\min\{\Phi(Z),1-\Phi(Z)\} \end{equation*}\]

Derivation 12.1. 直观上来讲，如果\(Q_{\pi}>q_0\)，那么\(K\)的实现值\(S^-\)应偏小，原假设的拒绝域应具有形式\(\{S^-\leqslant c\}\)。由性质 6.4.1(3)可知对于检验的\(p\)值可仅考虑\(Q_\pi=q_0\)，设\(K\)为\(S^-\)的随机变量形式，则应有\(K\sim\operatorname{Binom}(S^++S^-,\pi)\)，于是此时\(p\)值为：

\[\begin{equation*} p=P_{H_0}(K\leqslant S^-)=\sum_{i=0}^{S^-}\binom{S^++S^-}{i}\pi^{i}(1-\pi)^{S^++S^--i} \end{equation*}\]

同理可得另外两种情况。不对正态近似作解释。双边假设检验p值计算的争议问题

note 12.4. 符号检验仅对连续型随机变量成立，离散型随机变量无法满足\(K\sim\operatorname{Binom}(S^++S^-,\pi)\)（分布函数阶梯形，分位点唯一但小于等于该分位点的概率不一定等于对应的分位数，处在该阶梯上的所有点都小于等于该分位点）。

能否不删除等于\(q_0\)的数值？实际上可以，甚至我个人觉得从严谨的角度来讲应该不删除，考虑到符号检验的原理\(K\sim\operatorname{Binom}\)，二项分布的概率\(p\)实际上是\(X_i\)小于等于\(q_0\)的概率呀，但我们约定删除，删除也没影响（连续型随机变量在单点集上测度为\(0\)）。

sign.test <- function(x, q0, pi, exact = FALSE, 
                      alternative = c("two.sided", "greater", "less")) {
  alternative <- match.arg(alternative)
  s_pos <- sum(x > q0)
  STATISTIC <- sum(x < q0)
  n <- s_pos + STATISTIC
  STATISTIC <- setNames(STATISTIC, "s_neg")
  hypothesis_message <- switch(alternative,
                               "greater" = paste("Q0 >", pi),
                               "less" = paste("Q0 <", pi),
                               "two.sided" = paste("Q0 !=", pi))
  if (exact) {
    METHOD <- "Sign Test (Exact)"
    p.value <- switch(alternative,
                      "greater" = pbinom(q = STATISTIC, size = n, prob = pi),
                      "less" = pbinom(q = STATISTIC-1, size = n, prob = pi, lower.tail = FALSE),
                      "two.sided" = 2 * min(pbinom(q = STATISTIC-1, size = n, prob = pi, lower.tail = FALSE),
                                            pbinom(q = STATISTIC, size = n, prob = pi)))
  } else {
    METHOD <- "Sign Test (Approximate)"
    z_score <- (STATISTIC - n * pi) / sqrt(n * pi * (1 - pi)) 
    p.value <- switch(alternative,
                      "greater" = pnorm(z_score, lower.tail = TRUE),
                      "less" = pnorm(z_score, lower.tail = FALSE),
                      "two.sided" = 2 * min(pnorm(z_score), 1 - pnorm(z_score)))
  }
  
  RVAL <- list(statistic = STATISTIC, p.value = p.value,
               data.name = deparse(substitute(x)),
               method = METHOD, alternative = hypothesis_message)
  class(RVAL) <- "htest"
  RVAL
}

Method 12.6. （Wilcoxon Signed Rank Test）

设\((X,\mathscr{A},\mathscr{P})\)是统计结构，\(\mathbf{X}=(X_1, X_2, \dots, X_{n})\)为从总体\(F\)中抽取的简单样本，\(\mathscr{P}\mathbf{X}^{-1}\)为对称连续型随机变量的概率分布族，检验问题为：

\[\begin{equation*} H_0:M=M_0,\quad H_1: \begin{cases} M>M_0 \\ M<M_0 \\ M\ne M_0 \end{cases} \end{equation*}\]

记\(W^+,W^-\)为Wilcoxon秩统计量，\(W\)服从Wilcoxon秩分布，则上述检验问题的\(p\)值分别为：

\[\begin{equation*} P_{H_0}(W\geqslant W^+),\quad P_{H_0}(W\leqslant W^+) \end{equation*}\]

对\(i=1,2,\dots,n\)，\(d_i=|X_i-M_0|\)。
对\(d_i\)进行排序，求对应的秩\(R_i\)。
\(W^+=\sum\limits_{i\in\{i:X_i-M_0>0\}}R_i,\;W^-=\sum\limits_{i\in\{i:X_i-M_0<0\}}R_i\)。

Derivation 12.2. 若零假设成立，则\(W^+\)和\(W^-\)应相差不大，如果\(M>M_0\)，则\(W^+\)应较大，原假设的拒绝域应具有形式\(\{W^+\geqslant c\}\)，\(p\)值为\(p=P_{H_0}(W\geqslant W^+)\)；如果\(M<M_0\)，则\(W^+\)应较小，原假设的拒绝域应具有形式\(\{W^+\leqslant c\}\)，\(p\)值为\(p=P_{H_0}(W\leqslant W^+)\)；

12.2.2 趋势检验

Method 12.7. 设\(F(x)\)为某连续型随机变量的分布函数，\(X_1, X_2, \dots, X_{n}\)是按一定顺序收集到的样本，\(X_i\sim F(x-\theta_i)\)。检验问题为：

\[\begin{equation*} \begin{cases} H_0:\text{序列无上升趋势},\;H_1:\text{序列有上升趋势，即}\theta_1\leqslant\theta_2\leqslant\cdots\leqslant\theta_n\text{且至少成立一个严格不等式} \\ H_0:\text{序列无下降趋势},\;H_1:\text{序列有下降趋势，即}\theta_1\geqslant\theta_2\geqslant\cdots\geqslant\theta_n\text{且至少成立一个严格不等式} \\ H_0:\text{序列无趋势},\;H_1:\text{序列有上升或下降趋势} \end{cases} \end{equation*}\]

取数对\((X_i,X_{i+c})\)，其中：

\[\begin{equation*} c= \begin{cases} \dfrac{n}{2},&n\text{为偶数} \\ \dfrac{n+1}{2},&n\text{为奇数} \end{cases} \end{equation*}\]

记\(D_i=X_{i+c}-X_i\)，\(D_i\)中大于\(0\)的样本单元的个数为\(S^+\)，小于\(0\)的个数为\(S^-\)。令\(K\sim\operatorname{Binom}(S^++S^-,0.5)\)，则上述检验问题的\(p\)值分别为：

\[\begin{equation*} P_{H_0}(K\leqslant S^-),\quad P_{H_0}(K\geqslant S^-),\quad 2\min\{P_{H_0}(K\leqslant S^-),\;P_{H_0}(K\geqslant S^-)\} \end{equation*}\]

当\(n\to+\infty\)时，令\(Z=\dfrac{K-(S^++S^-)\pi}{\sqrt{(S^++S^-)\pi(1-\pi)}}\)，\(\Phi\)为标准正态分布的分布函数，则有如下对应的近似\(p\)值：

\[\begin{equation*} \Phi(Z),\quad1-\Phi(Z),\quad2\min\{\Phi(Z),1-\Phi(Z)\} \end{equation*}\]

Derivation 12.3. 原理与符号检验相同，不作解释。但这里的情况其实与符号检验还是有区别的，势函数在\(H_0\)上的上确界真的是由上述\(p\)值给出的吗？还是说在这里我们只是考虑一个近似呢？

cox.stuart.test <- function(x, exact = FALSE, 
alternative = c("two.sided", "greater", "less")) {
    alternative <- match.arg(alternative)
    n <- length(x)
    c <- ifelse(n %% 2 == 0, n / 2, (n + 1) / 2)
    D <- x[(c + 1):n] - x[1:c]
    D_nonzero <- D[D != 0]
    if (length(D_nonzero) == 0) {
        stop("All paired differences are zero; cannot perform Cox–Stuart test.")
    }
    RVAL <- sign.test(x = D_nonzero, q0 = 0, pi = 0.5, 
    exact = exact, alternative = alternative)
    RVA$method <- if (exact) {
        "Cox–Stuart Trend Test (Exact via Sign Test)"
    } else {
        "Cox–Stuart Trend Test (Approximate via Sign Test)"
    }
    RVA$data.name <- deparse(substitute(x))
    RVA
}

12.2.3 随机性检验

Definition 12.8. 设\(X_1, X_2, \dots, X_{n}\)是一个序列，若\(X_{i},X_{i+1},\dots,X_{j}\)满足某条件，而\(X_{i-1},X_{j+1}\)不满足该条件，则称\(X_{i},X_{i+1},\dots,X_{j}\)是一个游程（run）。

example 12.1. 将满足条件的数转换为\(1\)，不满足的转换为\(0\)（下面我们将延续该约定），则数据\(0,0,0,0,1,1,0,1,1,1,0,0,0\)中，一共有\(5\)个游程（\(3\)个\(0\)游程，\(2\)个\(1\)游程）。

Property 12.2.1. 设\(R\)为游程数，则有空证明：

\[\begin{gather*} P(R=2k)=\frac{2\binom{m-1}{k-1}\binom{n-1}{k-1}}{\binom{N}{n}} \\ P(R=2k+1)=\frac{\binom{m-1}{k-1}\binom{n-1}{k}+\binom{m-1}{k}\binom{n-1}{k-1}}{\binom{N}{n}} \\ E(R) = \frac{2mn}{m+1}+1\quad Var(R)=\dfrac{2mn(2mn-m-n)}{(m+n)^2(m+n-1)} \end{gather*}\]

Method 12.8. 设\(X_1, X_2, \dots, X_{n}\)独立同分布，将满足某条件的样本值转换为\(1\)，不满足的转换为\(0\)。