2013-08-11

BH法とlocal FDR法

FDRでは複数の検定統計量が与えられたときに、個々の検定について、「帰無仮説に合致しているか」「対立仮説に合致しているか」をfdr : 1-fdrという値で示す
fdrが小さいほど、対立仮説により強く合致することを示す(検定p値も小さいほど、対立仮説らしさが高い)
以下に２法示す
BH法は複数のp値から、それぞれのfdrを出す
local fdr法は統計量から、それぞれのfdrを出す

library(locfdr)
# 標準正規乱数に偏った値も加えてFDR用のz値分布を作る
zex <- rnorm(1000)
zex <- c(zex,rnorm(200,3))
# 両側で正規分布に照らしてp値にする
ps <- 1-2*abs(pnorm(zex)-0.5)
par(mfcol=c(2,3))
# 出したp値をBH法でFDR
p.bh <- p.adjust(ps,method="BH")
# 補正前p値とBH補正FDR値を比べる
plot(ps,p.bh,xlim=c(0,1),ylim=c(0,1))

# 正規分布を仮定してlocal fdr
w0 <- locfdr(zex,nulltype=0)
# z値のヒストグラム
# 帰無仮説合致分の標準正規分布(青破線)
# 推定した混合分布(緑実線)
# ヒストグラムのbinごとに対立仮説が真であると推定される分(紫)
# fdrが0.2未満となる閾値(黄色三角)
# 補正前p値とBH補正FDR値を比べる
plot(ps,w0$fdr,xlim=c(0,1),ylim=c(0,1))
# BH法とlocal FDR法(正規分布仮定)とを比べる
# シミュレーションした対立仮説はz値が正であるので小さいp値はz>0に集中しており、
# 左右非対称であるので、BH法とFDR法とのfdr値を比べたプロットはz 正負の２本になっている(対立仮説を正負の両方に置くと、よくわかる。fdr < 0.2 の閾値三角形も左右に現れる→後述)
plot(w0$fdr,p.bh,xlim=c(0,1),ylim=c(0,1))
# fdr値の分布をBH法、local fdr法で示す
hist(p.bh)
hist(w0$fdr)
par(mfcol=c(1,1))

local fdrの閾値の設定が両側で行われていることを明示するために、対立仮説を両側に立ててやってみる

# 正負のzに対立仮説を置く
zex <- rnorm(1000)
zex <- c(zex,rnorm(200,3))
zex <- c(zex,rnorm(100,-4))
ps <- 1-2*abs(pnorm(zex)-0.5)
par(mfcol=c(2,3))
p.bh <- p.adjust(ps,method="BH")
plot(ps,p.bh,xlim=c(0,1),ylim=c(0,1))
w0 <- locfdr(zex,nulltype=0)
# fdrが0.2未満となる閾値(黄色三角)が２か所にある
plot(ps,w0$fdr,xlim=c(0,1),ylim=c(0,1))
plot(w0$fdr,p.bh,xlim=c(0,1),ylim=c(0,1))
hist(p.bh)
hist(w0$fdr)
par(mfcol=c(1,1))

統計量が正規分布型でないと推定がうまくいかないが、たとえばカイ二乗統計量で得られていたら、それを対応するz値にして、local fdrを実施することは可能

# 自由度dfのカイ二乗乱数とその対立仮説分を発生
df <- 2
chi2<- rchisq(1000,df)
chi2 <- c(chi2,rchisq(200,df,10))
# 一度、p値にしてから
ps <- pchisq(chi2,df,lower.tail=FALSE)
# 正規分布対応値にする
zex <- qnorm(ps/2)
# 正負の偏りをランダムな正負割り付けで解消する
zex <- zex*sample(c(-1,1),length(zex),replace=TRUE)
# 念のため確認する
# plot(chi2,zex)
# 後の処理は同じ
par(mfcol=c(2,3))
p.bh <- p.adjust(ps,method="BH")
plot(ps,p.bh,xlim=c(0,1),ylim=c(0,1))
w1 <- locfdr(zex,nulltype=0)
# fdrが0.2未満となる閾値(黄色三角)が２か所にある
plot(ps,w1$fdr,xlim=c(0,1),ylim=c(0,1))
plot(w1$fdr,p.bh,xlim=c(0,1),ylim=c(0,1))
hist(p.bh)
hist(w1$fdr)
par(mfcol=c(1,1))

やや詳しい話はこちらとそこからのリンクで