em for recessive allele

The following code will estimate the recessive allele frequency from the counts of the dominant phenotype assuming Hardy Weinberg (see slides 80-81).

n<-1000
naa<-23
dat<-c(n-naa,naa)
n<-sum(dat)
(pa1<-(dat[2]/n)^.5)

## [1] 0.1516575

em.ra<-function(dat,ini=NULL,print.it=TRUE,eps=0.0001){
  if (is.null(ini)) pa2<-runif(1) else pa2<-ini 
  dif<-1
  acc<-0
  while (dif>=eps){
    acc<-acc+1
    pa2n<-2*(n*pa2+dat[2])/(2*n*(1+pa2))
    dif<-abs(pa2n-pa2)
    if (print.it) cat(paste(acc," : ",round(pa2,digits=3)," -> ",round(pa2n,digits=3),"\n"))
    pa2<-pa2n
  }
  pa2
}

em.ra(dat)

## 1  :  0.196  ->  0.183 
## 2  :  0.183  ->  0.174 
## 3  :  0.174  ->  0.168 
## 4  :  0.168  ->  0.163 
## 5  :  0.163  ->  0.16 
## 6  :  0.16  ->  0.158 
## 7  :  0.158  ->  0.156 
## 8  :  0.156  ->  0.155 
## 9  :  0.155  ->  0.154 
## 10  :  0.154  ->  0.153 
## 11  :  0.153  ->  0.153 
## 12  :  0.153  ->  0.153 
## 13  :  0.153  ->  0.152 
## 14  :  0.152  ->  0.152 
## 15  :  0.152  ->  0.152 
## 16  :  0.152  ->  0.152 
## 17  :  0.152  ->  0.152

## [1] 0.1518697