/*
Use the saved results file to calculate the summary statistics.
*/
program multivrs_calc_summary_stats, rclass
syntax [, saveas(string) replace noinfluencecalcs inf_means margins]
local model_namelist `"`r(model)'"'
local model_idlist `"`r(model_idlist)'"'
local nmodeltypes : word count `model_namelist'
local depvarlist `"`r(depvar)'"'
local intvarlist `"`r(intvar)'"'
local intvarlist_noline : subinstr local intvarlist "|" " ", all
local nintvar : word count `intvarlist_noline'
local depvarlist_noline : subinstr local depvarlist "|" " ", all
local ndepvar : word count `depvarlist_noline'
local bs_type `"`r(bs_type)'"'
local bs_opts `"`r(bs_opts)'"'
local nterms `"`r(nterms)'"'
local N `"`r(N)'"'
local nmodels `"`r(nmodels)'"'
local prefb `"`r(prefb)'"'
local prefse `"`r(prefse)'"'
local alpha `"`r(alpha)'"'
local opts_command `"`r(opts_command)'"'
local odds_ratio = r(odds_ratio)
local weights `"`r(weights)'"'
forvalues im = 1/`nmodeltypes' {
local model`im'_name `"`r(model`im'_name)'"'
local model`im'_opts `"`r(model`im'_opts)'"'
}
return add
local results_to_return totalParLB totalParUB extremeLB extremeUB ///
modeling95LB modeling95UB ///
p1 p5 p10 p90 p95 p99 kurtosis skewness bic_r2_cor ///
rratio totalSE samplingSE modelingSE /*modelingSE_biased modelingSE_unbiased_check*/ meanb meanr2
if "`bs_type'" != "" local results_to_return total95LB total95UB `results_to_return'
if "`prefb'" != "" local results_to_return prefLB prefUB prefpctile preftotalSE ///
`results_to_return'
tempname `results_to_return'
qui gen model_id = _n
qui generate se = sqrt(variance)
qui gen s2 = se^2
* CALCULATE INFLUENCE STATISTICS *
if "`influencecalcs'" != "noinfluencecalcs" {
GenerateInfStatsMatrix, modellist("`model_idlist'") depvarlist("`depvarlist_noline'") ///
intvarlist("`intvarlist_noline'") `inf_means'
matrix inf_coef = r(infstats_default)
return add
}
* CALCULATE WEIGHTS *
/* if weight = "no" then don't calculate any of them, but if it's
a specific option (or it's not specified, the default equal weighting),
then go ahead and calculate all of them */
if ("`weights'" == "no" | "`weights'" == "uniform") local weights_option ""
else local weights_option "[aweight=wt_`weights']"
if "`weights'" != "no" {
if "`influencecalcs'" != "noinfluencecalcs" {
* Influence Weights: Proportional to absolute value of influence regression coefficient.
* Note that the coefficients measure the effect of binary indicators for the presence
* of the specified variable in the model. So they're on the same scale.
/* Each rotated variable has a weight, and the model weight is the sum of
the weights of the included variables. */
local num = r(nterms)
mata: st_matrix("inf_sums", colsum(st_matrix("inf_coef")))
mata: st_matrix("weight_matrix", st_matrix("inf_coef")[,1]/st_matrix("inf_sums")[1,1])
local names : rownames inf_coef
matrix rownames weight_matrix = `names'
* Sum the weights of variables included in the model to arrive at the raw model weight.
local n = 0
foreach variable of local names {
local ++n
local weight = weight_matrix[`n',1]
qui gen wt_inf_`variable' = `variable'*`weight'
}
egen wt_inf_raw = rowtotal(wt_inf_*)
* Divide by the sum of all raw model weights to calculate weights that sum to 1 (wt_inf).
qui sum wt_inf_raw if i_bs == 1
gen wt_inf = wt_inf_raw/`r(sum)'
drop wt_inf_raw
}
* BIC Weights: Approximation to posterior probabilities.
/* Proportional to exp((-1/2)*BICdelta) where BICdelta is the difference between the model's BIC
and the minimum BIC over all the models considered */
qui sum bic if i_bs == 1
qui gen bic_min = r(min)
qui gen bic_delta = bic - bic_min if i_bs == 1
qui gen wt_bic_raw = exp(-0.5*bic_delta) if i_bs == 1
qui egen wt_bic_raw_total = total(wt_bic_raw) if i_bs == 1
qui gen wt_bic = wt_bic_raw / wt_bic_raw_total
* R2 Weights: Simply proportional to the R2 value.
qui sum r2 if i_bs == 1
qui gen wt_r2_total = r(sum)
qui gen wt_r2 = r2 / wt_r2_total
drop bic_min bic_delta wt_bic_raw* wt_r2_total
}
* (weighted) SUMMARY OF POINT ESTIMATES: MODELING DISTRIBUTION
quietly sum b_intvar `weights_option' if i_bs == 1, detail
* Return list at this point is (weighted) summary of b_intvar.
* The following calculations all rely on this summary.
* (weighted) MEAN & VARIANCE
scalar `meanb' = r(mean)
if ("`weights'" == "no" | "`weights'" == "uniform") {
scalar `modelingSE' = r(sd)*sqrt((_N-1)/_N)
}
else {
* Calculate unbiased (weighted) observed variance across estimates
* from https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#cite_note-3
qui gen wt_dif2 = wt_`weights'*(b_intvar - r(mean))^2
qui sum wt_`weights', meanonly
local v1 = r(sum)
qui gen wt_`weights'2 = wt_`weights'^2
qui sum wt_`weights'2
local v2 = r(sum)
qui sum wt_dif2
scalar `modelingSE' = sqrt(r(sum)/(`v1' - (`v2'/`v1')))
drop wt_dif2 wt_`weights'2
}
quietly sum b_intvar `weights_option' if i_bs == 1, detail
* (weighted) SKEW & KURTOSIS
scalar `skewness' = r(skewness)
scalar `kurtosis' = r(kurtosis)
* (unaffected by weights) MIN and MAX
scalar `extremeLB' = r(min)
scalar `extremeUB' = r(max)
* (unweighted) PERCENTILES, including EMPIRICAL 95% INTERVAL
local percentiles 1 2.5 5 10 90 95 97.5 99
quietly centile b_intvar if i_bs == 1, centile(`percentiles')
local npercentiles : word count `percentiles'
forvalues i = 1/`npercentiles' {
local p : word `i' of `percentiles'
if `p' == 2.5 scalar `modeling95LB' = r(c_`i')
else if `p' == 97.5 scalar `modeling95UB' = r(c_`i')
else scalar `p`p'' = r(c_`i')
}
* (weighted) AVERAGE R2
quietly summarize r2 `weights_option' if i_bs == 1
scalar `meanr2' = r(mean)
* (weighted) AVERAGE SAMPLING VARIANCE
quietly sum variance `weights_option' if i_bs == 1, meanonly
scalar `samplingSE' = sqrt(r(mean))
* DERIVED SUMMARY STATISTICS
if "`bs_type'" != "" {
quietly summarize b_bs
scalar `totalSE' = r(sd)
quietly centile b_bs, centile(2.5 97.5)
scalar `total95LB' = r(c_1)
scalar `total95UB' = r(c_2)
}
else scalar `totalSE' = sqrt(`modelingSE'^2 + `samplingSE'^2)
scalar `rratio' = `meanb'/`totalSE'
scalar `totalParLB' = `meanb' - 2*`totalSE'
scalar `totalParUB' = `meanb' + 2*`totalSE'
* (unweighted) BIC, R2 CORRELATION
qui sum bic if i_bs == 1
if `r(N)' > 0 {
capture cor bic r2 if i_bs == 1
scalar `bic_r2_cor' = r(rho)
}
else {
scalar `bic_r2_cor' = .
}
* STATISTICS FOR PREFERRED ESTIMATE
/* percentile of modeling distribution,
total standard error (= sqrt of sum of squares of preferred SE and modeling SE),
semi-parametric confidence interval
Note: because the raw unweighted percentile is calculated here, it may be
inconsistent with the weighted percentiles calculated above */
if "`prefb'" != "" {
if `odds_ratio' == 1 local prefb = log(`prefb')
quietly count if b_intvar <= `prefb' & i_bs == 1
scalar `prefpctile' = 100 * (r(N)/`nmodels')
if `prefpctile' == 0 | `prefpctile' >= 100 scalar `prefpctile' = .
scalar `preftotalSE' = sqrt(`prefse'^2 + `modelingSE'^2)
scalar `prefLB' = `prefb' - 2*`preftotalSE'
scalar `prefUB' = `prefb' + 2*`preftotalSE'
scalar `rratio' = `prefb'/`preftotalSE'
}
GenerateSigRatesMatrix, model_namelist("`model_idlist'") depvarlist("`depvarlist_noline'") ///
intvarlist("`intvarlist_noline'") alpha(`alpha')
return add
drop variance
// Case: everything uses the odds ratio (not log-OR coefficient) -->
// Exponentiate the relevant values and multiply the SEs
// See https://www.stata.com/manuals/rlogistic.pdf
quietly replace b_intvar= exp(b_intvar) if odds_ratio_ind == 1 & "`margins'" != "margins"
quietly replace se = b_intvar*se if odds_ratio_ind == 1 & "`margins'" != "margins"
quietly count if odds_ratio_ind == 0
if r(N) == 0 {
local to_exp meanb totalParLB totalParUB modeling95LB modeling95UB extremeLB extremeUB
if "`bs_type'" != "" local to_exp `"`to_exp' total95LB total95UB "'
if "`prefb'" != "" local to_exp `"`to_exp' prefLB prefUB"'
foreach t of local to_exp {
scalar ``t'' = exp(``t'')
}
local to_multiply samplingSE modelingSE totalSE
foreach t of local to_multiply {
scalar ``t'' = `meanb'*``t''
}
}
if "`bs_type'" == "" {
//quietly drop if i_bs != 1
quietly drop i_model i_bs
}
if `"`saveas'"' != "" quietly save `"`saveas'"' , `replace'
return local saveas `"`saveas'"'
//Return the regression type for display
local reg_type ""
forvalues im = 1/`nmodeltypes' {
local m = "`model`im'_name'"
if `im' > 1 local reg_type "`reg_type'/"
if ("`m'" == "logit") local reg_type "`reg_type'Logit"
else if ("`m'" == "logistic") local reg_type "`reg_type'Logistic"
else if ("`m'" == "probit") local reg_type "`reg_type'Probit"
else if ("`m'" == "poisson") local reg_type "`reg_type'Poisson"
else if ("`m'" == "nbreg") local reg_type "`reg_type'Negative Binomial"
else if ("`m'" == "areg") local reg_type "`reg_type'Absorbing"
else if ("`m'" == "rreg") local reg_type "`reg_type'Robust"
else if ("`m'" == "xtreg") {
local m_opts = `"`opts_command' `opts`im'' "'
local fe fe
local be be
if `: list fe in m_opts' == 1 local reg_type "`reg_type'Fixed-effects (within)"
else if `: list be in m_opts' == 1 local reg_type "`reg_type'Between (group means)"
else local reg_type "`reg_type'Random-effects GLS"
}
else local reg_type "`reg_type'Linear"
}
local reg_type "`reg_type' regression"
local reg_type : list retok reg_type
return local title `"`reg_type'"'
foreach r of local results_to_return {
return scalar `r' = ``r''
}
end
// End program CalcSummaryStats
/*
GenerateInfStatsMatrix
*/
program GenerateInfStatsMatrix, rclass
syntax , modellist(string) depvarlist(string) intvarlist(string) [inf_means]
quietly sum b_intvar if i_bs == 1
local meanb = r(mean)
/*Build the varlist for the influence regression, including the rotated control
variables and dummies for the functional form, depvar, and intvar.
*/
local inf_varnames ""
local inf_references ""
local groupleads ""
capture unab rotatedvars : r_*
local typelist model depvar intvar
foreach type of local typelist {
local n`type' : word count ``type'list'
if `n`type'' > 1 {
forvalues i = 1/`n`type'' {
pause
local curr : word `i' of ``type'list'
quietly gen byte `curr' = (`type' == "`curr'")
if `i' == 1 {
local inf_references `"`inf_references' `curr'"'
}
if `i' >= 2 {
local inf_varnames `"`inf_varnames' `curr'"'
}
}
}
}
// Always-In variables
unab alwaysinlist0 : __*
local alwaysinlist ""
local nalwaysin : word count `alwaysinlist0'
local dep_and_intvarlist `depvarlist' `intvarlist'
if `nalwaysin' > 1 {
forvalues i = 1/`nalwaysin' {
local curr : word `i' of `alwaysinlist0'
local b_curr = regexr("`curr'", "^__", "b_")
local curr_varname = regexr("`curr'", "^__", "")
if !`: list curr_varname in dep_and_intvarlist' {
capture confirm variable `curr'
if _rc == 0 {
local alwaysinlist `"`alwaysinlist' `curr'"'
}
}
}
}
local reference_grouplead `: word 1 of `alwaysinlist''
local inf_references "`inf_references' `reference_grouplead'"
foreach var1 of local alwaysinlist {
local part_of_group1 0
if "`var1'" != "`reference_grouplead'" {
quietly count if `var1' != `reference_grouplead'
if r(N) == 0 {
local part_of_group1 1
continue, break
}
if !`part_of_group1' {
local inf_varnames `"`inf_varnames' `var1'"'
}
}
}
local inf_varnames `"`inf_varnames' `rotatedvars'"'
/*
Run the regression for the marginal effect of variable inclusion
on the estimate. Store results.
*/
quietly reg b_intvar `inf_varnames' if i_bs == 1
return scalar infr2_b = e(r2)
return scalar infcons_b = _b[_cons]
local ngroups = e(df_m)
/*
Create the 2 matrices: infstats_default, which stores the coefficients from the
OLS influence regression, indicating the effect of inclusion of each variable
or functional form on the magnitude of the estimate, and infstats_sig_only, which
stores the odds ratios from the logistic influence regression, indicating the
effect of inclusion of each variable or functional form on the
likelihood of a positive or significant estimate.
*/
tempname inf_mem_default inf_mem_sig inf_mem_pos
tempfile inf_file_default inf_file_sig inf_file_pos
local postnames_default str32 var abs_coef coef pct_chg mean_b row_id
quietly postfile `inf_mem_default' `postnames_default' using `inf_file_default', replace
local postnames_sig str32 var abs_coef_sig coef_sig row_id_sig
quietly postfile `inf_mem_sig' `postnames_sig' using `inf_file_sig', replace
local postnames_pos str32 var coef_pos
quietly postfile `inf_mem_pos' `postnames_pos' using `inf_file_pos', replace
tempname infstats_default infstats_sig_only
if "`inf_means'" == "inf_means" {
local colnames_default abs_coef coef pct_chg mean_b row_id
}
else {
local colnames_default abs_coef coef pct_chg row_id
}
local colnames_sig_only abs_coef_sig coef_sig coef_pos row_id_sig
local ncols_default : word count `colnames_default'
local ncols_sig_only : word count `colnames_sig_only'
matrix `infstats_default' = J(`ngroups', `ncols_default', .)
matrix `infstats_sig_only' = J(`ngroups', `ncols_sig_only', .)
/*If the term has meaningful influence coefficient then add it to the matrix.
otherwise, check which variable it is grouped with and it to that group
*/
local lists colnames_default colnames_sig_only
foreach namelist of local lists {
local i_row 0
foreach name of local `namelist' {
local ++i_row
local `name' = `i_row'
}
}
local i_row 0
local b_colnames : colnames e(b)
qui gen double b_with_var = .
foreach var of local inf_varnames {
local omitted = regexm("`b_colnames'", "(^| )o\.`var'($| )")
local part_of_group 0
if "`omitted'" == "0" { // variable not omitted
local groupleads `"`groupleads' `var'"'
local group_`var' ""
qui replace b_with_var = b_intvar if `var' == 1
sum b_with_var, meanonly
qui replace b_with_var = .
quietly post `inf_mem_default' ("`var'") (abs(_b[`var'])) (_b[`var']) (100*(_b[`var']/`meanb')) (`r(mean)') (`i_row')
}
else if "`omitted'" == "1" { //variable omitted due to collinearity; it was either part of a group or a reference category
foreach groupvar of local groupleads {
quietly count if `var' != `groupvar'
if r(N) == 0 { // this variable is part of a group with some other variable recorded in groupleads
local part_of_group 1
mata: st_local("displayname", asarray(multivrs_varnames, "`var'"))
local group_`groupvar' `"`group_`groupvar'' `displayname'"'
continue, break
}
}
}
}
if "`inf_means'" == "inf_means" {
foreach var of local inf_references {
qui sum b_intvar if `var' == 1
qui replace b_with_var = b_intvar if `var' == 1
qui sum b_with_var, meanonly
qui replace b_with_var = .
quietly post `inf_mem_default' ("`var'") (.) (.) (.) (`r(mean)') (.)
}
}
postclose `inf_mem_default'
/*
Run the regressions for marginal effect of variable inclusion
on the probability of significant or positive estimate.
Store results in the matrix.
*/
local dummies sig pos
foreach d of local dummies {
capture reg `d' `inf_varnames' if i_bs == 1
local i_row 0
foreach var of local groupleads {
local ++i_row
//matrix `infstats_sig_only'[`i_row', `coef_`d''] = _b[`var']
if "`d'" == "sig" {
//matrix `infstats_sig_only'[`i_row', `abs_coef_sig'] = abs(`infstats_sig_only'[`i_row', `coef_sig'])
//matrix `infstats_sig_only'[`i_row',`row_id_sig'] = `i_row'
if _rc == 0 qui post `inf_mem_sig' ("`var'") (abs(_b[`var'])) (_b[`var']) (`i_row')
else qui post `inf_mem_sig' ("`var'") (.) (.) (`i_row')
}
else if "`d'" == "pos" {
if _rc == 0 qui post `inf_mem_pos' ("`var'") (_b[`var'])
else qui post `inf_mem_pos' ("`var'") (.)
}
}
if _rc == 0 return scalar infr2_`d' = e(r2)
else return scalar infr2_`d' = .
if _rc == 0 return scalar infcons_`d' = _b[_cons]
else return scalar infcons_`d' = .
return scalar `d'_rc = _rc
}
postclose `inf_mem_sig'
postclose `inf_mem_pos'
/*
Sort the rows of each infstats matrix in descending order of absolute value of the
influence coefficient. Iterate through the sorted row_ids and attach the rownames to matrices
in the proper order
*/
preserve
use `inf_file_default', clear
if "`inf_means'" == "inf_means" {
qui gen abs_mean_b = abs(mean_b)
gsort -abs_mean_b
qui drop abs_mean_b
}
else {
gsort -abs_coef
}
mata: infstats_default = st_data(.,(2,3,4,5,6))
mata: st_matrix("`infstats_default'", infstats_default)
qui valuesof var
matrix rownames `infstats_default' = `r(values)'
qui ds, has(type numeric)
matrix colnames `infstats_default' = `r(varlist)'
return matrix infstats_default = `infstats_default'
restore
// sig_only
preserve
use `inf_file_pos', clear
merge 1:1 var using `inf_file_sig', nogenerate noreport
gsort -abs_coef_sig
mata: infstats_sig_only = st_data(.,("abs_coef_sig", "coef_sig", "coef_pos"))
mata: st_matrix("`infstats_sig_only'", infstats_sig_only)
qui valuesof var
matrix rownames `infstats_sig_only' = `r(values)'
qui ds, has(type numeric)
matrix colnames `infstats_sig_only' = abs_coef_sig coef_sig coef_pos
return matrix infstats_sig_only = `infstats_sig_only'
restore
foreach var of local groupleads {
return local group_`var' `"`group_`var''"'
}
end
/*
GenerateSigRatesMatrix calculates the rates of significant and
positive coefficients for each functional form, depvar, and intvar, as well as
the pooled rates, stores these numbers in the matrix sigrates and returns the
matrix sigrates.
*/
program GenerateSigRatesMatrix, rclass
syntax , model_namelist(string) depvarlist(string) intvarlist(string) alpha(real)
//Create sigrates matrix
local lists model_namelist depvarlist intvarlist
local ncols 0
local colnames "all_models"
foreach l of local lists {
if `: word count ``l''' > 1 {
foreach entry of local `l' {
local colnames `colnames' `entry'
local ++ncols
}
}
}
local ++ncols
local rownames NModels SignStability SignificanceRate Positive PositiveandSig Negative NegativeandSig
local nrows : word count `rownames'
matrix sigrates = J(`nrows', `ncols', 0)
matrix colnames sigrates = `colnames'
matrix rownames sigrates = `rownames'
local i_row = 0
foreach rowname of local rownames {
local ++i_row
local `rowname' = `i_row'
}
//Calculate significance rates and store in matrix sigrates
forvalues i_col = 1/`ncols' {
local curr_colname : word `i_col' of `colnames'
if `:list curr_colname in model_namelist' == 1 local type model
else if `:list curr_colname in depvarlist' == 1 local type depvar
else local type intvar
local ifcond ""
if `i_col' > 1 local ifcond `"& `type' == "`curr_colname'""'
quietly count if i_bs == 1 `ifcond'
local curr_nmodels = r(N)
matrix sigrates[`NModels', `i_col'] = `curr_nmodels'
quietly count if pvalue <= `alpha' & i_bs == 1 `ifcond'
matrix sigrates[`SignificanceRate', `i_col'] = 100*(r(N)/`curr_nmodels')
quietly count if b_intvar> 0 & i_bs == 1 `ifcond'
matrix sigrates[`Positive',`i_col'] = 100*(r(N)/`curr_nmodels')
quietly count if b_intvar> 0 & pvalue <= `alpha' & i_bs == 1 `ifcond'
matrix sigrates[`PositiveandSig',`i_col'] = 100*(r(N)/`curr_nmodels')
quietly count if b_intvar< 0 & i_bs == 1`ifcond'
matrix sigrates[`Negative', `i_col'] = 100*(r(N)/`curr_nmodels')
quietly count if b_intvar < 0 & pvalue <= `alpha' & i_bs == 1 `ifcond'
matrix sigrates[`NegativeandSig',`i_col'] = 100*(r(N)/`curr_nmodels')
matrix sigrates[`SignStability',`i_col'] = max(sigrates[`Positive',`i_col'], sigrates[`Negative', `i_col'])
}
return matrix sigrates = sigrates
end
/*
* from ssc install valuesof
*! version 1.0.2 Ben Jann 20may2008
program valuesof
version 9.1, born(20Jan2006)
capt mata mata which mm_invtokens()
if _rc {
di as error "mm_invtokens() from -moremata- is required; type {stata ssc install moremata}"
error 499
}
syntax varname [if] [in] [, Format(str) MISSing ]
mata: st_rclear()
if `"`missing'"'=="" marksample touse, strok
else marksample touse, strok nov
capt confirm string var `varlist'
if _rc {
local dp = c(dp)
set dp period
if `"`format'"'=="" local format "%18.0g"
else confirm format `format'
mata: st_global("r(values)", ///
mm_invtokens(strofreal(st_data(.,"`varlist'","`touse'")', `"`format'"')))
set dp `dp'
}
else {
mata: st_global("r(values)", ///
mm_invtokens(st_sdata(.,"`varlist'","`touse'")'))
}
di `"`r(values)'"'
end
*/