/*
Parse the user input, estimate all of the possible models, store the results in
the results_tempfile.
*/
program multivrs_estimate, rclass
syntax anything (name = model_and_varlist id = "varlist") [if] [in] ///
[,  fe re be mle pa irr offset(passthru) exposure(passthru) absorb(passthru) vce(passthru) ///
noplot plotbs savelist(string) saveas(string) weights(string) alpha(passthru) ///
sample(passthru) size(passthru) bs(passthru) pref(passthru)  ///
estonly compare normal nozero nolistwise ///
 other(passthru) ///
noinfluencecalcs inf_means ///
noinfluence replace sig_only nosig bs_reps(integer 50) results_tempfile(string) intervals ///
margins marginsopts(string)]
//set up data file to store estimation results
tempname memhold
tempname mem_var_combos
tempfile temp_nonpar_bs_results
local results `"`results_tempfile'"'

tempfile tempfile_var_combos

//local tempfile_var_combos "temp\var_combos"

gettoken model_namelist varlist : model_and_varlist, bind match(par)
multivrs_parse_models `model_namelist'
local nmodeltypes `"`r(nmodeltypes)'"'
local model_namelist_uniq `"`r(model_namelist_uniq)'"'
local model_idlist `"`r(model_idlist)'"'
local model_namelist `"`r(model_namelist)'"'
forvalues i = 1/`nmodeltypes' {
	local model`i'_name `"`r(model`i'_name)'"'
	local model`i'_opts `"`r(model`i'_opts)'"'
	local model`i'_id `"`r(model`i'_id)'"'
}
return local model `model_namelist'
return local model_idlist `model_idlist'

return local varlist `varlist'
return scalar nmodeltypes = `nmodeltypes'
//If the user has not specified a sig_only option or marginal effects then assign it based on the model choices
if "`sig_only'" == "" & "`sig'" != "nosig" & "`margins'" =="" local sig_only `"`r(sig_only_models)'"'
//Otherwise it will stay as "sig_only" or not as the user specified
if "`weights'" == "" local weights "uniform"
// BIC does not work for rreg
if "`weights'" == "bic" & strmatch("`model_namelist_uniq'", "*rreg*") {
	di as error "BIC cannot be calculated for rreg; please select another type of weights (inf, r2, no, or uniform)."
	exit 198
}

//If the weighting option is not specified then set it to uniform

