* The Problem
/* Nested logit models are not directly comparable - a model which is
expressed as a subset of the variables in another model has coefficients
biased toward zero if the omitted variables have any discriminatory
power.*/

postfile coefs rep x1r x1f x2f cons using "khbproblem.dta", replace

quietly forvalues i=1/100 {
	clear
	set obs 200
	generate x1 = runiform(-3,3) // centered, for power
	generate x2 = runiform(-3,3)
	generate y = rlogistic() <(x1+2*x2)

	logit y x1
	scalar b_x1 = _b[x1]
	logit y x1 x2
	post coefs (`i') (b_x1) (_b[x1]) (_b[x2]) (_b[_cons])
}
postclose coefs

use khbproblem.dta, clear
generate bias = x1r-x1f
scatter bias x1r

summarize x1f x1r bias

/* So if we have a logit model with a mediator (an indirect effect), 
we cannot decompose total, direct, and indirect effects like we do
in linear models. */

clear
set obs 200
generate x1 = runiform(-3,3)
generate x2 = x1 + rnormal()
generate y = rlogistic() < (x1 + x2)

logit y x1 x2  // here, _b[x1] is the direct effect
logit y x1 // but in this model _b[x1] is NOT the total effect

khb logit y x1 || x2, summary
/* Here, "Reduced" means total effects, "Full" means direct effect,
and "Diff" is the indirect effect */

gsem (y <- x1 x2, logit) (x2 <- x1)

nlcom _b[y:x2]*_b[x2:x1]  // this is the indirect (difference) effect
nlcom (_b[y:x2]*_b[x2:x1]+_b[y:x1]) // total effect
nlcom _b[y:x2]*_b[x2:x1]/(_b[y:x2]*_b[x2:x1]+_b[y:x1]) // Conf_pct=indirect/total
nlcom (_b[y:x2]*_b[x2:x1]+_b[y:x1])/_b[y:x1] // Conf_ratio
// No rescaling ratio from gsem

/* We can implement the KHB approach in the gsem context */
quietly regress x2 x1
predict x2res, res
gsem (y <- x1 x2res, logit) (x2res <- x1)
/* Here _b[y:x1] is the total effect, _b[x2res:x1] is zero, and
_b[y:x2res] is the same as _b[y:x2] from the previous gsem model. */

* Using the KHB Example
clear
use dlsy_khb
summarize univ fses abil intact boy

khb logit univ fses || abil intact boy, summary disentangle

gsem (univ <- fses abil intact boy, logit) (abil intact boy <- fses)
estimates store base
nlcom _b[univ:abil]*_b[abil:fses]
nlcom _b[univ:intact]*_b[intact:fses]

/* To use the KHB method, we orthoganalize all the mediators
with respect to the independent variable(s). */
foreach var of varlist abil intact boy {
	quietly regress `var' fses
	predict `var'res, res
	}
	
gsem (univ <- fses abilres intactres boyres, logit) (abilres intactres boyres <- fses)
estimates store orthog
estimates table base orthog

/* compare with dropped mediators model */
gsem (univ <- fses, logit)
estimates table base orthog .

khb logit univ fses || abil i.intact i.boy, summary disentangle

gsem (univ <- fses abil i.intact i.boy, logit) (abil <- fses) ///
	 (intact boy <- fses, logit)
estimates table base .
nlcom _b[univ:1.intact]*_b[intact:fses]