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

version 8.1

/*
Stata do-file for some of the analysis of Section 9.5 of
Skrondal, A. and Rabe-Hesketh, S (2004). Generalized
Latent Variable Modeling. Multilevel, Longitudinal
and Structural Equation Models. Chapman & Hall/CRC
*/

* Read data
insheet using gum.dat, clear
rename study studynam
gen study=_n
list, clean

* Prepare data for analysis, creating variables
* gum=1 if in treatment arm, 0 otherwise
* quit=1 if quit smoking, 0 otherwise
* wt1: frequency weight

reshape long d n, i(study) j(gum)

gen y0 = n-d
rename d y1
gen i=_n
reshape long y, i(i) j(quit)

rename y wt1
drop i

sort study gum quit
list in 1/12, clean

* Continuous random intercept and slope model using gllamm
gen treat = gum -0.5
gen cons=1
eq int: cons
eq slop: treat

gllamm quit treat, i(study) nrf(2) eqs(int slop) l(logit) /*
  */ f(binom) weight(wt) adapt nocor trace allc


* Empirical Bayes predictions of study-specific effect sizes
gllapred eff, u

* 95% intervals
gen eff = effm2 + _b[treat]
gen lower = effm2 + _b[treat] - 1.96*effs2
gen upper = effm2 + _b[treat] + 1.96*effs2

* list results
sort study
qui by study: gen last=_n==_N
list studynam lower eff upper if last==1, clean


* plot results using forest plot

encode studynam, gen(s)
twoway (rcap lower upper s, horizontal blpattern(dash)) (scatter s eff),   /*
  */ ylabel(1(1)27,valuelabels angle(horizontal) labs(small)) legend(off)  /*
  */ xscale(range(-.5 1.5)) ysize(7) xlab(-.5(.5)1.5) xline(0) /*
  */ xtitle(Log odds ratio) saving(gumf, replace) 


* NPML

* two masses
gllamm quit treat, i(study) nrf(2) eqs(int slop) l(logit) /*
  */ f(binom) weight(wt) ip(f) nip(2) trace

* add one mass at a time
local n=2
while $HG_error==0&`n'<10{
    matrix a = e(b)
    local ll = e(ll)
    local k = e(k)
    local n = `n' + 1
    gllamm quit treat , i(study) nrf(2) eqs(int slop) l(logit) /*
       */ f(binom) weight(wt) ip(f) nip(`n') lf0(`k' `ll') /*
       */ gateaux(-5 5 40) from(a) trace
 }

* store locations and probabilities as variables

matrix a=e(zlc2)'
svmat a
matrix b=e(zlc3)'
svmat b
matrix prob = e(zps2)'
svmat prob
replace prob1=exp(prob1)
gen percent=round(prob1*100,1)

* scatterplot with symbol sizes determined by percent

twoway (scatter b1 a1 [w=percent], msymbol(Oh))   /*
  */ (scatter b1 a1, mlabel(percent) mlabpos(0) msymbol(none)), /*
  */ xtitle(Intercept) ytitle(Slope) legend(off) saving(blobs, replace)


log close
exit