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