*!chaid version 2.2 2/12/2015 - Joseph N. Luchman										
program define chaid, eclass	//version history at end of file								//program declarations

version 12.1																					//version declaration

if (replay()) {																					//replay syntax

	if ("`e(cmd)'" != "chaid") error 301
	
	if (_by()) error 190
	
	Display
	
	exit 0
	
}																					

syntax varlist(numeric min=1 max=1) [in] [if] [fw], ///											//only one dv and max 1 dv - allows if and in
[MINNode(integer 100) ///																		//minimum number of observations in a terminal node
MINSplit(integer 200) ///																		//minimum number of observations to allow splitting a node
MAXBranch(integer -1) ///																		//maximum number of splits/interactions stemming from the root node
UNOrdered(varlist numeric) ///																	//variable list for unordered splitting treatment in CHAID
ORDered(varlist numeric) ///																	//variable list for ordered splitting treatment in CHAID
MISsing ///																						//treat missing values as floating category
MERGAlpha(real .95) ///																			//alpha level for merging 2 allowable categories									
RESPAlpha(real -1) ///																			//alpha level for respliting 2 allowable, merged categories
SPLTAlpha(real .05) ///																			//alpha level for implementing a split in the data
NOIsily ///																						//shows a lot of output - if the user is interested
DVOrdered ///																					//declares that the outcome is ordered and uses ologit instead of mlogit
noadj ///																						//no Bonferroni adjustment - more permissive splitting
noDisplay	///																					//don't display the results
PREDicted ///																					//obtain predicted values - modal value for category
IMPortance ///																					//obtain permutation importance for input splitting variables
Permute ///																						//use permute to obtain p-values
svy ///																							//incorporate svy design to obtain p-values
///mi(string)	///																				//incorporate multiple imputation to obatain p-values (next relaease?)
XTile(string) ///																				//convenience command for obtaining ordered quantiles of predictors
EXhaust]																						//use exhaustive CHAID (binary split) instead of default multi-way splits

/*exit and warning conditions*/
if ((`:list sizeof unordered' == 0) & (`:list sizeof ordered' == 0) & ///
(`:list sizeof xtile' == 0)) {																	//exit if no splitting vars entered

	display "{err}No independent variables entered."
	
	exit 111
	
}

if (`minsplit' < `minnode') {																	//exit if you tell chaid that the minnode/cluster size is larger than the min size to make a split

	display "{err}{opt minsplit()} cannot be smaller than " ///
	"{opt minnode()}."

	exit 198
	
}

if ((`mergalpha' <= 0) | (`mergalpha' >= 1) | (`spltalpha' <= 0) | ///
(`spltalpha' >= 1)) {																			//ensure that the alphas are possible

	display "{err}{opt mergalpha()} and {opt spltalpha()} must range " ///
	"betweeen (but not include) 0 and 1."
	
	exit 198

} 

if (`minnode' <= 0) {																			//ensure that min cluster size is a possible number

	display "{err}{opt minnode()} must be a positive integer value."
	
	exit 198

}

if ((`respalpha' >= 1) | (`respalpha' > 1 - `mergalpha') | (`respalpha' == 0)) {				//ensure that alphas are possible values that don't induce splitting forever

	display "{err}{opt resplapha()}, when invoked, must range between (but not include)" _newline ///
	"0 and 1 and must be less than 1 - {opt mergalpha()} to prevent infinite merging and resplitting."
	
	exit 198

}

if (("`weight'" != "") & (("`permute'" != "") | ("`svy'" != ""))) {

	display "{err}{cmd:fweight}s not allowed with {opt permute} or {opt svy}."

	exit 198

}

if (("`permute'" != "") & ("`svy'" != "")) {

	display "{err}{opt permute} and {opt svy} not allowed together."

	exit 198

}

//1) no fw and permute, 2) no fw and svy, 3) no permute and svy

capture which lmoremata.mlib																	//is moremata present?

if (_rc != 0) {																					//if moremata cannot be found, tell user to install it.

	display "{err}Module {cmd:moremata} not found.  Install {cmd:moremata} here {stata ssc install moremata}."
	
	exit 198

}

if strlen("`permute'") local noadj "noadj"

/*begin data processing*/
//mark the estimation sample
tempvar keep touse stop cluster useob dv														//tempvar declarations - begin using "touse"; add "keep" markout variables in unordered() and ordered()

tempname cellcount fit																			//tempname declaration for matrix

quietly generate `keep' = 1 `if' `in'															//allows keeping of missings as valid category, but still adjusting for "if"s and "in"s

if strlen("`xtile'") {

	gettoken xtile xtile_opt: xtile, parse(",")

	foreach x of varlist `xtile' {																//convenience command for creating binned continuous predictors in CHAID

		capture inspect xt`x'
	
		if !_rc quietly drop xt`x'

		quietly xtile xt`x' = `x' if `keep' [`weight'`exp'] `xtile_opt'
	
		local ordered "`ordered' xt`x'"

	}
	
}

mark `touse'																					//mark estimation sample

if ("`missing'" == "") markout `touse' `varlist' `unordered' `ordered' `keep'					//markout all variables when missings as category is not desired

else if ("`missing'" != "") & ("`dvordered'" == "") markout `touse' `keep'						//markout only "if" and "in" when missings are a valid category

else if ("`missing'" != "") & ("`dvordered'" != "") markout `touse' `varlist' `keep'			//markout only "if" and "in" and on dv when missings are valid category and dv is declared ordered

drop `keep'																						//"keep" tempvar is dropped so that it does not find its way into mata

foreach x of varlist `varlist' `unordered' `ordered' {											//exit if too many categories (20 is max) or categories are negative integers

	quietly inspect `x' if `touse'
	
	if (r(N_unique) > 20) {
	
		display "{err}Number of distinct values in {cmd:`x'} is larger than allowed. " _newline ///
		"Consider collapsing across similar categories for unordered or using " _newline ///
		"{opt xtile()} option/command for ordered variables to reduce the number of levels."
		
		exit 198
	
	}
	
	else if (r(N) != r(N_0) + r(N_posint)) {
	
		display "{err}Levels of each variable must be non-negative integer values." _newline ///
		" Variable {cmd:`x'} contains negative or decimal values."
		
		exit 198
	
	}

}

//identify and process variables to use
quietly tab `varlist' [`weight'`exp'], matrow(`cellcount')										//identify levels of dependent/outcome variable to create inidicator varaibles

mata: st_local("dvlevs", invtokens(strofreal(st_matrix("`cellcount'")')))						//put levels into a macro for use
		
local dvlist "`varlist'"																		//produce local macro that's updated with tempvar names for use in mata

local dvdisp "`varlist'"																		//produce local macro that's updated with displayed names of outcome and levels
	
local type "`type' d"																			//note "type" of variable for dv

quietly generate `dv' = `varlist'																//local macro to denote variable type for splitting purposes - here "d" denotes dependent variable that will not be split

if ("`missing'" != "") quietly count if missing(`dv')											//discern whether outcome is missing values if missings are of interest

if (("`missing'" != "") & (r(N) > 0) & ("`dvordered'" == "")) {									//if "missing" invoked and dv is not ordered, turn missing into last category
	
	quietly summarize `dv' [`weight'`exp']														//find out what the larget category is
	
	quietly recode `dv' (missing = `=r(max)+1')													//make all missings one bigger than the biggest category													 
	
}

