Type S error: ベイズ的多重検定 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2016-08-24

Bayesian inference completely solves the multiple comparisons problem « Statistical Modeling, Causal Inference, and Social Science にType S error (Computational Statistics (2000) 15: 373.) という論文がある。
古典的な仮説検定では $H_0:\theta_1=\theta_2$ を帰無仮説とし、 $H_1:\theta_1\not{=}\theta_2$ の対立仮説を導く。 $\theta_1-\theta_2$ が信頼区間で0をまたぐかまたがないかが興味の焦点である。
ここで、Type S error (sign, 符号のS) と多重検定を考えようという話。
　
いま、 $J$ 個の(独立な)実験系があるとする。 $j$ 番目の実験に標本数 $n_j$ 個のデータが確認されている。ここで、 $n_j$ の統計量(適当に平均値とする)を $y_j$ で表されているものとする。
ここで、 $y_j$ は実験から得たデータなので、適当な真の効果量 $\theta_j$ と分散 $\sigma^2$ を用いて正規分布
$y_j|\theta_j, \sigma\sim N(\theta_j, \sigma^2)$
からサンプリングされている。また、 $\theta_j$ については適当な平均 $\mu$ と分散 $\tau^2$ を用いて
$\theta_j|\mu, \tau\sim N(\mu, \tau^2)$
とサンプリングされているものとする。以下のシミュレーションでは $\mu=0$ である。
$\theta_j$ は $j$ 番目の実験の真の効果量を適当にランダムサンプリングしている。例えばマイクロアレイで言えば、0ならば全く発現の変動がなく、正ならば高発現で逆もしかり、といったところだろうか。
　
古典的な仮説検定では、観測されたデータ $y_j$ がある幅(信頼区間)をはみ出ていれば有意に差がある、というロジックになる。このときの閾値は
$|y_j|>1.96\sqrt{2}\sigma$
である。1.96 は有意水準5% のときの正規分布の確率なので適当である。
ベイズ的な検定においては $\theta_j$ は
$\theta_j|y,\sigma,\mu,\tau\sim N(\hat{\theta}_j, V_j)$
ただし $\hat{\theta}_j=\frac{\frac{1}{\sigma^2}y_j + \frac{1}{\tau^2}\mu}{\frac{1}{\sigma^2}+\frac{1}{\tau^2}}$ , $V_j=\frac{1}{\frac{1}{\sigma^2}+\frac{1}{\tau^2}}$ .
から得られて、閾値は
$|y_j|>1.96\sqrt{2}\sigma\sqrt{1+\frac{\sigma^2}{\tau^2}}$
となって、古典的な閾値より大きい。
　
論文の図をそのまま再現するとこうなる。左の古典的検定ではy軸の水平線(閾値)が $\sigma$ , $\tau$ に依存せずに固定( $1.96\sigma$ あたり)されているのに対して、ベイズ的なのは $\frac{\tau}{\sigma}$ 比が大きくなるに従って狭くなっている(動かしているのは $\tau$ )。

　
$\frac{\tau}{\sigma}$ 比をいろいろ変えてシミュレーションしたときに、有意と言ってしまう確率(というか、 $J$ 個の実験のうち有意と言ってしまう実験の数)を再現してみる。古典的なほうは、 $\tau=0$ のときも0.05の確率で有意と言ってしまう。これは帰無仮説検定の定義から正しい。
$\tau=0$ のとき、Rでやればすぐにわかるがすべての $\theta$ は0である。というのも、 $\theta_j|\mu, \tau\sim N(\mu=0, \tau^2=0)$ だからである。
すべての実験において、 $\theta_j=0$ であることがわかっている。ここで $y_j$ は、実験をして少しランダムエラーのはいった $y_j\sim N(0, \sigma^2)$ であるから、差がないという帰無仮説が真であっても、仮説検定上、0.05の確率で有意になる。
一方のベイズは $|y_j|>1.96\sqrt{2}\sigma\sqrt{1+\frac{\sigma^2}{\tau^2}}\to\infty$ であり、 $\theta=0$ ならば棄却しようがない。
ということでconservative と言っている。
$\frac{\tau}{\sigma}$ が大きくなれば古典的もベイズ的もどちらも $J$ 個の実験が有意と言ってしまう確率は1に近づく。というのも、 $\frac{\tau}{\sigma}$ が大きくなれば $\frac{\sigma}{\tau}\to 0$ なので、ベイズの閾値が古典的な閾値と同じになる。ただし、 $\frac{\tau}{\sigma}$ がものすごく大きくなるような実験系、というのは $\tau$ が大きいデータだが、そうなれば $\theta_j$ がすごい大きかったり小さかったりとばらつきが大きいことになるので、そもそも検定を持ちださなくても(統計学的には)有意であることが明らかだと思う。

　
本当はもっと書きたかったけど書いている最中にこんがらがってきて忘れてしまった…
Type S error そのものについて書いてなかったけど、符号が違う場合を検出しよう、という話。
ほかにもType M error というのも提唱していて、それはMagnitude のM。

