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