if (`:list sizeof unordered' > 0) {																//process unordered splitting variables when present

	foreach x of varlist `unordered' {

		quietly tabulate `x' if `touse' [`weight'`exp'], matrow(`cellcount') `missing'			//discern levels of each unordered variable for creating binary variables
	
		mata: st_local("`x'levs", invtokens(strofreal(st_matrix("`cellcount'")')))				//put levels of unordered variable into local macro
		
		foreach y of numlist ``x'levs' {
		
			if (`:list posof "`y'" in `x'levs' == 1) {
			
				local type "`type' ("
				
				local unorddisp "`unorddisp' ("
			
			}
		
			if ("`y'" != ".") {
			
				tempvar `x'`y'																	//create tempvars with name of level of unordered, binary variable
				
				quietly generate byte ``x'`y'' = `x' == `y'										//generate binaries for each level of unordered variable
	
				label variable ``x'`y'' "`x' @ `y'"												//label the tempvar for later processing
	
				local unorderedlist "`unorderedlist' ``x'`y''"									//include name of unordered binary variable in list of tempvar names
			
				local unorddisp "`unorddisp' `x'`y'"											//include displayed name of unordered binary variable in list of tempvar names
			
				local type "`type' u"															//include variable type in type macro
	
			} 

			else {																				//if missing option invoked count missings to see whether following syntax needs to be run
			
				tempvar `x'ms																	//assuming missings are present, create tempvar for unordered varible indicating missing
				
				quietly generate byte ``x'ms' = cond(missing(`x'), 1, 0)						//generate indicator for missing category
				
				label variable ``x'ms' "`x' @ missing"											//label the tempvar for later processing
				
				local unorderedlist "`unorderedlist' ``x'ms'"									//include missing tempvar name in unordered list
			
				local unorddisp "`unorddisp' `x'ms"												//include missing display name in unordered display list
				
				local type "`type' f"															//include type "floating"
					
			}
			
			if (`: list posof "`y'" in `x'levs' == `: list sizeof `x'levs') {
			
				local type "`type')"
				
				local unorddisp "`unorddisp')"
			
			}
				
		}
			
	}
			
}
	
if (`:list sizeof ordered' > 0) {

	foreach x of varlist `ordered' {

		quietly tab `x' if `touse' [`weight'`exp'], matrow(`cellcount') `missing'				//if ordered variables present, discern their levels/values
	
		mata: st_local("`x'levs", invtokens(strofreal(st_matrix("`cellcount'")')))				//put levels of ordered variable into local macro
	
		foreach y of numlist ``x'levs' {
		
			if (`: list posof "`y'" in `x'levs' == 1) {
			
				local type "`type' ("
				
				local orddisp "`orddisp' ("
			
			}
		
			if ("`y'" != ".") {
			
				tempvar `x'`y'																	//declare tempvar for each value of ordered variable
				
				quietly generate byte ``x'`y'' = `x' == `y'										//generate binary for each level of ordered variable
		
				label variable ``x'`y'' "`x' @ `y'"												//label the tempvar for later processing
		
				local orderedlist "`orderedlist' ``x'`y''"										//include ordered tempvar in orderedlist 
			
				local orddisp "`orddisp' `x'`y'"												//include display name of variable in ordered display macro
			
				local type "`type' o"															//include variable type in type macro
				
			}
	
			else {																				//if missing option invoked count missings to see whether following syntax needs to be run
			
				tempvar `x'ms																	//assuming missings are present, create tempvar for ordered varible indicating missing
				
				quietly generate byte ``x'ms' = cond(missing(`x'), 1, 0)						//generate indicator for missing category
					
				label variable ``x'ms' "`x' @ missing"											//label the tempvar for later processing					
					
				local orderedlist "`orderedlist' ``x'ms'"										//include missing tempvar name in ordered list
			
				local orddisp "`orddisp' `x'ms"													//include missing display name in unordered display list
						
				local type "`type' f"															//include type "floating"
			
			}
			
			if (`: list posof "`y'" in `x'levs' == `: list sizeof `x'levs') {
			
				local type "`type')"
				
				local orddisp "`orddisp')"
			
			}
					
		}
		
	}

}

/*process estimation options*/
if ("`permute'" != "") local permute "permute `dv' e(chi2) e(ll) e(df_m) e(N):"

local test "e(chi2)"

if ("`svy'" != "") local test "e(F)"

local quietly "quietly"																			//is the chaid process going to display everything?

if ("`noisily'" != "") local quietly ""															//make it non-quiet

quietly generate byte `stop' = 1 if `touse'														//indicator telling chaid to not look for splits on this observation

quietly generate byte `cluster' = 1 if `touse'													//cluster membership - everyone starts in same cluster #1

quietly generate byte `useob' = 0 if `touse'													//indicator for "use this ob" in current analysis

`quietly' summarize `stop' `cluster' [`weight'`exp']											//look at data if non-quiet

local reg "mlogit"																				//default analysis is mlogit

if ("`dvordered'" ! = "") local reg "ologit"													//if dvordered then ologit

/*begin chaid algorithm in Mata*/
`quietly' mata: CHAIDsplit = chaid("`unorderedlist'", "`orderedlist'", "`type'", ///
"`dvdisp' `unorddisp' `orddisp'", "`dv' `unordered' `ordered'", `mergalpha', ///
`spltalpha', `minsplit', `minnode', `maxbranch', `respalpha', "`adj'", "`svy'", ///
"`exhaust'")																					//pass arguments to Mata

/*process chaid results*/
capture xi: count																				//nix xi: used in chaid - cleans up dataset														

capture inspect _CHAID																			//is there already a _CHAID variable?

if (_rc == 0) drop _CHAID																		//if there is a _CHAID variable, drop it

rename `cluster' _CHAID																			//rename the current tempvar "cluster" _CHAID

label variable _CHAID "CHAID-defined cluster"													//label the clustering variable

perm_import if `touse', dv("`varlist'") list(0) wexp("`exp'") `svy'

scalar `fit' = r(fit)

if ("`importance'" != "") mata: ifstmnts = J(rows(CHAIDsplit)-1, cols(CHAIDsplit), "")

mata: st_local("clusters", strofreal(rows(CHAIDsplit)-1))

if (`clusters' > 0) {

	forvalues x = 1/`clusters' {

		if ("`noisily'" == "noisily") mata: CHAIDsplit[`=`x'+1', .]
		
		local path`x' ""
	
		mata: st_local("search", strofreal(cols(CHAIDsplit)))
	
		forvalues z = 2/`search' {
	
			mata: st_local("currentname", CHAIDsplit[`=`x'+1', `z'])
	
			local vallist ""
		
			if ("`currentname'" != "") {
	
				foreach y of local currentname {
		
					local name "`:variable label ``y'''"
			
					gettoken varnam level: name, parse("@")
			
					gettoken atsign level: level, parse("@")
			
					if ("`level'" != "") local vallist "`vallist'`level'"
		
				}
			
				local varnam = trim("`varnam'")
			
				local vallist = trim("`vallist'")
		
				local path`x' "`path`x''`varnam'@`vallist';"
		
				if ("`importance'" != "") {
		
					if regexm("`vallist'", "missing") local vallist = subinword("`vallist'", "missing", ".", .)
		
					local ifcomma = subinstr("`vallist'", " ", ",", .)
				
					if (`z' == 2) mata: ifstmnts[`x', `=`z'-1'] = " inlist(" + st_local("varnam") + "," + st_local("ifcomma") + ")"
				
					else mata: ifstmnts[`x', `=`z'-1'] = " & inlist(" + st_local("varnam") + "," + st_local("ifcomma") + ")"
			
				}
			
			}
	
		}
	
	}
	
}


if ("`importance'" != "") & (`clusters' > 0) {
	
	tempname imp_v
	
	mata: imp_n = ""
	
	mata: imp_v = .
	
	foreach x of varlist `ordered' `unordered' {
	
		perm_import if `touse', permvar("`x'") dv("`varlist'") list(`clusters') ///
		predlist(`ordered' `unordered') wexp("`exp'") `missing' `svy'
		
		mata: imp_n = (imp_n, "`x'")
		
		mata: imp_v = (imp_v, st_numscalar("r(fit)"))
	
	}
	
	mata: imp_v = st_numscalar("`fit'"):-imp_v
	
	mata: imp_v = imp_v:/rowsum(imp_v)
	
	mata: imp_v = (imp_v \ (., mm_ranks((imp_v[1, 2..cols(imp_v)]:*-1)',1,3)'))
	
	mata: st_local("imp_n", invtokens(imp_n))
	
	mata: st_matrix("`imp_v'", imp_v[., 2..cols(imp_v)])
	
	matrix colnames `imp_v' = `imp_n'
	
	matrix rownames `imp_v' = "raw" "rank"

}

ereturn post , depname(`varlist') esample(`touse')												//clear ereturned values, post dv and obs

mata: st_local("search", strofreal(cols(CHAIDsplit)))

local splitcount = 1

forvalues x = 2/`search' {
	
	mata: t = tokeninit("", pchars=(","))														//set up processing CHAIDsplit matrix

	mata: tokenset(t, CHAIDsplit[1, `x'])

	mata: nameslist = tokengetall(t)
	
	mata: st_local("search_name", strofreal(cols(nameslist)))
	
	forvalues y = 1/`search_name' {
	
		local levels ""
	
		mata: st_local("currentname", nameslist[`y'])
		
		if !regexm("`currentname'", ",") {
		
			foreach z of local currentname {
		
				local name "`:variable label ``z'''"
			
				gettoken name level : name, parse("@")
			
				local level = subinstr("`level'", "@", "", .)
			
				local levels = rtrim("`levels' `level'")
				
			}
			
			local alllevels "`alllevels' (`levels')"
			
			if ("`quietly'" == "") {
		
				display "`name'"	//
			
				display "`alllevels'"	//
			
			}
			
		}
		
	}
			
	ereturn local split`splitcount++' = itrim("`name' `alllevels'")
			
	local alllevels ""	
	
}

ereturn matrix branches = `branchmat'

quietly summarize _CHAID if e(sample) [`weight'`exp'], meanonly

ereturn scalar N = r(N)

quietly inspect _CHAID

ereturn scalar N_clusters = r(N_unique)

forvalues x = 1/`=e(N_clusters)' {

	ereturn local path`x' "`path`x''"

}

quietly tabulate _CHAID, matcell(`cellcount')

matrix `cellcount' = `cellcount''

ereturn matrix sizes = `cellcount'

if ("`importance'" != "") & (`clusters' > 0) ereturn matrix importance = `imp_v'

if ("`predicted'" == "predicted") {

	tempvar order
	
	tempname maxcell maxval
	
	capture inspect CHAID_predict

	if (_rc == 0) drop CHAID_predict

	quietly inspect _CHAID
	
	quietly generate long CHAID_predict = .
	
	forvalues x = 1/`=r(N_unique)' {
	
		quietly tabulate `varlist' [`weight'`exp'] if _CHAID == `x', matcell(`maxcell') matrow(`maxval') `missing'
		
		mata: predicted = select(st_matrix("`maxval'"), st_matrix("`maxcell'"):==max(st_matrix("`maxcell'")))
		
		mata: st_numscalar("`maxcell'", rows(predicted))
		
		if (`maxcell' == 1) mata: st_numscalar("`maxval'", predicted)
		
		else mata: st_numscalar("`maxval'", min(predicted))
		
		quietly replace CHAID_predict = `maxval' if _CHAID == `x'
	
	}
	
	label variable CHAID_predict "predicted `varlist' values by CHAID"
	
}

ereturn scalar fit = `fit'

ereturn local title = "Chi-Square Automated Interaction Detection (CHAID)"

ereturn local cmd = "chaid"

ereturn local cmdline = "chaid `0'"


if ("`display'" != "nodisplay") Display

end

/*model fit program*/
program perm_import, rclass

syntax if, dv(string) list(integer) [permvar(varlist) predlist(string) wexp(string) missing svy]

tempvar temp_chaid rand_perm

tempfile chaid_perm

preserve
		
if (!strlen("`svy'")) quietly keep `if'

if strlen("`wexp'") quietly expand `wexp'

quietly generate int `temp_chaid' = .

if (`list' == 0) {

	quietly replace `temp_chaid' = _CHAID

}

else {

	local notpermvar: list predlist-permvar
	
	if (!strlen("`svy'")) quietly keep `predlist' `dv' `temp_chaid'
		
	quietly save `chaid_perm'
		
	quietly keep `permvar'
		
	quietly generate `rand_perm' = runiform()
		
	quietly sort `rand_perm'
	
	if (strlen("`svy'")) quietly merge 1:1 _n using `chaid_perm', nogenerate
		
	else quietly merge 1:1 _n using `chaid_perm', nogenerate keepusing(`notpermvar' `dv' `temp_chaid')

	forvalues x = 1/`list' {

		mata: stata("quietly replace \`temp_chaid' = \`x' if" + invtokens(ifstmnts[`x', .]))

	}
	
}

if (strlen("`svy'")) {

	tempvar svy
	
	tempname best_F
	
	quietly clonevar `svy' = `dv'
	
	quietly svy: tabulate `dv' `svy', `missing'
	
	scalar `best_F' = e(F_Full)
	
	quietly svy: tabulate `dv' `temp_chaid', `missing'
	
	capture scalar `best_F' = e(F_Full)/`best_F'
	
	if (!_rc) return scalar fit = `best_F'
	
	else return scalar fit = 0

}

else {

	quietly tabulate `dv' `temp_chaid', V `missing'

	capture assert r(CramersV) != .
	
	if (!_rc) return scalar fit = abs(r(CramersV))

	else return scalar fit = 0
	
}

restore

end

/*replayable display*/
program Display

display _newline "{res}Chi-Square Automated Interaction Detection (CHAID)" ///
" Tree Branching Results" _newline "{txt}{hline 80}" _newline

if (e(N_clusters) == 1) display "{res}No clusters uncovered.  Cluster #1 is null."

else {

	preserve
	
	clear
	
	mata: st_addobs(st_numscalar("e(N_clusters)"))
	
	quietly mata: labels = st_tempname(max(st_matrix("e(branches)")) + 1)
	
	quietly mata: st_addvar("str100", labels)
	
	forvalues x = 1/`=e(N_clusters)' {
		
		local parse "`e(path`x')'"
		
		gettoken first second: parse, parse(";")
		
		local count = 0
		
		while ("`second'" != "") {
		
			mata: st_local("labvar", labels[1, `++count'])
		
			quietly replace `labvar' = "`first'" in `x'
			
			gettoken first second: second, parse(";")
			
			gettoken first second: second, parse(";")
		
		}
		
		mata: st_local("labvar", labels[1, cols(labels)])
		
		quietly replace `labvar' = "Cluster #`x'" in `x'
	
	}
	
	mata: st_local("all", invtokens(labels))
	
	sort `all'
	
	quietly putmata namemat = (`all'), replace
	
	mata: namemat'
	
	mata: namemat =  namemat[., 1..cols(namemat)-1]'

	mata: namemat = (J(1, cols(namemat), "root") \ namemat)
	
	clear
	
	mata: count = 0
	
	mata: column = 1

	local column = 1

	mata: st_local("totcol", strofreal(cols(namemat)))
	
	quietly {

		generate long x = .
	
		generate long y = .

		generate int col = .

		generate byte down = .
	
		generate str100 label = ""

	}
	
	local justup = 0
	
	forvalues y = 1/`totcol' {

		mata: currentvec = namemat[., column]
	
		if (`y' == 1) {
	
			mata: currentvec = select(currentvec, currentvec:!="")
		
			mata: st_local("row", strofreal(rows(currentvec)))
		
			mata: row = 1
		
			forvalues x = 1/`row' {
		
				mata: st_addobs(1)
			
				mata: st_store(++count, st_varindex("x"), column)
			
				mata: st_store(count, st_varindex("y"), row)
			
				mata: st_sstore(count, st_varindex("label"), namemat[row++, column])
		
			}
		
			quietly replace down = 1 if missing(down)
		
			mata: countback = st_data(2::count-1, (st_varindex("x"), st_varindex("y")))
		
			mata: countback = sort(((rows(countback)::1), countback), 1)
		
			mata: st_addobs(rows(countback))
		
			mata: st_store(count+1::count+rows(countback), (st_varindex("x"), st_varindex("y")), countback[., 2..cols(countback)])
		
			mata: count = count + rows(countback)
		
			quietly replace down = 0 if missing(down)
		
			quietly replace col = `y' if missing(col)
		
			mata: column++
	
		}
	
		else if (`y' == `totcol') {
	
			mata: prevvec = namemat[., column-1]
		
			mata: retrace1 = currentvec:==prevvec
		
			mata: retrace2 = currentvec:!=""
		
			mata: retrace = retrace1:*retrace2
		
			mata: obsvd = st_data(., (st_varindex("x"), st_varindex("y"), st_varindex("col"), st_varindex("down")))
		
			mata: obsvd = select(obsvd[., (1..2,4)], obsvd[.,3]:==column-1)
		
			mata: obsvd = select(obsvd[., 1..2], obsvd[.,3]:==1)
		
			mata: obsvd = select((obsvd \ J(rows(retrace)-rows(obsvd), 2, .)), retrace)
		
			mata: st_addobs(rows(obsvd))
		
			mata: st_store(count+1::count+rows(obsvd), (st_varindex("x"), st_varindex("y")), obsvd)
		
			mata: count = count + rows(obsvd)
		
			mata: currentvec = select(currentvec, abs(retrace:-1))
		
			mata: currentvec = select(currentvec, currentvec:!="")
		
			mata: st_local("row", strofreal(rows(currentvec)))
		
			mata: row = rows(obsvd)+1
		
			forvalues x = 1/`row' {
		
				mata: st_addobs(1)
			
				mata: st_store(++count, st_varindex("x"), column)
			
				mata: st_store(count, st_varindex("y"), row)
			
				mata: st_sstore(count, st_varindex("label"), namemat[row++, column])
		
			}	
		
			quietly replace col = `y' if missing(col)
		
			quietly replace down = 1 if missing(down)
	
		}
	
		else {
	
			mata: prevvec = namemat[., column-1]
		
			mata: retrace1 = currentvec:==prevvec
		
			mata: retrace2 = currentvec:!=""
		
			mata: retrace = retrace1:*retrace2
		
			mata: obsvd = st_data(., (st_varindex("x"), st_varindex("y"), st_varindex("col"), st_varindex("down")))
		
			mata: obsvd = select(obsvd[., (1..2,4)], obsvd[.,3]:==column-1)
		
			mata: obsvd = select(obsvd[., 1..2], obsvd[.,3]:==1)
		
			mata: obsvd = select((obsvd \ J(rows(retrace)-rows(obsvd), 2, .)), retrace)
		
			mata: st_addobs(rows(obsvd))
		
			mata: st_store(count+1::count+rows(obsvd), (st_varindex("x"), st_varindex("y")), obsvd)
		
			mata: count = count + rows(obsvd)
		
			mata: currentvec = select(currentvec, abs(retrace:-1))
		
			mata: currentvec = select(currentvec, currentvec:!="")
		
			mata: st_local("row", strofreal(rows(currentvec)))
		
			mata: row = rows(obsvd)+1
		
			forvalues x = 1/`row' {
		
				mata: st_addobs(1)
			
				mata: st_store(++count, st_varindex("x"), column)
			
				mata: st_store(count, st_varindex("y"), row)
			
				mata: st_sstore(count, st_varindex("label"), namemat[row++, column])
		
			}
		
			quietly replace col = `y' if missing(col)
		
			quietly replace down = 1 if missing(down)
		
			mata: countback = st_data(., (st_varindex("x"), st_varindex("y"), st_varindex("col"), st_varindex("down")))
		
			mata: countback = select(countback[., (1..2,4)], countback[.,3]:==column)
		
			mata: countback = select(countback[., 1..2], countback[.,3]:==1)
		
			mata: countback = sort(((rows(countback)::1), countback), 1)
		
			mata: st_addobs(rows(countback))
		
			mata: st_store(count+1::count+rows(countback), (st_varindex("x"), st_varindex("y")), countback[., 2..cols(countback)])
		
			mata: count = count + rows(countback)
		
			quietly replace down = 0 if missing(down)
		
			quietly replace col = `y' if missing(col)
		
			mata: column++
	
		}
	
	}

	quietly summarize y

	quietly replace y = (y - r(min))*-1 + r(max)

	tempname max

	scalar `max' = r(max)

	quietly summarize x

	quietly replace x = `=r(max)/2+r(min)' if y == `max'

	quietly generate int group = .

	forvalues x = `=`max'-1'(-1)1 {

		quietly replace group = y == `x' & label != ""
	
		quietly replace group = group + group[_n-1] in 2/`=_N'
	
		quietly summarize group
	
		forvalues y = `=r(min)'/`=r(max)' {
	
			quietly summarize x if group == `y' & y <= `x'
		
			quietly replace x = `=r(max)/2+r(min)' if group == `y' & y == `x'
	
		}
	
		quietly replace group = .

	}
	
	graph twoway scatter y x, mlabel(lab) connect(direct)

	capture matrix list e(importance)

	if !_rc {

		display _newline "{res}Splitting variable permutation importance"

		matrix list e(importance), noheader
	
	}
	
	restore
	
}

end

/*mata function to recursively partition data based on CHAID algorithm*/
version 12.1

mata:

mata set matastrict on

string matrix chaid(string rowvector unordvar, string rowvector ordvar, string scalar type, ///
string scalar names, string scalar vars, real scalar alpha1, real scalar alpha3, ///
real scalar minsplit, real scalar minnode, real scalar maxbranch, real scalar alpha2, ///
string scalar bfadj, string scalar svy, string scalar exhaust)
{
	real matrix selectvar, select2var, countcheck, comparecount
	
	real rowvector var1, var2, merger_ps, merger_aic, withinselect, randresolve, ///
	branches
	
	real colvector unordpcomp
	
	real scalar position, moresplit, compare_p, current_p, compare_aic, current_aic, ///
	moremerge, clustercount, currentcluster, adj_p, countmerge, temploc1, temploc2, ///
	mergeruns, floating, match_p
	
	string matrix splits, CHAIDsplit
	
	string rowvector varsvec, namesvec, typevec, selecttype, collapsevec, selectnames, ///
	clusnames, comparenames, select2names, select2vec
	
	string scalar temploc1s, temploc2s
	
	transmorphic basis																				
	
	tokens(vars)		//																																						
	
	position = 2																						//updatable position of data to find in varsvec
	
	varsvec = tokens(vars)																				//vectorize the tempvars																		
	
	varsvec[position]	//
	
	moresplit = 1																						//scalar to indicate continue splitting process between-variable process
	
	compare_p = 1																						//currently lowest p-value for deciding which variable will be splitting variable
	
	compare_aic = .																						//currently lowesr AIC value for deciding which variable will be splitting variable
	
	clustercount = 1
	
	currentcluster = 1
	
	comparenames = ("")
	
	branches = (0)
	
	splits = ("")
	
	q = tokeninit(wchars = (" "), pchars = (" "), qchars = ("()"))										//vectorize the display names
		
	tokenset(q, names)
	
	namesvec = tokengetall(q)
	
	namesvec	//
	
	q = tokeninit(wchars = (" "), pchars = (" "), qchars = ("()"))										//vectorize the type inputs
		
	tokenset(q, type)
	
	typevec = tokengetall(q)
	
	typevec	//
	
	CHAIDsplit = ("Splits")
	
	while (moresplit == 1) {																			//do following process at least once - continue if between-variable splitting process is to continue (i.e., moresplit == 1)
	
		stata("replace \`useob' = \`stop' & \`cluster' == " + strofreal(currentcluster) + " if \`touse'")	
		
		selectnames = tokens(subinstr(subinstr(namesvec[position], ")", ""), "(", ""))					//pull out names for first independent variable (i.e., in position = 2] binaries using same procedure as with dependent variable - select() to pull out data - strpos() to create	vector of 1's, substr() to truncate names, strlen() to truncate substr() statement at length of focal independent variable to match												
		
		selectnames	//
		
		selecttype = tokens(subinstr(subinstr(typevec[position], ")", ""), "(", ""))					//pull out types for first independent variable using same procedure as with dependent variable - select() to pull out data - strpos() to create	vector of 1's, substr() to truncate names, strlen() to truncate substr() statement at length of focal independent variable to match												
		
		selecttype	//
		
		moremerge = 1																					//updatable scalar to indicate continute merging within-variable process
		
		stata("tabulate \`dv' \`useob' if \`touse' [\`weight'\`exp'], matcell(\`cellcount')") 			//tabulate dv with the current cluster set, keep the resulting matrix
		
		stata("ereturn matrix cellcount = \`cellcount'")												//ereturn the matrix to pull it into Mata - probably unnecessary, but it works
		
		countcheck = st_matrix("e(cellcount)")'															//pull matrix of counts into Mata, transposed
		
		countcheck	//
		
		if (rows(countcheck) > 1) countcheck = countcheck[2,]											//keep only the rows corresponding to the "useob" values if there are some that got a 0 on useob
		
		countcheck = select(countcheck, sign(countcheck))												//keep only the rows with some value (omit thosw with 0's)
		
		countcheck	//
		
		stata("summarize \`stop' [\`weight'\`exp']")													//quick summarize on stopping indicator
		
		if (st_numscalar("r(mean)") == 0) {																//if no-one has a 1 on stop, then end chaid as the process is now over
		
			moresplit = 0																				//don't split anymore - chaid is done
			
			moremerge = 0																				//don't merge anymore on this variable
			
		}
		
		else if ((rowsum(sign(countcheck)) == 1) | (rowsum(countcheck) <= minsplit)) {					//if there is only one value left in a cluster (i.e., it's pure) or it's size is below the minimum splitting size
			
			moremerge = 0																				//stop merging on this variable
			
			stata("replace \`stop' = 0 if \`useob' & \`touse'")											//update cluster to stop using this variable
			
			currentcluster++																			//increment current cluster - start using the next cluster up, look for splits there
			
			position = 2																				//restart "position" used in merging
			
		}
																										
		/*begin splitting variable category merging*/
		else {
		
			stata("tabulate " + varsvec[position] + ///
			" \`useob' if \`touse' [\`weight'\`exp'], \`missing' matcell(\`cellcount')") 
		
			stata("ereturn matrix cellcount = \`cellcount'")
		
			countcheck = st_matrix("e(cellcount)")'
		
			countcheck	//
		
			if (rows(countcheck) > 1) countcheck = countcheck[2,]
		
			selectnames = select(selectnames, sign(countcheck))
		
			selectnames	//
		
			while (moremerge == 1) {
		
				collapsevec = ("")																			//replace/generate collapsevec including names of variables to be collapsed into single variable
			
				merger_ps = (.)
				
				merger_aic = (.)
				
				floating = 0
			
				if (cols(selectnames) > 2)  {																//if there are more than 2 categories and the type of variabel is "unordered" then proceed
		
					/*Unordered merging*/
					if (selecttype[1] == "u") {
										
						base = (J(2, 1, 1) \ J(cols(selectnames) - 2, 1, 0))								//setup vector for use in cvpermute() to select all sets of 2 variables
			
						base	//
						
						basis = cvpermutesetup(base)														//register "base" vector with cvpermute()
				
						mergeruns = comb(cols(selectnames), 2)
				
					}
					
					else if (
					(selecttype[cols(selecttype)] == "f") & ///
					(colsum(regexm(selectnames[cols(selectnames)], "ms$")) == 1) & ///
					(colsum(cols(tokens(selectnames[cols(selectnames)]))) == 1)) {
					
						mergeruns = (cols(selectnames) - 1)*2 - 1
						
						floating = 1
						
						base = (J(1, 1, 1) \ J(cols(selectnames) - 2, 1, 0))								//setup vector for use in cvpermute() to select all sets of 2 variables
			
						base	//
						
						basis = cvpermutesetup(base)
					
					}
					
					else {
					
						mergeruns = cols(selectnames) - 1
					
					}
				
					countmerge = 1
				
					while (countmerge <= mergeruns) {														//go through all different combinations of predictor levels to find best levels to collapse
			
						if (selecttype[1] == "u") select2var = cvpermute(basis)'							//invoke cvpermute for selector vector
				
						else if (((selecttype[1] == "o") & (floating == 0)) | ///
						((floating == 1) & (sign(countmerge*2 - mergeruns) != 1))) {
						
							select2var = (J(countmerge - 1, 1, 0) \ J(2, 1, 1) \ ///
							J(cols(selectnames) - countmerge - 1, 1, 0))'
						
						}
						
						else if ((floating == 1) & (sign(countmerge*2 - mergeruns) == 1)) {
						
							select2var = (cvpermute(basis)', 1)
									
						}
				
						select2var	//
				
						collapsevec = (collapsevec, invtokens(strofreal(select2var)))						//create string vector indicating position of variabels to test
				
						collapsevec	//
				
						select2names = select(selectnames, select2var)										//proceed to select 2 columns of predictor to evaluate for potentially collapsing
				
						select2names	//
					
						select2names[1]	//
					
						select2vec = tokens(select2names[1])
					
						temploc1 = st_addvar("byte", st_tempname())
					
						temploc1s = st_varname(temploc1)
						
						stata("replace " + temploc1s + "=0")
					
						for (y = 1; y <= cols(tokens(select2names[1])); y++) {
						
							stata("replace " + temploc1s + "=`" ///
							+ select2vec[y] + "' + " + temploc1s)
								
						}
					
						select2names[2]	//
					
						select2vec = tokens(select2names[2])
					
						temploc2 = st_addvar("byte", st_tempname())
					
						temploc2s = st_varname(temploc2)
						
						stata("replace " + temploc2s + "=0")
					
						for (y = 1; y <= cols(tokens(select2names[2])); y++) {
						
							stata("replace " + temploc2s + "=`" ///
							+ select2vec[y] + "' + " + temploc2s)
								
						}
					
						stata("summarize " + temploc1s + " " + temploc2s + " if \`touse'" ///
						+ " & \`useob' [\`weight'\`exp']")	//
						 
						st_macroexpand("\`permute'version 10: \`reg' \`dv' " + temploc1s + ///
						" if \`touse' & \`useob' & (" + temploc1s + " | " + temploc2s + ") [\`weight'\`exp']")	//
						
						if (strlen(svy)) stata("capture version 10: svy, subpop(if \`touse' & \`useob' & (" + temploc1s + ///
						" | " + temploc2s + ")): \`reg' \`dv' " + temploc1s)
						
						else stata("capture \`permute'version 10: \`reg' \`dv' " + temploc1s + ///
							" if \`touse' & \`useob' & (" + temploc1s + " | " + temploc2s + ") [\`weight'\`exp']")
						
						if ((st_numscalar("c(rc)") == 0)) {
						
							if ((strlen(svy)) & !((st_numscalar(st_local("test")) == .) | (st_numscalar("df_m") == 0))) {
						
								stata("\`reg'")	//
																
								merger_ps = (merger_ps, Ftail(st_numscalar("e(df_m)"), ///
								st_numscalar("e(df_r)"), st_numscalar(st_local("test"))))
								
								stata("getic, type(aicw) ll(1) parm(1) obs(1)")
								
								merger_aic = (merger_aic, st_numscalar("r(ic)"))

						
							}
							
							else if (!(strlen(svy)) & !((st_numscalar(st_local("test")) == .) | (st_numscalar("df_m") == 0)) & ///
							!(rows(st_matrix("r(p)")))) {
					
								stata("\`reg'")	//
								
								merger_ps = (merger_ps, chi2tail(st_numscalar("e(df_m)"), ///
								st_numscalar(st_local("test"))))
								
								stata("getic, type(aicc) ll(\`=e(ll)') parm(\`=e(df_m)') obs(\`=e(N)')")
								
								merger_aic = (merger_aic, st_numscalar("r(ic)"))
						
							}
							
														
							else if (!(strlen(svy)) & (rows(st_matrix("r(p)")))) {
							
								merger_ps = (merger_ps, st_matrix("r(p)")[1,1])
								
								stata("getic, type(aicc) ll(" + strofreal(st_matrix("r(b)")[1, 2]) + ///
								") parm(" + strofreal(st_matrix("r(b)")[1, 3]) + ") obs(" ///
								+ strofreal(st_matrix("r(b)")[1, 4]) + ")")
								
								merger_aic = (merger_aic, st_numscalar("r(ic)"))
								
							}
					
							else { 
								
								merger_ps = (merger_ps, 1)
								
								merger_aic = (merger_aic, .)
								
							}
							
						}
						
						else {
						
							merger_ps = (merger_ps, 1)
							
							merger_aic = (merger_aic, .)
							
						}
						
						countmerge++
					
						countmerge	//
						
						st_dropvar((temploc1s, temploc2s))
				
					}
			
					if ((max(merger_ps) > 1 - alpha1) | (strlen(exhaust))) {								//if at least one p-value doesn't the "respalpha" limit offered by the user, collapse the levels as they appear sufficiently homogeneous - also ensures that, when the DV is a constant, the algorithm stops
					
						merger_ps = strmatch(strofreal(merger_ps), strofreal(max(merger_ps)))				//find the position of the largest p-value that passed the spltalpha criterion
				
						merger_ps	//
					
						if (rowsum(merger_ps) > 1) {
						
							merger_ps = (0, strmatch(strofreal(merger_aic[2..cols(merger_aic)]), ///
							strofreal(max(merger_aic[2..cols(merger_aic)]))))								//find the position of the largest aic value
				
							merger_ps	//
											
							if (rowsum(merger_ps) > 1) {
							
								randresolve = runiform(1, cols(merger_ps))
					
								merger_ps = strmatch(strofreal(merger_ps), strofreal(max(merger_ps)))
						
								merger_ps = randresolve:*merger_ps
						
								merger_ps = strmatch(strofreal(merger_ps), strofreal(max(merger_ps)))
								
								merger_ps	//
								
							}
					
						}
				
						withinselect = strtoreal(tokens(select(collapsevec, merger_ps)))					//pull out variables associated with largest p-value
						
						withinselect	//
						
						if (selecttype[1] == "u") selectnames = (invtokens(select(selectnames, withinselect)), ///
						select(selectnames, abs(withinselect:-1)))
						
						else if (selecttype[1] == "o") {
						
							countmerge = (strpos(invtokens(strofreal(withinselect)), "1") + 1)/2
							
							countmerge	//
							
							if (withinselect[countmerge] == withinselect[countmerge + 1]) { 
						
								if (countmerge == 1) selectnames = (invtokens(select(selectnames, withinselect)), ///
								select(selectnames[countmerge + 2..cols(withinselect)], ///
								abs(withinselect[countmerge + 2..cols(withinselect)]:-1)))
							
								else if (countmerge == cols(withinselect) - 1) selectnames = ///
								(select(selectnames[1..countmerge - 1], abs(withinselect[1..countmerge - 1]:-1)), ///
								invtokens(select(selectnames, withinselect)))
							
								else selectnames = (select(selectnames[1..countmerge - 1], abs(withinselect[1..countmerge - 1]:-1)), ///
								invtokens(select(selectnames, withinselect)), ///
								select(selectnames[countmerge + 2..cols(withinselect)], ///
								abs(withinselect[countmerge + 2..cols(withinselect)]:-1)))
							
							}
							
							else {
							
								if (countmerge == 1) selectnames = (invtokens(select(selectnames, withinselect)), ///
								select(selectnames[countmerge + 1..cols(withinselect) - 1], ///
								abs(withinselect[countmerge + 1..cols(withinselect) - 1]:-1)))
							
								else if (countmerge == cols(withinselect) - 1) selectnames = ///
								(select(selectnames[1..countmerge - 1], abs(withinselect[1..countmerge - 1]:-1)), ///
								invtokens(select(selectnames, withinselect)))
							
								else selectnames = (select(selectnames[1..countmerge - 1], abs(withinselect[1..countmerge - 1]:-1)), ///
								invtokens(select(selectnames, withinselect)), ///
								select(selectnames[countmerge + 1..cols(withinselect) - 1], ///
								abs(withinselect[countmerge + 1..cols(withinselect) - 1]:-1)))
							
							}
						
						}
					
						selectnames	//
						
						if (alpha2 > 0) {
						
							if (strlen(exhaust)) match_p = 1 - max(merger_ps)
							
							else match_p = -1
						
							selectnames = chaidsplit(selectnames, withinselect, selecttype[1], alpha2, svy, match_p)
							
							selectnames	//
							
						}
				
					}
			
					else moremerge = 0
				
				}
			
				else moremerge = 0
				
				adj_p = 0
						
				if (!strlen(bfadj)) {
					
					if (selecttype[1] == "u") {
						
						if (strlen(exhaust)) adj_p = cols(selectnames)*(cols(selectnames) - 1)/2
								
						else {
						
							unordpcomp = (0::cols(selectnames))
					
							for (y = 1; y <= cols(selectnames) - 1; y++) {
					
								adj_p = adj_p + (-1)^unordpcomp[y]*((cols(selectnames) - ///
								unordpcomp[y])^cols(selecttype)/(factorial(unordpcomp[y])* ///
								factorial(cols(selectnames) - unordpcomp[y])))
										
							}
								
						}
																						
					}
						
					else if (selecttype[1] == "o") {
						
						if (strlen(exhaust)) adj_p = cols(selectnames)*(cols(selectnames)^2 - 1)/2
						
						else adj_p = comb(cols(selecttype) - 1, cols(selectnames) - 1)
						
					}
					
					else if (selecttype[cols(selecttype)] == "f") {
						
						if (strlen(exhaust)) adj_p = cols(selectnames)*(cols(selectnames) - 1)/2
							
						else {
						
							comb(cols(selecttype) - 1, cols(selectnames) - 1)* ///
							((cols(selectnames) - 1) + cols(selectnames)*(cols(selecttype) - ///
							cols(selectnames)))/(cols(selecttype) - 1)
									
						}
							
					}
							
				}
						
				else adj_p = 1
						
				adj_p	//
			
				"Completed One Merging Session"	//
			
			}
		
			/*Compute chi-square used for splitting*/
		
			countmerge = 1
		
			temploc1 = st_addvar("byte", st_tempname())
		
			temploc1s = st_varname(temploc1)
						
			stata("replace " + temploc1s + "=0 if \`touse' & \`useob'")
			
			for (x = 1; x <= cols(selectnames); x++) {
			
				select2vec = tokens(selectnames[x])					
					
				for (y = 1; y <= cols(select2vec); y++) {
					
					st_macroexpand("replace " + temploc1s + "=`" ///
					+ select2vec[y] + "'*" + strofreal(countmerge) + " if `" ///
					+ select2vec[y] + "' & \`touse' & \`useob'")
						
					stata("replace " + temploc1s + "=`" ///
					+ select2vec[y] + "'*" + strofreal(countmerge) + " if `" ///
					+ select2vec[y] + "' & \`touse' & \`useob'")
						
				}
				
				countmerge++
				
				stata("summarize " + temploc1s + " if \`touse' & \`useob' [\`weight'\`exp']")	//
				
			} 
				
			stata("tabulate " + temploc1s + ///
			" \`useob' if \`touse' [\`weight'\`exp'], \`missing' matcell(\`cellcount')") 
		
			stata("ereturn matrix cellcount = \`cellcount'")
		
			countcheck = st_matrix("e(cellcount)")'
		
			countcheck	//
			
			if (selecttype[1] == "u") {
			
				st_local("xi", "xi:")
				
				st_local("i", "i.")
			
			}
			
			if (selecttype[1] == "o") {
			
				st_local("xi", "")
				
				st_local("i", "")
			
			}
			
			if (strlen(svy)) stata("capture version 10:svy, subpop(if \`touse' & \`useob'):\`xi' \`reg' \`dv' \`i'" + temploc1s)
			
			else stata("capture \`permute'version 10:\`xi' \`reg' \`dv' \`i'" + temploc1s + " if \`touse' & \`useob' [\`weight'\`exp']")
			
			if ((st_numscalar("c(rc)") == 0) & (cols(countcheck) > 1)) {
			
				if ((strlen(svy)) & !((st_numscalar(st_local("test")) == .) | (st_numscalar("df_m") == 0))) {
			
					stata("\`reg'")	//
				
					current_p = Ftail(st_numscalar("e(df_m)"), st_numscalar("e(df_r)"), st_numscalar(st_local("test")))
					
					stata("getic, type(aicw) ll(1) parm(1) obs(1)")
								
					current_aic = st_numscalar("r(ic)")
			
				}
			
				else if (!(strlen(svy)) & !((st_numscalar(st_local("test")) == .) | (st_numscalar("df_m") == 0)) ///
				& !(rows(st_matrix("r(p)")))) {
			
					stata("\`reg'")	//
			
					current_p = chi2tail(st_numscalar("e(df_m)"), st_numscalar(st_local("test")))
					
					stata("getic, type(aicc) ll(\`=e(ll)') parm(\`=e(df_m)') obs(\`=e(N)')")
								
					current_aic = st_numscalar("r(ic)")
				
				}
			
				else if (!(strlen(svy)) & (rows(st_matrix("r(p)")))) {
				
					current_p = st_matrix("r(p)")[1,1]
					
					stata("getic, type(aicc) ll(" + strofreal(st_matrix("r(b)")[1, 2]) + ///
					") parm(" + strofreal(st_matrix("r(b)")[1, 3]) + ") obs(" ///
					+ strofreal(st_matrix("r(b)")[1, 4]) + ")")
								
					current_aic = st_numscalar("r(ic)")
					
				}
			
				else {
				
					current_p = 1
					
					current_aic = .
					
				}
				
			}
			
			else {
			
				current_p = 1
				
				current_aic = .
				
			}
				
			st_dropvar(temploc1s)
			
			current_p = current_p*adj_p
				
			current_p	//
			
			current_aic	//
		
			compare_p	//
			
			compare_aic	//
		
			if (current_p < compare_p) {
		
				"better!"	//
		
				compare_p = current_p
		
				compare_p	//
				
				compare_aic = current_aic
				
				current_aic	//
			
				comparenames = selectnames 
			
				comparenames	//
				
				comparecount = countcheck
				
				comparecount	//
			
			}
		
			else if (current_p == compare_p) {	
		
				if (current_aic < compare_aic) {
				
					"better!"	//
		
					compare_p = current_p
		
					compare_p	//
					
					compare_aic = current_aic
				
					current_aic	//
			
					comparenames = selectnames 
			
					comparenames	//
				
					comparecount = countcheck
				
					comparecount	//
					
				}
					
				else {	
				
					randresolve = runiform(1, 2)
			
					if (randresolve[1] < randresolve[2]) {
		
						"better!"	//
		
						compare_p = current_p
		
						compare_p	//
						
						compare_aic = current_aic
				
						current_aic	//
			
						comparenames = selectnames 
			
						comparenames	//
						
						comparecount = countcheck
						
						comparecount	//
						
					}
						
				}
		
			}
		
			position =  position + 1
		
			position	//
		
			"Completed Merging Run"	//
		
			comparenames	//
	
			if (position == cols(varsvec) + 1) {
		
				if ((cols(tokens(invtokens(comparenames))) > 1) & (compare_p < alpha3) & ///
				(min(comparecount) >= minnode) & ///
				((branches[currentcluster] < maxbranch) | (maxbranch < 0))){
		
					"split!"	//
					
					branches[currentcluster] = branches[currentcluster] + 1
					
					branches //
					
					selectnames = (" ")
					
					for (x = 1; x <= cols(comparenames); x++) {
					
						comparenames[x]	// 
					
						selectnames = (selectnames, comparenames[x], ",")
					
					}
					
					selectnames	//
					
					if (branches[currentcluster] > cols(splits)) splits = (splits, J(rows(splits), 1, ""))
					
					splits[currentcluster, branches[currentcluster]] = comparenames[1]
					
					splits	//
					
					CHAIDsplit = (CHAIDsplit, invtokens(selectnames))
					
					for (x = 2; x <= cols(comparenames); x++) {
					
						clustercount++
						
						clustercount	//
						
						branches = (branches, branches[currentcluster])
						
						branches	//
						
						splits = (splits \ splits[currentcluster, .])
						
						splits[clustercount, branches[clustercount]] = comparenames[x]
						
						splits	//
																		
						selectnames = tokens(comparenames[x])
					
						for (y = 1; y <= cols(selectnames); y++) {
							
							st_macroexpand("replace \`cluster' = " + strofreal(clustercount)+ " if `" + ///
							selectnames[y] + "' & \`touse' & \`useob'")	//
				
							stata("replace \`cluster' = " + strofreal(clustercount)+ " if `" + ///
							selectnames[y] + "' & \`touse' & \`useob'")	//	
					
						}
					
					}
					
					position = 2
					
					comparenames = ("")
					
					compare_p = 1
					
					compare_aic = .
				
				}
				
				else {
				
					position = 2
				
					stata("replace \`stop'=0 if \`cluster' == " + strofreal(currentcluster))
					
					currentcluster++
				
					stata("summarize \`cluster'")
				
					if (st_numscalar("r(max)") < currentcluster) moresplit = 0
				
					currentcluster	//
				
					compare_p = 1
					
					compare_aic = .
					
					comparenames = ("")
						
				}
		
			}	
		
		}
	
	}
	
	currentcluster = 0
	
	for (x = 1; x <= rows(splits); x++) {
	
		currentcluster++
		
		if (cols(CHAIDsplit) > 1) CHAIDsplit = (CHAIDsplit \ ("path" + strofreal(currentcluster), splits[x, .], J(1, cols(CHAIDsplit)-cols(splits)-1, "")))
	
	}
	
	branchmat = st_tempname()
	
	st_matrix(branchmat, branches)
	
	st_local("branchmat", branchmat)
	
	"Finished Execution"	//
	
	return(CHAIDsplit)

}  

end

version 12.1

mata:

mata set matastrict on

string rowvector chaidsplit(string rowvector selectnames, real rowvector withinselect, string scalar type, ///
real scalar cutalpha, string scalar svy, real scalar match_p)
{
	string scalar trysplit, tploc1s, tploc2s
	
	real scalar count, float_miss
	
	real rowvector selectsplit, tploc1, tploc2, mg_ps, mg_aic, repeats
	
	real colvector  base, combin
	
	string rowvector splitvec, namesplit, comparevec
	
	real matrix unordmat
	
	transmorphic basis
	
	"chaidsplit begins"	//

	if (type == "o") {
	
		trysplit = select(selectnames, withinselect[1..cols(withinselect) - 1])
		
		trysplit = trysplit[1]
		
		float_miss = rowsum(strmatch(substr(tokens(trysplit), -2), "ms"))
		
	}
	
	else if (type == "u") trysplit = selectnames[1]
		
	trysplit	//
	
	type	//
	
	splitvec = ("")
	
	mg_ps = (.)
	
	mg_aic = (.)
		
	if (cols(tokens(trysplit)) > 2) {
	
		namesplit = tokens(invtokens(trysplit))
		
		namesplit	//
		
		if (float_miss == 1) {
		
			namesplit = (select(namesplit, abs(strmatch(substr(namesplit, -2), "ms")):-1), ///
			select(namesplit, strmatch(substr(namesplit, -2), "ms")))
			
			namesplit	//
		
		}
		
		if (type == "u") {
		
			repeats = 2^(cols(tokens(trysplit))) - 2
			
			unordmat = J(1, cols(tokens(trysplit)), .)
			
			for (x = 1; x <= cols(tokens(trysplit)); x++) {

				base = J(x, 1, 1)
	
				base = (base \ J(cols(tokens(trysplit)) - x, 1, 0))
	
				basis = cvpermutesetup(base)
	
				for (y = 1; y <= comb(cols(tokens(trysplit)), x); y++) {
	
					combin = cvpermute(basis)
		
					unordmat = (unordmat \ combin')
		
				}
	
			}
			
			unordmat = unordmat[2..rows(unordmat)-1, .]
			
			unordmat	//
			
		}
		
		else if (type == "o") {
		
			repeats = cols(tokens(trysplit)) - 1
			
			if (float_miss == 1) repeats = repeats*2 - 1
			
		}
		
		repeats	//
		
		count = 1
		
		while (count <= repeats) {
			
			if (type == "u") {
				
				selectsplit = unordmat[count, .]
				
			}
			
			else if (type == "o") {
			
				if (count <= cols(namesplit) - 1) selectsplit = (J(1, count, 1), J(1, cols(namesplit) - count, 0))
				
				else if (count > cols(namesplit) - 1) {
				
					selectsplit = (J(1, count - (cols(namesplit) - 1), 1), ///
					J(1, cols(namesplit) - 1 - (count - (cols(namesplit) - 1)), 0), 1)
				
				}
			
			}
	
			selectsplit	//
				
			splitvec = (splitvec, invtokens(strofreal(selectsplit)))						//create string vector indicating position of variables to test
				
			splitvec	//
				
			comparevec = select(namesplit, selectsplit)										//proceed to select 2 columns of predictor to evaluate for potentially collapsing
				
			comparevec	//
				
			tploc1 = st_addvar("byte", st_tempname())
					
			tploc1s = st_varname(tploc1)
						
			stata("replace " + tploc1s + "=0")
					
			for (y = 1; y <= cols(comparevec); y++) {
						
				stata("replace " + tploc1s + "=`" ///
				+ comparevec[y] + "' + " + tploc1s)
								
			}
					
			comparevec = select(namesplit, abs(selectsplit:-1))
			
			comparevec	//
					
			tploc2 = st_addvar("byte", st_tempname())
					
			tploc2s = st_varname(tploc2)
						
			stata("replace " + tploc2s + "=0")
					
			for (y = 1; y <= cols(comparevec); y++) {
						
				stata("replace " + tploc2s + "=`" ///
				+ comparevec[y] + "' + " + tploc2s)
								
			}
					
			stata("summarize " + tploc1s + " " + tploc2s + " if \`touse'" ///
			+ " & \`useob' [\`weight'\`exp']")	//
						 
			st_macroexpand("\`permute'version 10: \`reg' \`dv' " + tploc1s + ///
			" if \`touse' & \`useob' & (" + tploc1s + " | " + tploc2s + ") [\`weight'\`exp']")	//
			
			if (strlen(svy)) stata("capture version 10: svy, subpop(if \`touse' & \`useob' & (" + tploc1s + " | " + ///
			tploc2s + "): \`reg' \`dv' " + tploc1s)
			
			else stata("capture \`permute'version 10: \`reg' \`dv' " + tploc1s + ///
			" if \`touse' & \`useob' & (" + tploc1s + " | " + tploc2s + ") [\`weight'\`exp']")
						
			if ((st_numscalar("c(rc)") == 0)) {
						
				if ((strlen(svy)) & !((st_numscalar(st_local("test")) == .) | (st_numscalar("df_m") == 0))) {
						
					stata("\`reg'")	//
							
					mg_ps = (mg_ps, Ftail(st_numscalar("e(df_m)"), ///
					st_numscalar("e(df_r)"), st_numscalar(st_local("test"))))
					
					stata("getic, type(aicw) ll(1) parm(1) obs(1)")
								
					mg_aic = (mg_aic, st_numscalar("r(ic)"))
						
				}
							
				else if (!(strlen(svy)) & !((st_numscalar(st_local("test")) == .) | (st_numscalar("df_m") == 0)) & ///
				!(rows(st_matrix("r(p)")))) {
					
					stata("\`reg'")	//
						
					mg_ps = (mg_ps, chi2tail(st_numscalar("e(df_m)"), ///
					st_numscalar(st_local("test"))))
					
					stata("getic, type(aicc) ll(\`=e(ll)') parm(\`=e(df_m)') obs(\`=e(N)')")
								
					mg_aic = (mg_aic, st_numscalar("r(ic)"))
						
				}
							
											
				else if (!(strlen(svy)) & (rows(st_matrix("r(p)")))) {
				
					mg_ps = (mg_ps, st_matrix("r(p)")[1,1])
					
					stata("getic, type(aicc) ll(" + strofreal(st_matrix("r(b)")[1, 2]) + ///
					") parm(" + strofreal(st_matrix("r(b)")[1, 3]) + ") obs(" ///
					+ strofreal(st_matrix("r(b)")[1, 4]) + ")")
								
					mg_aic = (mg_aic, st_numscalar("r(ic)"))
					
				}
					
				else {
				
					mg_ps = (mg_ps, 1)
					
					mg_aic = (mg_aic, .)
					
				}
							
			}
						
			else {
			
				mg_ps = (mg_ps, 1)
				
				mg_aic = (mg_aic, .)
				
			}
						
			st_dropvar((tploc1s, tploc2s))
			
			count++
				
		}
		
		if ((min(mg_ps) < cutalpha) & (min(mg_ps) < match_p)) {
			
			mg_ps = strpos(strofreal(min(mg_ps[2..cols(mg_ps)])), strofreal(mg_ps[2..cols(mg_ps)]))
			
			mg_ps	//
				
			if (rowsum(mg_ps) > 1) {
			
				mg_ps = strpos(strofreal(min(mg_aic[2..cols(mg_aic)])), strofreal(mg_aic[2..cols(mg_aic)]))
			
				mg_ps	//
				
				if (rowsum(mg_ps) > 1) {
			
					randresolve = runiform(1, cols(mg_ps))
				
					mg_ps = strmatch(strofreal(mg_ps), strofreal(max(mg_ps)))
						
					mg_ps = randresolve:*mg_ps
						
					mg_ps = strmatch(strofreal(mg_ps), strofreal(max(mg_ps)))
								
					mg_ps	//
					
				}
	
			}
				
			trysplit = (invtokens(select(namesplit, strtoreal(tokens(select(splitvec[2..cols(splitvec)], mg_ps))))), ///
			invtokens(select(namesplit, abs(strtoreal(tokens(select(splitvec[2..cols(splitvec)], mg_ps))):-1))))
			
			trysplit	//
		
			if (type == "u") selectnames = (trysplit, selectnames[2..cols(selectnames)])
		
			else if (type == "o") {
		
				if (withinselect[1] == 1) selectnames = (trysplit, selectnames[2..cols(selectnames)])
			
				else if (withinselect[cols(withinselect)-1] == 1) selectnames = (selectnames[1..cols(selectnames)-1], trysplit)
			
				else {
			
					count = (strpos(invtokens(strofreal(withinselect)), "1") + 1)/2
					
					count	//
				
					selectnames = (selectnames[1..count-1], trysplit, selectnames[count..cols(selectnames)])
			
				}
				
				if (rowsum(strmatch(substr(selectnames, -2), "ms")) == 1) {
				
					if (cols(tokens((select(selectnames, strmatch(substr(selectnames, -2), "ms")))))) {
					
						selectnames = (select(selectnames, abs(strmatch(substr(selectnames, -2), "ms")):-1), ///
						select(selectnames, strmatch(substr(selectnames, -2), "ms")))
					
					}
					
				}				
				
			}
			
		}
		
	}
	
	"end chaidsplit"	//
	
	return(selectnames)

}

end

//program to obtain information criteria
program define getic, rclass

syntax, type(string) ll(real) parm(integer) obs(integer)

if ("`type'" == "aicc") return scalar ic = -2*`ll'+2*`parm'+(2*`parm'*`parm' + 1)/(`obs' - `parm' - 1)

else if ("`type'" == "aicw") {

	tempname ilog gdeff
	
	matrix `ilog' = e(ilog)
	
	mata: ilog = select(st_matrix("`ilog'"), abs(sign(st_matrix("`ilog'"))))
	
	mata: st_matrix("`ilog'", ilog[1, cols(ilog)])
	
	matrix `gdeff' = trace(invsym(e(V_srs))*e(V))
	
	return scalar ic = -2*`ilog'[1,1]*(e(N)/e(N_pop))+2*`gdeff'[1,1]

}

else exit 198

end

/* programming notes and history

- chaid version 1.0 - date - November 18, 2013

Basic version

-----

- chaid version 1.1 - date - March 20, 2014

//notable changes\\
a] added frequency weights
b] updated display to highlight hierarchical tree-structure using Mata matrix and fix errors with original display
c] implemented the remerge process as an option
d] allowed a no Bonferroni adjustment option
e] fixed merging inaccuracy that excluded categories arbitrarily
f] error returning an extra "root" fixed
g] predicted value offered based on mode of node
h] additional returned values
i] updated helpfile - going over collapse option to expedite chaid fitting

-----

- chaid version 2.0 - date - August 22, 2014

//notable changes\\
a] added permutation-based p-values
b] added svy prefix compatability
c] added graphical depiction to display output
d] added accessible mata matrices
e] now can "replay" chaid results again as with other Stata programs
f] added permutation importance metric
g] fixed N-obs issue w/ fweights
h] added fit metric based on cramer's V (and svy compatible version with similar properties)
i] fixed Q1 Q11 nested variable selection issue - also allows to predict itself (dv and iv, though non-sensical, no more error)
j] added exhaustive CHAID algorithm
k] convenience command allows user to "xtile" for continuous predictors

-----

- chaid version 2.1 - date - December 16, 2014

//notable changes\\
a] fixed missing colon for svy with respalpha
b] moved location of bonferroni adjustment (doesn't require "exhaust" or there to have been a merge implemented to compute it) 
	to ensure no missing adjusted p-value
c] use aic for between model comparisons to decide on splitting (p-value has numerical precision issues... just use it to meet 
	the p<.05 threshold and that's it)
d] fix importance w/ missing (extra space added to "ifstmnts" for missings erroneously that was turned into comma
e] fix DV levels check that occurs before markout

-----

- chaid version 2.2 - date - February 12, 2015

//notable changes\\
a] fixed Bonferroni adjustment (was not invoked unless noadj option was used) - had opposite behavior as intended
b] fixed but which prevented splits on splitting variables with a single level

v 3?
- allow custom xtile() for each variable
- save info for all splits (as obj) 
- postestimation for predict & importance
- make interactive fitting - graphing?
- add mi?