clear all
set obs 1000
forvalues i=1/4 {
	generate x`i' = runiform(-3,3) // x's independent and balanced
	}
gen ylo = x1+2*x2+3*x3+4*x4
//histogram ylo // normal distribution
generate ypr = invlogit(ylo)
//histogram ypr  // beta distribution
generate y = runiform() < ypr

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")
A4c = st_matrixcolstripe("B4")
A3 = st_matrix("B3")
A3c = st_matrixcolstripe("B3")
A2 = st_matrix("B2")
A2c = st_matrixcolstripe("B2")
A1 = st_matrix("B1")
A1c = st_matrixcolstripe("B1")
end

mata // a function to match varnames across coefficient vectors
mata drop varmatchpos()
real vector function varmatchpos(from, to)
{
	matches = J(1,rows(from),0)
	for (i=1; i<=rows(to); i++) {
		matchi = (from[.,2]:==to[i,2])':*i
		matches = matches + matchi
	}
	return(matches)
}
end

mata // build conforming coefficient-ratio vectors
S4 = J(1,cols(A4),.)
S3 = J(1,cols(A4),.)
S2 = J(1,cols(A4),.)
S1 = J(1,cols(A4),.)

S4 = A4
S3[1,varmatchpos(A3c, A4c)] = A3
S2[1,varmatchpos(A2c, A4c)] = A2
S1[1,varmatchpos(A1c, A4c)] = A1
S4 = A4

S4 = S4:/A4
S3 = S3:/A4
S2 = S2:/A4
S1 = S1:/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