FDR: Benjamini-Hockberg

• 昨日の記事はKnockoff 変数を用いたFDRの制御の話
• FDRといえば、Benjamini-Hochbergもある
• これは、「ある閾値で変数の取捨選択をするとする」ときに、すべの変数が帰無仮説ＯＫ変数だったとした場合に、何個の変数がFalselyに帰無仮説を棄却するかの期待個数を算出する
• その数を、実際に帰無仮説を棄却した変数の数で割ってやる。この中には、真に関連のある変数が入っているので、全部帰無の仮定の下での期待個数より大きいと考えられる
• この割合をある値q(たとえば0.05)以下にしてくれるような閾値のうち、もっとも小さい値をとろう、という戦略
• そうすることで、その割合(FDR)をq以下にコントロールしつつ、閾値を越える変数の数を最大にできる、という作戦
• p個の変数が独立であることを前提に考えている
• 以下、Rでやってみる
  • p個の変数があり、k個は真に関連があり、残りは帰無
  • 統計量が正規乱数で得られるとする(ほかの統計量でももちろんよい)。p値でやっても、もちろんよい
  • 閾値をどうするかを探索するためのプロットが黒と赤水平線
  • 実際に偽陽性の割合を緑でプロット
  • 黒カーブが赤水平線と交わるところを閾値にすることになり、そこより右側では、緑は赤水平線より下に来る。これがFDRをコントロールしてある、という意味

p <- 10000
k <- 200

z <- rnorm(p)
for(i in 1:k){
	z[i] <- rnorm(1,rnorm(5,2),1)*sample(c(-1,1),1)
}

t <- seq(from=0,to=5,length=1000)

q <- 0.05

x <- rep(0,length(t))
y <- rep(0,length(t)) # no. selected features
w <- rep(0,length(t)) # no. true positive
for(i in 1:length(x)){
	x[i] <- p * pnorm(t[i],lower.tail=FALSE)*2/max(sum(abs(z)>=t[i]),1)
	
	y[i] <- sum(abs(z)>=t[i])
	tmp <- which(abs(z)>=t[i])
	w[i] <- y[i]-sum(tmp %in% 1:k)
	
}

plot(t,x)
abline(h=q,col=2)
y. <- y
y.[which(y==0)] <- 1
points(t,w/y.,col=3)