* Version 1.0; 10/10/97   STB-42 snp14
* Written by Andy Perkins (a.perkins@uea.ac.uk)
* Address: School of Health Policy & Practice, Univ. East Anglia, Norwich, NR4 7TJ, UK.
program define mv2snp
	version 5.0
	local varlist "required existing min(2)"	/* Minimum of 2 variables required */
	local options "BY(string) VERBose"		/* option for grouping variable and */
							/* verbose output */
	parse "`*'"
	parse "`options'", parse(" ")
	if "`by'" == "" {						/* Grouping variable name missing */
		display in red "You must specify a grouping variable using " in yellow ",by(varname)"
		exit 102
	}

	preserve
	
	* Use touse (and variations) for implementing case-wise deletion
	tempvar touse S_touse ngtouse poshold neghold
	tempname vct u v u1 u2 u3 nv nv1
	mark `touse'
	markout `touse' `varlist'
	generate S_touse = `touse'
	generate `ngtouse' = -`touse'

quietly {		/* Suppress output */

	encode `by', generate(dum)
	generate poshold = dum			/* Create 2 variables to sort by */
	generate neghold = -dum			/* so that sub-totals can be calculated */
	drop dum

	scalar vct = 0				/* Counter for number of variables */
	parse "`varlist'", parse(" ")

	while "`1'" ~= "" {
		scalar vct = vct + 1

			local stub=string(vct)
			local vhold "amp`stub'"		/* Create variable name to hold rankings */
			local w1hold "n1`stub'"		/* To carry size of 1st group */
			local x1hold "s1`stub'"		/* To carry sum of ranks of 1st group */
			local w2hold "n2`stub'"		/* To carry size of 2nd group */
			local x2hold "s2`stub'"		/* To carry sum of ranks of 2nd group */

			egen `vhold' = rank(`1') if `touse'

			sort poshold		
			quietly by poshold: /*
				*/ summarize `vhold' if `touse'
			scalar `w1hold' = _result(1)
			scalar `x1hold' = _result(3) * _result(1)
			
			sort neghold		
			quietly by neghold: /*
				*/ summarize `vhold' if `touse'
			scalar `w2hold' = _result(1)
			scalar `x2hold' = _result(3) * _result(1)

		macro shift
	}

	* Variables amp1 to ampZZ now hold ranks of 1st to ZZth variable
	* Size of 1st group on each variable held in n1x where x = 1 to vct
	* Size of 2nd group on each variable held in n2x where x = 1 to vct
	* Sum of ranks of 1st group on each variable held in s1x where x = 1 to vct	
	* Sum of ranks of 2nd group on each variable held in s2x where x = 1 to vct
	* We don't actually need all n1x and n2x values since case-wise deletion will remove
	* cases with incomplete data.
	* Therefore find smallest of n11 and n21. This becomes nsmall
	* ntotal is n11 + n21
	* For each variable, R values are sum of ranks from smallest group

	if n11 == 0 | n21 ==0 {
		display in red "One of the variables does not have any data."
		display in red "This may be due to case-wise deletion."
		exit 2000
	}

	if n11 < n21 {
		global S_nsmall = n11
		global S_small = 1
		global S_nlarge = n21
		global S_large = 2
	}
	else if n11 > n21 {
		global S_nsmall = n21
		global S_small = 2
		global S_nlarge = n11
		global S_large = 1
	}
	else if n11 == n21 {
		global S_vct = vct
		global S_ntotal = n11 + n21

		* Send both combinations of small/large to compute two Ustars
		global S_nsmall = n11
		global S_small = 1
		global S_nlarge = n21
		global S_large = 2
		mvcalcs
		local uuhold1 = trace(u3)

		global S_nsmall = n21
		global S_small = 2
		global S_nlarge = n11
		global S_large = 1
		mvcalcs
		local uuhold2 = trace(u3)

		if `uuhold1' <= `uuhold2' {	/* We want the smaller (worse) U */
			global S_nsmall = n11
			global S_small = 1
			global S_nlarge = n21
			global S_large = 2
		}
		else {
			global S_nsmall = n21
			global S_small = 2
			global S_nlarge = n11
			global S_large = 1
		}
		* I think that the two Ustars are always equal, so this comparison part
		* may not be necessary. However, leave it in to be safe
	}

	global S_vct = vct
	global S_ntotal = n11 + n21

}	/* End of quietly */

	* Perform vector and matrix calculations
	mvcalcs

	display " "
	display in white "Multi-variate non-parametric test (Leach, 1991)"
	display " "
	if "`verbose'" == "verbose" {
		display in blue "Vector of test statistics, U"
		matrix list u ,nohalf noheader nonames format(%9.5f)
		display " "
		display in blue "Variance-Covariance matrix, V"
		matrix list v ,nohalf noheader nonames format(%9.5f)
		display " "
		display in blue "Inverse of matrix NV"
		matrix list nv1 ,nohalf noheader nonames format(%9.5f)
		display " "
	}

	sort `ngtouse'				/* Find number of dropped cases */
	quietly by `ngtouse': summarize `touse' if `touse' == 0
	local drppd = _result(1)

	local w = trace(u3)
	local vv = vct 
	local p = chiprob(`vv',`w')
	display _dup(79) "-"
	display in white "Test on variables: " in yellow "`varlist'"
	parse "`varlist'", parse(" ")

	while "`1'" ~= "" {
		local varlab: variable label `1'
		display in yellow _col(20) "`1' " _col(29) in blue "- `varlab'"
		macro shift
	}
