capture log close
log using lips, replace text
set more off

use lips, clear
gen smr = 100*o/e

* random intercept model
gen lne = ln(e)
gllamm o, i(county) offset(lne) f(poiss) adapt
estimates store model1

* add covariate
matrix a=e(b)
gllamm o x, i(county) offset(lne) f(poiss) adapt from(a)
estimates restore model1

* predicted SMRs
gllapred mu, mu nooffset
gen thet = 100*mu
sort county
list county thet in 1/10, clean noobs

twoway (scatter thet smr, msymbol(none) mlabpos(0) mlabel(county)) ///
 (function y=x, range(0 600)), xline(108) yline(108) ///
 xtitle(Crude SMR) ytitle(Empirical Bayes SMR) legend(off)

 
/* NPML */
gllamm o, i(county) offset(lne) f(poiss) nip(2) ip(f)

matrix a=e(b)
local ll=e(ll)
local k=e(k)

gllamm o, i(county) offset(lne) f(poiss) nip(3) ip(f) ///
  gateaux(-5 5 100) from(a) lf0(`k' `ll')

matrix a=e(b)
local ll=e(ll)
local k=e(k)

gllamm o, i(county) offset(lne) f(poiss) nip(4) ip(f) ///
  gateaux(-5 5 100) from(a) lf0(`k' `ll')
  
matrix a=e(b)
local ll=e(ll)
local k=e(k)

gllamm o, i(county) offset(lne) f(poiss) nip(5) ip(f) ///
  gateaux(-5 5 100) from(a) lf0(`k' `ll')
  
gllamm o, i(county) offset(lne) f(poiss) nip(5) ip(f) ///
  gateaux(-5 5 1000) from(a) lf0(`k' `ll')

* plot discrete distribution
matrix locs = e(zlc2)'
matrix lp = e(zps2)'
svmat locs
svmat lp

gen p=exp(lp1)
gen smrloc = 100*exp(_b[_cons]+locs1)
twoway (dropline p smrloc), xtitle(Location) ytitle(Probability)

* predicted SMRs
gllapred mu2, mu nooffset
gen thet2 = 100*mu2
sort county
list county thet thet2 in 1/10, clean noobs

log close
exit