multivrs_parse_varlist `varlist'

local nSets `"`r(nSets)'"'
local n_var_combinations `"`r(n_var_combinations)'"'
local allvarnames `"`r(allvarnames)'"'
local allvarnames_clean `"`r(allvarnames_clean)'"'
local b_always_in_clean `"`r(b_always_in_clean)'"'
local always_in_varnames `"`r(always_in_varnames)'"'
local intvarlist `"`r(intvarlist)'"'
local depvarlist `"`r(depvarlist)'"'
local intvarlist_noline : subinstr local intvarlist "|" " ", all
local depvarlist_noline : subinstr local depvarlist "|" " ", all
forvalues is = 1/`nSets' {
	local s`is'NSX `"`r(s`is'NSX)'"'
	forvalues isx = 1/`s`is'NSX' {
		local s`is'sx`isx'NT = `"`r(s`is'sx`isx'NT)'"'
		forvalues it = 1/`s`is'sx`isx'NT' {
		 local s`is'sx`isx't`it'NTX = `"`r(s`is'sx`isx't`it'NTX)'"'
			forvalues itx = 1/`s`is'sx`isx't`it'NTX' {
				local s`is'sx`isx't`it'tx`itx' `"`r(s`is'sx`isx't`it'tx`itx')'"'
			}
		}
	}
}
if "`sig_only'" == "" {
	local nwarnings 0
	local lists_to_check intvarlist depvarlist
	foreach list of local lists_to_check {
		if `: word count ``list'_noline'' > 1 {
			local list_with_commas : subinstr local `list'  "|" " ", all
			//local list_with_commas = subinstr("``list''", "|",", ",.)
			local ++nwarnings
			di as text " Warning:  Check that " as result `"`list_with_commas'"' as text " are on a comparable scale, "
			di as text "           so that the meaning of a one unit change is consistent."
		}
	}
	if `nwarnings' > 0 	di as text "           If not, consider standardizing the variables or use the " as result "sig_only" as text " option."
}

multivrs_parse_est_options `model_namelist_uniq' , `fe' `re' `be' `mle' `pa' `irr' ///
`vce' `offset' `exposure' `absorb' `other'
local opts_command `"`r(opts_command)'"'
multivrs_parse_multivrs_options, nSets(`nSets') `sample' `size' `bs' `pref' `alpha' weights(`weights')
local sample_pct `"`r(sample_pct)'"'
local sizemin `"`r(sizemin)'"'
local sizemax `"`r(sizemax)'"'
local bs_type `"`r(bs_type)'"'
if "`intervals'" == "intervals" & "`bs_type'" == "" local bs_type par
local bs_options_user `"`r(bs_options_user)'"'
local bs_opts_all `"`r(bs_opts_all)'"'
local prefb `"`r(prefb)'"'
local prefse `"`r(prefse)'"'
local alpha `"`r(alpha)'"'
local weights `"`r(weights)'"'
local opts_multivrs `"`r(opts_multivrs)'"'


//identify observations meeting "if" and "in" criteria
tempvar touse
mark `touse' `if' `in'
quietly count if `touse'
local total_N = r(N)
if `total_N' == 0 {
	error 2000
}

/*Listwise deletion:  mark out observations with any missing values.
Alert the user to how many observations have been marked out.
TODO: "Sample size varies across model specifications" may not be strictly true, e.g. if the only var
with missing values is always-in.
*/
if "`listwise'" == "nolistwise" {
	tempvar touse_temp
	gen `touse_temp' = `touse'
}

markout `touse' `allvarnames'
quietly count if `touse'
local N = r(N)
if `N' == 0 {
	error 2000
}
local nmissing = `total_N' - `N'
if (`nmissing' > 0) di as text "note:  sample size varies across model specifications."
if "`listwise'" != "nolistwise" {
	if (`nmissing' > 0) di  "Listwise deletion:  " `nmissing' as text " out of "  `total_N' as text " observations will not be used."
}
if "`listwise'" == "nolistwise" {
	local N = `total_N'
	drop `touse'
	gen `touse' = `touse_temp'
	drop `touse_temp'
}

//save the list of control term combinations in a text file
if "`savelist'" != "" {
	quietly file open listfile using `savelist', write append
	* if there are options specified to drop any variables, then included this section.
	if !("`if'" == "" & "`in'" == "" & "`listwise'" == "nolistwise") {
		file write listfile _n _n "* Generate the indicator variable for observations to use (incorporating any listwise delete, [if] and/or [in] statements)."
		file write listfile _n "tempvar touse"
		file write listfile _n "mark \`touse' `if' `in'"
		if "`listwise'" != "nolistwise"	{
			file write listfile _n "markout \`touse' `allvarnames'"
		}
	}
	file write listfile _n _n "* The number preceding each model corresponds to the variable model_id in the results dataset."
}


*check on variables of interest for xtreg Fixed Effects
forvalues im = 1/`nmodeltypes' {
	local fe fe
	local estimation_options `"`model`im'_opts' `opts_command'"'
	if "`model`im'_name'" == "xtreg" & `: list fe in estimation_options' == 1 {
		foreach intvar of local intvarlist_noline {
			quietly xtsum `intvar' if `touse'
			if r(sd_w)==0 {
				di as err "Variable of interest `intvar0' has 0 within standard deviation."
				exit 198
			}
		}
	}
}

//set margins option
// default for 1 model type : no margins.
// but it can be overriden by specifying margins or marginopts.