sigma <- 1
tau <- .5
N <- 1e5
alpha <- 0.05

par(mfrow=c(3, 2), las=1, cex.axis=1.5)
for(tau in c(0.5, 1, 2)){

theta <- rnorm(N, 0, tau)
y <- theta + rnorm(N, 0, sigma)
signif_classical <- abs(y) > qnorm(1-alpha/2)*sigma # 頻度主義的に有意なもののindex

mean(abs(y)[signif_classical]) # 有意なものの y 統計量平均
mean(abs(theta)[signif_classical]) # 有意なもののtheta 平均
mean((sign(theta)!=sign(y))[signif_classical]) # 有意なものの符号違い Type S error

xl <- yl <- c(-1, 1)*10
plot(theta, y, xlab=expression(theta), col=signif_classical+2, pch=16, xlim=xl, ylim=yl,
     main=substitute("Classical "*over(tau, sigma)==list(x), list(x=round(tau/sigma, 2))))
abline(0, 1, lty=3, h=c(-1, 1)*qnorm(1-alpha/2)*sigma, v=0)
led(cex=1.5)

theta_hat_bayes <- y * (1/sigma^2) / (1/sigma^2 + 1/tau^2)
theta_se_bayes <- sqrt(1 / (1/sigma^2 + 1/tau^2))
signif_bayes <- abs(theta_hat_bayes) > qnorm(1-alpha/2)*sqrt(2)*theta_se_bayes

mean(abs(theta_hat_bayes)[signif_bayes])
mean(abs(theta)[signif_bayes])
mean((sign(theta)!=sign(theta_hat_bayes))[signif_bayes])

plot(theta, y, xlab=expression(theta), col=signif_bayes+2, pch=16, xlim=xl, ylim=yl,
     main=substitute("Bayesian "*over(tau, sigma)==list(x), list(x=round(tau/sigma, 2))))
abline(0, 1, lty=3, h=c(-1, 1)*qnorm(1-alpha/2)*sqrt(2)*sigma*sqrt(1 + (sigma/tau)^2), v=0)
led(cex=1.5)

}

# legend プロット用
led <- function(cex=1){
  legend("topleft", legend="Type S error", bty="n", cex=cex)
  legend("topright", legend="No Type S error", bty="n", cex=cex)
  legend("right", legend="No claim with\nconfidence", bty="n", cex=cex)
  legend("left", legend="No claim with\nconfidence", bty="n", cex=cex)
  legend("bottomright", legend="Type S error", bty="n", cex=cex)
  legend("bottomleft", legend="No Type S error", bty="n", cex=cex)
}

# tau / sigma 比のシミュレーション
sim <- function(N, tau, sigma, alpha=0.05){
  theta <- rnorm(N, 0, tau)
  y <- theta + rnorm(N, 0, sigma)
  signif_classical <- abs(y) > qnorm(1-alpha/2)*sigma # 頻度主義的に有意なもののindex
  theta_hat_bayes <- y * (1/sigma^2) / (1/sigma^2 + 1/tau^2)
  theta_se_bayes <- sqrt(1 / (1/sigma^2 + 1/tau^2))
  signif_bayes <- abs(theta_hat_bayes) > qnorm(1-alpha/2)*sqrt(2)*theta_se_bayes
  sapply(list(Freq=signif_classical, Bayes=signif_bayes), mean)
}

taus <- seq(0, 5, by=0.01)
s <- mapply(sim, N=N, sigma=sigma, tau=taus)
matplot(taus/sigma, t(s), xlab=expression(over(tau, sigma)), ylab="Probability of claiming significance", type="l", lty=1, lwd=3)
legend("bottomright", legend=rownames(s), lty=1, lwd=3, col=1:2)