/* mtest - returns the observed M statistic and the p-value (Monte Carlo resampling) of M under the null hypothesis that Cases and Controls follow the same spatial distribution
Pietro Tebaldi , Harvard School of Public Health and Bocconi University
July 2010 */

program define mtest , rclass
version 10.1
syntax , x(varname) y(varname) g(varname) [ bins(integer 20) iter(integer 100) level(real 95) SCatter DENsity ] ///
		[scolor0(string asis)] 	///
		[scolor1(string asis)] 	///
		[smarker0(string)] 		///
		[smarker1(string)] 		///
		[ssize0(string asis)] 		///
		[ssize1(string asis)] 		///
		[slabel0(string)] 		///
		[slabel1(string)] 		///
		[stitle(string)] 		///
		[sytitle(string)] 		///
		[sxtitle(string)] 		///
		[dcolor0(string)] 		///
		[dcolor1(string)]		///
		[dpattern0(string)] 	///
		[dpattern1(string)] 	///
		[dwidth0(string)] 		///
		[dwidth1(string)]		///
		[dlabel0(string)]		///
		[dlabel1(string)]		///
		[dtitle(string)]

tempname Sinv d

capture count
local N = r(N)
forvalues i = 1/`N' {
if `g'[`i']!=1 & `g'[`i']!=0 {
display as error "The group variable g must be bivariate 0,1"
error
}
}

