*! emh.ado Version 1.0 JRC 2010-10-12
program define emh, rclass byable(onecall)
	version 11.1
	syntax varlist(min=2 max=2 numeric) [if] [in] [fweight], [Strata(varlist) Correlation ///
		Anova General Transformation(string)]

	// Validate input
	local model `correlation' `anova' `general'
	if (1 < `: word count `model'') {
		display in smcl as error "Only one of " as input "`model'" as error " may be specified"
		error 178
	}

	local transformation = cond("`transformation'" == "", "table", "`transformation'")
	if !inlist("`transformation'", "table", "rank", "ridit", "modridit") &  ///
		!inlist("`transformation'", "logrank", "savage", "median", "mood", "klotz", "vdw") {
		display in smcl as error "`transformation' scores not recognized"
		error 178
	}

	marksample touse
	markout `touse' `strata', strok

	// Generate column of ones for -summarize-
	tempname ones
	generate byte `ones' = 1

	// Generate unique identifier to strata variable combinations
	tempvar strata_var
	local strata_tally : word count `strata'
	if (`strata_tally' == 1) {
		local strata_type : type `strata'
		if (substr("`strata_type'", 1, 3) == "str") {
			createUniqueStrataVariable `strata_var' if `touse', strata(`strata') one string
		}
		else if inlist("`strata_var'", "float", "double") {
			createUniqueStrataVariable `strata_var' if `touse', strata(`strata')
		}
		else {
			createUniqueStrataVariable `strata_var' if `touse', strata(`strata') one
		}
	}
	else if (1 < `strata_tally') {
		createUniqueStrataVariable `strata_var' if `touse', strata(`strata')
	}
	else {
		quietly generate byte `strata_var' = 1 if `touse'
	}

	// Process weight option
	local weight = cond("`weight'" == "", "", "[`weight' `exp']")

	// Set up model
	tempname G V R C
	if inlist("`model'", "correlation", "") { // Default is Correlation
		local model correlation
		quietly correlate `varlist' if `touse'
		local df = !missing(r(rho))
	}
	else {
		// ANOVA and General Association need the list of levels of row variable (union set over all strata)
		gettoken Y X : varlist
		quietly tabulate `X' if `touse', matrow(`R')
		local df = r(r) - 1
		// General Association needs also the list of levels of column variable (union set over all strata)
		if ("`model'" == "general") {
			quietly tabulate `Y' if `touse', matrow(`C')
			local df = `df' * (r(r) - 1)
		}
	}

	if (`df' == 0) {
		display in smcl as error "Insufficient variation in data"
		error 178
	}

	matrix define `G' = J(`df', 1, 0)
	matrix define `V' = J(`df', `df', 0)

	// Process strata, cumulating G and V
	tempname Gh Vh
	summarize `strata_var', meanonly
	forvalues stratum = 1/`r(max)' {
		`model' `varlist' if `touse' & `strata_var' == `stratum',  ///
			weight("`weight'") transformation(`transformation') ///
			ones(`ones') r(`R') c(`C') g(`Gh') v(`Vh')
		matrix define `G' = `G' + `Gh'
		matrix define `V' = `V' + `Vh'
	}

	// Invert V
	tempname invV
	matrix define `invV' = invsym(`V')

	// Compute Q and p
	tempname Q
	if 0 < diag0cnt(`invV') {
		local singularity_warning (Covariance matrix is singular)
		return scalar chi2 = .
	}
	else {
		matrix define `Q' = `G'' * `invV' * `G'
		return scalar chi2 = `Q'[1, 1]
	}

	return scalar df = `df'
	return scalar p = chi2tail(`df', return(chi2))

	// Get formal name of transformation score type
	tempname ScoreType
	if "`transformation'" == "table" {
		scalar define `ScoreType' = "Table Scores (Untransformed Data)"
	}
	else if "`transformation'" == "rank" {
		scalar define `ScoreType' = "Ranks"
	}
	else if "`transformation'" == "ridit" {
		scalar define `ScoreType' = "Ridits"
	}
	else if "`transformation'" == "modridit" {
		scalar define `ScoreType' = "Standardized Midranks (Modified Ridits)"
	}
	else if inlist("`transformation'", "logrank", "savage") {
		scalar define `ScoreType' = "Savage (log-rank) Scores"
	}
	else if "`transformation'" == "median" {
		scalar define `ScoreType' = "Median Scores"
	}
	else if "`transformation'" == "mood" {
		scalar define `ScoreType' = "Mood Scores"
	}
	else if "`transformation'" == "klotz" {
		scalar define `ScoreType' = "Klotz Scores"
	}
	else if "`transformation'" == "vdw" {
		scalar define `ScoreType' = "van der Waerden Scores"
	}
	else {
		scalar define `ScoreType' = ""
	}
	return local scoretype = "`transformation'"
	return local ScoreType = `ScoreType'

	// Display results
	display in smcl as text _newline(1) ///
		"Extended Mantel-Haenszel (Cochran-Mantel-Haenszel) Stratified Test of Association"
	if ("`model'" == "general") {
		local statistic General Association
	}
	else if ("`model'" == "anova") {
		local statistic ANOVA (Row Mean Scores)
	}
	else {
		local statistic Correlation
	}
	display in smcl as text _newline(1) "`statistic' Statistic:"
	display in smcl as text "Q (" as result `df' as text ") = " as result %7.4f return(chi2) _continue
	display in smcl as text ", P = " as result %06.4f return(p) as text " `singularity_warning'"
	display in smcl as text "Transformation: " `ScoreType' _newline(1)
end

program define createUniqueStrataVariable
	version 11.1
	syntax anything [if], strata(varlist) [one string]

	if ("`one'" == "") {
		preserve
		quietly keep `if'
		contract `strata'
		drop _freq
		generate long `anything' = _n

		tempfile tmpfil0
		quietly save "`tmpfil0'"

		restore
		quietly merge m:1 `strata' using "`tmpfil0'", assert(match master) nogenerate
		erase "`tmpfil0'"
	}
	else {
		tempname Label
		if ("`string'" == "") {
			tempvar tmpvar0
			generate str `tmpvar0' = string(`strata')
			encode `tmpvar0' `if', generate(`anything') label(`Label')
			drop `tmpvar0'
		}
		else {
			encode `strata' `if', generate(`anything') label(`Label')
		}
		label drop `Label'
	}
end

program define correlation
	version 11.1
	syntax varlist [if], [weight(string)] transformation(string) ///
			[ones(varname) r(name) c(name)] g(name) v(name)

	quietly count `if'
	if r(N) < 2 {
		matrix define `g' = J(1, 1, 0)
		matrix define `v' = J(1, 1, 0)
		exit
	}

	summarize `ones' `if' `weight', meanonly
	local nh `r(N)'

	gettoken response groupvar : varlist

	// Get boldnh
	tempname boldnh
	quietly tabulate `response' `groupvar' `if' `weight', matcell(`boldnh')
	matrix define `boldnh' = vec(`boldnh')

	// Get Phstardot
	tempname Phstardot Rh
	quietly tabulate `groupvar' `if' `weight', matcell(`Phstardot') matrow(`Rh')
	getTransformation , scores(`transformation') marginaltotals(`Phstardot') values(`Rh')
	matrix define `Phstardot' = `Phstardot' / `nh'

	// Get Phdotstar
	tempname Phdotstar Ch
	quietly tabulate `response' `if' `weight', matcell(`Phdotstar') matrow(`Ch')
	getTransformation , scores(`transformation') marginaltotals(`Phdotstar') values(`Ch')
	matrix define `Phdotstar' = `Phdotstar' / `nh'

	computeGV , boldnh(`boldnh') phdotstar(`Phdotstar') phstardot(`Phstardot') ///
		rh(`Rh') ch(`Ch') nh(`nh') g(`g') v(`v')
end

program define anova
	version 11.1
	syntax varlist [if], [weight(string)] transformation(string) ///
			ones(varname) r(name) [c(name)] g(name) v(name)

	local R = rowsof(`r')
	local df = `R' - 1

	quietly count `if'
	if r(N) < 2 {
		matrix define `g' = J(`df', 1, 0)
		matrix define `v' = J(`df', `df', 0)
		exit
	}

	summarize `ones' `if' `weight', meanonly
	local nh `r(N)'

	gettoken response groupvar : varlist

	// Get boldnh
	tempname C boldnh
	quietly tabulate `response' `if', matrow(`C')
	forvalues row = 1/`R' {
		forvalues index = 1/`=rowsof(`C')' {
			summarize `ones' ///
				`if' & `groupvar' == `r'[`row', 1] & `response' == `C'[`index', 1] ///
				`weight', meanonly
			matrix define `boldnh' = (nullmat(`boldnh') \ r(N))
		}
	}

	// Get Phstardot
	tempname Phstardot
	matrix define `Phstardot' = J(`R', 1, 0)
	forvalues row = 1/`R' {
		summarize `ones' ///
			`if' & `groupvar' == `r'[`row', 1] ///
			`weight', meanonly
		matrix define `Phstardot'[`row', 1] = r(N) / `nh'
	}

	// Get Rh
	tempname Rh
	matrix define `Rh' = ( I(`df') , J(`df', 1, -1) )

	// Get Phdotstar
	tempname Phdotstar Ch
	quietly tabulate `response' `if' `weight', matcell(`Phdotstar') matrow(`Ch')
	getTransformation , scores(`transformation') marginaltotals(`Phdotstar') values(`Ch')
	matrix define `Phdotstar' = `Phdotstar' / `nh'

	computeGV , boldnh(`boldnh') phdotstar(`Phdotstar') phstardot(`Phstardot') ///
		rh(`Rh') ch(`Ch') nh(`nh') g(`g') v(`v')
end

program define general
	version 11.1
	syntax varlist [if], [weight(string)] [transformation(string)] ///
			ones(varname) r(name) c(name) g(name) v(name)

	local R = rowsof(`r')
	local C = rowsof(`c')
	local dfR = `R' - 1
	local dfC = `C' - 1

	quietly count `if'
	if r(N) < 2 {
		local df = `dfR' * `dfC'
		matrix define `g' = J(`df', 1, 0)
		matrix define `v' = J(`df', `df', 0)
		exit
	}

	summarize `ones' `if' `weight', meanonly
	local nh `r(N)'

	gettoken response groupvar : varlist

	// Get boldnh
	tempname boldnh
	matrix define `boldnh' = J(`=`R' * `C'', 1, 0)
	forvalues row = 1/`R' {
		local RC = (`row' - 1) * `C'
		forvalues index = 1/`C' {
			summarize `ones' ///
				`if' & `groupvar' == `r'[`row', 1] & `response' == `c'[`index', 1] ///
				`weight', meanonly
			local ++RC
			matrix define `boldnh'[`RC', 1] = r(N)
		}
	}

	// Get Phstardot
	tempname Phstardot
	matrix define `Phstardot' = J(`R', 1, 0)
	forvalues row = 1/`R' {
		summarize `ones' ///
			`if' & `groupvar' == `r'[`row', 1] ///
			`weight', meanonly
		matrix define `Phstardot'[`row', 1] = r(N) / `nh'
	}

	// Get Bh
	tempname Rh
	matrix define `Rh' = ( I(`dfR') , J(`dfR', 1, -1) )

	// Get Phdotstar
	tempname Phdotstar
	matrix define `Phdotstar' = J(`C', 1, 0)
	forvalues index = 1/`C' {
		summarize `ones' ///
			`if' & `response' == `c'[`index', 1] ///
			`weight', meanonly
		matrix define `Phdotstar'[`index', 1] = r(N) / `nh'
	}

	// Get Ch
	tempname Ch
	matrix define `Ch' = ( I(`dfC') , J(`dfC', 1, -1) )

	computeGV , boldnh(`boldnh') phdotstar(`Phdotstar') phstardot(`Phstardot') ///
		rh(`Rh') ch(`Ch') nh(`nh') g(`g') v(`v')
end

program define computeGV
	version 11.1
	syntax , boldnh(name) phdotstar(name) phstardot(name) rh(name) ch(name) nh(integer) g(name) v(name)

	// Get expected row marginal frequencies under the null hypothesis
	tempname mh
	matrix define `mh' = `nh' * (`phstardot' # `phdotstar')

	// Get variance of row marginal frequencies under the null hypothesis
	tempname var
	matrix define `var' = `nh' * `nh' / (`nh' - 1) * ///
		(diag(`phstardot') - `phstardot' * `phstardot'') # ///
		(diag(`phdotstar') - `phdotstar' * `phdotstar'')

	// Get Bh
	tempname Bh
	matrix define `Bh' = `rh' # `ch'

	// Compute Gh, VGh
	matrix define `g' = `Bh' * (`boldnh' - `mh')
	matrix define `v' = `Bh' * `var' * `Bh''
end

program define getTransformation
	version 11.1
	syntax , scores(string) marginaltotals(name) values(name)
	if ("`scores'" == "table") {
		// Do nothing
	}
	else {
		local row_tally = rowsof(`values')

		if ("`scores'" == "integer") {
			forvalues index = 1/`row_tally' {
				matrix define `values'[`index', 1] = `index'
			}
		}
		else {
			matrix define `values' = J(`row_tally', 1, 0)
			tempname cumulator

			if ("`scores'" == "savage") {
				tempname n
				scalar define `n' = 0
				forvalues i = 1/`row_tally' {
					scalar define `n' = `n' + `marginaltotals'[`i', 1]
				}
				tempname A
				matrix define `A' = J(`n', 1, 0)
				forvalues rank = 1/`=`n'' {
					matrix define `A'[`rank', 1] = 1 / (`n' - `rank' + 1)
				}
				forvalues rank = 2/`=`n'' {
					matrix define `A'[`rank', 1] = `A'[`=`rank'-1', 1] + `A'[`rank', 1]
				}

				tempname B start_value cumulator marginaltotal
				matrix define `B' = `values'
				scalar define `start_value' = 1

				forvalues Bindex = 1/`=rowsof(`B')' {
					scalar define `cumulator' = 0
					scalar define `marginaltotal' = `marginaltotals'[`Bindex', 1]
					forvalues Aindex = `=`start_value''/`=`start_value' + `marginaltotal' - 1' {
						scalar define `cumulator' = `cumulator' + `A'[`Aindex', 1]
					}
					matrix define `B'[`Bindex', 1] = `cumulator' / `marginaltotal' - 1
					scalar define `start_value' = `start_value' + `marginaltotal'
				}

				matrix define `values' = `B'
			}
			else if ("`scores'" == "logrank") {
			// Taken from Stokes, Davis & Koch, Page 73, which differs from other algorithms

				forvalues j = 1/`row_tally' {
					forvalues k = 1/`j' {
						scalar define `cumulator' = 0
						forvalues m = `k'/`row_tally' {
							scalar define `cumulator' = `cumulator' + `marginaltotals'[`m', 1]
						}
					matrix define `values'[`j', 1] = `values'[`j', 1] + `marginaltotals'[`k', 1] / `cumulator'
					}
				}
				matrix define `values' = J(`row_tally', 1, 1) - `values'
			}
			else {
				scalar define `cumulator' = 0
				forvalues index = 1/`row_tally' {
					matrix define `values'[`index', 1] = `cumulator' + (`marginaltotals'[`index', 1] + 1) / 2
					scalar define `cumulator' = `cumulator' + `marginaltotals'[`index', 1]
				}

				if ("`scores'" == "rank") {
					// Do nothing
				}
				else if ("`scores'" == "ridit") {
					matrix define `values' = `values' / `cumulator'
				}
				else if inlist("`scores'", "modridit", "vdw", "klotz") {
					matrix define `values' = `values' / (`cumulator' + 1)
					if ("`scores'" == "modridit") {
						// Do nothing
					}
					else {
						forvalues index = 1/`row_tally' {
							matrix define `values'[`index', 1] = invnormal(`values'[`index', 1])
						}
						if ("`scores'" == "vdw") {
							// Do nothing
						}
						else {
							tempname XtX
							if  `row_tally' < 2500 { // < 50 megabyte matrix
								matrix define `XtX' = `values' * `values''
								matrix define `values' = vecdiag(`XtX')'
							}
							else {
								forvalues index = 1/`row_tally' {
									`values'[`index', 1] = `values'[`index', 1] * `values'[`index', 1]
								}
							}
						}
					}
				}
				else if inlist("`scores'", "median", "mood") {
					tempname average_rank // Location of median
					scalar define `average_rank' = (`cumulator' + 1) / 2
					if ("`scores'" == "median") {
						forvalues index = 1/`row_tally' {
							matrix define `values'[`index', 1] = (`values'[`index', 1] > `average_rank')
						}
					}
					else {
						forvalues index = 1/`row_tally' {
							matrix define `values'[`index', 1] = (`values'[`index', 1] - `average_rank')^2
						}
					}
				}
				else {
					display in smcl as error "`scores' transformation not yet implemented"
					exit 178
				}
			}
		}
	}
	matrix define `values' = `values''
end