父権否定確率

DNA鑑定の続き。
父権否定確率なるものを求める。
母子関係が確定した母子を考える。
子について、母から受け継がれたアレルと父から受け継がれたアレルがあるが、父から受け継がれるアレルであり得ないパターンを考える。
例えば、母PQ、子PPとしよう。母からは必ずPを受け継ぐ。父からはPを受け継がないといけないが、Pが受け継がれるというのを否定しようとすると、アレル頻度というプールでないところから取ってくればいい、というのは余事象で、染色体が2つであることを考慮すると

複数のローカスのマーカーを考えると、こういうときは連鎖しない独立なところを選んであるので

library(MCMCpack)
marker<- 20                                    # 調べるローカス数。各ローカスは独立と仮定する。
allele<- 10                                    # 存在するアレル数。
A<- rdirichlet(marker,rep(1,allele))           # アレル存在頻度。行がローカス、列がアレル番号。
dimnames(A)<- list(1:marker,1:allele)
names(dimnames(A))<- c("No.marker","allele probability")
F<- function(marker,allele,A){
al<- function(vec){
  sample(1:length(vec),size=1,prob=vec)
}
MO1<- apply(A,1,al)
MO2<- apply(A,1,al)
AF<- apply(A,1,al)
data<- rbind(MO1,MO2,MO1,AF)
negative.prob<- rep(0,marker)
for(d in 1:ncol(data)){
  if(length(unique(data[,d]))==3){
    negative.prob[d]<- A[d,data[4,d]]
  }else{
     p<- (data[3,d]==data[4,d])+1
     switch(p,
            switch((data[1,d]==data[2,d])+1,
                   negative.prob[d]<- A[d,data[1,d]]+A[d,data[2,d]],
                   negative.prob[d]<- A[d,data[4,d]]
                   ),
            negative.prob[d]<- A[d,data[1,d]] 
            )
   }
}
#negative.prob
PE<-  1-prod((1-(1-negative.prob)^2)^2)
#PE
result<- matrix(c(marker,allele,PE),nr=1)
colnames(result)<- c("marker","allele","PE")
return(result)
}

allele<- 50
marker<- 50
R<- matrix(0,marker/5,allele/5)
for(a in 1:(allele/5)){
  for(m in 1:(marker/5)){
    Data<- rep(0,1000)
    for(k in 1:1000){
      A<- rdirichlet(m,rep(1,a)) 
      Data[k]<- F(m,a,A)[3]
    }
    R[m,a]<- mean(Data)
  }
}
dimnames(R)<- list((1:(marker/5))*5,(1:(allele/5)*5))
names(dimnames(R))<- c("No.marker","No.allele")
persp(R,col=rainbow(9),
      xlim=c(0,1),ylim=c(0,1),zlim=c(0,1),theta=-60,phi=30,
      xlab="No.locus",
      ylab="No.allele",
      zlab="Negative probability")
R
         No.allele
No.marker 5        10        15        20        25        30        35        40        45        50
       5  1 0.8709114 0.7507584 0.6627568 0.5590630 0.5077778 0.4545591 0.4069119 0.3777021 0.3555796
       10 1 0.9827929 0.9330957 0.8774903 0.8140328 0.7554255 0.7040656 0.6601192 0.6003302 0.5723280
       15 1 0.9985821 0.9846393 0.9588770 0.9208368 0.8803012 0.8392827 0.8032682 0.7667504 0.7167236
       20 1 0.9997170 0.9954594 0.9861300 0.9678163 0.9405994 0.9095452 0.8827655 0.8471607 0.8244041
       25 1 0.9999476 0.9987473 0.9949791 0.9855154 0.9715579 0.9552347 0.9294924 0.9027719 0.8799162
       30 1 0.9999986 0.9998100 0.9984151 0.9935364 0.9865050 0.9751875 0.9597047 0.9408409 0.9223335
       35 1 0.9999997 0.9999288 0.9993800 0.9974552 0.9929459 0.9858177 0.9757530 0.9648577 0.9477928
       40 1 1.0000000 0.9999871 0.9997760 0.9991394 0.9965889 0.9924593 0.9865758 0.9772950 0.9672214
       45 1 1.0000000 0.9999979 0.9999376 0.9995749 0.9982885 0.9956014 0.9915084 0.9850804 0.9774683
       50 1 1.0000000 0.9999992 0.9999729 0.9997571 0.9990701 0.9976638 0.9955583 0.9916851 0.9856949

アレルが多すぎると否定確率は小さくなるのか…
rdirichletで振りすぎると、1アレルあたりの頻度が小さくなって、余事象を取ると大きくなって、かけ合わせまくってもなかなか小さくならないからか…
dataというのが母子関係。

 data
     1 2  3 4 5  6 7  8 9 10 11 12 13 14 15 16 17 18 19 20
MO1  8 2 10 3 2 10 4  2 7  3  8  7  3  8  1  2 10  1  9  1
MO2 10 7  6 9 1  1 4  9 8  3  6 10  4  3  8  2 10  2  3  4
MO1  8 2 10 3 2 10 4  2 7  3  8  7  3  8  1  2 10  1  9  1
AF   7 2 10 1 3  1 7 10 7  3 10  7  3  7  5  5  7  1  8  8

MO1とMO2というのが母のジェノタイプ。
MO1とAFというのが子のジェノタイプ。
MO1は母から子に伝えられた母方染色体上のローカスとアレル。
AFは擬父からの父方染色体。
母と子のアレルパターンは、ホモ-ホモ、ホモ-ヘテロ、相同なヘテロ-ヘテロ、一方だけ相同なヘテロ-ヘテロで、アレルP、Q、Rを用いて下図の通り