clear all
set obs 1000
forvalues i=1/4 {
	generate x`i' = runiform(-3,3) // x's independent and balanced
	}
summarize x*
pwcorr x*, sig // no substantial correlation

egen ylo = rowtotal(x*)
//histogram ylo // normal distribution
generate ypr = invlogit(ylo)
//histogram ypr  // beta distribution
generate y = runiform() > ypr
tabulate y  // should be roughly "fair"

quietly logit y x*
estimates store four
quietly logit y x1-x3
estimates store three
quietly logit y x1-x2
estimates store two
quietly logit y x1
estimates store one

estimates table four three two one, se stats(bic r2_p)
// Note se's become smaller
coefplot (four)(three)(two)(one), drop(_cons)

khb logit y x1-x2 || x3-x4