// SCATTER OPTIONS
if "`scatter'" == "scatter" {
// color
	foreach j in 0 1 {
		if "`scolor`j''" != "" {
			local scolor`j' = "`scolor`j''"
		}
		else {
			local scolor`j' = "black"
		}
	}
// marker
	if "`smarker0'" != "" {
		local smarker0 = "`smarker0'"
	}
	else {
		local smarker0 = "oh"
	}
	if "`smarker1'" != "" {
		local smarker1 = "`smarker1'"
	}
	else {
		local smarker1 = "x"
	}

// size
	foreach j in 0 1 {
		if "`ssize`j''" != "" {
			local ssize`j' = "`ssize`j''"
		}
		else {
			local ssize`j' = "medium"
		}
	}

// label
	foreach j in 0 1 {
		if "`slabel`j''" != "" {
			local slabel`j' = "`slabel`j''"
		}
		else {
			local slabel`j' = "Group `j'"
		}
	}

// title
	if "`stitle'" != "" {
		local stitle = "`stitle'"
	}
	else {
		local stitle = "Spatial Distribution of the two groups"
	}

// axis title
	if "`sytitle'" != "" {
		local sytitle = "`sytitle'"
	}
	else {
		local sytitle = "`y'"
	}
	if "`sxtitle'" != "" {
		local sxtitle = "`sxtitle'"
	}
	else {
		local sxtitle = "`x'"
	}


	tempvar X1 X2 Y1 Y2
	quietly {
		gen `X1' = `x' if `g'==0
		gen `Y1' = `y' if `g'==0
		gen `X2' = `x' if `g'==1
		gen `Y2' = `y' if `g'==1
		label var `X1' "`slabel0'"
		label var `Y1' "`slabel0'"
		label var `X2' "`slabel1'"
		label var `Y2' "`slabel1'"
	}
	capture scatter `Y1' `X1'  , msymbol(`smarker0') mcolor(`scolor0') msize(`ssize0') ytitle(`sytitle') xtitle(`sxtitle') title(`stitle') || scatter `Y2' `X2' , msymbol(`smarker1') mcolor(`scolor1') msize(`ssize1') , nodraw saving(scatter)
}
// DENSITY OPTIONS
if "`density'" == "density" {

// color
	foreach j in 0 1 {
		if "`dcolor`j''" != "" {
			local dcolor`j' = "`dcolor`j''"
		}
		else {
			local dcolor`j' = "black"
		}
	}

// pattern
	if "`dpattern0'" != "" {
		local dpattern0 = "`dpattern0'"
	}
	else {
		local dpattern0 = "dash"
	}
	if "`dpattern1'" != "" {
		local dpattern1 = "`dpattern1'"
	}
	else {
		local dpattern1 = "solid"
	}

// width
	foreach j in 0 1 {
		if "`dwidth`j''" != "" {
			local dwidth`j' = "`dwidth`j''"
		}
		else {
			local dwidth`j' = "medium"
		}
	}

// label
	foreach j in 0 1 {
		if "`dlabel`j''" != "" {
			local dlabel`j' = "`dlabel`j''"
		}
		else {
			local dlabel`j' = "Group `j'"
		}
	}


// title
	if "`dtitle'" != "" {
		local dtitle = "`dtitle'"
	}
	else {
		local dtitle = "IDD Kernel Densities"
	}


	preserve
	capture count
	local N = r(N)
	local Ndist = comb(`N',2)
	quietly {
		gen X1 = `x' if `g'==0
		gen Y1 = `y' if `g'==0
		drop if X1 == .
		eucldist X1 Y1
		mata: g0 = st_data(.,1)
		clear
		restore
		preserve
		gen X2 = `x' if `g'==1
		gen Y2 = `y' if `g'==1
		drop if X2 == .
		eucldist X2 Y2
		mata: g1 = st_data(.,1)
		clear
	}
	mata: Nca = rows(g1)
	mata: Nco = rows(g0)
	quietly {
		set obs `Ndist'
		mata: ca = st_addvar("float","dlabel1")
		mata: co = st_addvar("float","dlabel0")
		mata: st_store((1,Nca),ca,g1)
		mata: st_store((1,Nco),co,g0)
		label var dlabel1 "`dlabel1'"
		label var dlabel0 "`dlabel0'"
		kdensity dlabel1 , gen(x1 d1) nograph
		kdensity dlabel0 , gen(x2 d2) nograph
		line d2 x2 , lcolor(`dcolor1') lpattern(`dpattern1') lwidth(`dwidth1') ytitle("density") xtitle("interpoint distance") title(`dtitle') || line d1 x1 , lcolor(`dcolor0') lpattern(`dpattern0') lwidth(`dwidth0')  , nodraw saving(density)
		clear
		restore
	}
}
if "`scatter'" == "scatter" {
	graph use scatter.gph
	erase scatter.gph
}
if "`density'" == "density" {
	graph use density.gph
	erase density.gph
}

capture mstat , x(`x') y(`y') g(`g') bins(`bins')
matrix `d' = r(d)
matrix `Sinv' = r(Sinv)
tempname b c reps p se ci
capture permute `g' M=r(M) , reps(`iter') right level(`level') nodots nowarn noheader: Mst_iter , x(`x') y(`y') g(`g') sinv(`Sinv') d(`d')
capture return list
matrix `b' = r(b)
matrix `c' = r(c)
matrix `reps' = r(reps)
matrix `p' = r(p)
matrix `se' = r(se)
matrix `ci' = r(ci)
display "{txt}{hline 65}"
display as text "M statistic"
display as text "Monte Carlo permutation results"
display as text "H0: The two groups have the same spatial distribution"
display as text "Number of bins = `bins'"
display as text "Number of permutations = `iter'"

tempname Tab
	.`Tab' = ._tab.new, col(7) lmargin(0) ignore(.b)
	ret list
	// column           1      2     3     4     5     6     7
	.`Tab'.width	12     |8     8     8     8    10    10
	.`Tab'.titlefmt  .     .     .     .     .  %20s     .
	.`Tab'.pad	     2     0     0     0     0     0     1
	.`Tab'.numfmt   %9.0g     .     . %7.4f %7.4f     .     .
// begin display
	.`Tab'.sep, top
	.`Tab'.titles "M(obs)" "c" "n" "p=c/n" "SE(p)" ///
		"[`level'% Conf. Interval]" ""
	.`Tab'.sep
	.`Tab'.row 	`b'[1,1] `c'[1,1] `reps'[1,1] `p'[1,1] `se'[1,1] `ci'[1,1]	`ci'[2,1]
    .`Tab'.sep, bottom
dis as text "note: c = #{M>=M(obs)}"
dis as text "note: exact binomial confindece interval with respect to p=c/n"
return matrix M = `b'
return matrix c = `c'
return scalar N = `N'
return matrix d = `d'
return matrix Sinv = `Sinv'
return matrix reps = `reps'
return matrix p = `p'
return matrix se = `se'
return matrix ci = `ci'

end