最大尤度をもたらす変数の値を見る

  • SNPのディプロタイプ頻度をg(1:3)=c(p^2+p(1-p)f,2p(1-p)(1-f),(1-p)^2+p(1-p)f)として表したときに、観察3ディプロタイプ人数G(1:3)を用いてその対数尤度はLL=\sum_{i=1}^3 G(i)*log(g(i))と表される
  • HWEを仮定すると変数fは0に固定することになるので、変数はpのみ
  • HWEの仮定とHWDの仮定では変数の個数が2,1と、その差が1。このためHWE検定の自由度は1
  • 今適当に、p,fを与えて、Gを作り、LLを最大にするp,fをグリッド計算で出してみよう
# 「真実データ」を作る
p<-0.3
f<-0.1
N<-100
g<-c(p^2+p*(1-p)*f,2*p*(1-p)*(1-f),(1-p)^2+p*(1-p)*f)*N

# グリッド探索用にp,fの変数を取る
Ps<-seq(from=0.01,to=0.99,by=0.01)
fs<-seq(from=-0.99,to=0.99,by=0.01)

Pfs<-expand.grid(Ps,fs)

Ls2<-rep(-Inf,length(Pfs[,1]))
for(i in 1:length(Pfs[,1])){
	tmpp<-Pfs[i,1]
	tmpf<-Pfs[i,2]
	tmpG<-c(tmpp^2+tmpp*(1-tmpp)*tmpf,2*tmpp*(1-tmpp)*(1-tmpf),(1-tmpp)^2+tmpp*(1-tmpp)*tmpf)
	if(min(tmpG>0)){
		Ls2[i]<-sum(g*log(tmpG))
	}
	

}
library(rgl)
# plotする
Ls2[which(Ls2==-Inf)]<-min(Ls2[which(Ls2!=-Inf)])
plot3d(Pfs[,1],Pfs[,2],Ls2)
# 最大尤度のポイントの変数の値を確認する
maxpoint<-which(Ls2==max(Ls2))
Pfs[maxpoint,]

# こうすると、最大尤度を出すp,fは、初期値で与えたアレル頻度とHWDの係数になる