use schools, clear
gen id = _n


** Produce table of estimates
eq het: sd              /* SD = alpha*sd */

cons def 1 [s1]sd = 1   /* alpha = 1 */


* ML for sigma
gllamm est, i(id) s(het) nats constr(1) nip(20) adapt allc
gllamm, robust


* Bayes modal, sigma, alpha=2
gllamm est, i(id) s(het) nats constr(1) prior(gamma, scale(10000) shape(2)) nip(20) adapt allc
gllamm, robust

* Bayes modal, sigma, alpha=3
gllamm est, i(id) s(het) nats constr(1) prior(gamma, scale(10000) shape(3)) nip(20) adapt allc
gllamm, robust


* Bayes modal, sigma^2, alpha=1.5
gllamm est, i(id) s(het) nats constr(1) prior(gamma, scale(10000) shape(1.5) variance) nip(20) adapt allc
gllamm, robust

* Bayes modal, sigma^2, alpha=2
gllamm est, i(id) s(het) nats constr(1) prior(gamma, scale(10000) shape(2) variance) nip(20) adapt allc
gllamm, robust