ABFの関数 - ryamadaの遺伝学・遺伝統計学メモ

背景はこちらの記事とこちらの記事で。
~~注！もろもろ、未検証のβ版です。~~ おそらく大丈夫
まとめなおした記載はこちら
添付は、このベータ版関数の出力について次の３相関をプロットしたもの。上段が、６カウントについて、ロジスティックのｐとトレンド検定のｐとの相関。下段左側が、ロジスティックｐとABFの相関。下段右側が、トレンド検定ｐとABFの相関。
カウントデータファイルから漸近近似ベイズ因子計算の関数
- 第１引数：入力ファイル名
  - ファイルの例

rs1	10	20	30	40	20	10
rs2	11	21	31	41	21	11

- 第２引数：出力ファイル名
  - 出力ファイルの例

"g1"	"g2"	"g3"	"g4"	"g5"	"g6"	"abf"	"theta"	"sd"	"z"	"p"
"10"	"20"	"30"	"40"	"20"	"10"	"0.00142233962067752"	"-1.24514568888400"	"0.254021868403497"	"-4.90172636202394"	"9.49981016296989e-07"
"11"	"21"	"31"	"41"	"21"	"11"	"0.00158057066870834"	"-1.17841828507336"	"0.244047808377548"	"-4.82863703184878"	"1.37470719730343e-06"

- 第３引数：入力ファイルの６カウントのカラム番号
- 第４引数：リスク $\theta$ について正規分布仮定を行い、そのときの、指定リスク(1/g〜g)の範囲にfの割合のリスクローカスが入るとしたときの、g(g>0)
- 第５引数：割合 f
コマンド
- "ABFfromFile("infile.txt","outfile.txt",c(2,3,4,5,6,7),2,0.95)
出力
- 第１〜６カラムにカウント
- 第７カラムにABF
- 第８カラムに $\hat{\theta}$
- 第９カラムに $\sqrt{V}$
- 第10カラムに $Z$
- 第11カラムに $p$

ABFfromFile<-function(file,outfile,columns,g,f){
  ABF<-function(counts,g,f){
   getW<-function(g,f){
    ret<-log(g)/qnorm((1+f)/2,mean=0,sd=1)
    return(ret)
   }
   W<-getW(g,f)
   P<-c(rep(0,counts[1]+counts[2]+counts[3]),rep(1,counts[4]+counts[5]+counts[6]))
   G<-c(rep(0,counts[1]),rep(1,counts[2]),rep(2,counts[3]),rep(0,counts[4]),rep(1,counts[5]),rep(2,counts[6]))
   result<-glm(P~G,family=binomial)
   summary<-summary(result)
   V<-summary$coefficients[4]
   r<-W/(W+V)
   Z<-summary$coefficients[6]
   ABF<-1/sqrt(1-r)*exp(-Z^2/2*r)
   ret<-c(ABF,summary$coefficients[2],summary$coefficients[4],summary$coefficients[6],summary$coefficients[8])
   return(ret)
  }

 data<-read.table(file)
 res<-c()
 label<-c("g1","g2","g3","g4","g5","g6","abf","theta","sd","z","p")
 counter<-0
 res<-c(res,label)
 for(i in 1:length(data$V1)){
  counts<-c(data[i,columns[1]],data[i,columns[2]],data[i,columns[3]],data[i,columns[4]],data[i,columns[5]],data[i,columns[6]])
  result<-ABF(counts,g,f)
  res<-c(res,counts,result)
 }
 mat<-matrix(res,nrow=length(label))
 #mat<-t(mat)
 write.table(t(mat),outfile,sep="\t",col.names=FALSE,row.names=FALSE)
 
 return(mat)
}

ファイルからではなく、６カウントからのABF(とその他の関連統計量の)算出関数(上記の関数 "ABFfromFile"の中にある

ABF<-function(counts,g,f){
 getW<-function(g,f){
 ret<-log(g)/qnorm((1+f)/2,mean=0,sd=1)
 return(ret)
 }

６カウントから、ロジスティック回帰ベースでWを求める関数。関数 ABF 内にある

 getW<-function(g,f){
 ret<-log(g)/qnorm((1+f)/2,mean=0,sd=1)
 return(ret)
 }