Fisher 精确检验

一句话总结

给定了零假设 (sharp null hypothesis) 后,所有样本的潜在结果 \(Y_i(0), Y_i(1)\) 都已经确定,统计量 \(T\) 的随机性完全由分配机制 \(W_i\) 引入。此时我们可以穷举所有的分配方式并计算对应的概率,然后统计出现比当前观测数据下的统计量更“极端”的概率 (Fisher Exact P-values, FEPs),并以此作为判据接受/拒绝原假设。

使用的数据

这里使用的数据是 Paul et al.(2007) 里关于对上呼吸道感染的儿童的治疗手段的效果的检验。这里我们选择了 \(N_t=35\) 名接受荞麦蜂蜜的儿童作为实验组,\(N_c=37\) 名不治疗的儿童作为对照组。观察的变量中我们选择咳嗽频率和咳嗽强度,这些量都被调整到 0(几乎没有)——6(十分严重)。

数据的概览

变量 均值 标准差 控制组均值 实验组均值
处理前咳嗽频率 (cfp) 3.86 0.92 3.73 4.00
处理后咳嗽频率 (cfa) 2.47 1.61 2.81 2.11
处理前咳嗽强度 (csp) 3.99 1.03 3.97 4.00
处理后咳嗽强度 (csa) 2.54 1.73 2.86 2.20

试验后结果的分层累计分布

控制组 cfa 实验组 cfa 控制组 cfa 实验组 cfa
0 0.14 0.14 0.16 0.17
1 0.19 0.40 0.22 0.46
2 0.32 0.63 0.35 0.54
3 0.73 0.83 0.59 0.77
4 0.89 0.91 0.86 0.91
5 0.92 0.97 0.95 0.94
6 1.00 1.00 1.00 1.00

六个具体的样本

样本 cfa \(Y_i(0)\) cfa \(Y_i(1)\) W_i cfp \(X_i\) cfa \(Y_i^{\mathrm{obs}}\)
1 ? 3 1 4 3
2 ? 5 1 6 5
3 ? 0 1 4 0
4 4 ? 0 4 4
5 0 ? 0 1 0
6 1 ? 0 5 1

我们先使用这六个样本作为例子

检验步骤

确定零假设(Fisher Sharp Null Hypothesis)

Fisher 当时提出的原假设为 \(H_0: Y_i(0)=Y_i(1)\),但很自然的我们可以将它扩展,只要最后能确定所有的潜在结果即可,比如

  • \(H_0: Y_i(1)=Y_i(0)+C\)
  • \(H_0: Y_i(1)=Y_i(0)*C\)
  • \(H_0: Y_i(1)=Y_i(0)+C_i\)

我们使用最经典的原假设 \(H_0: Y_i(0)=Y_i(1)\),也就是说处理完全没有作用,填入上表得

样本 cfa \(Y_i(0)\) cfa \(Y_i(1)\) W_i cfp \(X_i\) cfa \(Y_i^{\mathrm{obs}}\) rank(\(Y_i^{\mathrm{obs}}\))
1 (3) 3 1 4 3 4
2 (5) 5 1 6 5 6
3 (0) 0 1 4 0 1.5
4 4 (4) 0 4 4 5
5 0 (0) 0 1 0 1.5
6 1 (1) 0 5 1 3

(注: rank 的计算方式即直接排序,然后分配123456。如果值相等则取平均数作为rank。如果将其都减去 \((N+1)/2\)后,则其和为0。我们将平移后的结果记为 \(R_i\))

选择统计量 \(T\)

我们先使用 \[T_1(\boldsymbol{W},\boldsymbol{Y}^\mathrm{obs})=\lvert \bar{Y}_t^\mathrm{obs} - \bar{Y}_c^\mathrm{obs}\rvert=\lvert \frac83-\frac53\rvert=1\]\[T_2(\boldsymbol{W},\boldsymbol{Y}^\mathrm{obs})=\lvert \bar{R}_t^\mathrm{obs} - \bar{R}_c^\mathrm{obs}\rvert=\lvert \frac{11.5}3-\frac{9.5}3\rvert=\frac{2}{3}\]

下表展示了\(\binom{6}{3}=20\)中情况下的 \(T\)

\(W_1\) \(W_2\) \(W_3\) \(W_4\) \(W_5\) \(W_6\) \(T_1\) \(T_2\)
0 0 0 1 1 1 -1.00 -0.67
0 0 1 0 1 1 -3.67 -3.00
0 0 1 1 0 1 -1.00 -0.67
0 0 1 1 1 0 -1.67 -1.67
0 1 0 0 1 1 -0.33 0.00
0 1 0 1 0 1 2.33 2.33
0 1 0 1 1 0 1.67 1.33
0 1 1 0 0 1 -0.33 0.00
0 1 1 0 1 0 -1.00 -1.00
0 1 1 1 0 0 1.67 1.33
1 0 0 0 1 1 -1.67 -1.33
1 0 0 1 0 1 1.00 1.00
1 0 0 1 1 0 0.33 0.00
1 0 1 0 0 1 -1.67 -1.33
1 0 1 0 1 0 -2.33 -2.33
1 1 1 1 0 0 0.33 0.00
1 1 0 0 0 1 1.67 1.67
1 1 0 0 1 0 1.00 0.67
1 1 0 1 0 0 3.67 3.00
1 1 1 0 0 0 1.00 0.67

