ヘテロ型の数

• アレル数 naのとき、ディプロタイプ数は
• ディプロタイプ数 ng のとき、観察数は ng
• 総人数はN
• ここで、na個のアレルのアレル頻度を p1,,,pnaとして推定します。
• ここから期待ディプロタイプ人数がng個の数として推定される。
• 観察数 ng があって、その内訳を自由に変えてもよいときには、総人数は固定し

ているので 自由度がng-1。

• アレル頻度を自由に変えてよいとすると、na-1の自由度がある。
• (ng-1)-(na-1)=ng-na がHWE検定の自由度。
• を使えば
• が自由度
• na=1のときは1,na=3のときは3。
• はヘテロなディプロタイプの数。
• もしくは、アレルのペアの取り方の場合の数。
• ヘテロの取り分け方の場合の数が自由度、というように言い換えることもできる。
• 一応、これを任意のnaについて確かめる。

library(MCMCpack)
# 観測試行数
N<-300
# 観測サンプル数
N.sample<-100
# 試すアレル数
nas<-2:7
# hwe-chisq 統計量の格納用行列
hwe.chi2s<-matrix(0,length(nas),N)
# hwe-chisq テスト用の自由度格納行列
# シミュレーションで想定したアレルのすべてが必ずしも観測されるわけではないので
# 観測されたアレル数を基に、自由度を定めるのが「本当」だろう
hwe.dfs<-hwe.chi2s
# アレル数を変えてループ
for(naloop in 1:length(nas)){
	# このループでのアレル数
	na<-nas[naloop]
	# ジェノタイプ数
	ng<-na*(na+1)/2
	# アレル頻度を発生
	ps<-rdirichlet(1,rep(1,na))
	# このループでのhwe-chisq 統計量の格納ベクトル
	hwe.chi2<-rep(0,N)
	# このループでのhwe.df 格納ベクトル
	hwe.df<-rep(0,N)
	# N回試行
	for(i in 1:N){
		# HWEを仮定してランダムサンプリング
		a1<-sample(0:(na-1),N.sample,prob=ps,replace=TRUE)
		a2<-sample(0:(na-1),N.sample,prob=ps,replace=TRUE)
		# ひとまずa1, a2の順序を区別してジェノタイプを数える
		g.count.mat<-matrix(0,na,na)
		# アレルの観測数を数える
		a.count<-rep(0,na)
		for(j in 1:N.sample){
			g.count.mat[a1[j]+1,a2[j]+1]<-g.count.mat[a1[j]+1,a2[j]+1]+1
			a.count[a1[j]+1]<-a.count[a1[j]+1]+1
			a.count[a2[j]+1]<-a.count[a2[j]+1]+1
		}
		# 観測アレル頻度
		a.freq<-a.count/sum(a.count)
		# a1,a2の順序を区別したジェノタイプの期待数を計算する
		g.exp.mat<-outer(a.freq,a.freq,FUN="*")*N.sample
		# heteros to be jointed
		# 行列の対称なセルの表すヘテロジェノタイプは同一視する
		g.count.mat2<-g.count.mat+t(g.count.mat)
		# 対角成分(ホモ)の増えすぎを是正する
		diag(g.count.mat2)<-diag(g.count.mat)
		# 長さngのベクトルが観測ジェノタイプテーブル
		g.table<-c(diag(g.count.mat2),g.count.mat2[upper.tri(g.count.mat2)])
		# 期待ジェノタイプ行列でも同じことをする
		g.exp.mat2<-g.exp.mat+t(g.exp.mat)
		diag(g.exp.mat2)<-diag(g.exp.mat)
		e.table<-c(diag(g.exp.mat2),g.exp.mat2[upper.tri(g.exp.mat2)])
		# 期待値が非０のジェノタイプのみでHWEする
		nonzeros<-which(e.table!=0)
		hwe.chi2[i]<-sum((g.table[nonzeros]-e.table[nonzeros])^2/e.table[nonzeros])
		tmp.na<-length(which(a.freq!=0))
		hwe.df[i]<-tmp.na*(tmp.na-1)/2
	}
	hwe.chi2s[naloop,]<-hwe.chi2
	hwe.dfs[naloop,]<-hwe.df
	hwe.dfs.x[naloop,]<-hwe.dfx
	#plot(sort(hwe.chi2))
	par(mfcol=c(1,2))
	# 自由度を「あるべきアレル数」とする
	plot(sort(pchisq(hwe.chi2,df=na*(na-1)/2)))
	plot(sort(pchisq(hwe.chi2,df=hwe.df)))

}