*! chaidforest_pr - version 2.0 - 9/14/2015 - Joseph N. Luchman
program define chaidforest_pr //program history and notes at end of file
version 12.1
syntax anything, [pr mode ctpr ctmode useboot]
/*set-up for prediction*/
tempvar id touse temp_dv //temp variable declarations
quietly generate long `id' = _n //id variable made for future merging
quietly generate byte `touse' = e(sample) //restrict to estimation sample based on e(sample) which should be in memory still
if !strlen("`pr'`mode'") local rftype "pr" //if there is no option specified, make empirical probabilities the default for random forest predictions
else if (strlen("`pr'`mode'") == 6) { //...or if both are specitied, indicate problem...
display "{err}Only one of {opt pr} and {opt mode} can be selected."
exit 198
}
else local rftype "`pr'`mode'" //...otherwise user specified prediction for random forest predictions
if !strlen("`ctpr'`ctmode'") local cttype "ctpr" //make empirical probabilities the default for individual trees
else if (strlen("`ctpr'`ctmode'") == 10) { //problem - conflicting options for individual trees
display "{err}Only one of {opt ctpr} and {opt ctmode} can be selected."
exit 198
}
else local cttype "`ctpr'`ctmode'" //otherwise user specified prediction for individual trees
quietly clonevar `temp_dv' = `e(depvar)' //generate temporary response variable to recode missings if needed
if `=e(validmiss)' {
local missing "missing" //can missings be predicted?
quietly replace `temp_dv' = . if missing(`temp_dv') //recode any extended missings to "." - necessary to get right number of predicted probabilities
}
quietly levelsof `temp_dv' if e(sample), local(levs) `missing' //obtain levels of response variable
capture mata: length(results) //is the needed "results" chaidforest object in memory?
if _rc {
display "{err}Chaidtree object {cmd:results} not found." _newline ///
"{cmd:chaidforest} will have to be re-executed."
exit 198
}
if ("`rftype'" == "pr") capture mata: stata("describe " + ///
invtokens(st_local("anything"):+strofreal((0..cols(tokens(st_local("levs")))-1))) ", varlist") //are there already a set of variables in the dataset with the user-specified stem?
else capture describe `anything', varlist //is there already a variable in the dataset with the user-specified name?
if (!_rc & regexm("`r(varlist)'", "`anything'[1]")) { //notify user of conflict
display "{err}`anything' already defined."
exit 110
}
if ("`rftype'" == "pr") local levels = `:list sizeof levs' //if multiple variables will be constructed, note number of levels of response variable...
else local levels = 1 //...otherwise, there's just 1 level/variable
mata: predict_cf(results, `=e(ntree)', "`rftype'", `=e(validmiss)', "`levs'", /// //implement predicted value computation
"`temp_dv'", "`touse'", "`anything'", "`cttype'", "`e(depvar)'", strlen("`useboot'"), ///
strtrim(regexr("`e(wexp)'", "=", "")))
end
/*Re-invoke object to use for prediction*/
version 12.1
mata:
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
/*Predicted values function*/
version 12.1
mata:
mata set matastrict on
real matrix predict_cf(class chaidtree results, real scalar ntree, ///
string scalar type, real scalar missing, string scalar levs, string scalar dvname, ///
string scalar touse, string scalar predname, string scalar cttype, string scalar displayname, ///
real scalar bootstrap_use, string scalar weight_var)
{
/*declarations*/
real matrix preds, predicted, matchup
real colvector current_preds, dv, weights, rand_resolve, tlevs
real scalar countup, varnum, x, y
string colvector storednames
/*begin set-up*/
countup = -1 //placeholder which is incremented when naming variables
weights = 1 //make default weight value 1
if (sign(strlen(weight_var))) weights = st_data(., st_varindex(weight_var), st_varindex(touse)) //if there is an 'fweight', use it
if (bootstrap_use) weights = results.bootstrap_weight //if the user would prefer making the results depend on the bootstrapped results, replace weight with the bootstrapped weights
tlevs = strtoreal(tokens(levs)') //tokenizes levels of response variable
dv = st_data(., dvname, st_varindex(touse)) //pull in the response variable from Stata
if (missing) { //remove missings from levels of response variable - calls them -1's, pull them into frequncy computations
levs = subinstr(levs, ".", "-1")
tlevs = editmissing(tlevs, -1)
dv = editmissing(dv, -1)
}
/*process individual trees - ready them for incorporation to individual-level predictions*/
for (x = 1; x <= ntree; x++) { //for each tree in the chaidforest...
if (cttype == "ctmode") { //if the desired predicted value per observation for this specific tree is the modal value...
matchup = mm_collapse(dv, weights, results[x].clusters, &mode2()) //find the modal value of the response for this tree by cluster/node
current_preds = rowsum((J(1, rows(matchup), ///
results[x].clusters):==(matchup[., 1]')):*matchup[., 2]') //"spread" the appropriate modal response variable value to members of the appropriate cluster/node
current_preds = current_preds:*exp(ln(results[x].clusters:!=.)) //replace non-sampled values with missing
}
else { //...otherwise the desired prediction is the empirical proportion falling into each response category per observation for this specific tree
for (y = 1; y <= rows(tlevs); y++) { //for each level of the response variable...
matchup = mm_collapse(dv:==tlevs[y], weights, /// //obtain mean/proportion in focal level of response variable by cluster for this tree
results[x].clusters)
matchup = rowsum((J(1, rows(matchup), /// //"spread" the appropriate proportion value to members of the appropriate cluster/node
results[x].clusters):==(matchup[., 1]')):*matchup[., 2]')
if (y == 1) current_preds = matchup:*exp(ln(results[x].clusters:!=.)) //if at the first level, make it the leading column vector, all non-sampled observations are assigned as missing...
else current_preds = current_preds, matchup:*exp(ln(results[x].clusters:!=.)) //...otherwise, make the focal column vector follow the lead, all non-sampled observations are assigned as missing...
}
}
if (x == 1) preds = current_preds //if at the first tree, make it the leading column vector for the forest's predictions...
else preds = (preds, current_preds) //...otherwise, make the focal tree's predictions follow the lead
}
/*process across trees to obtain individual-level predictions*/
for (x = 1; x <= rows(preds); x++) { //loop over observations...
if (x == 1) { //if at the first observation...
if (cttype == "ctpr") predicted = /// //if tree-level empirical probabilities asked for, average proportions in each level of the response across all trees in the forest for the focal individual...
mm_collapse(preds[x, .]', 1, J(cols(preds)/rows(tlevs), 1, (1::rows(tlevs))))[., 2]'
else if ((cttype == "ctmode") & (type == "pr")) predicted = /// //...or if tree-level modal values were asked for, yet the predicted probabilities in each category across all trees are desired, find the proportion of all trees predicting each level of the response variable for the focal individual...
(mm_freq(preds[x, .]', 1, tlevs)'):/sum(mm_freq(preds[x, .]', 1, tlevs)')
else predicted = mode(preds[x, .]', levs, 1) //...otherwise, find the modal value across all trees in the forest given the modal values at each individual tree
}
else { //...otherwise, for all other observations
if (cttype == "ctpr") predicted = /// //if tree-level empirical probabilities asked for, average proportions in each level of the response across all trees in the forest for the focal individual...
predicted \ mm_collapse(preds[x, .]', 1, J(cols(preds)/rows(tlevs), 1, (1::rows(tlevs))))[., 2]'
else if ((cttype == "ctmode") & (type == "pr")) predicted = (predicted \ /// //...or if tree-level modal values were asked for, yet the predicted probabilities in each category across all trees are desired, find the proportion of all trees predicting each level of the response variable for the focal individual...
(mm_freq(preds[x, .]', 1, tlevs)'):/sum(mm_freq(preds[x, .]', 1, tlevs)'))
else predicted = (predicted \ mode(preds[x, .]', levs, 1)) //...otherwise, find the modal value across all trees in the forest given the modal values at each individual tree
}
}
/*finalize predictions*/
if ((cttype == "ctpr") & (type == "mode")) { //if tree-level empirical probabilities asked for, yet the user wants modal/most likely values based on those averaged probabilities...
matchup = predicted:==rowmax(predicted) //first, identify which row has the highest probability for each observation
if (sum(matchup):>1 != rows(matchup)) { //in the case there are ties...
rand_resolve = /// //randomly resolve which is predicted
((rowsum(matchup):>1):*runiform(rows(matchup), rows(tlevs))):*matchup
matchup = matchup:+((rand_resolve:==rowmax(rand_resolve)):*(rowsum(matchup):>1)) //merge the randomly chosen values in to replace the ties
}
predicted = rowsum(matchup:*tlevs'):*editvalue((rowmax(predicted):!=.), 0, .) //apply the levels to the binaries indiciating column associated with each level - basically, this assigns the predicted value and ensures missings persist
}
if (missing) predicted = editvalue(predicted, -1, .p) //if there are missings, replace -1's with "predicted missings" or ".p"
if (type == "mode") { //if a modal value only is desired...
varnum = st_addvar("double", predname) //generate a double with the desired name
mata: st_varlabel(varnum, "Modal value: " + displayname) //give the new variable an appropriate label
storednames = predname //storednames is more flexible than predname here, just use storednames (shortens syntax)
}
else { //...otherwise the user asked for empirical probabilities
for (y = 1; y <= cols(tokens(levs)); y++) { //for each level of the response...
varnum = st_addvar("double", predname + strofreal(++countup)) //generate a double as the predict function would
if (tokens(levs)[y] == "-1") /// //if the value of the level is -1, replace it with ".p" - give the variable a label...
st_varlabel(varnum, "Pr(" + displayname + "==.p)")
else st_varlabel(varnum, "Pr(" + displayname + "==" + tokens(levs)[y] + ")") //otherwise give the variable a label as predict would
if (y == 1) storednames = st_varname(varnum) //if the first level lead the names...
else storednames = (storednames, st_varname(varnum)) //...otherwise follow the leader
}
}
st_store(., storednames, st_varindex(touse), predicted) //put the predictions into Stata as requested
}
end
/*Mata function to return a modal value across all trees in the forest (ties are randomly decided)*/
version 12.1
mata:
mata set matastrict on
function mode(real colvector data, string scalar levs, real scalar nouse) {
real scalar rand_resolve
transmorphic result
result = mm_freq(data, 1, strtoreal(tokens(levs)')) //get frequencies on focal variable
if (!sum(result)) result = . //if frequencies are 0, assign missing
else {
result = strtoreal(select(tokens(levs)', result:==max(result))) //retain the level(s) associated with the highest frequency
if (rows(result) > 1) { //if there's a tie...
rand_resolve = mm_sample(1, rows(result)) //randomly sample a value from "result"
result = result[rand_resolve, 1] //assign the randomly sampled value to "result"
}
}
return(result) //return the modal value
}
end
/*Mata function to return a modal value by cluster within a tree (ties are randomly decided)*/
version 12.1
mata:
mata set matastrict on
function mode2(real colvector data, real scalar nouse) {
real matrix result, rand_resolve
result = mm_freq(data, 1, uniqrows(data)) //get frequencies on focal variable
result = select(uniqrows(data), result:==max(result)) //retain the level(s) associated with the highest frequency
if (rows(result) > 1) { //if there's a tie..
rand_resolve = runiform(rows(result), 1) //obtain random numbers
result = select(result, rand_resolve:==max(rand_resolve)) //select the max of the random numbers ans assign it to "result"
}
return(result) //return the mode
}
end
/* programming notes and history
- chaidforest_pr version 1.0 - date - October 27, 2014
Basic version
-----
- chaidforest_pr version 2.0 - date - September 14, 2015
a] more prediction options (empirical probabilities at individual tree-level) - moves all prediction to this function
b] nixes lasting "predicted" matrix
c] allows user to "use" bootstrap weights to form predicted values (for fun!)
d] compatible with, and checks for, fweights