/*
if `nmodeltypes' == 1 {
	if "`margins'" == "margins" | "`marginsopts'" != "" local margins "margins"
	else local margins "nomargins"
}
else if `nmodeltypes' > 1 {
	/*if "`margins'" == "nomargins" {
		di "Note: Coefficient estimates from different model types may not be directly comparable."
		di "The -margins- option is recommended to compare estimated effects across the specified models."
	}
	else {*/
		local margins "margins"
	//}
}
*/

// prepare margin options : make sure there is exactly 1 "post" included.
if "`margins'" == "margins" {
	local marginsopts_nopost : subinstr local marginsopts "post" "", all
	local marginsopts_post "`marginsopts_nopost' post"
}

//save the user's dataset as a tempfile with touse
tempfile orig_touse
quietly save `orig_touse'

local nrotated = `nSets' - 2
local bs_reps = 50

/*
Estimate the time required for the multivrs command by timing 3 regressions,
averaging the time required and multiplying by the number of models.
*/
local timer_2 99
local timesum 0
local ntests 5
local time_multiplier 1.3
local curr_varlist ""
forvalues is = 1/`nSets' {
	forvalues it = 1/`s`is'sx1NT' {
		local curr_varlist `"`curr_varlist' `s`is'sx1t`it'tx1'"'
	}
}
forvalues im = 1/`nmodeltypes' {
	local model `"`model`im'_name'"'
	forvalues iv = 1/`ntests' {
		timer clear `timer_2'
		timer on `timer_2'

		if "`bs_type'" == "nonpar" quietly bootstrap `bs_opts_all' : ///
			`model' `curr_varlist' if `touse', `opts_command' `opts`im''
		else {
			quietly `model' `curr_varlist' if `touse', `opts_command' `opts`im''
			forvalues ib = 1/`bs_reps' {
				local x = rnormal(0,1)
			}
		}
		timer off `timer_2'
		quietly timer list `timer_2'
		local timesum = `timesum' + r(t`timer_2')
	}
}
local time_avgreg = `time_multiplier'*`timesum'/(`ntests'*`nmodeltypes')
di _newline _continue
local nmodels = `n_var_combinations' * `nmodeltypes'
local large_nmodels 0
local large_threshold 1000

