【Note】 因果推断(二)——Fisher 精确检验

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(十分严重)。

数据的概览

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

六个具体的样本

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

检验步骤

确定零假设(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)\)，也就是说处理完全没有作用，填入上表得

(注: 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\) 值

以此计算 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 提出的重复采样方法中得到了解决。