*! version 5.0.1  20jan1999     STB-48 sg110
program def genhw, rclass
version	 6
	syntax varlist(min=2 max=2) [if] [in] [fweight] /*
		*/ [, Label(string) Binvar Exact REPS(integer 1000) ]
	marksample touse, strok
qui {
		tokenize `varlist'
		local type1: type `1'
		local type2: type `2'
		tempvar a1 a2 both grp  al1 al2 
		if substr("`type1'",1,3)!= "str" {
		        gen str8 `al1'=string(`1')
		}
		else gen str8 `al1'= `1'
		if substr("`type2'",1,3)!= "str" {
		        gen str8 `al2'= string(`2')
		}
		else gen str8 `al2'= `2'
		
		gen str8 `a1'= `al1' if `al1'<=`al2'
		gen str8 `a2'= `al2' if `al1'<=`al2'
		replace  `a1'= `al2' if `a1'==""
		replace  `a2'= `al1' if `a2'==""
	
		tempvar cnt
		if "`weight'" != "" { 
			gen double `cnt' `exp'
			compress `cnt'
		}
		else	gen byte `cnt' = 1
	
		preserve
		
		drop `1' `2'
		rename `a1' `1'
		rename `a2' `2'
		drop if `cnt'==0
		expand `cnt'
		keep if `touse'
		stack `1' `2' , into(`2') clear
		tempname allele palle
		tab `2', matcell(`allele')
		local numall=r(r)           /* number of different alleles */
		if `numall'==2 {
			restore, preserve
			drop if `cnt'==0
			expand `cnt'
			keep if `touse'
			if `"`binvar'"'=="" {
				local binvar="0"
			}
			if `"`label'"'=="" {
				local label=""
			} 
			noi HW2 `1' `2' `binvar' `"`label'"'
			ret add
			exit
		}
			
		mat `palle' = (1/r(N))*`allele'
		keep `2' 
		sort `2'
		qui by `2': keep if _n==1 
		local allname=`2'[1]
		local i 2
		while `i'<=_N {
			local allname="`allname'"+" "+ `2'[`i']
			local i=`i'+1
		}
		tempfile tt
		save `"`tt'"'
		rename `2' `1'
		Mcross `"`tt'"'
		sort `1' `2'
		save `"`tt'"', replace
		
		restore, pres
		keep if `touse'
		drop `1' `2'
		rename `a1' `1'
		rename `a2' `2'
		keep `1' `2' `cnt'
	
		sort `1' `2'
		merge `1' `2' using `"`tt'"'
		drop _merge
		replace `cnt'=0 if `cnt'==.
		/* ROUTINE FOR DISPLAYING CROSS TABULATION" */
		/* 
			if "`notable'"~="notable" {
		 	noi table `1' `2',c(sum `cnt') row col center f(%5.0f)
			}
		*/
		sort `1' `2'
		collapse (sum) `cnt', by(`1' `2')
		
		tempvar a1 a2 
		/* order alleles */
		gen str12 `a1'=`1' if `1'<=`2'
		gen str12 `a2'=`2' if `1'<=`2'
		replace `a1'=`2' if `2'<=`1'
		replace `a2'=`1' if `2'<=`1'
		replace `1'=`a1'
		replace `2'=`a2'
		drop `a1' `a2'
		sort `1' `2'
		collapse (sum) `cnt', by(`1' `2')
		sum `cnt' 
		local obs=r(sum)
		/* genotype counts(geno), p(geno)=pgeno, falle=freq. allele,
		   palle=prob. allele*/
		gen double pgeno=`cnt'/`obs'
		gen double fall1=0
		gen double fall2=0
		gen double pall1=0
		gen double pall2=0
		tempname all1 all2
		encode `1', gen(`all1')
		encode `2', gen(`all2')
		replace fall1=`allele'[`all1',1]
		replace fall2=`allele'[`all2',1]
		replace pall1=`palle'[`all1',1]
		replace pall2=`palle'[`all2',1]
		gen int n=`obs'
		replace n=2*n if `1'~=`2'
		gen double exp=n*pall1*pall2
		drop n
		replace `1'=trim(`1')
		replace `2'=trim(`2')
		gen str16 geno=trim(`1'+`2')
		gen double D=pall1*pall2 - (0.5*pgeno) if `1'~=`2'
		/* noi di _n in gr "     Genotype counts" */
		noi di in gr /*
		*/ "                                          Disequilibrium"
		noi di in gr /*
		*/ "        Genotype | Observed   Expected   Coefficient (D)"
		noi di in gr /*
		*/ "     ------------+--------------------------------------"
		tempname hom
		mat `hom'=J(`numall',1,0)
		local i 1
		local j 1
		while `i'<=_N {
			noi di in gr %16s /*
			*/ geno[`i'] " |" in ye _col(20) %8.0f `cnt'[`i'] /*
			*/   _col(31) %8.2f exp[`i'] _c
			if `1'[`i']~=`2'[`i'] {
				noi di in ye _col(11) %8.4f D[`i']
			}
			else noi di ""
			if `1'[`i']==`2'[`i'] {
				mat `hom'[`j',1]=pgeno[`i']		
				local j = `j' + 1
	
			}
			local i = `i' + 1
		}
		noi di in gr /*
		*/ "     ------------+--------------------------------------"
		gen tg=sum(`cnt')
		gen te=sum(exp)
		noi di in gr _col(12) "Total" " |" /*
		*/  in ye _col(20) %8.0f tg[_N]  _col(31) %8.2f te[_N] 
		tempname x2 llo lle lr
		gen double sx=sum( ((`cnt'-exp)^2)/(exp))	
		scalar `x2'=sx[_N]	             /* chi2 test */
		gen double llo= sum(`cnt'*ln(`cnt'/`obs')) 
		gen double lle= sum(`cnt'*ln(exp/te[_N]))
		scalar `lr'=2*(llo[_N]-lle[_N])
		local df=`numall'*(`numall'-1)/2
		drop sx llo lle
		noi di ""
		noi di in gr /*
		*/ "          Allele | Observed     Frequency      Std. Err."
		noi di in gr /*
		*/ "     ------------+--------------------------------------"
		tokenize `allname'
		local i 1
		while `i'<=`numall' {
			noi di in gr %16s /*
			*/ "``i''" " |" in ye _col(20) %8.0f `allele'[`i',1] /*
			*/  _col(34) %8.4f `palle'[`i',1] _c
			if "`binvar'"!="" {
	   			local vartype " (binomial)"
				noi di in ye _col(8) %8.4f /*
		 	 	*/ sqrt(`palle'[`i',1]* /*
				*/ (1-`palle'[`i',1])/(2*`obs')) /*
				*/  in gr "`vartype'"
			}
			else {
	   			local vartype 
				noi di in ye _col(8) %8.4f /*
				*/ sqrt((1/(2*`obs'))*(`palle'[`i',1]+ /*
				*/ `hom'[`i',1]-2* /*
				*/ (`palle'[`i',1]*`palle'[`i',1])))
			}
			local i=`i'+1
		}
		noi di in gr /*
		*/ "     ------------+--------------------------------------"

		noi di _n in gr "     Hardy-Weinberg Equilibrium Test: "
		noi di in gr /*
		*/ "              Pearson chi2 (" in ye `df' in gr ") = " /*
		*/ in ye %8.3f `x2' /*
		*/  in gr "  Pr= " in ye %5.4f  chiprob(`df',`x2')
		noi di in gr /*
		*/ "     likelihood-ratio chi2 (" in ye `df' in gr ") = " /*
		*/ in ye %8.3f `lr' /*
		*/  in gr "  Pr= " in ye %5.4f  chiprob(`df',`lr')
        	ret scalar chi2_p = chiprob(`df',`x2')
        	ret scalar df = `df'
		ret scalar chi2 = `x2'
        	ret scalar lr_p = chiprob(`df',`lr')
		ret scalar lr_chi2 = `lr'
		if "`exact'"~="" {
			drop if `cnt'==0
			keep `all1' `all2' `cnt'
			expand `cnt'
			keep `all1' `all2' 
			tempvar exapv
			EXAct `all1' `all2' `reps'  `exapv'
			noi di in gr  /*
			*/ "      Exact significance prob  =               " /*
			*/ in ye %5.4g `exapv'
		}

	}
end
program define Mcross          /* using `tt' */
        args using
        local nob = _N
        tempfile cross2
        tempvar order midx
        preserve
        quietly use `"`using'"', clear
        quietly {
                gen long `order'=_n
                expand `nob', clear
                sort `order'
                by `order': gen long `midx' = _n
                sort `midx' `order'
                drop `order'
                save `"`cross2'"', replace
                restore, preserve
                gen long `midx' = _n
                sort `midx'
                merge `midx' using "`cross2'"
                drop `midx' _merge
                restore, not
        }
end

program def HW2, rclass
qui {
	local type1: type `1'
	local type2: type `2'
	local binvar="`3'"
	local label="`4'"
	if "`binvar'"=="0" { local binvar }
	tempvar a1 a2 both grp  al1 al2
	if substr("`type1'",1,3)!= "str" {
		gen str8 `al1'=string(`1')
	}
	else gen str8 `al1'= `1'
	if substr("`type2'",1,3)!= "str" { 
		gen str8 `al2'= string(`2')
	}
	else gen str8 `al2'= `2'

		gen str8 `a1'= `al1' if `al1'<=`al2'
		gen str8 `a2'= `al2' if `al1'<=`al2'
		replace  `a1'= `al2' if `a1'==""
		replace  `a2'= `al1' if `a2'==""
		replace `a1'=trim(`a1')
		replace `a2'=trim(`a2')
		gen str8 `both'=trim(`a1'+`a2')
		tab `both'
		local ngrp=_result(2)
		sort `both'
		by `both': replace `both'="zzzzzzz" if _n~=1
		egen `grp'=group(`a1' `a2')
		local vars=" "
		local i=1
		while `i'<=`ngrp' {
			count if `grp'==`i'
			local na`i'=_result(1)
			local vars="`vars'"+"  "+"`na`i''"
			local i=`i'+1
		}
		sort `both'
		local lvars=" "
		local i=1
		while `i'<=`ngrp' {
			local la`i'=`both'[`i']
			local lvars="`lvars'"+"  "+"`la`i''"
			local i=`i'+1
		}
		if "`label'"~="" {
			local labopt="label(`label')"
		}
		else local labopt="label(`lvars')"
        }
	genhwi `vars', `labopt' `binvar'
	ret add
end

program def EXAct
	gen long a1=`1'
	gen long a2=`2'
	local reps=`3'
	local sub = _N
	preserve
	/* begin calculate observed prob */
	tempname obsprob 
	PRob `obsprob' 
	/* end calculate observed prob */
	restore, preserve
	stack a1 a2, into(a1) clear
	drop _stack
	tempfile stack
	qui save `stack', replace
	tempname  rand allele k1 k2
	tempfile myfile
	local count = 0
	local i 1
	while `i'<=`reps' {
		qui postfile `allele' a1 a2 using `myfile', double replace
			gen `rand' = uniform()
			sort `rand'
		local k 1
		local j 1
		while `j' <= `sub' {
			scalar `k1' = a1[`k']
			scalar `k2' = a1[`k'+1]
			post `allele' `k1' `k2'
			local k=`k'+2
			local j = `j' + 1
		} 	
		postclose `allele'
		use `myfile', clear
		/* begin calculate observed prob */
		tempname prob 
		PRob `prob' 
		/* end calculate observed prob */
		if `prob'<=`obsprob' {
			local count = `count'+1
		}

		use `stack', clear	
		local i = `i' +1
	}			
*noi di "count= " `count'
		scalar `4'= `count'/`reps'
end	
prog def PRob
	tempvar na1 na2 cnt num den
	tempname H
	qui gen int het=1 if a1~=a2
	qui count if het==1
	scalar `H' = (r(N)) 
	qui gen `na1'=a1 if a1<=a2
	qui replace `na1'=a2 if a2<a1
	qui gen `na2'=a2 if a2>=a1   
	qui replace `na2'=a1 if a1>a2
	contract `na1' `na2', freq(`cnt')
	qui gen double `den' = lngamma(`cnt'+1)
	qui replace `den'= sum(`den')
	scalar `1' = (2^`H')/exp(`den'[_N])
end