******* Do-file for UK Stata Users Group Meeting, 2008 
* Prediction in multilevel logistic regression
* Sophia Rabe-Hesketh and Anders Skrondal

* Run this file interactively to produce all the graphs
* of the talk

*** install latest version of gllamm (must be after 15 Sept 2008)
ssc install gllamm, replace

*** install ci_marg_mu (requires Stata version 10.0 or later)
ssc install ci_marg_mu, replace

*** read and prepare data
use antibiotics, clear
gen id = _n
recode abuse 1=0 2=1 3=1

*** fit model
gllamm abuse WHO age temp Paymed Selfmed Wrdiag DRed, ///
   i(doc hosp) link(logit) family(binom) adapt

gllamm, eform
estimates store full

*** predictions of random intercepts and confidence intervals (CIs)
(note that CIs ignore parameter uncertainty)

gllapred eb, u fsample

egen pickone = tag(hosp)
gsort +ebm2 -pickone
gen rank = sum(pickone)

*gen labpos = ebm2 + 1.96*ebs2 + .1

gen lower = ebm2-1.96*ebs2
gen upper = ebm2+1.96*ebs2

gen labpos = upper+.1
gen rank_e = rank-.1
gen rankp_e = rank+.1


serrbar ebm2 ebs2 rank if pickone==1,                             ///
  addplot(scatter labpos rank, mlabel(hosp) msymbol(none)         ///
     mlabpos(0) mlabcol(black))                                   ///
  scale(1.96) xtitle(Rank of predicted hospital random intercept) ///
  ytitle(Hospital random intercept with 95% CI) legend(off)       ///
  xscale(range(1 36)) xlabel(1/36, alternate) 

* alternative graph command  
set scheme s1mono
twoway (rcap lower upper rank, blpatt(solid) lcol(black)) ///
  (pcspike ebm2 rank_e ebm2 rankp_e,lwidth(thick) lcol(black) ) ///
  (scatter labpos rank, mlabel(hosp) msymbol(none) mlabpos(0) ///
    mlabcol(black) mlabsiz(medium)),  ///
  xtitle(Rank of predicted hospital random intercept) ///
  ytitle(Hospital random intercept with 95% CI) legend(off)       ///
  xscale(range(1 36)) xlabel(1/36) ysize(1.4)
  


*** predicted posterior mean probabilities
use postprob, clear  /* prediction data */

gllapred postp, fsample mu

sort hosp DRed
twoway (line postp DRed if hosppred==1, connect(ascending)), ///
  xtitle("Doctor's qualification")                           ///
  ytitle(Predicted posterior mean probability)

egen interesting = anymatch(hosp), values(1 8 10 4 6 3 22 28 25 20 19 13)
gen DRedjit=0
replace DRedjit = DRed + .5*(uniform()-.5)
twoway (line postp DRed if hosppred==1, sort)                           ///
   (scatter postp DRedjit if docpred==1, msym(Oh) msize(large))         ///
      if interesting==1, by(hosp, compact legend(off))                  ///
   xscale(range(1 6)) xlabel(1(1)6) xtitle("Doctor's qualification")    ///
   ytitle(Predicted posterior mean probability)


*** predicted marginal probabilities
use margprob, clear  /* prediction data */

gllapred prob, mu marg fsample

twoway ///
  (line prob DRed if WHO==0, sort lpatt(dash) lcolor(navy) ) ///
  (line prob DRed if WHO==1, sort lcolor(maroon)),           ///
  legend(order(1 "no intervention" 2 "intervention"))        ///
  xtitle("Doctor's qualification")                           ///
  ytitle(Predicted population averaged probability)          ///
  yscale(range(0 1)) ylabel(0(.2)1,format(%2.1f))
  

*** Approximate 95% CIs for predicted marginal probabilities
ci_marg_mu lower upper, level(95) dots  

twoway ///
(rarea lower upper DRed if WHO==0, sort fintensity(inten10) lwidth(none)) ///
(rarea lower upper DRed if WHO==1, sort fintensity(inten10) lwidth(none)) ///
(rline lower upper DRed if WHO==0, lcolor(navy))           ///
(rline lower upper DRed if WHO==1, lcolor(maroon))         ///
(line prob DRed if WHO==0, sort lpatt(dash) lcolor(navy) ) ///
(line prob DRed if WHO==1, sort lcolor(maroon)),           ///
 legend(order(5 "no intervention" 6 "intervention"))       ///
 xtitle("Doctor's qualification")                          ///
 ytitle(Population averaged probability with 95% CI)       ///
 yscale(range(0 1)) ylabel(0(.2)1,format(%2.1f))