以此计算 P 值,以 \(T_1\) 为例,出现比观测更极端的概率为 \(16/20=0.8\), 显然无法拒绝原假设。使用 \(T_2\) 可以得到同样的结论。

一般的,统计量 \(T(\boldsymbol{W},\boldsymbol{Y}^\mathrm{obs}\boldsymbol{X})\) 需要是一个实值的函数,且只和分配 \(\boldsymbol{W}\) 、观测到的数据 \(\boldsymbol{Y}^\mathrm{obs}\) (注意它是 \(\boldsymbol{W},\boldsymbol{Y}(0),\boldsymbol{Y}(1)\) 的函数)以及实验前得到的数据 \(\boldsymbol{X}\) 有关。下面列举一下常用的统计量

  • \(T^\mathrm{dif}=\lvert \bar{Y}_t^\mathrm{obs} - \bar{Y}_c^\mathrm{obs}\rvert=\lvert\frac{\sum_{i:W_i=1}Y_i^\mathrm{obs}}{N_t}-\frac{\sum_{i:W_i=0}Y_i^\mathrm{obs}}{N_c}\rvert\)
  • \(T^\mathrm{log}=\lvert\frac{\sum_{i:W_i=1}\ln(Y_i^\mathrm{obs})}{N_t}-\frac{\sum_{i:W_i=0}\ln(Y_i^\mathrm{obs})}{N_c}\rvert\)
  • \(T^\mathrm{median}=\lvert \mathrm{med}_t(Y_i^{\mathrm{obs}})-\mathrm{med}_c(Y_i^{\mathrm{obs}})\rvert\)
  • \(T^\mathrm{t-stat}=\lvert \frac{\bar{Y}_t^\mathrm{obs} - \bar{Y}_c^\mathrm{obs}}{\sqrt{s_c^2/N_c+s_t^2/N_t}}\rvert,\,s_t^2=\sum_{i:W_i=1}(Y_i^\mathrm{obs}-\bar{Y}_t^\mathrm{obs})^2/(N_t-1)\)\(s_c^2\) 类似
  • \(T^\mathrm{quant}=\lvert q_{\delta,t}(Y_i^{\mathrm{obs}})-q_{\delta,t}(Y_i^{\mathrm{obs}})\rvert\), \(\delta_{q,c}\) 是经验累计分布 (e.c.d.f) 的 \(\delta\) 分位数
  • \(T^\mathrm{rank}=\lvert \bar{R}_t - \bar{R}_c\rvert=\lvert\frac{\sum_{i:W_i=1}\bar{R}_i}{N_t}-\frac{\sum_{i:W_i=0}\bar{R}_i}{N_c}\rvert\rvert\)

除此之外还有一些基于模型的统计量,如

  • \(Y_i(0)\sim\mathcal{N}(\mu_c,\sigma_c^2),\,Y_i(1)\sim\mathcal{N}(\mu_t,\sigma_t^2)\)

    \(T^{\mathrm{model}}=\lvert\widehat{\mu}_t-\widehat{\mu}_c\rvert=\lvert\bar{Y}_t^\mathrm{obs} - \bar{Y}_c^\mathrm{obs}\rvert=T^\mathrm{dif}\)

  • \(\log Y_i(0)\sim\mathcal{N}(\mu_c,\sigma_c^2),\,\log Y_i(1)\sim\mathcal{N}(\mu_t,\sigma_t^2)\)

    \(T^{\mathrm{model}}=\lvert\widehat{\mu}_{mle,t}-\widehat{\mu}_{mle,c}\rvert\)

  • 还可以用非参数的模型,如 Kolmogorov-Smirnov Statistic

    \(T^{\mathrm{ks}}=\sup\limits_y\lvert\widehat{F}_t(y)-\widehat{F}_c(y)\rvert=\max\limits_{i=1,\ldots,N}\lvert\widehat{F}_t(Y_i^\mathrm{obs})-\widehat{F}_c(Y_i^\mathrm{obs})\rvert\)

    其中 \(\widehat{F}_c(y)=\frac{1}{N_c}\sum_{i:W_i=0}\boldsymbol{1}_{Y_i^\mathrm{obs}\leq y}\) ,即经验累计分布函数 (e.c.d.f)。\(\widehat{F}_t(y)\) 类似