/*Publish number of models to calculate and calculate estimated time.
This is relevant only when the size option is not specified because
have not figured out how to calculate the number of models of each size
taking into account all of the potential for either/or sets.
However, it could be useful to publish an estimate not taking into account
either/or, just as an estimate, would definitely be closer than the full model space
and more useful than no output? */
local extra_digits 3
local two_minutes 120
local sec_per_min 60
local sec_per_hr 3600
//local sample_pct = `sample_pct'
local dot_interval 1000
if "`size'" == "" {
	if `sample_pct' == 0 {
		local ndigits = ceil(log10(`nmodels')) + `extra_digits'
		di "Calculating " %-`ndigits'.0fc `nmodels' "models..."
		local time_est = `time_avgreg' * `nmodels'
		if `nmodels' > `large_threshold' local large_nmodels 1
	}
	else {
		local nmodels_est = `nmodels'*(`sample_pct'/100)
		local ndigits = ceil(log10(`nmodels_est')) + `extra_digits'
		local time_est = `time_avgreg' * `nmodels_est'
		local nmodels_est = round(`nmodels_est')
		di "Calculating " %-`ndigits'.0fc `nmodels_est' "models..."
		if `nmodels_est' > `large_threshold' local large_nmodels 1
	}
if `time_est' < 10 local time_est 10
else local time_est = `time_est'*3

if `time_est' > `two_minutes' {
		local timeunit min
		local ndigits = ceil(log10(`time_est'/`sec_per_min')) + `extra_digits'
		di "Estimated time is " %-`ndigits'.0fc round(`time_est'/`sec_per_min', 1)  " minutes (" round(`time_est'/`sec_per_hr', .1) " hours)."
}

else {
	local timeunit sec
	di "Estimated time is " round(`time_est', 1)  " seconds (" round(`time_est'/`sec_per_min', .1) " minutes)."
}
}
if `large_nmodels' == 1 di "Each dot represents `dot_interval' models calculated"

// Setup to save the results in a dataset to return to the user.
local postnames i_model str40 (model depvar intvar  )  `allvarnames_clean' str40 (opts) ///
	n r2 df b_intvar ///
	variance pvalue sig pos bic odds_ratio_ind  ///
	`b_always_in_clean'  i_bs

if "`bs_type'" != "" local postnames `"`postnames' b_bs sig_bs pos_bs "'
quietly postfile `memhold' `postnames' using `results', replace

// Setup to save the var combinations in a dataset just for internal use.
quietly postfile `mem_var_combos' `allvarnames_clean' str40 (depvar intvar model) using `tempfile_var_combos', replace

//indicator vector for which sets are present in the model, sets 1 and 2 (depvar and intvar) always in
matrix curr_Sets_ind = J(1, `nSets', 0)
local Sets_always_in 1 2
foreach Set of local Sets_always_in {
	matrix curr_Sets_ind[1, `Set'] = 1
}
local first_rotated_set 3
local nrounds = 2^(`nSets'-2)
local nmodels_actual 0

*********************************************************
// Loop to list out and save the variable combinations
********************************************************
// First, loop over the number of rounds: number of combinations of presence/absence of Sets
forvalues i_round = 1/`nrounds' {
	/* Select the Sets present in this series of models by iterating through
	the binary numbers in the vector curr_Sets_ind */
	forvalues i = `first_rotated_set'/`nSets' {
		if curr_Sets_ind[1, `i'] == 0 {
			matrix curr_Sets_ind[1, `i'] = 1
			continue, break
		}
		else {
			matrix curr_Sets_ind[1, `i'] = 0
		}
	}
	local curr_n_control_terms = 1
	local curr_Setlist "`Sets_always_in'"

	//Track the number of SetX within each of the included sets
	local NSX_by_Set "`s1NSX' `s2NSX'"
	forvalues i_Set = `first_rotated_set'/`nSets' {
		if curr_Sets_ind[1, `i_Set'] == 1 {
			local ++curr_n_control_terms
			local curr_Setlist `"`curr_Setlist' `i_Set'"'
			local NSX_by_Set `NSX_by_Set' `s`i_Set'NSX'
		}
	}

	// Check whether the specified model meets the user-input model size criteria
	if `curr_n_control_terms' >= `sizemin' & `curr_n_control_terms' <= `sizemax' {
		local curr_n_sets = `curr_n_control_terms' + 1
		/* Subroutine to generate the list of all possible combinations of
		either/or subsets of the included sets.*/
		GenerateSetXList,  nSX_by_Set(`NSX_by_Set') setXlist("") count_Xlists(0) ///
			setlist(`curr_Setlist') unit(Set) curr_length(0)
		local count_Xlists = r(count_Xlists)
		forvalues i = 1/`count_Xlists' {
			local SetXlist`i' = r(SetXlist`i')
		}
		clear results
		/* Iterate through each of the SetXlists and generate the list of all
		the terms of each of the included SetX */
		forvalues i_SetXlist = 1/`count_Xlists' {
			local curr_Termlist ""
			forvalues i_Set = 1/`curr_n_sets' {
				local Set : word `i_Set' of `curr_Setlist'
				local SetX : word `i_Set' of `SetXlist`i_SetXlist''
				forvalues i_Term = 1/`s`Set'sx`SetX'NT' {
					local curr_Termlist "`curr_Termlist' s`Set'sx`SetX't`i_Term'"
				}
			}
			//track the number of Termx within each of the included Terms
			local NTX_by_Term ""
			foreach term of local curr_Termlist {
				local NTX_by_Term "`NTX_by_Term' ``term'NTX'"
			}
			/*Generate the list of all possible combinations of
			either/or variables of the included Terms */
			GenerateSetXList, curr_length(0) setlist(`curr_Termlist') ///
			 nSX_by_Set(`NTX_by_Term') setXlist("") count_Xlists(0) unit(Term)
			local count_TermXlists = r(count_Xlists)
			forvalues i = 1/`count_TermXlists' {
				local TermXlist`i' = r(TermXlist`i')
			}
			clear results

			// Save the variable combination to a dataset.
			forvalues i_TermXlist = 1/`count_TermXlists' {
				local curr_varlist ""
				foreach TermX of local TermXlist`i_TermXlist' {
					local curr_varlist "`curr_varlist' ``TermX''"
				}
				local depvar : word 1 of `curr_varlist'
				local intvar : word 2 of `curr_varlist'
				local curr_varlist_ind ""
				foreach varname of local allvarnames {
					local IsIncludedVariable : list varname & curr_varlist
					if "`IsIncludedVariable'" != "" {
						local curr_varlist_ind `curr_varlist_ind' (1)
					}
					else local curr_varlist_ind `curr_varlist_ind' (0)
				}
				forvalues i_model_type = 1/`nmodeltypes' {
					quietly post `mem_var_combos' `curr_varlist_ind' ("`depvar'") ("`intvar'") ("`model`i_model_type'_id'")
				}
			}
		}
	}
}
postclose `mem_var_combos'

// Translate var combos to a stata matrix to iterate in the next loop
preserve
use `tempfile_var_combos', clear
quietly duplicates drop
quietly sample `sample_pct'
quietly gen i_model = _n
qui save `tempfile_var_combos', replace
local n_vars_considered : word count `allvarnames'
local n_var_combos = `c(N)'

