******* Do-file for UK Stata Users Group Meeting, 2008 * Prediction for multilevel models * Sophia Rabe-Hesketh * 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 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) *** 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))