原书里还介绍了其他的统计量,以及统计量间的组合,不一而足,感兴趣的读者可以自行查阅。

关于统计量的选择没有固定的方式。一方面你需要考虑零假设和备择假设下数据的分布,然后尝试选择功效较大的统计量。比如说,你期望处理后增大了结果的分散程度,但没有改变结果的均值,那么选择衡量分散程度的统计量,如方差或者 \(q\) 分位数区间的长度就会是一个具有更高功效的统计量。另一方面你可能需要考察已有数据的信息。比如你发现了观测数据中存在某些离群点,那么选择中位数而非均值作为统计量可能会得到更高的功效。

下面的模拟实验说明了一般来说,基于秩(rank)的统计量是比较好的。他在保持较好的功效时有更好的稳健性,其他的统计量遇到不适合的数据可能出现较差的情况。

  • 实验一:\(Y_i(0)\sim\mathcal{N}(0,1),\,Y_i(1)=Y_i(0)+\tau\)
  • 实验二:在上一个实验的基础上,给 \(20\%\) 的数据加上了 \(5\) 的偏置使其变为离群值
  • 实验三:在第一个实验的基础上,都取\(\exp\),即 \(Y_i(0)\sim\log-\mathcal{N}(0,1),\,Y_i(1)=Y_i(0)\cdot\exp{\tau}\)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
gen.threshold <- function(method, q=0.9, trail=20000, Nc=1000, Nt=1000, seed=42) {
set.seed(seed)
controls = rnorm(Nc+Nt, 0, 1)
if (method == 'outlier') {
w = sample(1:(Nc+Nt), floor((Nc+Nt)/5))
controls[w] = controls[w] + 5
}
treatment = controls
if (method == 'multiplicative') {
controls = exp(controls)
treatment = exp(treatment)
}
T.difs = rep(NA, trail)
T.meds = rep(NA, trail)
T.raks = rep(NA, trail)
for (i in 1:trail) {
w = sample(1:(Nc+Nt), Nt)
T.difs[i] = abs(mean(treatment[w]) - mean(controls[-w]))
T.meds[i] = abs(median(treatment[w]) - median(controls[-w]))
ranks = rank(c(treatment[w], controls[-w]))
T.raks[i] = abs(mean(ranks[1:Nt]) - mean(ranks[(Nt+1):(Nt+Nc)]))
}
return(c(quantile(T.difs, q), quantile(T.meds, q), quantile(T.raks, q)))
}
get.power <- function(tau, thres, method, trail=500, Nc=1000, Nt=1000, seed=42) {
set.seed(seed)
controls = rnorm(Nc+Nt, 0, 1)
if (method == 'outlier') {
w = sample(1:(Nc+Nt), floor((Nc+Nt)/5))
controls[w] = controls[w] + 5
}
treatment = controls + tau
if (method == 'multiplicative') {
controls = exp(controls)
treatment = exp(treatment)
}
T.difs = rep(NA, trail)
T.meds = rep(NA, trail)
T.raks = rep(NA, trail)
for (i in 1:trail) {
w = sample(1:(Nc+Nt), Nt)
T.difs[i] = abs(mean(treatment[w]) - mean(controls[-w]))
T.meds[i] = abs(median(treatment[w]) - median(controls[-w]))
ranks = rank(c(treatment[w], controls[-w]))
T.raks[i] = abs(mean(ranks[1:Nt]) - mean(ranks[(Nt+1):(Nt+Nc)]))
}
return(c(sum(T.difs>thres[1])/trail,sum(T.meds>thres[2])/trail,sum(T.raks>thres[3])/trail))
}
thres1 <- gen.threshold('normal')
thres2 <- gen.threshold('outlier')
thres3 <- gen.threshold('multiplicative')
dat1 <- matrix(NA, nrow=3, ncol=200)
dat2 <- matrix(NA, nrow=3, ncol=200)
dat3 <- matrix(NA, nrow=3, ncol=200)
for (i in 1:200) {
dat1[,i] <- get.power(i/1000, thres1, 'normal')
dat2[,i] <- get.power(i/1000, thres2, 'outlier')
dat3[,i] <- get.power(i/1000, thres3, 'multiplicative')
}
plot(x=1:200/1000, y=dat1[1,], col='red', type='l', ylim=c(0,1), main='normal', xlab='tau', ylab='Power')
lines(x=1:200/1000, y=dat1[2,], col='green', ylim=c(0,1))
lines(x=1:200/1000, y=dat1[3,], col='blue', ylim=c(0,1))
legend('bottomright', legend=c('dif', 'median', 'rank'), col=c('red', 'green', 'blue'), lty=1)

