2011-03-23

自由度が高くなると、関係が薄い代用マーカーの方がパワーをより高く上げる

昨日まで、検出力と非心カイ二乗分布のことを書いている

こんな場合を考えよう
あるマーカーが真の関係があるとする
その真のマーカーだけで関連検定をするときのパワーは簡単に計算できる
非心カイ二乗分布を持ち出せば、任意の自由度での検定にあっても、そのパワーが簡単に算出できる
さて、あるマーカーが真であって、それとある程度の相関のあるマーカーがあるとする
真のマーカーとその関連マーカーのどれであってもいいから、どれかで、関連がでないかなー、というタイプの「検定繰り返し」を行っているものとする
これは、多型と形質の関連検定について、複数の遺伝形式で検定する場合にあてはまる
これは、あるDNAセグメント上にある複数の多型で検定を行って、そのうちのどれでもいいから、「相加モデル(ある一つの遺伝形式)」で検定してやる場合にもあてはまる
これは、あるDNAセグメント上の複数の多型について、複数の遺伝形式で検定するという２重の「検定繰り返し」の場合にもあてはまる
さて。
これを考えていくために、単純化しよう
真のマーカーにおける、真の関連の強さを設定する
真のマーカーと相関のあるマーカーについて、その相関の強さをパラメタ化する
真のマーカーと相関の強さが同じマーカーがあるとき、データの次元に応じて異なる方向(相互に直交する異なる方向)に別のマーカーを置くこととして、これを「多次元空間における、代用マーカーの配置」のモデルとする
考えるべき要素は
- 真の関連の強さ
- 検定閾値
- 代用マーカーと真のマーカーの相関の強さ
- 空間の次元
代用マーカーと真のマーカーの相関の強さは、「強すぎても弱すぎても寄与は小さく、ちょうどいい値」がある
次元が大きくなると、代用マーカーがパワーを押し上げる程度は強くなる
代用マーカーと真のマーカーの相関の強さのうち、「もっとも影響が大きくなる相関の強さ」は、次元の増大によって、「より弱い相関の代用マーカーが有用」になる(らしい)
もちろん、真の関連の強さが、検定閾値よりはるかに大きいとき(逆にはるかに小さいとき)は、その他の因子を振っても影響は小さく、パワーは、それぞれ１と０

# a 検定半径 a^2:検定カイ二乗値
# b 対立仮説半径 b^2:対立仮説最頻カイ二乗値
# df 自由度・次元
# t 0<= theta <= pi/2

as<-sqrt(c(2^(2)))
#as[1]<-0.000001
bs<-sqrt(c(1.1))
#bs[1]<-as[1]
#dfs<-c(2^(0:3))
dfs<-c(10)
ts<-seq(from=0,to=1,length.out=30)*pi/2
nperm<-10000
out0323.7<-array(0,c(length(dfs),length(as),length(bs)))
out0323.8<-array(0,c(length(dfs),length(as),length(bs),length(ts)))

for(i in 1:length(dfs)){
	for(j in 1:length(as)){
		for(k in 1:length(bs)){
			out0323.7[i,j,k]<-pchisq(as[j]^2,dfs[i],as[j]^2*bs[k]^2,lower.tail=FALSE)
			for(i2 in 1:length(ts)){
				df<-dfs[i]
				Kv<-c(1,rep(0,df-1))
				Tiis<-matrix(c(1,rep(0,df-1)),nrow=1)
				#c1<-qchisq(p,df=dfs[i],lower.tail=FALSE)
				Ks<-as[j]*bs[k]
				CheckDists<-as[j]

				if(df>1){
					for(i3 in 2:df){
						tmpt1<-rep(0,df)
						tmpt1[1]<-cos(ts[i2])
						tmpt1[i3]<-sin(ts[i2])
						tmpt2<-tmpt1
						tmpt2[i3]<--tmpt1[i3]
						Tiis<-rbind(Tiis,tmpt1,tmpt2)
					}
				}

				SPout<-SpherePower(df=df,Kv=Kv,Ks=Ks,CheckDists=CheckDists,Tiis=Tiis,nperm=nperm)
				#matplot(t(SPout$p.out),type="l")
				#SPouts[i,,]<-SPout$p.out
				out0323.8[i,j,k,i2]<-SPout$p.out
			}
			plot(out0323.8[i,j,k,])
		}
	}
}

plot(out0323.4,type="l")
plot(c(out0323.4,out0323.6[1,,,],out0323.6[2,,,],out0323.8),type="l")

which(out0323.4==max(out0323.4))
which(out0323.6[1,,,]==max(out0323.6[1,,,]))
which(out0323.6[2,,,]==max(out0323.6[2,,,]))
which(out0323.8==max(out0323.8))