*! 3.0.1 Ariel Linden 20Mar2024 // added z-distribution option (default is t-distribution with noncentrality parameter)
*! 3.0.0 Ariel Linden 07Mar2024 // added Hedges g option
								// -esizereg- now uses variance produced by margins to compute pooled std. dev.
*! 2.0.2 Ariel Linden 26Feb2024 // added -intreg-, -meintreg- and -metobit- models
*! 2.0.1 Ariel Linden 27Oct2021 // changed 'est' to a scalar to avoid issues with squaring negative values (which happens with local)
*! 2.0.0 Ariel Linden 29May2019 // made esizereg a postestimation command and converted version 1.0.0 to an immediate command (esizeregi)
*! 1.0.0 Ariel Linden 02Feb2019

capture program drop esizereg
program define esizereg, rclass
version 11.0

			syntax anything [, COHensd HEDgesg Zdistribution LEVel(cilevel) ]

			gettoken treat : 0
			// * store model estimates * //
			if "`e(cmd)'" == "margins" {
				di as err "You must re-estimate the regression model"
				exit 198
			}
			else estimates store results

			if !inlist("`e(cmd)'", "regress", "tobit", "truncreg", "hetregress", "xtreg", "intreg") & !inlist("`e(cmd2)'", "meintreg", "metobit") {
				di as err `"`e(cmd)' is not supported by {bf:esizereg}"'
				exit 198
			}

			// * verify that the previous model estimates are available * //
			capture assert matrix list r(table)
				if _rc {
					qui estimates restore results
					qui estimates replay results
				}

			// * save table of estimates as matrix * //
			qui matrix b = r(table)


			// * generate error if treat(name) is not in regression table * //
			local colnames : colnames b
			if !`: list treat in colnames' {
				di as err "The regressor {bf:`treat'} is not found." ///
				" Use the model's option -coeflegend- to display coefficient names"
				exit 498
			}

			tempname est
			scalar `est' = b[1,colnumb(matrix(b),"`treat'")]

			local weightexp `e(wtype)' `e(wexp)'

			// extract treat varname from that produced in -regress- using factor variables (e.g. 1.treat)
			local treatvar = substr("`treat'", strpos("`treat'", ".") + 1, .)

			// * ensure binary variable *
			qui tabulate `treatvar' if e(sample)
				if r(r) != 2 {
					di as err "With a binary treatment, {bf:`treatvar'} must have exactly two values (coded 0 or 1)."
					exit 420
				}
			else if r(r) == 2 {
				capture assert inlist(`treat', 0, 1) if e(sample)
				if _rc {
					di as err "With a binary treatment, {bf:`treatvar'} must be coded as either 0 or 1."
					exit 450
				}
			}

			// * get N per group * //
			qui count if `treatvar' == 1 & e(sample)
			local n1 = r(N)
			qui count if `treatvar' == 0 & e(sample)
			local n0 = r(N)

			// * get variance of model using margins * //
			tempname V  N sdpooled CohensD v se iz CohensD_Lower CohensD_Upper

			if "`e(cmd)'" == "meglm" {
				qui margins, post predict(mu fixed)
			}
			else {
				qui margins, post
			}
			matrix `V' = r(V)
			scalar `V' = `V'[1,1]
			scalar `N' = `n1' + `n0'

			// CALCULATE COHEN'S D
			// ==============================================================
			scalar `sdpooled' = sqrt(`V' * `N')
			scalar `CohensD' = `est' / `sdpooled'
			scalar `v' = (`n1' + `n0') / (`n1' * `n0') + (`CohensD'^2) / (2 *(`n1' + `n0'))
			scalar `se' = sqrt(`v')

			if "`zdistribution'" != "" {
				scalar `iz' = invnorm(1-(1-`level'/100)/2)
				scalar `CohensD_Lower' = `CohensD' - `iz' * sqrt(`v')
				scalar `CohensD_Upper' = `CohensD' + `iz' * sqrt(`v')
			}
			else {
				tempname alpha AlphaLower AlphaUpper ns df LowerLambda UpperLambda
				scalar `alpha' = 1-(`level'/100)
				scalar `AlphaLower' = `alpha'/2
				scalar `AlphaUpper' = 1 - (`alpha'/2)
				scalar `ns' = sqrt((`n1'*`n0')/(`n1'+`n0'))
				scalar `df' = `N' - 2
				scalar `LowerLambda' = npnt(`df',`CohensD'*`ns',`AlphaUpper')
				scalar `UpperLambda' = npnt(`df',`CohensD'*`ns',`AlphaLower')
				scalar `CohensD_Lower' = `LowerLambda' * sqrt((`n1'+`n0')/(`n1'*`n0'))
				scalar `CohensD_Upper' = `UpperLambda' * sqrt((`n1'+`n0')/(`n1'*`n0'))
			}

			// CALCULATE HEDGE'S G
			// =================================================================
			// EXACT BIAS CORRECTION: Hedges (1981) pg 111, Equation 6e
			tempname m BiasCorrectionFactor HedgesG HedgesG_Lower HedgesG_Upper
			scalar `m' = (`n1'+`n0'- 2)
			scalar `BiasCorrectionFactor' = exp(lngamma(`m'/2) - 1/2 * ln(`m'/2) - lngamma((`m'-1)/2))
			// Turner & Bernard (2006) , Eq 4
			scalar `HedgesG' = `CohensD' * `BiasCorrectionFactor'
			scalar `HedgesG_Lower' = `CohensD_Lower' * `BiasCorrectionFactor'
			scalar `HedgesG_Upper' = `CohensD_Upper' * `BiasCorrectionFactor'


			// DISPLAY OUTPUT
			// ====================================================================
			// SET DEFAULT OUTPUT
			if "`cohensd'"== "" & "`hedgesg'"== "" {
				local cohensd "cohensd"
				local hedgesg "hedgesg"
			}

			// Display Title (weighted or unweighted)
			if "`weightexp'" == "" {
				disp _newline as text "Effect size based on the regression coefficient of the treatment (exposure) variable"
            }
			else {
				disp _newline as text "{bf:Weighted} effect size based on the regression coefficient of the treatment (exposure) variable"
			}

			// Display table header information
			disp _newline %45s "Obs per group:"
			disp %47s "Group 1 = " %10.0fc `n1'
			disp %47s "Group 2 = " %10.0fc `n0'

			// Display output table
			tempname mytab
			.`mytab' = ._tab.new, col(5) lmargin(0)
			.`mytab'.width    20   |11  12  12    12
			.`mytab'.titlefmt  .     .   . %24s   .
			.`mytab'.pad       .     1   1  3     3
			.`mytab'.numfmt    . %9.0g %9.0g %9.0g %9.0g
			.`mytab'.strcolor result  .  .  .  .
			.`mytab'.strfmt    %19s  .  .  .  .
			.`mytab'.strcolor   text  .  .  .  .
			.`mytab'.sep, top
			.`mytab'.titles "Effect Size"							/// 1
							"Estimate"								/// 2
							"Std. Err."								/// 3
							"[`level'% Conf. Interval]" ""          //  4 5
			.`mytab'.sep, middle
                .`mytab'.strfmt    %24s  .  .  .  .
			if "`cohensd'" != "" {
                .`mytab'.row    "Cohen's {it:d}"        ///
                        `CohensD' 	                   	///
                        `se'							///
						`CohensD_Lower'                 ///
                        `CohensD_Upper'
			}
			if "`hedgesg'" != "" {
                .`mytab'.row    "Hedges's {it:g}"       ///
                        `HedgesG'                      	///
                        `se'							///
						`HedgesG_Lower'                 ///
                        `HedgesG_Upper'
			}
				.`mytab'.sep, bottom

			// Return results
			if "`hedgesg'" != "" {
                return scalar ub_g = `HedgesG_Upper'
                return scalar lb_g = `HedgesG_Lower'
                return scalar g = `HedgesG'
			}
			if "`cohensd'" != "" {
                return scalar ub_d = `CohensD_Upper'
                return scalar lb_d = `CohensD_Lower'
                return scalar d = `CohensD'
			}

			return scalar se = `se'
			return scalar n0 = `n0'
			return scalar n1 = `n1'
			return scalar sdpooled = `sdpooled'
			return scalar V = `V'
			return scalar est = `est'

			// Make a c_local macro of d, g and se
			c_local d = `CohensD'
			c_local g = `HedgesG'
			c_local se = `se'

end