cap prog drop baltable
*! program -baltable- constructs a balance table by assembling the results of multiple -ttest- executions
*! version 1.4 25Apr2024 by Mead Over
program define baltable, rclass

	version 17
	syntax varlist(fv min=1 numeric) [if] [in] [using/] ,  ///
		BY(varname numeric)                                               /// Treatment variable, 2 values
		[EXPortopts(string)                                               /// Options for exporting to `using'
		UNEQual HETerogeneous nojoint                                     /// Unequal variances
		noLIst FVExpand                                                   /// Factor variables
		noTOTals nostddev nostderr nonobs nostddiff notstat ponly         /// Output simplifying options
		SPACEr split verbose extraspace(passthru) title(string)           /// Reporting options
		noGRPLabels noVARLabels                                           /// Group value label & Variable label options
		WRAPLength(int 25) wrap(int 10) wrapon(string)                    /// Variable label stacking options
		format(string) nformat(string) sformat(string) pformat(string)    /// Numerical formatting options
		p_upper p_lower level(passthru)                                   /// One-sided tests
		TRUNCGrplbl(int 15) saving(passthru) ]                            //  Group label length, collection saving

	marksample touse
	markout `by'

	preserve

	if "`verbose'"~="verbose" {
		local shhh quietly
	}
	else {
		local noisily noisily
	}

	if "`unequal'"=="unequal" {
		local heterogeneous heterogeneous
	}
	if "`heterogeneous'"=="heterogeneous" {
		local unequal unequal
	}

	`shhh' keep if `touse'

	if "`split'"=="" {
		local stack stack
	}

	if "`ponly'"=="ponly" {
		local stddev  nostddev
		local stderr  nostderr
		local nobs    nonobs
		local stddiff nostddiff
		local tstat   notstat
	}
	//	Pull r() from -summarize-
	if "`totals'"~="nototals" &   ///  nostddev nostderr nonobs
	  ("`ponly'"~="ponly" &  ///
	  ~("`stddev'"=="nostddev" & "`stderr'"=="nostderr" & "`nobs'"=="nonobs"))  {
		local StdDevTot StdDev = r(sd)
		local StdErrTot StdErr = sqrt(r(Var)/r(N))
	}
	//	From here on, pull r() from -ttest-
	if "`stddev'"~="nostddev" {
		local StdDev1 StdDev = r(sd_1)
	    local StdDev2 StdDev = r(sd_2)
	    local StdDev  StdDev = r(sd)
	}
	if "`stderr'"~="nostderr" {
		local StdErr1 StdErr = r(se_1)
		local StdErr2 StdErr = r(se_2)
		local StdErr  StdErr = r(se)
	}
	if "`nobs'"~="nonobs" {
		local Nobs1  Nobs = r(N_1)
		local Nobs2  Nobs = r(N_2)
		local Nobs   Nobs = r(N)
		local Nobst  Nobs = r(df_t)
	}

	tempname balance  // matrices

	collect clear
	collect rename default baltable , replace
	collect style autolevels result , clear

	*	Check the -by- variable to be sure it has only two levels
	qui levelsof `by'
		if `r(r)'>2 {
			di as err "The -by- variable, `by', must have exactly 2 integer values."  ///
				_n "Here it has the following `r(r)' levels: `r(levels)'"
			exit 420
		}
		local i = 0
		foreach lvl of numlist `r(levels)' {
			if "`grplabels'"~="nogrplabels" {
				local lbl`i' : label (`by') `lvl'
				local lbl`i' = abbrev(`"`lbl`i''"',`truncgrplbl')
			}
			else {
				local lbl`i' `lvl'
			}
			local i = `i' + 1
		}

	// Since -ttest- does not accept factor variables, by default
	// -baltable- specifies the -fvrvar- option -list- to strip factor-variable
	// prefixes from `varlist' and use the base variables.
	// The -baltable- option -nolist- (synonym -fvexpand-) expands
	// the factor variables into temporary variables, but -fvrevar-
	// unfortunately does not provide useful variable names or labels.
	// Use Newson's -fprevar- if it is installed.
	// N.B.: -baltable- option -fvexpand- is a synonym for option -nolist-
	if "`list'" ~= "nolist" & "`fvexpand'"=="" {
		local list list
	}
	else {
		local list
	}
	cap which fvprevar
	tempname tmp
	if _rc==0 {
		fvprevar `varlist' , `list' gen(`tmp')
	}
	else {
		di as txt "Suggestion: To improve the row labels for expanded factor variables,"  ///
			_n "install Roger Newson's utility {help fvprevar} from {stata ssc describe fvprevar:SSC}."
		fvrevar	`varlist' , `list'
	}
	local varlist `r(varlist)'


*	At this point we need to give a full list of variables to the loop that -ttest- can handle.
*	Applying -fvrevar- to the entire list creates a list of temporary variables
*	numbered from 1 to the total number of levels of all factor variables.
*	We need to expand each factor variable to LSDV dummies individually,
*	and then give the variables names correspondin to their levels.
*	Newson's -fvprevar- does part of the job by converting each of the
*	temporary variables created by -fvrevar- to a permanent variable
*	with a name onsisting of the user's designated "stub" suffixed
*	by the factor variable level.

	foreach v of varlist `varlist'  {

		//	In the "collect" command system, this could probably be done with
		//	-collect label levels dim level "label", modify-
		if "`varlabels'"=="novarlabels" {
			label var `v'
		}

		if "`totals'"~="nototals" {
			qui summarize `v'
			collect get Mean=r(mean)  `StdDevTot' `StdErrTot' `Nobs'  , tags(var[`v'] colmn[Total])
			local totstats `r(mean)' , `r(sd)' , sqrt(`r(Var)'/`r(N)')  ,  `r(N)' ,
			local totnms  tot_mean tot_sd tot_se tot_N
		}

		`shhh' _ttest_plus `v', by(`by') `unequal' `level'

		scalar level = `r(level)'

		collect get Mean=r(mu_2)  `StdDev2' `StdErr2' `Nobs2'  , tags(var[`v'] colmn[`lbl0'])  //  Control
		collect get Mean=r(mu_1)  `StdDev1' `StdErr1' `Nobs1'  , tags(var[`v'] colmn[`lbl1'])  //  Treatment

		collect get Mean=r(mndif) `StdDev'  `StdErr'  `Nobs'   , tags(var[`v'] colmn[Difference])
		if "`stddiff'"~="nostddiff" {
			collect get Mean=r(sdiffpct)                          , tags(var[`v'] colmn[StdDiff])
		}
		if "`tstat'"~="notstat" {
			collect get Mean=r(t)                     `Nobst'  , tags(var[`v'] colmn[t-stat])
		}
		//  Report only one of the p-values
		opts_exclusive "`p_lower' `p_upper'"
		if "`p_upper'"~="" {
			collect get Mean=r(p_u)                            , tags(var[`v'] colmn[p-value*])
			local p_note "*p-values for one-sided tests of Ho: mu_2<=mu_1"
		}
		else if "`p_lower'"~="" {
			collect get Mean=r(p_l)                            , tags(var[`v'] colmn[p-value*])
			local p_note "*p-values for one-sided tests of Ho: mu_2>=mu_1"
		}
		else {  //  Default is a two-sided p-value
			collect get Mean=r(p)                              , tags(var[`v'] colmn[p-value*])
			local p_note "*p-values for two-sided tests of Ho: mu_2==mu_1"
		}

		//	Assemble the results in a matrix to be returnd to the user
		mat define `balance' = nullmat(`balance') \ ( ///
			`totstats'  ///
			`r(mu_1)' , `r(N_1)' , `r(sd_1)' , `r(se_1)' ,   ///
			`r(mu_2)' , `r(N_2)' , `r(sd_2)' , `r(se_2)' ,  ///
			`r(mndif)' , `r(sd)' , `r(se)' ,  `r(mndif_lb)' , `r(mndif_ub)' ,  ///
			`r(sdiffpct)' , `r(sdiffpct_lb)' , `r(sdiffpct_ub)' ,  ///
			`r(sdiffabs)' , `r(sdiffabs_lb)' , `r(sdiffabs_ub)' ,  ///
			`r(t)' , `r(p)' , `r(p_u)', `r(p_l)' )

		local rownames `rownames' `v'
	}
	collect notes "`p_note'"

	*	Name the rows and columns of the matrix -balance-
	mat colnames `balance' = `totnms' mu_1 N_1 StdDev1 StdErr1 mu_2 N_2   ///
		StdDev2 StdErr2 mndif StdDev StdErr mndif_lb mndif_ub   ///
		StdDiff    StdDiff_lb    StdDiff_ub  ///
		StdDiffAbs StdDiffAbs_lb StdDiffAbs_ub  ///
		t-stat p-value p_upper p_lower
	mat rownames `balance' = `rownames'

	*	Remap the column variables into Group Variables and Balance Statistics
	`shhh' collect remap colmn[Total `lbl0' `lbl1'] = grpvar
		local varlbl : variable label `by'
		if `"`varlbl'"'~="" {
			local grpvarname `"`varlbl'"'
		}
		else {
			local grpvarname "`by'"
		}
		collect label dim grpvar `"`grpvarname'"' , modify
	`shhh' collect remap colmn[Difference StdDiff t-stat p-value*] = balstats
		collect label dim balstats "Balance Statistics" , modify

	*	Label the -by- (i.e. -grpvar-) super-column with value labels
	collect style header grpvar, title(label)

	*	Label the balance statistics super-column
	collect style header balstats, title(label)

	*	Base style formatting for all cells	(option format)
	if "`format'"==""   local format %7.2f
	collect style cell result, nformat(`format')

	*	Style for row of Nobs & d.f.'s (option nformat)
	if "`nobs'"~="nonobs"{
		if "`nformat'"==""   local nformat %7.0f
		collect style cell result[Nobs], nformat(`nformat')
		collect label levels result Nobs "N, df", modify
	}

	*	Style for column of standardized differences (option sformat)
	if "`stddiff'"~="nostddiff" {
		if "`sformat'"==""   local sformat %7.1f
		collect style cell balstats[StdDiff] , nformat(`sformat') sformat("%s%%")
	}

	*	Style for column of p-values (option pformat)
	if "`pformat'"==""   local pformat %7.4f
	collect style cell balstats[p-value*] , nformat(`pformat') minimum(.0001)

	*	Add space between columns
	collect style column , dups(center) `extraspace'

	*	Format rows
	if "`stack'"~="" {
		if "`wrapon'"==""   local wrapon word
		if "`abbreviate'"==""   local abbreviate abbreviate
		local stackopts length(`wraplength') wrapon(`wrapon') wrap(`wrap')
	}
	collect style row `stack' `split', `spacer' `stackopts'

	if "`ponly'"=="ponly" | ("`stddev'"=="nostddev" & "`stderr'"=="nostderr" & "`nobs'"=="nonobs")  {
		collect style header result, level(hide)
	}

	if "`stddiff'"~="nostddiff" {
		collect notes "The standardized difference is defined as 100*(mu_1 - mu_2)/[(sd_1^2 + sd_2^2)/2]^0.5"
	}

	//  Perform the joint test of means equality only on the non-collinear subset of variables
	if "`joint'"~="nojoint"  {
		qui _rmcoll `varlist'
		foreach v in `r(varlist)' {
			local prefix = substr("`v'",1,strpos("`v'","."))
			if strpos("`prefix'","o")==0 & strpos("`prefix'","b")==0  {
				local cleanvarlist `cleanvarlist' `v'
			}
			else {
				local vv = substr("`v'",strpos("`v'",".")+1,.)
				local vvv : char `vv'[fvrevar]
				//  If there is no -fvrevar- characteristic, use the variable name, stripped of the -fv- prefix
				if "`vvv'"=="" {
					local collvars `collvars' `vv'
				}
				//  If there is a -fvrevar- characteristic, use the name stored in the characteristic
				else {
					local collvars `collvars' `vvv'
				}
			}
		}
		if "`collvars'"~="" {
			di as txt "The following " as res  `r(k_omitted)' as txt " collinear variable(s) is(are) omitted"  ///
			_n "from the joint test that all means are equal:"
		di as res "    `collvars'"
		}
		cap `noisily' mvtest means `cleanvarlist', by(`by') `heterogeneous'
		if _rc==0 {
			local mvdf1   : di %6.0f `r(df_m)'
			local mvdf2   : di %6.1f `r(df_r)'
			local mvFstat : di %6.3f `r(F)'
			local mvp_F   : di `pformat' `r(p_F)'
		}
		else {
			di as err "Stata's {help mvtest} has failed. Suppress it by specifying the -baltable- option -nojoint-"
			exit _rc
		}

		collect notes "Joint test that all means are equal: F(`mvdf1',`mvdf2') = `mvFstat'; p-value = `mvp_F'"
	}
	if "`heterogeneous'"=="heterogeneous" | "`unequal'"=="unequal"  {
		collect notes "Assuming unequal (co)variances between groups consumes degrees of freedom."
	}
	else {
		collect notes "Assuming equal (co)variances between groups conserves the degrees of freedom."
	}

	if "`title'" ~= "" {
		collect title `"`title'"'
	}

	`shhh' collect layout (var#result) ( grpvar balstats)

	collect preview

	if "`saving'"~="" {
		local saveopts = substr("`saving'", strpos("`saving'","(")+1 , strpos("`saving'",")") - 1)
		collect save `saveopts'
	}

	return scalar level = scalar(level)
	return local unequal `unequal'
	return local heterogeneous `heterogeneous'

	return local groupvar `by'
	return local varlist `varlist'
	local coefplotcmd
	return matrix balance `balance'

	if "`using'"~="" {
		collect export `using'	, `exportopts'
	}

	restore

end   /*  End of main program  */

cap prog drop _ttest_plus
program define _ttest_plus, rclass

	syntax varname , BY(varname) [unequal level(passthru) ]

	//  Option -reverse- computes the difference as mu_1 - mu_2
	//	Assumes group mu_1 is for treatment group, mu_2 is for control
	//	mu_1 is the mean when the -by- variable ==1, mu_2 when -by- == 0.
	ttest `varlist' , by(`by') reverse `unequal' `level'
		local level `r(level)'

	tempname se_1 se_2 mndif sdiffpct stdmu_1 stdmu_2  //  scalars
	scalar `se_1' = sqrt(`r(sd_1)'^2/`r(N_1)')
		return scalar se_1 = scalar(`se_1')
	scalar `se_2' = sqrt(`r(sd_2)'^2/`r(N_2)')
		return scalar se_2 = scalar(`se_2')
	scalar `mndif' = `r(mu_1)' - `r(mu_2)'

	//	Rosenbaum & Rubin's "standardized distance"
	scalar `sdiffpct' = 100*scalar(`mndif')/((scalar(`r(sd_1)')^2 + scalar(`r(sd_2)')^2)/2)^(1/2)

	//	Define "standardized means" as the mean divided by the geometric average standard error
	scalar `stdmu_1' = 100*scalar(`r(mu_1)')/((scalar(`se_1')^2 + scalar(`se_2')^2)/2)^(1/2)
	scalar `stdmu_2' = 100*scalar(`r(mu_1)')/((scalar(`se_1')^2 + scalar(`se_2')^2)/2)^(1/2)

	//	Critical values of t
	local tcrit_lb = invt(`r(df_t)',     (1-`level'/100)/2)
	local tcrit_ub = invt(`r(df_t)', 1 - (1-`level'/100)/2)

	//  Confidence bounds on -mndif-
	return scalar mndif_lb = `mndif' + `tcrit_lb' * `r(se)'
	return scalar mndif_ub = `mndif' + `tcrit_ub' * `r(se)'

	return scalar level = `r(level)'
	return scalar mndif = `mndif'
	return scalar sdiffpct = scalar(`sdiffpct')
	return scalar sdiff    = scalar(`sdiffpct')/100
	return scalar sdiffabs = abs(scalar(`sdiffpct')/100)
	return scalar N=`r(N_1)'+`r(N_2)'

	//  Confidence bounds on -sdiffpct- & -sdiff-
	return scalar sdiffpct_lb = `sdiffpct'     + `tcrit_lb' * 100* `r(se)'/((scalar(`r(sd_1)')^2 + scalar(`r(sd_2)')^2)/2)^(1/2)
	return scalar sdiffpct_ub = `sdiffpct'     + `tcrit_ub' * 100* `r(se)'/((scalar(`r(sd_1)')^2 + scalar(`r(sd_2)')^2)/2)^(1/2)
	return scalar sdiff_lb    = `sdiffpct'/100 + `tcrit_lb' *      `r(se)'/((scalar(`r(sd_1)')^2 + scalar(`r(sd_2)')^2)/2)^(1/2)
	return scalar sdiff_ub    = `sdiffpct'/100 + `tcrit_ub' *      `r(se)'/((scalar(`r(sd_1)')^2 + scalar(`r(sd_2)')^2)/2)^(1/2)
	return scalar sdiffabs_lb = abs(`sdiffpct'/100) + `tcrit_lb' * `r(se)'/((scalar(`r(sd_1)')^2 + scalar(`r(sd_2)')^2)/2)^(1/2)
	return scalar sdiffabs_ub = abs(`sdiffpct'/100) + `tcrit_ub' * `r(se)'/((scalar(`r(sd_1)')^2 + scalar(`r(sd_2)')^2)/2)^(1/2)

	//	From -ttest-, returning the scalar r(se) for the Std Err of the difference.
	return add

end  /* End of program -_ttest_plus- */


/*
See the references in -balance_plot- here:
	https://github.com/Chris-Larkin/balance-plot

"The standardized difference in percent is the mean difference as a percentage of the average standard deviation:
100*(mu_1 - mu_2)/[(sd_1^2 + sd_2^2)/2]^0.5, where for each covariate, mu_1 and mu_2 are the sample means
in the treated group and the control reservoir and sd_2^2 and sd_1^2 are the corresponding sample variances."
Source: Rosenbaum & Rubin (1985), Note to Table 1.  Notation altered because this program uses option -reverse- on -ttest-.

*/
* Version
*	1.0 21Feb2024
*	1.1 2Mar2024	Add the check for collinearity before the command -mvtest-
*			Use -fvrevar- or Roger Newson's -fprevar- to expand factor variable lists
*	1.2 6Mar2024	Assure that dummy variable 0 is control and 1 is treatment
*	1.3 13Mar2024	Added variable label stacking options for wrapping variable labels in row headers
*			Add the upper and lower bounds of the standardized difference
*	1.4 25Apr2024   -markout- the treatment variable.
*			Add the option -truncgrplbl(15)- to control size of Control & Treatment group labels, if any.
*			Add the absolute value of the standardized difference and its bounds to -r(balance)-
*			Use -abbev()- with option -truncgrplbl- to truncate group value labels
*TO DO:
*	Replace -fvprevar- with my -fvprevar2- once its debugged.