hwe.test <- function (loc,loc.2=NULL,idx=1:length(loc)) {
 if (!is.null(loc.2)) loc=paste(loc,loc.2,sep="/")
 counts<-table(loc[idx])
 gtp<-unlist(dimnames(counts))
 alleles<-sort(unique(unlist(strsplit(gtp,"/"))))
 all.genos<-outer(alleles,alleles,paste,sep="/")
 all.genos<-all.genos[!lower.tri(all.genos)]
 g<-t(data.frame(strsplit(all.genos,"/")))
 fullcount<-integer(length(all.genos))
 names(fullcount)<-all.genos
 fullcount[all.genos %in% gtp]<-counts
 n<-length(alleles)
 dum<-matrix(0,nc=n, nr=length(fullcount))
 for(i in 1:n) {
   dum[,i]<-as.integer(alleles[i]==g[,1])+as.integer(alleles[i]==g[,2])
 } 
 colnames(dum)<-as.character(alleles)
 offset<-double(length(fullcount))
 offset[g[,1]!=g[,2]]<-log(2)
 offset<-offset+log(sum(fullcount))
 m0<-glm(fullcount ~ dum-1, offset=offset, family=poisson())
 return(alleles,fullcount,dum,m0)
}
summary.hwe.test <- function(res) {
 x<-summary(res$m0)$coef
 m<-x[,1]
 se<-x[,2]
 df<-summary(res$m0)$df.residual
 print(data.frame(Count=c(res$fullcount),
                  Expected=fitted(res$m0)))
 cat("\nLRTS=",round(res$m0$deviance,digits=3),
         " df=",df,
         " (P=",round(pchisq(res$m0$deviance,df,lower=F),digits=3),")\n\n",
         sep="")
 x<-round(exp(cbind(m,m-1.96*se,m+1.96*se)),4)
 x[x[,2]<0,2]<-0
 x[x[,3]>1,3]<-1
 colnames(x)<-c("Allele Freq","lower 95%CL", "upper 95%Cl")
 rownames(x)<-res$alleles
 print(x)
}

print.hwe.test <- function(res) {
 x<-c(summary(res$m0)$deviance,
      summary(res$m0)$df.residual,
      exp(coef(res$m0)))
 names(x)<-c("LRTS","df",paste("Allele",res$alleles))
 x
}

#
# Exact test for HWE: 2 alleles
#
# See eg Emigh TH.  Comparison of tests for Hardy-Weinberg Equilibrium.  
#  Biometrics 1980; 36: 627-642
#
dhwe2 <- function(n11, n12, n22) {
  f <- function(x) lgamma(x+1)
  n <- n11+n12+n22
  n1 <- 2*n11+n12
  n2 <- 2*n22+n12
  exp(log(2)*(n12) + f(n) - f(n11) - f(n12) - f(n22) - f(2*n) + f(n1) + f(n2))
}

hwe2.test <- function(n11, n12, n22) {
  n1 <- 2*n11+n12
  n2 <- 2*n22+n12
  x12 <- seq(n1 %% 2,min(n1,n2),2)
  x11 <- (n1-x12)/2
  x22 <- (n2-x12)/2
  dist <- data.frame(n11=x11,n12=x12,n22=x22,density=dhwe2(x11,x12,x22))
  dist <- dist[order(dist$density),]
  data.frame(n11,n12,n22,
             freq=n1/(n1+n2),
             prob=cumsum(dist$density)[dist$n11==n11 & 
                                       dist$n12==n12 & 
                                       dist$n22==n22])
 }