*! dominance0.mata version 0.0.0 8/11/2023 Joseph N. Luchman
version 15
**# Container object for domin specs
mata:
mata set matastrict on
struct domin_specs {
string scalar mi, reg, dv, all, weight, exp, touse, regopts
}
end
**# Mata function to compute all combinations of predictors or predictor sets run all subsets regression, and compute all dominance criteria
mata:
mata set matastrict on
void dominance0(struct domin_specs scalar model_specs, pointer scalar model_call, ///
string scalar cdlcompu, string scalar cptcompu, string scalar iv_string, ///
real scalar all_subsets_fitstat, real scalar constant_model_fitstat)
{
/*# object declarations*/
real matrix IV_antiindicator_matrix, conditional_dominance, ///
weighted_order_fitstats, weighted_order1less_fitstats, ///
weight_matrix, complete_dominance, select2IVfits, ///
IV_indicator_matrix, sorted_two_IVs, compare_two_IVs
string matrix IV_name_matrix
real rowvector fitstat_vector, orders, combin_at_order, ///
combin_at_order_1less, general_dominance, standardized, ///
focal_row_numbers, nonfocal_row_numbers, cpt_desig, ///
indicator_weight, antiindicator_weight
real colvector binary_pattern, cdl3, cpt_setup_vec, select2IVs
string colvector IVs
real scalar number_of_IVs, number_of_regressions, display, fitstat, ///
var1, var2, x, y
string scalar IVs_in_model
transmorphic cpt_permute_container, t, wchars, pchars, qchars
/*parse the predictor inputs*/
t = tokeninit(wchars = (" "), pchars = (" "), qchars = ("<>")) //set up parsing rules
tokenset(t, iv_string) //register the "iv_string" string scalar for parsing/tokenizing
IVs = tokengetall(t)' //obtain all IV sets and IVs as a vector
/*remove characters binding sets together (i.e., "<>")*/
for (x = 1; x <= rows(IVs); x++) {
if (substr(IVs[x], 1, 1) == "<") { //if any entry begins with "<" (will be the case for any/all sets)...
IVs[x] = substr(IVs[x], 1, strlen(IVs[x]) - 1) //first character removed ("<")
IVs[x] = substr(IVs[x], 2, strlen(IVs[x])) //last character removed (">")
}
}
/*set-up and compute all combinations of predictors and predictor sets*/
number_of_IVs = rows(IVs) //compute total number of IV sets and IVs
number_of_regressions = 2^number_of_IVs - 1 //compute total number of regressions
printf("\n{txt}Total of {res}%f {txt}sub-models\n", number_of_regressions)
if (number_of_IVs > 12) printf("\n{txt}Computing all independent variable combination sub-models\n")
IV_indicator_matrix = J(number_of_IVs, 2^number_of_IVs, .) //indicating the IV by it's sequence in the rows (each row is an IV) and presence in a model by the columns (each column is a model); all models and the IVs in those models will be represented in this matrix; the matrix starts off empty/as all missings
for (x = 1; x <= rows(IV_indicator_matrix); x++) { //fills in the IV indicators matrix - for each row...
binary_pattern = J(1, 2^(x-1), 0), J(1, 2^(x-1), 1) //...make a binary matrix that grows exponentially; (0,1), then (0,0,1,1), then (0,0,0,0,1,1,1,1) growing in size until it fills the entire set of columns with sequences of 0's and 1's...
IV_indicator_matrix[x, .] = J(1, 2^(number_of_IVs-x), binary_pattern) //...spread the binary pattern across all rows forcing it to repeat when not long enough to fill all columns - net effect is staggering all binary patters across rows to obtain all subsets in the final matrix
}
IV_indicator_matrix = IV_indicator_matrix[., 2..cols(IV_indicator_matrix)] //omit the first column that indicates a model with no IVs
IV_indicator_matrix = (colsum(IV_indicator_matrix) \ IV_indicator_matrix)' //transpose the indicator matrix and create a column indicating the "order"/number of IVs in the model - used to sort all the models
IV_indicator_matrix = sort(IV_indicator_matrix, (1..cols(IV_indicator_matrix))) //sort beginning with the order of the model indicator column followed by all other columns's binary values - net effect results in same sort order as cvpermute() as desired/originally designed in domin 3.0
orders = (IV_indicator_matrix[.,1]')[rows(IV_indicator_matrix)..1] //keep the orders of each model - sort order is reversed below which is also implemented here
IV_indicator_matrix = IV_indicator_matrix[., 2..cols(IV_indicator_matrix)]' //omit orders from IV indicators matrix and transpose back to the models as rows
//IV_indicator_matrix = sort(((cols(IV_indicator_matrix)::1), IV_indicator_matrix'), 1)[., 2..rows(IV_indicator_matrix)+1]' //reverse sort order, functions below expect reversed order
IV_indicator_matrix = IV_indicator_matrix[., cols(IV_indicator_matrix)..1] // reverse the IV indicator matrix's order; functions below expect a reversed order as originally designed in domin 3.0
IV_name_matrix = IV_indicator_matrix:*IVs //apply string variable names to all subsets indicator matrix
/*all subsets regressions and progress bar syntax if predictors or sets of predictors is above 5*/
display = 1 //for the display of dots during estimation - keeps track of where the regressions are - every 5% of models complete another "." is added
if (number_of_IVs > 4) {
printf("\n{txt}Progress in running all sub-models\n{res}0%%{txt}{hline 6}{res}50%%{txt}{hline 6}{res}100%%\n")
printf(".")
displayflush()
}
fitstat_vector = J(1, cols(IV_indicator_matrix), .) //pre-allocate container vector that will contain fitstats across all models
for (x = 1; x <= number_of_regressions; x++) { //loop to obtain all possible regression subsets
if (number_of_IVs > 4) {
if (floor(x/number_of_regressions*20) > display) {
printf(".")
displayflush()
display++
}
}
IVs_in_model = invtokens(IV_name_matrix[., x]') //collect a distinct subset of IVs, then collpase names into single string separated by spaces
fitstat = (*model_call)(IVs_in_model, all_subsets_fitstat, ///
constant_model_fitstat, model_specs) //implement called model - will differ for domin vs. domme
fitstat_vector[1, x] = fitstat //add fitstat to vector of fitstats
}
/*define the incremental prediction matrices and combination rules*/
IV_antiindicator_matrix = (IV_indicator_matrix:-1) //matrix flagging which IVs are not included in each model and setting them up for a subtractive effect
combin_at_order = J(1, number_of_regressions, 1):*comb(number_of_IVs, orders) //vector indicating the number of model combinations at each "order"/# of predictors for each models - this is used as a "weight" to construct general dominance statistics
combin_at_order_1less = J(1, number_of_regressions, 1):*comb(number_of_IVs - 1, orders) //vector indicating the number of model combinations from the previous "order" - this is also used as a "weight" to construct general dominance statistics
combin_at_order_1less[1] = 0 //replace missing first value in vector with 0
combin_at_order = combin_at_order - combin_at_order_1less //remove # of combinations for the "order" less the value at "order" - 1; this gives the number of relevant combinations at each order that include the focal IV and thus produce the right weight for averaging
indicator_weight = IV_indicator_matrix:*combin_at_order //"spread" the vector indicating number of combinations at the current order involving the relevant IV to all models including that IV - to be used as a weight for summing the fitstat values in raw form (not increments)
antiindicator_weight = IV_antiindicator_matrix:*combin_at_order_1less //"spread" the vector indicating number of combinations at the previous order not involving the relevant IV to all models not including the relevant IV - these components assist to subtract out values to make the values "increments"
/*define the full design matrix - compute general dominance (average conditional dominance across orders)*/
weight_matrix = ((indicator_weight + antiindicator_weight):*number_of_IVs):^-1 //combine weight matrices (which reflect the conditional dominance weights/within order averages) with the number of ivs (between order averages) and invert cell-wise - now can be multiplied by fit stat vector and summed to obtain general dominance statistics
general_dominance = colsum((weight_matrix:*fitstat_vector)') //general dominance weights created by computing product of weights and fitstats and summing for each IV row-wise; in implementing the rows are transposed and column summed so it forms a row vector as will be needed to make it an "e(b)" vector
fitstat = rowsum(general_dominance) + all_subsets_fitstat + constant_model_fitstat //total fitstat is then sum of gen. dom. wgts replacing the constant-only model and the "all" subsets stat
st_matrix("r(domwgts)", general_dominance) //return the general dom. wgts as r-class matrix
standardized = general_dominance:*fitstat^-1 //generate the standardized gen. dom. wgts
st_matrix("r(sdomwgts)", standardized) //return the stdzd general dom. wgts as r-class matrix
st_matrix("r(ranks)", mm_ranks(general_dominance'*-1, 1, 1)') //return the ranks of the general dom. wgts as r-class matrix
st_numscalar("r(fs)", fitstat) //return overall fit statistic in r-class scalar
/*compute conditional dominance*/
if (strlen(cdlcompu) == 0) {
if (number_of_IVs > 5) printf("\n\n{txt}Computing conditional dominance\n")
conditional_dominance = J(number_of_IVs, number_of_IVs, 0) //pre-allocate contrainer matrix to hold conditional dominance stats
weighted_order_fitstats = ((IV_indicator_matrix:*combin_at_order):^-1):*fitstat_vector // create matrix fit stats weighted by within-order counts - to be summed to create the conditional dominance statistics; these fit stats are "raw"/not incrments without the antiindicators
weighted_order1less_fitstats = ((IV_antiindicator_matrix:*combin_at_order_1less):^-1):*fitstat_vector // create matrix fit stats weighted by within-order counts at the previous order - these are negative and also to be summed to create the conditional dominance statistics; they ensure that the values are increments
/*loop over orders/number of predictors to obtain average incremental prediction within order*/
for (x = 1; x <= number_of_IVs; x++) { //proceed order by order
conditional_dominance[., x] = rowsum(select(weighted_order_fitstats, orders:==x)):+rowsum(select(weighted_order1less_fitstats, orders:==x-1)) //sum the weighted fit statistics at the focal order and the weighted (negative) fit statistics from the previous order
}
st_matrix("r(cdldom)", conditional_dominance) //return r-class matrix "cdldom"
}
/*compute complete dominance*/
if (strlen(cptcompu) == 0) {
if (number_of_IVs > 5) printf("\n{txt}Computing complete dominance\n")
complete_dominance = J(number_of_IVs, number_of_IVs, 0) //pre-allocate matrix for complete dominance
cpt_setup_vec = (J(2, 1, 1) \ J(number_of_IVs - 2, 1, 0)) //set-up a vector the length of the number of IVs with two "1"s and the rest "0"s; used to compare each 2 IVs via 'cvpermute()'
cpt_permute_container = cvpermutesetup(cpt_setup_vec) //create a 'cvpermute' container to extract all permutations of two IVs
for (x = 1; x <= comb(number_of_IVs, 2); x++) {
select2IVs = cvpermute(cpt_permute_container) //invoke 'cvpermute' on the container - selects a unique combination of two IVs
select2IVfits = colsum(select(IV_indicator_matrix, select2IVs)):==1 //filter IV indicator matrix to include just those fit statistic columns that inlude focal IVs - then only keep the columns that have one (never both or neither) of the IVs
focal_row_numbers = (1..rows(IV_indicator_matrix)) //sequence of numbers for selecting specific rows indicating focal IVs
nonfocal_row_numbers = select(focal_row_numbers, (!select2IVs)') //select row numbers for all non-focal IVs; needed for sorting (below)
focal_row_numbers = select(focal_row_numbers, select2IVs') //select row numbers for all focal IVs; also needed for sorting
sorted_two_IVs = ///
sort((select((IV_indicator_matrix \ fitstat_vector), select2IVfits))', /// //combine the indicators for IVs and fit statistic matrix - selecting only those columns corresponding with fit statistics where the focal two IVs are located, then...
(nonfocal_row_numbers, focal_row_numbers)) // ...sort this matrix on the nonfocal IVs first then the focal IVs. This ensures that sequential rows are comparable ~ all odd numbered rows can be compared with its subsequent even numbered row
compare_two_IVs = ///
(select(sorted_two_IVs[,cols(sorted_two_IVs)], mod(1::rows(sorted_two_IVs), 2)), /// //select the odd rows for comparison (see constuction of 'sorted_two_IVs')
select(sorted_two_IVs[,cols(sorted_two_IVs)], mod((1::rows(sorted_two_IVs)):+1, 2))) //select the even rows for comparison
cpt_desig = ///
(all(compare_two_IVs[,1]:>compare_two_IVs[,2]), /// //are all fit statistics in odd/first variable larger than the even/second variable?
all(compare_two_IVs[,1]:<compare_two_IVs[,2])) //are all fit statistics in even variable larger than the odd variable?
if (cpt_desig[2] == 1) complete_dominance[focal_row_numbers[1], focal_row_numbers[2]] = 1 //if the even/second variables' results are all larger, record "1"...
else if (cpt_desig[1] == 1) complete_dominance[focal_row_numbers[1], focal_row_numbers[2]] = -1 //else if the odd/first variables' results are all larger, record "-1"...
else complete_dominance[focal_row_numbers[1], focal_row_numbers[2]] = 0 //otherwise record "0"...
}
complete_dominance = complete_dominance + complete_dominance'*-1 //make cptdom matrix symmetric
st_matrix("r(cptdom)", complete_dominance) //return r-class matrix "cptdom"
}
}
end