*! version 1.0.0  30apr2008
program haplologit, eclass
	version 10.0
	if replay() {
		if "`e(cmd)'" != "haplologit" {
			error 301
		}
		if `e(rc)' != 0 {
			di as err "`e(cmd)' exited with error `e(rc)'; " ///
				  "replay is not allowed"
			error 301
		}
		local depvar `e(depvar)'
		Display `0'
		eret local depvar `"`depvar'"'
		exit
	}
	cap noi Estimate `0'
	local rc = _rc
	eret local cmdline `"haplologit `0'"'
	eret local cmd "haplologit"
	eret scalar rc = `rc'
	sret clear
	global MY_panel
	global MY_hap1
	global MY_hap2
	global MY_inhermodel
	global MY_riskhapmat
	global MY_rhapconsmat
	exit `rc'
end

program Estimate, eclass sortpreserve
	version 10.0
	syntax [varlist] [if] [in] ,			///
				SNPvars(varlist) [	/// Required
				noCONstant		/// Model opts
				INHERitance(string)	///
				RISKHAP(string asis) 	///
				RISKHAP1(string asis) 	///
				HFThreshold(string)	///
				CONSTraints(passthru) 	///
				COLlinear 		///
				Level(cilevel)		/// Reporting
				OR			///
				ALLDOTs 		///
				noCOEf noFREq noHEADer	///
				HAPPREFIX(string)	///
				EMITERate(passthru) 	/// EM opts
				EMTOLerance(passthru)	///
				EMSAMPle(string)	///
				eminit(string)		///
				EMLOG EMDOTs 		///
				noEMSHOW noEMTable	///
				DIFficult		/// ML opts
				ITERate(passthru)	///
				noLOG TRace		///	
				GRADient HESSian	///
				showstep		///
				TOLerance(passthru)	///
				LTOLerance(passthru)	///
				GTOLerance(passthru)	///
				NRTOLerance(passthru)	///
				NONRTOLerance		///
				SHOWNRtolerance		///
				MLMethod(string) 	/// undocumented opts
				RARE			///
				PHASEINFO(passthru) 	///
				NRISKHAPlos(real 0) 	///
				TECHnique(passthru)	///
				SORT			///
				BASE10			///
				*			/// riskhap#() opts
					]
	if `"`happrefix'"' == "" {
		local happrefix "hap_"
	}
	if `"`rare'"' == "" {
		local rare rare
	}
	local maxopts `difficult' `iterate' `log' `trace' `gradient' `hessian'
	local maxopts `maxopts' `showstep' `tolerance' `ltolerance' `gtolerance'
	local maxopts `maxopts' `nrtolerance' `nonrtolerance' `shownrtolerance'
	local maxopts `maxopts' `technique'
	if `"`hfthreshold'"' != "" {
		cap confirm number `hfthreshold'
		if _rc {
			di as err "hfthreshold() must be between 0 and 1"
			exit 198
		}
		if `hfthreshold' <= 0 | `hfthreshold'>=1 {
			di as err "hfthreshold() must be between 0 and 1"
			exit 198
		}
		if `hfthreshold' < 1e-15 {
			di as err "hfthreshold(`hfthreshold') is too small"
			exit 498
		}
	}
	if "`emshow'" != "" {
		local emshow qui
	}
	if "`emtable'" != "" {
		local emtable nohaptable
	}
	if "`alldots'" != "" {
		if "`emlog'" != "" {
			di as err "alldots and emlog may not be combined"
			exit 198
		}
		local emdots emdots
	}
	marksample touse
	// main equation
	local ycons `constant'
	local constant
	if "`ycons'" == "" {
		local c = 1
	}
	else {
		local c = 0
		local note nonote
	}
	gettoken depvar xvars: varlist
	tempvar y
	qui gen byte `y' = (`depvar'!=0) if `touse'
	qui count if `depvar'==0
	if r(N) == 0 {
		di as txt "outcome does not vary; remember:"
		di as txt "                          0 = negative outcome,"
		di as txt "all other nonmissing values = positive outcome"
		exit 2000
	}
	if `"`xvars'"' == "" & "`ycons'" != "" {
		di as err "noconstant is not allowed without " ///
			  "independent variables"
		exit 198
	}
	_rmcoll `xvars' if `touse', `collinear' `ycons'
	local xvars `r(varlist)'
	local nxvars: word count `xvars'
	
	// check SNP vars to have values 0, 1, 2, or missing
	local nsnp: word count `snpvars'
	tokenize `snpvars'
	tempvar snpchk
	qui gen byte `snpchk' = .
	forvalues i = 1/`nsnp' {
		qui replace `snpchk'=(``i''==0)+(``i''==1)+(``i''==2)+(``i''==.)
		qui summ `snpchk', meanonly
		if r(min) == 0 {
			di as err `"snpvars(): variable ``i'' must "' ///
				   "contain values 0, 1, 2, and . only"
			exit 198
		}
	}
	drop `snpchk'
	// nriskhaplos()
	local sign = -1
	if `nriskhaplos' < 0 {
		local sign = 1
		local nriskhaplos = abs(`nriskhaplos')
	}
	if `sign' < 0 {
		local dihaptype most
	}
	else {
		local dihaptype less
	}
	
	local isriskhap `"`riskhap'`riskhap1'"'
	if `"`riskhap'"' != "" {
		if `"`riskhap1'"' != "" {
			di as err "riskhap() and riskhap1() may not be combined"
			exit 198
		}
		local options riskhap1(`riskhap') `options'
	}
	else if `"`riskhap1'"' != "" {
		local options riskhap1(`riskhap1') `options'
	}
	// parse riskhap#() options
	tempname K
	scalar `K' = 2^`nsnp'
	local nriskhap = 0
	// local options should contain only riskhap#() options
	if `"`options'"' != "" {
		if `nriskhaplos' > 0 {
			di as err "nriskhap() may not be specified " ///
				  "with risk haplotypes"
			exit 198
		}
		local 0 , `options'
		tokenize `"`options'"'
		local i = 1
		while `"``i''"' != "" {
			syntax [, RISKHAP`i'(string asis) * ]
			if `"`riskhap`i''"' != "" {
				local ++nriskhap
				gettoken riskhap`i' riskhap`i'vars: ///
						riskhap`i', parse(",") quotes
				// parse suboptions: interaction() noconstant
				local opts `options'
				local 0 `riskhap`i'vars'
				syntax [, INTERaction(varlist)	///
					  noCONstant 		///
					   *			]
				if `"`interaction'"'=="" & "`constant'"!="" {
				       di as err "riskhap`i'(): noconstant " ///
					  "is not allowed without interaction()"
					exit 198
				}
				markout `touse' `interaction'
				local riskhap`i'vars `interaction'
				local riskhap`i'cons `constant'
				local options `opts'
				cap confirm integer number `riskhap`i''
				qui if _rc { // riskhap#() contains bin. seq. 
local 0 , riskhap`i'(`riskhap`i'') `options'
syntax [, RISKHAP`i'(string) * ]
local riskhap`i'str `riskhap`i''
mata: st_local("riskhap`i'",strofreal(frombase(2,st_local("riskhap`i'"))+1))
				}
				else {      // riskhap#() contains hap. index
cap assert `riskhap`i'' > 0
if _rc {
	di as err "riskhap`i'(): invalid haplotype index `riskhap`i''"
	exit 198
}
mata: st_local("riskhap`i'str",inbase(2,strtoreal(st_local("riskhap`i'"))-1))
				}
				// riskhap`i' stores ith ref. hap. index
				// riskhap`i'str stores ith riskhap as bin. seq.
				local len = strlen(`"`riskhap`i'str'"')
				if `len' > `nsnp' {
di as err "{p 0 12 4}riskhap`i'(): length of the risk haplotype "	///
   `""`riskhap`i'str'" may not exceed the number of variables in snpvars()"' ///
   " or haplotype index `riskhap`i'' may not exceed 2^N_snp = `=`K''{p_end}"
exit 198
				}
				if `"`riskhap`i'str'"' == "" {
di as err `"riskhap`i'() is incorrectly specified"'
exit 198
				}
				local riskhap`i'str ///
				  `: di _dup(`=`nsnp'-`len'') 0'`riskhap`i'str'
				if `riskhap`i'' == . & `"`isriskhap'"' != "" {
di as err `"{p 0 16 4}riskhap`i'(): risk haplotype "`riskhap`i'str'" "' ///
	  "is incorrectly specified{p_end}"
exit 198
				}
				local riskhaplist `riskhaplist' `riskhap`i''
				local duphap: list dups riskhaplist
				if `"`duphap'"' != "" {
di as err "riskhap`i'(): duplicate risk haplotype "	///
	  `""`riskhap`i'str'" (`riskhap`i'') is specified "'
exit 198
				}
			}
			local 0 , `options'
			local ++i
		} /* end of riskhap#() parsing */
		if `"`options'"' != "" {
			di as err `"option `options' is not allowed"'
			exit 198
		}
	}
	else { // the default is main effect of the most frequent haplotype
		if `nriskhaplos' == 0 {
			local nriskhaplos = 1
		}
		local nriskhap = 0
	}
	
	// parse inheritance()
	if `"`inheritance'"' == "" {
		local inheritance additive
	}
	else {
		local 0 , `inheritance'
		cap syntax [, Additive Dominant Recessive Codomimant]
		if _rc {
			di as err "inheritance(): `inheritance' is not allowed"
			exit 198
		}
		local inheritance `additive' `dominant' `recessive' `codominant'
		if `:word count `inheritance'' > 1 {
			di as err "inheritance(): only one model is allowed"
			exit 198
		}
		if `"`inheritance'"' == "codominant" {
			di as err "codominant model is not available"
			exit 198
		}
		
	}

	if `"`emsample'"' == "" {
		local emsample controls
		local emsampti control sample
		local ifopt ifopt(`y'==0)
	}
	else {
		local emsamplen = length(`"`emsample'"')
		if `emsamplen' < 2 {
			di as err "emsample() is incorrectly specified"
			exit 198
		}
		else {
			if `"`emsample'"'==substr("controls",1,`emsamplen') {
				local emsample controls
				local emsampti control sample
				local ifopt ifopt(`y'==0)
			}
			else if `"`emsample'"'==substr("cases",1,`emsamplen') {
				local emsample cases
				local emsampti case sample
				local ifopt ifopt(`y'==1)
			}
			else if `"`emsample'"'==substr("all",1,`emsamplen') {
				local emsample all
				local emsampti case-control sample
				local ifopt ifopt(`y'==0|`y'==1)
			}
			else {
				di as err "emsample() is incorrectly specified"
				exit 198
			}
		}
	}
	// estimate haplotype frequencies from the control sample
	// retain only haplotypes with sample freq. less than max(2/N, 0.001)
	qui summ if `touse', meanonly
	local N = r(N)
	if `"`hfthreshold'"' == "" {
		local cutoff = max(2/`N', 0.001)
		if `cutoff' >= 0.1 {
			local cutoff = 0.001
		}
		local hfthreshold = `cutoff'
	}
	else {
		local cutoff = `hfthreshold'
	}
	tempfile fhaplo
	tempname id indhap1 indhap2
	// expand snp-based genotype data to diplotype data (hap. pairs)
	_gntp2dptp `snpvars' if `touse', 				 ///	
					usevar(`id' `indhap1' `indhap2') ///
					saving(`"`fhaplo'"')		 ///
					base10				 ///
					`ifopt'			 	 ///
					nohead				 ///
					`alldots' `phaseinfo'
	// obtain initial haplotype frequencies from emsample() (_ifcond)
	di
	di as txt "Obtaining initial haplotype frequency estimates"	///
		  " from the `emsampti':"
	local emopts `emtolerance' `emiterate' `base10' `emlog' `emdots'
	local emopts `emopts' `emtable' `sort' eminit(`eminit')
	`emshow' _haplo_em_phased `id' `indhap1' `indhap2' _ifcond	///
						using `"`fhaplo'"',   	///
						nsnp(`nsnp') 	  	///
						cutoff(`cutoff')  	///
						onlycutoff	  	///
						nohead		  	///
						`emopts'
	tempname em_ll em_N
	scalar `em_ll'  = e(ll)
	scalar `em_N'   = e(N_s)
	// redefine K to be the number of retained haplotypes
	tempname freqmat
	mat `freqmat' = e(freq)
	scalar `K' = colsof(`freqmat')
	if `K' == 1 {
		di as err "only one haplotype retained"
		exit 498
	}
	else if `K' < 1 {
		di as err "no haplotypes retained"
		exit 498
	}
	if `nriskhap' >= `K' {
		di "{p 0 6 4}{txt}Note: too many risk haplotypes are " ///
		   "specified; using {res}" `=`K'-1' "{txt} haplotype(s){p_end}"
		local nriskhap = `K'-1
	}
	di
	di as txt "Performing gradient-based optimization:"
	//di as txt "Note: using rare-disease approximation"
	tempfile fesamp
	// sid is used to save the estimation sample
	tempvar sid
	qui gen long `sid' = _n
	preserve
	qui keep if `touse'
	qui gen long `id' = _n if `touse'
	sort `id'
	tempvar merge
	// long form: consistent haplotype pairs are nested within subjects
	qui merge `id' using `"`fhaplo'"', uniqmaster _merge(`merge')
	sort `id'
	// macro ind contains retained haplotype indices
	mata: ///
st_local("ind",invtokens(substr(st_matrixcolstripe("e(freq)")[,2]',5,.)))
	local nlevels: word count `ind'
	// remove obs not compatible with retained haplotype pairs
	tempvar cnt subjdrop
	qui egen byte `cnt' = anycount(`indhap1' `indhap2'), values(`ind')
	// count how many subjects are removed
	qui by `id': gen long `subjdrop' = sum(`cnt'!=2)
	qui by `id': replace `subjdrop' = (`subjdrop'[_N]==_N & _n==_N)
	qui count if `subjdrop' & `touse'
	if r(N)!=0 {
		di as txt "{p 0 6 4}Note: removing " ///
		   as res `=string(r(N),"%9.0g")' ///
		   as txt " observation(s); constituent haplotype "	///
			  "frequencies are smaller than " ///
		   as res `"`=string(`hfthreshold',"%9.0g")'{p_end}"'
	}
	qui replace `touse' = 0 if `cnt'!=2
	drop `merge' `cnt' `subjdrop'
	qui keep if `touse'
	// save esample in a file
	qui outfile `sid' using `"`fesamp'"'
	// ifreqsort: matrix of hap. indices sorted from lowest to highest freq.
	tempname ifreqsort
	qui mata: ///
	st_matrix(st_local("ifreqsort"), sort((st_matrix("e(freq)")', st_matrix("e(index)")'),`sign')[.,2]')
	if `nriskhaplos' > 0 {
		if `nriskhaplos' >= `nlevels' {
			local nriskhaplos = `nlevels' - 1
			di as txt "{p 0 6 4}Note: nriskhap() must be less " ///
				  "than the number of retained haplotypes " ///
			   as res `nlevels' as txt "; using "		    ///
			   as res `nriskhaplos' as txt " `dihaptype' "	    ///
			   "frequent haplotypes from the `emsampti'{p_end}"
		}
		else if `nriskhaplos' == 1 {
			if `sign' < 0 {
di as txt "{p 0 6 4}Note: using the most frequent haplotype " ///
	  "from the `emsampti' as a risk haplotype{p_end}"
			}
			else {
di as txt "{p 0 6 4}Note: using the least frequent haplotype " ///
	  "from the `emsampti'{p_end}"
			}
		}
		else {
			di as txt "{p 0 6 4}Note: using "		///
			   as res `nriskhaplos' as txt " `dihaptype' "	///
			   "frequent haplotypes from the `emsampti'{p_end}"
		}
	}
	local concat X
	// map retained indices and risk haplotype indices
	// to a sequence of ordered integers
	// ind1, ind2 contain new indices
	tempvar ind1 ind2
	qui gen double `ind1' = `indhap1'
	qui gen double `ind2' = `indhap2'
	tokenize `ind'
	local nretained = 0
	local Nrhapvars = 0
	tempname riskhapmat rhapconsmat
	forvalues i=1/`nlevels' {
		if ``i'' != `i' {
			qui replace `ind1' = `i' if `ind1'==``i''
			qui replace `ind2' = `i' if `ind2'==``i''
		}
		if `nriskhaplos' > 0 { // if nriskhap() specified
			forvalues j=1/`nriskhaplos' {
if `=`ifreqsort'[1,`j']' == ``i'' {
	mat `riskhapmat' = (nullmat(`riskhapmat'), `i')
	local ++nretained
	// build haplotype ml and test equations
	mata: st_local("riskhap`nretained'str",inbase(2,strtoreal("``i''")-1))
	local len = `nsnp' - strlen(`"`riskhap`nretained'str'"')
	local riskhap`nretained'str `: di _dup(`len') 0'`riskhap`nretained'str'
	local eqriskhap `eqriskhap' (hap_`riskhap`nretained'str'`concat':)
	local testriskhap `testriskhap' ([hap_`riskhap`nretained'str'`concat']: _cons)
	mat `rhapconsmat' = (nullmat(`rhapconsmat'), 1)
	local ++Nrhapvars
}
			}
		}
		else { // if opts riskhap#() specified
			forvalues j=1/`nriskhap' {
if "`riskhap`j''" == "" {
	di as err "option riskhap`j'() is required with `nriskhap' " ///
		  "risk haplotype(s)"
	exit 198
}
if `riskhap`j'' == ``i'' {
	local riskhap`j' = `i'
	mat `riskhapmat' = (nullmat(`riskhapmat'), `riskhap`j'')
	local ++nretained
	// build haplotype ml and test equations
	local eqriskhap `eqriskhap' (hap_`riskhap`j'str'`concat': `riskhap`j'vars', ///
							`riskhap`j'cons')
	if "`riskhap`j'cons'" == "" {
		local testriskhap `testriskhap' ///
				([hap_`riskhap`j'str'`concat']: `riskhap`j'vars' _cons)
	}
	else {
		local testriskhap `testriskhap' ///
				([hap_`riskhap`j'str'`concat']: `riskhap`j'vars')
	}
	local riskhapxvars `riskhapxvars' `riskhap`j'vars'
	local Nrhapvars = `Nrhapvars' + `:word count `riskhap`j'vars''
	if "`riskhap`j'cons'" == "" {
		local ++Nrhapvars
		mat `rhapconsmat' = (nullmat(`rhapconsmat'), 1)
	}
	else {
		mat `rhapconsmat' = (nullmat(`rhapconsmat'), 0)
	}
				}
			}
		}
		// build equation names for frequency parameters and 
		// nlcom expression
		if `i' != `nlevels' {
			local eqnu `eqnu' /nu`i'
			local nl_denom `nl_denom'`plus'exp([nu`i']_b[_cons])
			local plus +
		}
	} 
	if `nretained' == 0 {
		di as txt  "{p 0 6 4}Note: " as res `nretained'	   ///
		   as txt  " out of " as res `nriskhap'			   ///
		   as txt  " specified risk haplotypes are retained"  /// 
			   " in the computation; using the most frequent"  ///
		           " haplotype from the `emsampti'{p_end}"
		local riskhap1 = `:list posof "`=`ifreqsort'[1,1]'" in ind'
		mat `riskhapmat' = `riskhap1'
		mata: st_local("riskhap1str",inbase(2,`=`ifreqsort'[1,1]-1'))
		local len = `nsnp' - strlen(`"`riskhap1str'"')
		local riskhap1str `: di _dup(`len') 0'`riskhap1str'
		local nretained = 1
		local eqriskhap (hap_`riskhap1str'`concat': )
		local testriskhap ([hap_`riskhap1str'`concat']: _cons)
		mat `rhapconsmat' = 1
		local Nrhapvars = 1 
	}
	else if `nretained'<`nriskhap' {
		di as txt  "{p 0 6 4}Note: " as res `nretained'	///
		   as txt  " out of " as res `nriskhap'			///
		   as txt  " specified risk haplotypes" 		///
			   " are used in the computation{p_end}"
	}

	if `"`rare'"' != "" { // algorithm for rare disease
		global MY_panel `id'
		global MY_hap1 `ind1'
		global MY_hap2 `ind2'
		global MY_inhermodel `inheritance'
		global MY_riskhapmat `riskhapmat'
		global MY_rhapconsmat `rhapconsmat'
		tempname mat initmat basefq hapfreq
		mat `hapfreq' = e(freq)
		// specify base haplotype frequency
		scalar `basefq' = `hapfreq'[1, `nlevels']
		if `basefq' < 1e-8 {
			di as err "base haplotype frequency is too small"
			exit 498
		}
		qui mata: ///
st_matrix(st_local("mat"),ln(st_matrix("e(freq)")/st_numscalar(st_local("basefq"))))
		sort `id' `ind1' `ind2'
		if `"`mlmethod'"' == "" {
			local mlmethod d1
		}
		local dtype = substr(`"`mlmethod'"',1,2)
		local noprofile = (`"`xvars'"'=="" & `"`riskhapxvars'"' == "")
					
		if `noprofile' { 
// no covariates in the model
local note nonote
mat `initmat' = J(1,`Nrhapvars',0), `mat'[1,1..`=`nlevels'-1']
global MY_depvar `depvar'
ml model `mlmethod' _haplologit_d1_rare 		///
		`eqriskhap'				///
		`eqnu'					///
		if `touse'				///
		, init(`initmat', copy) maximize	///
		  search(off) `constraints' `collinear' ///
		  crittype(Retrospective log likelihood) ///
		  `maxopts'
		}
		else {
mat `initmat' = ///
	J(1,`=`nxvars'+`c'',0), J(1,`Nrhapvars',0), `mat'[1,1..`=`nlevels'-1']
ml model `mlmethod' _haplologit_`dtype'_rarex 		///
		(`depvar': `y' = `xvars', `ycons')	///
		`eqriskhap'				///
		`eqnu'					///
		if `touse'				///
		, init(`initmat', copy) maximize	///
		  search(off) `constraints' `collinear' ///
		  crittype(Retrosp. profile log likelihood) ///
		  `maxopts'
		}
		// combine haplotype equations with depvar equation
		// and order the estimation results
		tempname b V
		mata: _est_reorder()
		eret repost b = `b' V=`V', rename
		qui test [`depvar']
		local df_m = r(df)
		local chi2 = r(chi2)
		local chi2_p = r(p)
		eret scalar k_eq_nu = `nlevels'-1
		eret scalar k_eq = e(k_eq_nu) + 1
		eret local snpvars `"`snpvars'"'
		eret local inheritance `"`inheritance'"'
		eret local depvar `"`depvar'"'
		eret local genepop "Hardy-Weinberg equilibrium"
		eret local title "Haplotype-effects logistic regression"
		eret local emsample `emsample'
		eret local happrefix `happrefix'
		// clear predict
		eret local predict ""
		eret matrix em_freq = `freqmat'
		qui infile `sid' using `"`fesamp'"', clear
		qui by `sid', sort: keep if _n==1
		qui save `"`fesamp'"', replace
		restore
		tempvar merge
		qui merge `sid' using `"`fesamp'"', _merge(`merge') sort
		qui replace `touse' = 0 if `merge'!=3
		// compute number of phased, unphased, and missing genotypes
		tempvar miss one
		qui egen byte `miss' = rowmiss(`snpvars')
		qui egen byte `one' = anycount(`snpvars'), values(1)
		qui count if `miss'!=0 & `touse'
		local N_miss = r(N)
		qui count if `one'<=1 & `miss'==0 & `touse'
		local N_phased = r(N)
		qui count if `one'>1 & `miss'==0 & `touse'
		local N_unphased = r(N)
		qui count if `touse'
		local N_s = r(N)
		ereturn repost, esample(`touse')
		eret scalar N		= `N_s'
		eret scalar N_phased	= `N_phased'
		eret scalar N_unphased	= `N_unphased'
		eret scalar N_miss	= `N_miss'
		eret scalar ll		= `e(ll)'
		eret scalar converged	= `e(converged)'
		eret scalar df_m	= `df_m'
		eret scalar chi2 	= `chi2'
		eret scalar p 		= `chi2_p'
	
		// display results
		Display, ind(`ind') level(`level') `header' `coef' `freq' ///
			`or' `note'
		eret local depvar `"`depvar'"'
		// e-scalars (specific order)
		eret scalar N		= `N_s'
		eret scalar N_phased	= `N_phased'
		eret scalar N_unphased	= `N_unphased'
		eret scalar N_miss	= `N_miss'
		eret scalar ll		= `e(ll)'
		eret scalar converged	= `e(converged)'
		eret scalar df_m	= `df_m'
		eret scalar chi2 	= `chi2'
		eret scalar p 		= `chi2_p'
		eret scalar em_N = `em_N'
		eret scalar em_ll = `em_ll'
		eret scalar cutoff = `cutoff'
	}
	else {
		di as err "algorithm without rare disease assumption " ///
			  "is not available; specify option rare"
		exit
	}
end

program Display
	syntax [, ind(numlist) noHEADer noCOEf noFREq Level(cilevel) OR ///
		  noNote ]
	if "`header'" == "" {
		DiHeader
	}
	if "`coef'" == "" {
		local neq = e(k_eq) - e(k_eq_nu)
		if `neq' == 1 {
			local eq first
		}
		else {
			local eq neq(`neq')
		}
		di
		cap noi ml display, `eq' noheader level(`level') `or' nofoot
		if "`note'" == "" & "`or'" == "" {
di as txt "Note: _cons = b0 + ln(N1/N0) - ln{Pr(D=1)/Pr(D=0)}"
		}
	}
	if "`freq'" == "" {
		di
		ClearDepvar
		DiFreqTable, ind(`ind') level(`level')
	}
	ml footnote
end

program ClearDepvar, eclass
	eret local depvar
end

program DiHeader
	local crittype `e(crittype)'
	local piece2 `e(snpvars)'
	local piece1: piece 1 27 of `"`piece2'"', nobreak
	local piece2 = abbrev(subinstr(`"`piece2'"', `"`piece1'"', "", 1), 27)
	if missing(e(chi2)) {
		local dichi2 ///
		  "{help j_robustsingular##|_new:`e(chi2type)' chi2(`e(df_m)')}"
	}
	else {
		local dichi2 "{txt}`e(chi2type)' chi2({res}`e(df_m)'{txt})"
	}
	di
	di 	   `"{txt}`e(title)'"' 
	di 	    "{txt}Mode of inheritance: {res}" `"`e(inheritance)'"' _c
	di _col(49) "{txt}Number of obs      = {res}" %9.0g e(N)
	di
	di	    "{txt}Genetic distribution: {res}" ///
			"`=subinstr("`e(genepop)'","rium",".",1)'" _c
	di _col(49) "{txt}Number phased      = {res}" %9.0g e(N_phased)
	di	    "{txt}Genotype: {res}" `"`piece1'"' _c
	di _col(49) "{txt}Number unphased    = {res}" %9.1g e(N_unphased)
	di	         "          {res}" `"`piece2'"' _c
	di _col(49) "{txt}Number missing     = {res}" %9.0g e(N_miss)
	di
	di _col(49) "`dichi2'" _col(68) "{txt}=" _col(70) as res %9.2f e(chi2)
        di          "{txt}`crittype' = {res}" %10.0g e(ll) _c
	di _col(49) "{txt}Prob > chi2        = {res}" ///
						%9.4f chiprob(e(df_m), e(chi2))
end

program DiFreqTable
	syntax [, ind(numlist) Level(cilevel) ]
	local happrefix `e(happrefix)'
	if `"`ind'"' == "" {
		mata: ///
st_local("ind",invtokens(substr(st_matrixcolstripe("e(em_freq)")[,2]',5,.)))
	}
	local nsnp: word count `e(snpvars)'
	mata: ///
st_local("indstr",invtokens(_inbase2fill(`nsnp',strtoreal(tokens(st_local("ind"))):-1)))
	// build nl expression
	local nlevels: word count `ind'
	local nlevels1 = `nlevels' - 1
	tokenize `indstr'
	forvalues i=1/`nlevels1' {
		local nldenom `nldenom'`plus'exp([nu`i']_b[_cons])
		local plus +
	}
	forvalues i=1/`nlevels1' {
		local nl_expr `nl_expr' ///
		      (`"`happrefix'``i''"':exp([nu`i']_b[_cons])/(1+`nldenom'))
	}
	local nl_expr `nl_expr' (`"`happrefix'``nlevels''"': 1/(1+`nldenom'))
	di as txt "Haplotype frequencies"
	nlcom `nl_expr', nohead level(`level')
end

version 10

mata:

/* 
modifies haplotype equation names
	form			to
	hap_"str"X:_cons	status:-prefix-_"str"
	hap_"str"X:xvar		status:-prefix-_"str"*xvar

and puts estimation results in the following order
	xvars
	haplotype constants
	haplotype interactions
	constant
*/
void _est_reorder()
{
	string matrix colname_hap, colname_xb, colname_nu, colstripe
	real colvector sel, sel_dv, sel_hap, sel_cons, p, i
	real matrix b, V, w
	string scalar depvar
	real scalar nrows, nparms, indcons

	b 		= st_matrix("e(b)")
	V 		= st_matrix("e(V)")
	colstripe	= st_matrixcolstripe("e(b)")
	nrows 		= rows(colstripe)
	depvar		= st_local("depvar")
	
	sel_dv		= (substr(colstripe[.,1],1,strlen(depvar)):==depvar)
	sel_hap		= (substr(colstripe[.,1],1,4):=="hap_")
	sel_cons	= strmatch(colstripe[.,2], "_cons")
	sel		= (substr(colstripe[.,1],1,2):=="nu")
	// change stripe names	
	colname_xb = select(colstripe, sel_dv)
	colname_hap = select(colstripe, sel_hap)
	colname_nu = select(colstripe, sel)
	colname_hap[.,2] = colname_hap[.,1]:+colname_hap[.,2]
	colname_hap[.,2] = subinstr(colname_hap[.,2], 
				    st_local("concat")+"_cons", "", .)
	colname_hap[.,2] = subinstr(colname_hap[.,2], 
				    st_local("concat"), "*", .)
	if (st_local("happrefix")!="hap_") {
		colname_hap[.,2] = subinstr(colname_hap[.,2], 
				    	    "hap_", st_local("happrefix"), .)
	}
	colname_hap[.,1] = J(rows(colname_hap), 1, depvar)
	colstripe = (colname_xb\colname_hap\colname_nu)

	nparms = colsum(sel_dv)
	if (nparms != 0 ) { // change the order of estimation results
		if (nparms == 1) {
			if (st_local("xvars") == "") {
				indcons = 1
				p = J(0,0,.)
			}
			else {
				indcons = 0
				p = 1
			}
		}
		else {
			if (st_local("ycons") == "") {
				p = range(1, nparms-1, 1)
				indcons = nparms
			}
			else {
				p = range(1, nparms, 1)
				indcons = 0
			}
		}
		nparms = colsum(sel_hap)
		if (nparms != 0) { 
// should always be true since always have haplo effects
			sel = sel_hap:*sel_cons
			if (colsum(sel) != 0) { // haplotype constants
				minindex(sel, -1, i, w)
				if (cols(p) == 0 ) p = i
				else p = p \ i
			}
			sel = sel_hap:*(1:-sel_cons)
			if (colsum(sel) != 0) { // haplotype interactions
				minindex(sel, -1, i, w)
				if (cols(p) == 0 ) p = i
				else p = p \ i
			}
		}
		if (indcons != 0 ) p = p \ indcons
		nparms = colsum(sel_nu)
		if (nparms != 0) { // should always be true
			p = p \ range(rows(p)+1, nrows, 1)
		}
		// reorder rows and columns
		colstripe = colstripe[p,.]
		b = b[.,p']
		V = V[.,p']
		V = V[p, .]
	}
	st_matrix(st_local("b"), b)
	st_matrix(st_local("V"), V)
	st_matrixcolstripe(st_local("b"), colstripe)
	st_matrixcolstripe(st_local("V"), colstripe)
	st_matrixrowstripe(st_local("V"), colstripe)
}

end