plot(x=1:200/1000, y=dat2[1,], col='red', type='l', ylim=c(0,1), main='outlier', xlab='tau', ylab='Power')
lines(x=1:200/1000, y=dat2[2,], col='green', ylim=c(0,1))
lines(x=1:200/1000, y=dat2[3,], col='blue', ylim=c(0,1))
legend('bottomright', legend=c('dif', 'median', 'rank'), col=c('red', 'green', 'blue'), lty=1)

plot(x=1:200/1000, y=dat3[1,], col='red', type='l', ylim=c(0,1), main='multiplicative', xlab='tau', ylab='Power')
lines(x=1:200/1000, y=dat3[2,], col='green', ylim=c(0,1))
lines(x=1:200/1000, y=dat3[3,], col='blue', ylim=c(0,1))
legend('bottomright', legend=c('dif', 'median', 'rank'), col=c('red', 'green', 'blue'), lty=1)

造“置信区间”

首先明确一点,该置信区间和频率学派的置信区间有差别。事实上 Fisher 方法无法给出传统意义下的置信区间。在这里,我们将其称为 Fisher 置信区间,因为其和 Fisher exact P value 密切相关。 简单说来,其通过不断调整零假设 \(Y_i(1)=Y_i(0)+C\) 来构造出 \(C\) 的一个区间。我们知道,对于每一个给定的 \(C\),所有的潜在结果都被确定,因此我们可以算出出现我们观测的 P-value. 而当 \(C\) 极大或极小时,计算所得的 P-value 都将成为0。因此我们可以取使得 P-value 大于0.05 的那些 \(C\) 值,然后组成一个 \(95\%\) 的“置信区间”,来作为 \(Y_i(1)-Y_i(0)\) 的一个区间估计。

计算 P 值

显然,计算 P 值需要穷举所有的可能性,一共\(\binom{N_t+N_c}{N_t}\)种,当 \(N_t,N_c\) 都不太小时显然是不可计算的。一个简单的近似,也是上面代码里所用到的,就是我随机的采 1000 或 10000 次样本做重复实验,然后根据这些样本算出的 \(T\) 取分位数来构造置信区间。可以证明,给定真实的 P值 \(p^\star\),如果我们抽取 \(K\) 个样本,那么使用这 \(K\) 个样本计算得到标准差为 \(\sqrt{p^\star(1-p^\star)/K}\leq1/\sqrt{4k}\)。即使我们让标准差小于 \(0.0001\),也只需要 \(K=250000\) 次采样,这在现在的电脑上是十分容易实现的。实际应用中,由于你已经采了足够的样本,因此使用 \(\widehat{p}\) 近似真值 \(p^\star\) 来估计方差可以得到更好的效果。

带协变量的检验

之前我们都没有使用实验前观测得到的数据,也就是协变量。加入协变量我们又可以构造出若干统计量。

一个常见的情况是协变量和实验变量是实验前后同一个量的观测,那么做一个减法,即 \(Y_i'(w)=Y_i(w)-X_i\) 是十分自然的,此时构造出一个新的统计量

\[T^\mathrm{gain}=\frac{\sum_{i:W_i=1}(Y_i^\mathrm{obs}-X_i)}{N_t}-\frac{\sum_{i:W_i=0}(Y_i^\mathrm{obs}-X_i)}{N_c}=\bar{Y}_t^{\mathrm{obs}}-\bar{Y}_c^{\mathrm{obs}}-(\bar{X}_t-\bar{X}_c)\]

当然你可以再进行一个”归一化“,简单地除以 \(X_i\) 本身,即\(Y_i'‘(w)=\frac{Y_i(w)-X_i}{X_i}\),此时得到新的统计量

\[T^\mathrm{prop-change}=\frac{1}{N_t}\sum\limits_{i:W_i=1}\frac{Y_i^\mathrm{obs}-X_i}{X_i}-\frac{1}{N_c}\sum\limits_{i:W_i=0}\frac{Y_i^\mathrm{obs}-X_i}{X_i}\]

还有一种可能就是做一个线性回归模型,即

\[(\widehat{\beta}_0,\widehat{\beta}_X,\widehat{\beta}_W)=\arg\min\limits_{\widehat{\beta}_0,\widehat{\beta}_X,\widehat{\beta}_W}\sum\limits_{i=1}^{N}(Y_i^\mathrm{obs}-\widehat{\beta}_0-\widehat{\beta}_X\cdot X_i-\widehat{\beta}_W\cdot W_i)^2\]

得到新统计量 \(T^{\mathrm{reg-coef}}=\widehat\beta_W\)

局限

Fisher Exact Test 是一个十分直观且简单的检验方法,但也带来了他的局限性,主要体现在两点

  • Sharp null hypothesis 带来了太多的限制
  • 拒绝了原假设后,对于组间效应不能给出一个预测量,连点估计都无法给出。

这些局限在下一个方法,也就是 Neyman 提出的重复采样方法中得到了解决。