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*
matrix B4 = e(b)
estimates store four
quietly logit y x1-x3
matrix B3 = e(b)
estimates store three
quietly logit y x1-x2
matrix B2 = e(b)
estimates store two
quietly logit y x1
matrix B1 = e(b)
estimates store one

estimates table four three two one, se stats(bic r2_p)


mata // bring coefficient matrices into Mata
A4 = st_matrix("B4")
A4r = st_matrixrowstripe("B4")
A4c = st_matrixcolstripe("B4")
_matrix_list(A4, A4r, A4c)
A3 = st_matrix("B3")
A3
A2 = st_matrix("B2")
A2
A1 = st_matrix("B1")
A1
end

mata // build conforming coefficient-ratio vectors
S4 = A4:/A4
S3 = (A3[1::3],J(1,1,.),A3[4::4]):/A4
S2 = (A2[1::2],J(1,2,.),A2[3::3]):/A4
S1 = (A1[1::1],J(1,3,.),A1[2::2]):/A4
S4\S3\S2\S1
exp(S4\S3\S2\S1)
end

mata // row means for each vector, w/o constants
S4mu = mean(S4[1::4]')
S3mu = mean(S3[1::4]')
S2mu = mean(S2[1::4]')
S1mu = mean(S1[1::4]')
S4mu, S3mu, S2mu, S1mu
end

mata // rescaled ratios
R4 = S4/S4mu
R3 = S3/S3mu
R2 = S2/S2mu
R1 = S1/S1mu
R4\R3\R2\R1
exp(R4\R3\R2\R1)
end