mata {
	n_vars_considered = strtoreal(st_local("n_vars_considered"))
	var_combos = st_data(.,1..n_vars_considered)
	key_vars = st_sdata(., ("depvar", "intvar", "model"))
}
restore

*********************************************************
// Main loop for estimation of all the models
*********************************************************
// Loop over the var combos
forvalues i_var_combo = 1/`n_var_combos' {

	/*
	//Sample barrier
	local PassedSampleThreshold 0
	if `threshold' == 0 local PassedSampleThreshold 1
	else {
		local r = runiform()
		if `r' <= `threshold' local PassedSampleThreshold 1
	}
	if `PassedSampleThreshold' == 1 {
	*/
		// identify the model type (reg, logit, etc)
		mata: st_local("model_id", key_vars[strtoreal(st_local("i_var_combo")), 3])
		forvalues i_model_type = 1/`nmodeltypes' {
			if "`model_id'" == "`model`i_model_type'_id'" {
				local model `model`i_model_type'_name'
				local opts_touse `"`opts_command' `model`i_model_type'_opts'"'
				local odds_ratio_ind 0
				local irr irr
				if "`model'" == "logistic" | `: list irr in opts_touse' local odds_ratio_ind 1
			}
		}

		//Identify the variable names included in this model.

		//mata: st_matrix("curr_row", var_combos[`i_var_combo',])
		//mata: st_local("curr_depvar", var_combos[`i_var_combo',
		//matlist curr_row

		local curr_varlist ""
		forvalues i_var = 1/`n_vars_considered' {
			//mata: var_combos[1,1] // column : this var name

			mata: st_local("var_in_model", /// assign to this local var
				 strofreal(var_combos[strtoreal(st_local("i_var_combo")), /// row : this var combo
						   strtoreal(st_local("i_var"))] )) // column : this var name
			if  "`var_in_model'" == "1" { // build up the list of var names in this model.
				local curr_varlist "`curr_varlist' `: word `i_var' of `allvarnames''"
			}
		}
		local ++nmodels_actual

		if "`savelist'" != "" {
			file write listfile _n "* `nmodels_actual'"
			file write listfile _n `"`model' `curr_varlist'"'
			// note: for some reason, would not write `touse' but would
			// evaluate `touse' to __0000004 if the following was all on one line
			if !("`listwise'" == "nolistwise" & "`if'" == "" & "`in'" == "") {
				file write listfile " if \`touse'"
			}
			if (trim("`opts_touse'") != "") file write listfile `", `opts_touse'"'

		}

		mata: st_local("depvar", key_vars[strtoreal(st_local("i_var_combo")), 1])
		mata: st_local("intvar", key_vars[strtoreal(st_local("i_var_combo")), 2])

		if strmatch("`intvar'", "i.*") {
			local intvar_for_results : subinstr local intvar "i." "1."
		}
		else {
			local intvar_for_results "`intvar'"
		}

		/*
		local depvar : word 1 of `curr_varlist'
		local intvar : word 2 of `curr_varlist'
		*/

		//Estimate the model
		if "`bs_type'" == "nonpar" quietly bootstrap `bs_opts_all' ///
			saving(`temp_nonpar_bs_results', replace): ///
			`model' `curr_varlist' if `touse', `opts_touse'

		else quietly `model' `curr_varlist' if `touse', `opts_touse'
		mat results_table = r(table)
		// Calculate information criteria including BIC
		// note that estat ic does not work for rreg
		capture estat ic
		if "`_rc'" == "0"  {
			local is_bic_ok 1
		}
		else {
			local is_bic_ok 0
		}

		//Save the results
		if "`model'" == "xtreg"	{
			if ("`be'" == "be"| "`fe'" == "fe"  | "`re'" == "re") local e_r2 = e(r2_o)
			else local e_r2 = .
		}
		else if inlist("`model'", "regress", "areg", "rreg") local e_r2 = e(r2)
		else local e_r2 = e(r2_p)

		if "`margins'" != "margins" {
			local b = _b[`intvar_for_results']
			local variance = (_se[`intvar_for_results'])^2
			local pval = results_table[rownumb(results_table, "pvalue"), colnumb(results_table, "`intvar_for_results'")]
		}
		else if "`margins'" == "margins" {

			capture margins , dydx(`intvar_for_results') `marginsopts_post'
			local marginsrc = _rc
			if `marginsrc' != 0 {
				di _newline
				di as error "Error with -margins- command."
				di as result "Last estimation result completed:"
				estimates replay
				di _newline
				di as result  "-margins- command that produced the error:"
				di as error `"margins, dydx(`intvar_for_results') `marginsopts_post' "'
				exit `marginsrc'
			}
			mat margins_table = r(table)
			local b = margins_table[rownumb(margins_table, "b"), colnumb(margins_table, "`intvar_for_results'")]
			local pval = margins_table[rownumb(margins_table, "pvalue"), colnumb(margins_table, "`intvar_for_results'")]
			local variance = margins_table[rownumb(margins_table, "se"), colnumb(margins_table, "`intvar_for_results'")]^2

			if "`savelist'" != "" file write listfile _n `"margins, dydx(`intvar_for_results') `marginsopts_post' "'

		}
		local always_in_estimates ""
		foreach var of local always_in_varnames {
			capture local b_var = _b[`var']
			if _rc == 0 local always_in_estimates `always_in_estimates' (`b_var')
			else local always_in_estimates `always_in_estimates' (.)
		}


		/*
		if inlist("`e(model)'", "ols", "fe") & "`bs_type'" != "nonpar" ///
			local pval = (2*ttail(e(df_r), abs(_b[`intvar_for_results']/_se[`intvar_for_results'])))
		else local pval = 2*(1-normal(abs(_b[`intvar_for_results']/_se[`intvar_for_results'])))
		*/
		local sig = (`pval' <= `alpha')
		local pos = (`b' > 0)
		//Save the BIC
		matrix ic_mat = r(S)
		local bic = ic_mat[1,6]

		// Keep track of sample size range
		// useful when listwise delete is off and sample size varies across models
		local sampsize = e(N)
		if "`cur_max_sampsize'" != "" {
			if `sampsize' > `cur_max_sampsize' {
				local cur_max_sampsize = `sampsize'
			}
		}
		else {
			local cur_max_sampsize = `sampsize'
		}
		if "`cur_min_sampsize'" != "" {
			if `sampsize' < `cur_min_sampsize' {
				local cur_min_sampsize = `sampsize'
			}
		}
		else {
			local cur_min_sampsize = `sampsize'
		}
		// Prepare to post results in the results dataset.
		local topost (`nmodels_actual') ("`model_id'") ("`depvar'") ("`intvar'") ///
			`curr_varlist_ind' ("`model`i_model_type'_opts'")  ///
			(`sampsize') (`e_r2') (e(df_m)) (`b') ///
			(`variance') (`pval')  (`sig') (`pos') (`bic') (`odds_ratio_ind') ///
			`always_in_estimates'


		/* If user specified nonparametric bootstrap then open the saved
		results from the bootstrap routine and store them in the model
		results file */
		if "`bs_type'" == "nonpar" {
			capture use `temp_nonpar_bs_results', clear
			if inlist("`model'", "probit", "logit", "logistic", "poisson") ///
				local _bvar `depvar'_b_`intvar'
			else local _bvar _b_`intvar'
			local n_bsreps = _N
			forvalues i_bs = 1/`n_bsreps' {
				if 2*(1-normal(abs(`_bvar'[`i_bs']/_se[`intvar']))) <= `alpha' ///
					local sig_bs 1
				else local sig_bs 0
				local pos_bs = ( `_bvar'[`i_bs'] > 0 )
				post `memhold' `topost' (`i_bs') (`_bvar'[`i_bs']) (`sig_bs') (`pos_bs')
			}
			capture use `orig_touse', clear
		}
		/* Otherwise resample the coefficients with parametric bootstrap
		according to the normal distribution */
		else if "`bs_type'" == "par"{
			forvalues i_bs = 1/`bs_reps' {
				local bs_sample = rnormal(_b[`intvar'], _se[`intvar'])
				local pos_bs = ( `bs_sample' > 0 )
				if inlist("`e(model)'", "ols", "fe") {
					if (2*ttail(e(df_r), abs(`bs_sample'/_se[`intvar']))) <= `alpha' ///
						local sig_bs 1
					else local sig_bs 0
				}
				else if 2*(1-normal(abs(`bs_sample'/_se[`intvar']))) <= `alpha' ///
					local sig_bs 1
				else local sig_bs 0
				post `memhold' `topost'  (`i_bs') (`bs_sample') (`sig_bs') (`pos_bs')
			}
		}
		else post `memhold' `topost'  (1)

		// Display progress for user.
		if `large_nmodels' == 1 {
			if mod(`nmodels_actual', `dot_interval') == 0 {
				di _continue "."
			}
		}

	//} // if sample threshold is met
} // loop over var combos

postclose `memhold'
if "`savelist'" != "" quietly file close listfile


//Merge the 2 datasets: results and var combos
preserve
use `tempfile_var_combos', clear
merge 1:1 i_model depvar intvar model using `results', nogenerate noreport
qui save `results', replace
restore

/*calculate multicollinearity of var of interest
if there are >1 options for intvar, take the mean of the MCs*/
local mcsum = 0
local nintvar : word count `intvarlist_noline'
local xvars : list allvarnames - depvarlist_noline
local xvars : list xvars - intvarlist_noline
foreach intvar of local intvarlist_noline {

	if strmatch("`intvar'", "i.*") {
		local intvar_for_multicol : subinstr local intvar "i." ""
	}
	else {
		local intvar_for_multicol "`intvar'"
	}

	quietly reg `intvar_for_multicol' `xvars' if `touse'
	local mcsum = `mcsum' + e(r2)
}
return scalar mc = `mcsum'/`nintvar'

/*calculate the SD of the variable of interest
if there are >1 options for intvar, take the mean of the SDs
*/
local sd2intvar 0
foreach intvar of local intvarlist_noline {
	quietly summarize `intvar' if `touse'
	local sd2intvar = `sd2intvar' + (r(sd)^2)
}
return scalar sdintvar = sqrt(`sd2intvar')
return local bs_type `"`bs_type'"'
return local bs_opts `"`bs_options_user'"'
return scalar bs_reps = `bs_reps'
return scalar alpha = `alpha'
if "`pref'" != "" {
	return scalar prefb = `prefb'
	return scalar prefse = `prefse'
}
return local sig_only "`sig_only'"
return local savelist "`savelist'"
return local opts_command `"`opts_command'"'
return local opts_multivrs `"`opts_multivrs'"'
forvalues im = 1/`nmodeltypes' {
	return local model`im'_name "`model`im'_name'"
	return local model`im'_opts "`model`im'_opts'"
	return local model`im'_id "`model`im'_id'"

}
return local intvar `"`intvarlist'"'
return local depvar `"`depvarlist'"'
return local model `"`model_namelist'"'
return local cmd multivrs
return scalar nmodels = `nmodels_actual'
return scalar nterms = `nrotated'
return scalar nvars = `: word count  `allvarnames'' - `: word count `depvarlist_noline''
return scalar N = `N'
return scalar Nmax = `cur_max_sampsize'
return scalar Nmin = `cur_min_sampsize'
return local weights `"`weights'"'
return local margins `margins'
return local marginsopts `marginsopts'

if "`listwise'" == "nolistwise" {
	return scalar listwise_delete = 0
}
else {
	return scalar listwise_delete = 1
}

end

// End program Multivrs_estimate

/*
Recursive branching subroutine to generate the list of all possible combinations of
of either/or SetX units of the included Sets.
We want exactly one SetX for each of the included n Sets.  The base case is we have
assembled a SetXlist to return containing n SetX ( 1 SetX for each Set ).
If the assembled SetXlist contains i SetX, i < n, then for each of the SetX j
in the next Set i+1, we add j to the SetXlist and pass this incremented SetXList
of length i+1 back into GenerateSetXList.
Arguments:
curr_length length of current SetXList to return
Setlist is the list of Sets included
NSX_by_set lists the number of setx in each of the included sets
SetXlist lists the current assembled setxlist to return
count_SetXlist holds the total number of generated setxlists so far
unit indicates whether we are dealing with Sets or Terms
Return values
r(count_Xlists) is the total number of distinct SetXlists generated
for values i = 1 to r(count_Xlists),
r(SetXlist`i') is a string indicating the components of the ith SetXlist
*/

//program GenerateUnitList

program GenerateSetXList, rclass
syntax [, setlist(string) nSX_by_Set(numlist) setXlist(string) count_Xlists(integer 0)  unit(string) curr_length(integer 0) ]

	if `curr_length' == `: word count `setlist'' {
		local ++count_Xlists
		return local count_Xlists `count_Xlists'
		return local `unit'Xlist`count_Xlists' `setXlist'
	}
	else {
		local ++curr_length
		local curr_set : word `curr_length' of `setlist'
		local nSX_this_set : word `curr_length' of `nSX_by_Set'
		forvalues i = 1/`nSX_this_set' {
			if "`unit'" == "Set" local curr_SetXlist "`setXlist' `i'"
			else local curr_SetXlist "`setXlist' `curr_set'tx`i'"
			local count = r(count_Xlists)
			if "`count'" == "."  local count = 0
			GenerateSetXList, curr_length(`curr_length') setlist("`setlist'") nSX_by_Set(`nSX_by_Set') setXlist(`curr_SetXlist') count_Xlists(`count') unit(`unit')
			foreach mac in `:r(macros)' {
				return local `mac' `=r(`mac')'
			}
		}
	}
end