*! chaidforest_estat - version 1.0 - 10/12/2015 - Joseph N. Luchman

program chaidforest_estat, eclass																			//program history and notes at end of file

version 12.1

syntax anything [, TRee(integer -1) INSample GRaph mata NOIsily]

if ("`anything'" == "prox") {																				//if a proximity matrix is desired...

	if strlen("`mata'") {																					//if the matrix should be saved to Mata...

		mata: prox = proximity(results, `=e(ntree)', sign(strlen("`mata'")))								//conduct proximity computation, keep the matrix in Mata

		display "{txt}Matrix {cmd:prox} created in Mata" _newline ///
		"see {matacmd mata describe prox}"

	}

	else {																									//...otherwise if the matrix should be added to ereturned results

		quietly mata: proximity(results, `=e(ntree)', sign(strlen("`mata'")))

		display "{txt}Matrix {cmd:e(prox)} added to estimation results" _newline ///
		"see {stata ereturn list}"

	}

}

else if ("`anything'" == "fit") {																			//...or if in or out of sample overall fit based on Cramer's V is desired...

	mata: V_fit(results, `=e(ntree)', "`e(depvar)' `e(splitvars)'", sign(strlen("`insample'")), ///			//conduct the fit computation in Mata
	`=e(dvorder)', "", sign(strlen("`noisily'")), `=e(validmiss)', ///
	strtrim(regexr("`e(wexp)'", "=", "")), `=e(w_o_replace)')

	ereturn scalar fit = r(fit)																				//V_fit() produces an "r(fit)" metric - "ereturn" it

	if strlen("`insample'") display "{txt}In-sample only fit: {res}{col 40}" %16.4f e(fit)					//if "in-sample" estimates are requested - note this is the nature of the metric...

	else display "{txt}Resample fit: {res}{col 40}" %16.4f e(fit)											//...otherwise "out-of-sample" estimates are requested - noted in the display

}

else if ("`anything'" == "import") {																		//...or if a permutation importance matrix is desired...

	mata: permute_import(results, `=e(ntree)', ///															//conduct the importance computations in Mata
	"`e(depvar)' `e(splitvars)'", sign(strlen("`insample'")), ///
	`=e(dvorder)', sign(strlen("`noisily'")), `=e(validmiss)', strtrim(regexr("`e(wexp)'", "=", "")), ///
	`=e(w_o_replace)')

	tempname import																							//define temporary matrix to pass on returned matrix to ereturn

	matrix `import' = r(import)																				//pass returned matrix to temp for further processing

	matrix colnames `import' = `e(splitvars)'																//provide the columns names of the variables

	matrix rownames `import' = rank raw																		//label the rows

	ereturn matrix import = `import'																		//ereturn the importance matrix

	if strlen("`insample'") display "{txt}In-sample only importance:"										//if "in-sample" estimates are requested - note this is the nature of the metric...

	else display "{txt}Resample permutation importance:"													//...otherwise "out-of-sample" estimates are requested - noted in the display

	matrix list e(import), noheader																			//show the obtained matrix

}



else if ("`anything'" == "gettree")	{ 																		//...or if the results from a tree within the chaidforest is desired...

	if (`tree' > e(ntree)) | (`tree' <= 0) {																//must be a valid tree

		display "{err}Tree requested in {opt tree()} not valid."

		exit 198

	}

	capture summarize _CHAID_`tree' bwgt_`tree', meanonly													//see if _CHAID_ and bwgt_ variables corresponding to focal tree already exist

	if _rc == 0 drop _CHAID_`tree' bwgt_`tree'																//if _CHAID_ and bwgt_ variables exist, drop them

	mata: get_rules(results, `tree', "`e(depvar)' `e(splitvars)'", `=e(dvorder)', `=e(validmiss)', ///		//conduct the "rule generation" based on the stored information from the chaidtree object with a specific tree
	strtrim(regexr("`e(wexp)'", "=", "")), `=e(w_o_replace)')

	label variable _CHAID_`tree' "CHAID-defined cluster for tree `tree'"

	label variable bwgt_`tree' "bootstrapped frequency weight for tree `tree'"

	if strlen("`graph'") Display, tree(`tree')																//if desired, generate a graph and results display as that from chaid

}

else { 																										//...othwerwise the estat command is unrecognized

	display "{err}{cmd:estat} command not recognized."

	exit 199

}

end

/*Program to allow returning r-class results*/
program return_local, rclass

version 12.1

syntax anything

gettoken one two: anything, bind																			//split the syntax into components to return it - bind parentheses; all components are intended to be parentheses bound

local one = strtrim(subinstr(subinstr("`one'", "(", "", .), ")", "", .))									//remove binding parentheses in syntax

gettoken three four: one																					//split components of parentheses bound syntax

return local `three' "`four'"																				//return given the name 'three' the string found in 'four' - basically always used to return 'rules'

while strlen("`two'") {																						//while there are other rules left...

	gettoken one two: two, bind																				//...again, split the syntax into components to return it - bind parentheses; all components are intended to be parentheses bound

	local one = strtrim(subinstr(subinstr("`one'", "(", "", .), ")", "", .))								//...again, remove binding parentheses in syntax

	gettoken three four: one																				//...again, split components of parentheses bound syntax

	return local `three' "`four'"																			//...again, return given the name 'three' the string found in 'four' - basically always used to return 'rules'

}

end

/*Re-invoke object for CHAID forests use*/
version 12.1

mata:

mata set matastrict on

class chaidtree {

	real colvector clusters																					//each tree retains clusters

	string matrix CHAID_rules																				//each tree retains if/then statements

	string scalar used_vars																					//each tree retains variables used

	real colvector bootstrap_weight																			//each tree retains frequency weights from bootstrap or e(sample)-like data from prop_oos

}

end

/*Mata function to process rules*/
version 12.1

mata:

mata set matastrict on

function get_rules(class chaidtree results, real scalar tree, string scalar fit_inputs, ///
real scalar dvordered, real scalar missing, string scalar weight_var, real scalar w_o_replace) {

	/*declarations*/
	real matrix variable_locations

	real scalar path, variable

	string matrix rules, parsed_rules, variable_matrix, variable_levels, ///
	paths

	string rowvector variables, pass_locals, fullpath

	string scalar esample

	/*process resuls from specific tree in a way similar to base CHAID*/
	esample = st_tempname(1)																				//establish tempname for e(sample) for use by Mata

	stata("generate byte " + esample + " = e(sample)")														//generate variable that can be used by Mata

	rules = results[tree].CHAID_rules																		//pull the rules matrix from CHAIDtree object for the tree in question to parse and return

	variables = results[tree].used_vars																		//pull the variables used in focal CHAID tree from CHAIDtree object

	if (rows(rules) > 1) {																					//if there were any results, parse rules to pass...

		variable_matrix = ///																				//first, generate a block diagonal matrix "spreading" splitting variable names across blocks for matching
		I(cols(variables[2..cols(variables)]))#J(rows(rules[2..rows(rules), 2..cols(rules)]), ///
		cols(rules[2..rows(rules), 2..cols(rules)]), 1):* ///
		vec(J(rows(rules[2..rows(rules), 2..cols(rules)]), 1, variables[2..cols(variables)]))

		parsed_rules = J(cols(variables[2..cols(variables)]), cols(variables[2..cols(variables)]), ///		//then generate a matrix conformable with the variable_matrix above with all rules contained repeated to allow for matching
		rules[2..rows(rules), 2..cols(rules)])

		parsed_rules = regexm(parsed_rules, variable_matrix):*variable_matrix								//make indicator matrix indicating matches of variable names in specific locations corresponding to the CHAID rules

		variable_locations = regexm(parsed_rules, variable_matrix):*I(cols( ///								//notes the locations of a 'match' from the variables_matrix in the parsed_rules matrix
		variables[2..cols(variables)]))#J(rows(rules[2..rows(rules), 2..cols(rules)]), ///
		cols(rules[2..rows(rules), 2..cols(rules)]), 1)

		variable_levels = subinstr(subinstr( ///															//cuts out the variable name and underscore from the rules to keep only the value levels (also changes 'ms' to '.')
		J(cols(variables[2..cols(variables)]), cols(variables[2..cols(variables)]), ///
		rules[2..rows(rules), 2..cols(rules)]), (parsed_rules:+"_"), "", .), "ms", ".")

		variable_levels = variable_levels:*variable_locations												//omit all entries in places where there are no valid rules in the paths

		paths = ((parsed_rules:+"@"):*variable_locations):+variable_levels									//pair up variable names with levels and add in '@' to complete the paths matrix

		for (path = 1; path <= rows(rules) - 1; path++) {													//for each row of the paths section of the rules matrix...

			for (variable = 1; variable <= cols(variables) - 1; variable++) { 								//...and for each column potentially containing a rule

				if (variable == 1) fullpath = paths[path, 1..cols(rules)-1]									//put components of paths matrix into single column to collapse cleanly - first column

				else fullpath = fullpath:+paths[path+(variable-1)*(rows(rules)-1), ///						//put components of paths matrix into single column to collapse cleanly - additional columns beyond the first
				(cols(rules)-1)*(variable-1)+1..(cols(rules)-1)*variable]

			}

			if (path == 1) pass_locals = "(path1 " + ///													//enclose the entire entry in parentheses to parse and return by return_local - first path
			invtokens(stritrim(strtrim(subinword((fullpath:+";"), ";", "", .)))) + ")"

			else pass_locals = pass_locals, "(path" + strofreal(path)+ " " ///								//enclose the entire entry in parentheses to parse and return by return_local - additional paths beyond the first
			+ invtokens(stritrim(strtrim(subinword((fullpath:+";"), ";", "", .)))) + ")"

		}

		pass_locals = pass_locals, "(splitvars " + invtokens(variables[2..cols(variables)]) + ")"			//add in splitting variables to be returned by return_locals

	}

	/*return results*/
	else pass_locals = "(splitvars " + invtokens(variables[2..cols(variables)]) + ")"						//...otherwise return the splitvars as the only data to pass

	stata("return_local " + invtokens(pass_locals))															//invoke return local program with contents of pass_locals

	V_fit(results[tree], 1, fit_inputs, 1, dvordered, "", 0, missing, weight_var, w_o_replace)				//conduct the fit computation in Mata

	st_numscalar("r(N_clusters)", rows(rules) - 1)															//return number of clusters in focal CHAID tree

	if (rows(rules) > 1) {																					//if there are rules to return...

		st_matrix("r(sizes)", mm_freq(select(results[tree].clusters, results[tree].clusters:!=.))')			//return the sizes of each cluster - not weighted by bootstraps

		st_matrix("r(branches)", (rowsum(sign(strlen(rules[2..rows(rules), 2..cols(rules)]))))')			//return the number of branches made in each cluster

	}

	else {																									//...otherwise return the non-informative results

		st_matrix("r(sizes)", sum(sign(results[tree].clusters)))											//return number of non-0 weighted obs in first cluster

		st_matrix("r(branches)", 1)																			//only one branch

	}

	variable_locations = st_addvar(("float", "int"), (("_CHAID_" + strofreal(tree)), ///					//add 2 variables, bwgt_ and _CHAID - corresponding to the CHAID clusters as well as the bootstrapped frequency weight
	("bwgt_" + strofreal(tree))))

	st_store(., (("_CHAID_" + strofreal(tree)), ("bwgt_" + strofreal(tree))), ///							//store the _CHAID cluster and bootstrap weight data in the new variables
	st_varindex(esample), (results[tree].clusters, results[tree].bootstrap_weight))

}

end


/*Mata function to generate a proximity matrix between observations*/
version 12.1

mata:

mata set matastrict on

function proximity(class chaidtree results, real scalar number_of_trees, real scalar mata) {

	/*declarations*/
	real matrix clusters, proximity, valid_comparisons

	real scalar tree, obs

	string scalar matname

	/*set up for processing*/
	for (tree = 1; tree <= number_of_trees; tree++) {														//for each tree in the forest...

		if (tree == 1) clusters = results[1].clusters														//if the focal tree is first, make it the leader...

		else clusters = (clusters, results[tree].clusters)													//...otherwise, the focal tree follows the leader

	}

	proximity = valid_comparisons = J(rows(clusters), rows(clusters), 0)									//make a "dummy" proximity and valid_comparisons matrix of 0's

	clusters = clusters:*(clusters:!=.) + runiform(rows(clusters), cols(clusters)):*(clusters:==.)			//replace all missing values with random unformly distributed value in 0-1 range to (basically) ensure no matches on missings

	/*obtain matches*/
	for (obs = 2; obs <= rows(clusters) - 1; obs++) {														//foreach obervation...

		proximity[obs..rows(proximity), obs-1] = ///														//obtain number of matches with each other observation as it relates to cluster number - lower triangle
		rowsum((clusters[obs..rows(proximity), .]:==J(rows(clusters)-(obs-1), 1, clusters[obs-1, .])))

		valid_comparisons[obs..rows(valid_comparisons), obs-1] = ///										//obtain number of valid comparisons (i.e., in which both observations are non-missing - lower triangle
		rowsum(((clusters[obs..rows(valid_comparisons), .]:!=.):==J(rows(clusters)-(obs-1), 1, ///
		(clusters[obs-1, .]:!=.))))

	}

	proximity = proximity + proximity'																		//addin the upper triangle to make symmetric - proximity

	valid_comparisons = valid_comparisons + valid_comparisons'												//addin the upper triangle to make symmetric - valid_comparisons

	proximity = (proximity + I(rows(clusters))):/(valid_comparisons + I(rows(clusters)))					//obtain proportion matching

	proximity = sqrt(1:-proximity)																			//Breiman suggests making a dissimilarity matrix by reflecting across 0 and square-rooting

	/*return results*/
	if (!mata) {																							//if a Stata matrix...

		matname = st_tempname(1)																			//get a tempname

		st_matrix(matname, proximity)																		//save the proximity matrix as tempname

		stata("ereturn matrix prox " + matname)																//add the proximity matrix to ereturned results

	}

	else return(proximity)																					//...otherwise return Mata proximity matrix - necessary for large numbers of observations

}

end

/*Mata function to estimate permuataion importance - calls V_fit()*/
version 12.1

mata:

mata set matastrict on

function permute_import(class chaidtree results, real scalar number_of_trees, ///
string scalar splitting_variables, real scalar insample, real scalar ordered_response, ///
real scalar display, real scalar missing, string scalar weight_var, real scalar w_o_replace) {

	/*declarations*/
	real colvector importance

	real scalar variable

	/*proceed through variables to obtain importance*/
	importance = J(1, cols(tokens(splitting_variables)) - 1, .)												//generate dummy vector for fit values

	for (variable = 1; variable <= cols(tokens(splitting_variables)) - 1; variable++) {						//for each variable in the chaidforest() run...

		if (display) tokens(splitting_variables)[variable]	/**/

		V_fit(results, number_of_trees, splitting_variables, insample, ordered_response, ///				//proceed by feeding the data to V_fit() and choosing a permutation variable
		tokens(splitting_variables)[variable + 1], display, missing, weight_var, w_o_replace)

		importance[variable] = st_numscalar("r(fit)")														//record the fit value as permuted

	}

	st_matrix("r(import)", (mm_ranks(importance')' \ importance))											//return the ranked and raw importance vectors

}

end

/*Postestimation function to obtain overall model fit in- or out-of-sample; based on Cramer's V*/
version 12.1

mata:

mata set matastrict on

function V_fit(class chaidtree results, real scalar number_of_trees, ///
string scalar splitting_variables, real scalar insample, ///
real scalar ordered_response, string scalar permute, real scalar display, real scalar missing, ///
string scalar weight_var, real scalar w_o_replace) {

	/*declarations*/
	real matrix data, data_permuted

	real colvector unused_ob, clusters_reproduced, build_rule, results_vec, weights, ///
	weights_permuted

	real scalar V, V_sum, reduce_sum, permute_proceed, splitting_var, tree, path, rule

	string matrix current_rules

	string colvector names, variable_list, variable_list_permuted

	string scalar esample

	/*set-up for fit computation*/
	esample = st_tempname(1)																				//establish tempname for e(sample) when used by Mata

	stata("generate byte " + esample + " = e(sample)")														//generate "touse" variable that can be referenced by Mata

	data = st_data(., tokens(splitting_variables), st_varindex(esample))									//pull data into Mata

	if (missing) data = editmissing(data, .)																//change extended missings to regular missing for all variables (response and splitters)

	if (strlen(weight_var)) weights = weights_permuted = ///												//if there is an 'fweight', use it...
	st_data(., st_varindex(weight_var), st_varindex(esample))

	else weights = weights_permuted = J(rows(data), 1, 1)													//...otherwise assume equal weights

	names = tokens(splitting_variables)																		//split up variable names list

	variable_list = names[1]																				//start with response variable in actual names used list

	V_sum = reduce_sum = 0																					//make summed fit metric (to eventually average) and reduced sum (i.e., reduces denominator of average) variable start at 0

	for (splitting_var = 2; splitting_var <= cols(tokens(splitting_variables)); splitting_var++) {			//syntax to expand the imported data to a form usable by CHAID_rules - creates binary variables indicating levels of the splitting variable

		names[splitting_var] = ///																			//take name of splitting variable and "spread" it to all "levelsof" that variable separated by "_" as this is how they're registered in CHAID_rules
		invtokens((tokens(splitting_variables)[splitting_var]:+"_"):+strofreal(uniqrows(data[., ///
		cols(data)+(splitting_var-cols(tokens(splitting_variables)))])'))

		variable_list = variable_list, tokens(names[splitting_var])											//add the newly created names to the names vector to be matched on later

		if (splitting_var < cols(tokens(splitting_variables))) ///											//if the end of the variable list has not been reached, expand the data into set of binaries reflecting levels of the splitting variable, must "fit" the binaries in where the variable was to match with the names created...
		data = data[., 1..cols(data)+(splitting_var-cols(tokens(splitting_variables)))-1], ///
		J(1, rows(uniqrows(data[., cols(data)+(splitting_var-cols(tokens(splitting_variables)))])), ///
		data[., cols(data)+(splitting_var-cols(tokens(splitting_variables)))]):==uniqrows(data[., ///
		cols(data)+(splitting_var-cols(tokens(splitting_variables)))])', ///
		data[., cols(data)+(splitting_var-cols(tokens(splitting_variables)))+1..cols(data)]

		else ///																							//...otherwise the end of the variable list has been reached, expand the data into set of binaries reflecting levels of the splitting variable which can be plugged in at the end
		data = data[., 1..cols(data)-1], ///
		J(1, rows(uniqrows(data[., cols(data)])), data[., cols(data)]):==uniqrows(data[., cols(data)])'

	}

	/*computing fit and executing permutations*/
	for (tree = 1; tree <= number_of_trees; tree++) {														//for each tree...

		if (sign(strlen(permute))) permute_proceed = sum(strmatch(results[tree].used_vars, permute))		//if this is a part of a call from "permute_import()," only proceed if the variable chosen is in "used_vars"...

		else permute_proceed = 1																			//...otherwise this is a "regular" V_fit() call and just proceed

		if (permute_proceed) {																				//if the variable is present or this is a regular run...

			if (sign(strlen(permute))) {																	//again, if this is a part of a call from "permute_import()" ...

				variable_list_permuted  = variable_list, ""													//add a dummy-space to the list of "levelsof"-d splitting variables

				data_permuted = data, J(rows(data), 1, .)													//add a column of missings to the data

				if (strlen(weight_var)) data_permuted = mm_expand(data_permuted, weights) 					//if there are frequency weights - expand the data to allow for an accurate permutation

				data_permuted = select(data_permuted, ///													//permute or jumble() the data associated with the variable to permute, put it at the end of the data
				!regexm(variable_list_permuted, permute + "_[0-9]+")), ///
				select(data_permuted, regexm(variable_list_permuted, permute + ///
				"_[0-9]+"))[jumble(1::rows(data_permuted)), .]

				if (strlen(weight_var)) {

					weights_permuted = mm_freq(rowsum(data_permuted:*(10:^(cols( ///						//if there are frequency weights, obtain the updated frequency weights based on the newly permuted data
					data_permuted)..1)):/10), 1, uniqrows(rowsum(data_permuted:*(10:^( ///
					cols(data_permuted)..1)):/10)))

					data_permuted = mm_collapse(data_permuted, 1, (rowsum(data_permuted:*( ///				//collapse data based on newly pemuted variables
					10:^(cols(data_permuted)..1)):/10)))[., 2..cols(data_permuted)+1]

				}

				variable_list_permuted = select(variable_list_permuted, ///									//take the names associated with the permuted variable and put them at the end where the data are now
				!regexm(variable_list_permuted, permute + "_[0-9]+")), ///
				select(variable_list_permuted, regexm(variable_list_permuted, permute + "_[0-9]+"))

			}

			else {																							//accomodate the possibility of permuted data

				variable_list_permuted = variable_list

				data_permuted = data

			}

			if ((insample) & !sign(strlen(permute))) unused_ob = results[tree].bootstrap_weight				//if in-sample and not permuted, just use boostrap weight...

			else if ((insample) & sign(strlen(permute))) unused_ob = weights_permuted						//...otherwise if in-sample and permuted, use the newly formed weights from above...

			else unused_ob = mm_sample(sum(weights_permuted), rows(weights_permuted), ., ///				//...otherwise out-of-sample, based on the weight vector - do a fresh resampling
			weights_permuted, w_o_replace, 1)

			if ((insample) & !(sign(strlen(permute)))) clusters_reproduced = ///							//if in-sample fit is desired and this is not part of permute_import() run, use all the originally obtained clusters, changing all missing clusters to cluster "0"...
			editmissing(results[tree].clusters, 0)

			else {																							//...otherwise out-of-sample fit and/or a permute_import() run

				current_rules = results[tree].CHAID_rules													//pull in CHAID rules to apply to new observations

				clusters_reproduced = J(rows(data_permuted), 1, 0)											//start with clusters at "0" or, basically, missing

				for (path = 2; path <= rows(current_rules); path++) {										//for each row/path in CHAID rules...

					build_rule = J(rows(data_permuted), 1, 1)												//refresh the rule "builder" - start with everyone in-cluster and let them drop out

					for (rule = 2; rule <= cols(current_rules); rule++) {									//for each specific rule in the rules that needs accounting for...

						if (cols(tokens(current_rules[path, rule])) > 0) ///								//There are variables to account for in the focal rule, multiply all the binaries from the applicable splitting variable levels (i.e., the ones that match the rules) with one another as well as the current state of build_rule - this is an "and" or intersection situation, must meet all the rules' criteria or the observation gets a 0...
						build_rule = build_rule:*sign(rowsum(select(data_permuted, ///
						colsum(J(cols(tokens(current_rules[path, rule])), 1, ///
						variable_list_permuted):==tokens(current_rules[path, rule])'))))


					}

					clusters_reproduced = clusters_reproduced:+(sign(unused_ob):*(build_rule:*(path-1)))	//build the new cluster with the set of rules just obtained and its own unique number

					if (display) mm_freq(clusters_reproduced), uniqrows(clusters_reproduced)	//

				}

			}

			if ((sum(clusters_reproduced) > 0) & (rows(uniqrows(select(clusters_reproduced, ///				//if there are any non-0 clusters, and the number of clusters and levels of the response is more than 1...
			sign(unused_ob)))) > 1) & (rows(uniqrows(select(data_permuted[., 1], sign(unused_ob)))) > 1)) {

				results_vec = goodman((data_permuted[., 1], clusters_reproduced), ///						//implement loglinear/poisson model to obtain expected counts
				ordered_response, display, unused_ob)


				V = sqrt((results_vec[1]/sum(unused_ob))/min( ///											//compute the chi-square
				(rows(uniqrows(data[., 1]))-1, rows(uniqrows(clusters_reproduced))-1)))


			}

			else {																							//...otherwise the variable is not included

				V = 0																						//don't add to the "fit" sum

				reduce_sum++																				//increment "reduce_sum" to adjust average for invalid comparison (not a "true" 0)

			}

		}

		else {																								//...otherwise the variable is not included

			V = 0																							//don't add to the "fit" sum

			reduce_sum++																					//increment "reduce_sum" to adjust average for invalid comparison (not a "true" 0)

		}

		if (display) V	//

		V_sum = V_sum + V																					//add the value in "V" to the "fit" sum

	}

	st_numscalar("r(fit)", V_sum/(number_of_trees - reduce_sum))											//average the fit sum by the valid number of V's and pass to Stata as r() scalar

}

end


/*replayable display - copied from chaid.ado*/
program Display

version 12.1

syntax, tree(integer)

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

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

else {

	preserve

	clear

	mata: st_addobs(st_numscalar("r(N_clusters)"))

	quietly mata: labels = st_tempname(max(st_matrix("r(branches)")) + 1)

	quietly mata: st_addvar("str100", labels)

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

		local parse "`r(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)

	restore

}

end

/*Goodman association model based on loglinear/Poisson regression*/
version 12.1

mata:

mata set matastrict on

real rowvector goodman(real matrix data, real scalar ordered_response, ///
real scalar noisily, real colvector weight) {

	/*declarations*/
	real matrix IVs, DV

	real colvector observed_count, expected_count, combined_row_col

	real rowvector b, rescaled_b, xb

	real scalar converged, categories_IV, chi2, df, mean_count, responses_missing, ///
	ll, number_of_estimates, base_ll

	transmorphic row, col, estimation_object

	/*process inputs to goodman to estimate as poisson/loglinear model*/
	if (noisily) "Goodman/loglinear modeling Begins"	//

	if (sum(rowmissing(uniqrows(data[., 1])))) {												//are there missing values on response? Change them to -1's and activate "responses_missing"

		data[., 1] = ///																		//changes all missings to -1's
		editmissing(data[., 1], -1)

		responses_missing = 1																	//responses_missing is activated

	}

	else responses_missing = 0																	//otherwise deactivate responses_missing

	row = rowsum((J(1, rows(uniqrows(data[., 1])), ///
	data[., 1]):==(uniqrows(data[., 1])')):*(1..rows(uniqrows(data[., 1]))))					//set up response into column vector

	col = data[., 2]:*100																		//set up binary splitting variable/IV into column vector that can be "added to" response

	combined_row_col = (row + col)																//combines unique row/column cross-tab combinations in a vector

	if (noisily) uniqrows(combined_row_col)'	//

	observed_count = mm_freq(combined_row_col, weight, J(rows(uniqrows(col)), 1, ///			//obtain observed counts for all unique combinations in the cross-tabulation
	uniqrows(row)):+(uniqrows(col)#J(rows(uniqrows(row)), 1, 1)))

	if (noisily) observed_count	//

	/*row and/or column Goodman model - needs ml estimation*/
	if ((ordered_response) & (rows(observed_count) > 3)) {										//if there is a ordered response variable and there is a contingency table worth of data (at least 4 rows)...

		col = uniqrows(data[., 2])#J(rows(uniqrows(data[., 1])), 1, 1)							//set up columns of the cross-tabulation approproately for analysis (i.e., spread them across levels of the response as binary dummies - kind of irrelevant as only binaries are ever applied in chaidforest...)

		categories_IV = rows(uniqrows(data[., 2]))												//how many unique splitting variable levels are there?

		row = J(categories_IV, 1, uniqrows(data[., 1]))											//ordered response - treat row variable as ordered IV "predicing" counts as a loglinear model

		if (responses_missing) row = (editvalue(row, -1, 1), row:==min(row))					//when there is an ordered response, separate out the floating missing category into its own dummy (missing effect) but make it equal to the smallest value to reduce collinearity in the ordered variable

		if (noisily) row	//

		if (noisily) col	//

		estimation_object = moptimize_init()													//initalize ml estimation object

		//obtain starting values - based on poisson.ado's routine
		mean_count = mean(observed_count)														//obtains a mean of the counts as an imput to starting values

		IVs = cross((row, col), 1 , (row, col), 1)												//cross-products of IV's (for a linear regression)

		DV = cross((row, col), 1 , observed_count:-mean_count, 0)								//cross-products of IV's with DV (for a linear regression)

		b  = invsym(IVs)*DV																		//linear regression estimation by least squares

		rescaled_b = b:/mean_count																//rescale all estimates by mean number of observations which are used as initial/start values for poisson model

		/*set up ml estimation*/
		moptimize_init_evaluator(estimation_object, &loglin())									//tell ml evaluator what function to optimize

		moptimize_init_evaluatortype(estimation_object, "lf2")									//tell ml evaluator the optimization scheme (likelihood function with known 1st and 2nd derivatives of likelihood)

		moptimize_init_depvar(estimation_object, 1, observed_count)								//tell ml evaluator what is the dependent variable is

		moptimize_init_eq_indepvars(estimation_object, 1, (row, col))							//tell ml evaluator what is the independent variables are

		moptimize_init_search_rescale(estimation_object, "off") 								//tell ml evaluator not to spend time rescaling estimates to improve estimation - takes time to do

		moptimize_init_search(estimation_object, "off")											//tell ml evaluator not to spend time searching for better start values to improve estimation - takes time to do

		moptimize_init_eq_coefs(estimation_object, 1, rescaled_b')								//tell ml evaluator what the starting values of the parameters are

		moptimize_init_conv_maxiter(estimation_object, 20)										//tell ml evaluator to stop after 20 iterations - really should converge or have pretty reasonable values by then; keeps speed managable

		if (noisily) moptimize_init_tracelevel(estimation_object, "value")						//tell ml evaluator to show iterations if "noisily" is asked for by user...

		else moptimize_init_tracelevel(estimation_object, "none") 								//...otherwise show nothing

		moptimize_init_conv_warning(estimation_object, "off")									//tell ml evaluator not to warn when convergence is not reached

		moptimize(estimation_object)															//tell ml evaluator to proceed with estimation

		xb = moptimize_result_coefs(estimation_object)											//pull estimates from results

		if (noisily) xb	//

		converged = moptimize_result_converged(estimation_object)								//is everyting well?  That is, did estimation converge?

		if (sum(rowmissing(xb)) == 0) converged = 1												//keep estimates if they're usable in computation, even if convergence wasn't reached - we're data mining here

		expected_count = exp((row, col, J(rows(observed_count), 1, 1))*xb')						//compute expected frequencies

		if (noisily) expected_count	//

		df = (rows(uniqrows(data[., 2])))*(rows(uniqrows(data[., 1]))) - cols(xb)				//obtain DF

		number_of_estimates = cols(xb)															//obtain number of parameters estimated

	}

	/*unstructured model - traditional contingency table/product of margins approach*/
	else if (rows(observed_count) > 3) {														//if a contingency table can be made...

		row = mm_freq(data[., 1], weight)														//row/response variable margins

		if (noisily) row	//

		col = mm_freq(rowsum(data[., 2..cols(data)]:*(1..cols(data)-1)), weight)				//column/splitting variable margins

		if (noisily) col	//

		expected_count = (col#row):/sum(row)													//compute expected frequencies

		if (noisily) expected_count	//

		df = (rows(uniqrows(data[., 2])) - 1)*(rows(uniqrows(data[., 1])) - 1)					//obtain DF

		number_of_estimates = (rows(uniqrows(data[., 2])))*(rows(uniqrows(data[., 1]))) - ///
		(rows(uniqrows(data[., 2])) - 1)*(rows(uniqrows(data[., 1])) - 1)

		converged = 1																			//is everything well? (answer is always yes here)

	}

	else converged = 0																			//...otherwise, no contingency table - all's not well...

	if (converged) {																			//if all's well...

		chi2 = sum(((observed_count - expected_count):^2):/expected_count)						//obtain chi-square value

		row = mm_freq(data[., 1], weight)														//row/response variable margins to infer the ll from the lr test stat

		expected_count = select(expected_count, observed_count:>0)								//remove expected counts that will affect a sum

		observed_count = select(observed_count, observed_count:>0)								//remove 0's that will affect a sum

		ll = sum(observed_count:*ln(observed_count:/expected_count)) + ///						//take the liklihood ratio chi-square and use the constant-only multinomial distribution to obtain model log-likelihood by working backward
		sum((ln((row:/sum(row)))):*row)

		base_ll = sum((ln((row:/sum(row)))):*row)												//baseline/constant-only log-likilihood

	}

	else {																						//...otherwise return nothing

		chi2 = .

		ll = .

	}

	return((chi2, df, ll, number_of_estimates, base_ll))										//chaid needs the chi-square, degrees of freedom, base and model log-likelihood, and number of estimates - pass it to chaid

}

end


/*Loglinear model for moptimize()*/
version 12.1

mata:

mata set matastrict on

function loglin(transmorphic M, todo, real rowvector b, ll, S, H) {

	real rowvector inter, y, xb

	y = moptimize_util_depvar(M, 1)

	xb  = moptimize_util_xb(M, b, 1)

	ll = -exp(xb) + xb:*y - lngamma(y:+1)

	S = -exp(xb) + y

	inter = -exp(xb)

	H = moptimize_util_matsum(M, 1, 1, inter, moptimize_util_sum(M, ll))

}

end

/*sum function for mm_collapse()*/
version 12.1

mata:

mata set matastrict on

function dosum(real colvector data, real scalar nouse) {

	real scalar result

	result = sum(data)

	return(result)

}

end

/* programming notes and history

- chaidforest_estat version 1.0 - date - October 12, 2015

Basic version

-----