*	 in yellow "`varlist'"

	local varlab: variable label `by'
	display in white "Grouping variable: " in yellow "`by'" _col(29) in blue "- `varlab'"
	if `drppd' ~= 0 {
		display in white "Dropped cases: " in yellow "`drppd'"
	}
	display " "
	display "Size of smaller group = $S_nsmall" _col(35) "Number of variables = " vct
	display "Size of larger group  = $S_nlarge" _col(35) "U* = " %9.2f `w'
	display "Combined sample size  = $S_ntotal" _col(35) "P = " %9.5f `p'
	display _dup(79) "-"

	* Save results for other programs
	global S_1 = `w'
	global S_2 = $S_ntotal
	global S_3 = $S_nsmall
	global S_4 = vct

end


program define mvcalcs
	version 5.0

quietly {

	local small = $S_small
	local ntotal = $S_ntotal
	local nsmall = $S_nsmall
	local large = $S_large
	local nlarge = $S_nlarge


	* Calculate vector U
	local rws = $S_vct
	matrix u = J(`rws',1,0)

	local lpi = 1
	while `lpi' <= `rws' {
		local r = s`small'`lpi'
		local k = (`r'/(`ntotal'+1))-(`nsmall'/2)
		local rw = `lpi'
		matrix u[`rw',1] = `k'
		local lpi = `lpi' + 1
	}

	* Calculate matrix V
	matrix v = J(`rws',`rws',0)

	* Off-diagonal entries
	local lpi = 1
	local lpj = 1
	while `lpi' <= `rws' {
		while `lpj' <= `rws' {
			gen amprod = amp`lpi' * amp`lpj' if S_touse  	/* Calculate cross-products on a pair of variables */
			quietly summarize amprod if S_touse
			local rr = _result(3) * _result(1)
			drop amprod
			local k1 = (`nlarge'*`nsmall')/((`ntotal'*`ntotal')*(`ntotal'-1)*(`ntotal'+1)*(`ntotal'+1))
			local k2 = `rr'-(`ntotal'*((`ntotal'+1)*(`ntotal'+1)))/4
			local k = `k1' * `k2'
			local rw = `lpi'
			local cl = `lpj'
			matrix v[`rw',`cl'] = `k'
			local lpj = `lpj' + 1
		}
		local lpj = 1			/* Reset lpj */
		local lpi = `lpi' + 1
		
	}

	* Add-in (replace) the diagonals as it is easier to program this way
	local lpi = 1
	while `lpi' <= `rws' {
		local k = (`nlarge'*`nsmall')/(12*`ntotal'*(`ntotal'+1))
		local rw = `lpi' 
		matrix v[`rw',`rw'] = `k'
		local lpi = `lpi' + 1
	}

	matrix u1 = u'
	matrix nv = `ntotal' * v
	matrix nv1 = inv(nv)
	matrix u2 = u1 * nv1
	matrix u3 = u2 * u		/* u3 is Ustar */

}	/* End of quietly */

end