/* do file for gllamm companion for: Rabe-Hesketh, S. and Skrondal, A. (2012). Multilevel and Longitudinal Modeling Using Stata (3rd Edition). Volume I: Continuous Responses. College Station, TX: Stata Press. */ * install gllamm ssc install gllamm, replace ******************* Chapter 2 *********************************** ** Data preparation: use http://www.stata-press.com/data/mlmus3/pefr, clear reshape long wp wm, i(id) j(occasion) *** 2.5 Estimation using Stata ** 2.5.4 Using gllamm gllamm wm, i(id) nip(12) adapt estimates store RI *** 2.6 Hypothesis tests and confidence intervals ** 2.6.2 for between-cluster variance * likelihood ratio test quietly gllamm wm, i(id) init lrtest RI . *** 2.9 Crossed versus nested effects * random intercept model with occasion 2 dummy generate occ2 = occasion==2 gllamm wm occ2, i(id) nip(12) adapt *** 2.11 Assigning values to the random intercepts ** 2.11.2 Empirical Bayes prediction * using gllamm estimates restore RI gllapred eb, u sort id format ebm1 %8.2f list id ebm1 if occasion==1, clean noobs *** 2.11.3 Empirical Bayes standard errors * comparative standard errors display ebs1[1] * diagnostic standard errors generate diag_SE = sqrt(([id1]_cons)^2 - ebs1^2) display diag_SE[1] ******************* Chapter 3 *********************************** ** Data preparation: use http://www.stata-press.com/data/mlmus3/smoking, clear *** 3.4 Estimation using Stata ** 3.4.3 Using gllamm gllamm birwt smoke male mage hsgrad somecoll collgrad /// married black kessner2 kessner3 novisit pretri2 pretri3, /// i(momid) adapt estimates store gllamm *compare with xtmixed display 370.6656^2 display 338.7669^2 ** 3.6.1 Hypothesis tests for regression coefficients * Wald test testparm kessner2 kessner3 *** 3.9 Assigning values to random effects: Residual diagnostics estimates restore gllamm gllapred lev2, ustd gllapred lev1, pearson histogram lev1, normal xtitle(Standardized level-1 residuals) histogram lev2m1 if idx==1, normal xtitle(Standardized level-2 residuals) ******************* Chapter 4 *********************************** ** Data preparation: use http://www.stata-press.com/data/mlmus3/gcse, clear *** 4.5 (g) Estimation using gllamm ** 4.5.1 (g) Random-intercept model gllamm gcse lrt, i(school) adapt estimates store RI *compare with xtmixed display 7.521481^2 display 3.035269^2 ** 4.5.2 (g) Random-coefficient model eq slope: lrt generate cons = 1 eq inter: cons gllamm gcse lrt, i(school) nrf(2) eqs(inter slope) ip(m) nip(15) /// adapt estimates store RC *compare with xtmixed display 7.440787^2 display 3.007444^2 display .1205646^2 * use starting values from random-intercept model estimates restore RI matrix a = e(b) matrix a = (a,0,0) gllamm gcse lrt, i(school) nrf(2) eqs(inter slope) ip(m) nip(15) /// adapt from(a) copy *** 4.6 Testing the slope variance lrtest RI RC display 0.5*chi2tail(1,40.37) + 0.5*chi2tail(2,40.37) *** 4.8 Assigning values to the random intercepts and slopes ** 4.8.2 Empirical Bayes prediction estimates restore RC gllapred eb, u egen pickone = tag(school) list school ebm1 ebm2 if pickone==1 & (school<10 | school==48), noobs ** 4.8.3 Model visualization gllapred muRC, linpred sort school lrt twoway (line muRC lrt, connect(ascending)), xtitle(LRT) /// ytitle(Empirical Bayes regression lines for model 2) estimates restore RI gllapred muRI, linpred sort school lrt twoway (line muRI lrt, connect(ascending)), xtitle(LRT) /// ytitle(Empirical Bayes regression lines for model 1) ******************* Chapter 6 *********************************** ** Data preparation: use http://www.stata-press.com/data/mlmus3/wagepan, clear generate educt = educ - 12 generate yeart = year - 1980 ** 6.4.2 Heteroskedastic level-1 residuals over occasions tabulate yeart, generate(yr) eq het: yr1-yr8 generate one = 1 eq inter: one eq slope: yeart gllamm lwage black hisp union married exper yeart educt, i(nr) /// nrf(2) eqs(inter slope) s(het) ip(m) nip(11) adapt *compare with xtmixed display .3804544^2 display .053767^2 display exp([lns1]yr1) display exp([lns1]yr2) display exp([lns1]yr3) display exp([lns1]yr4) display exp([lns1]yr5) display exp([lns1]yr6) display exp([lns1]yr7) display exp([lns1]yr8) ** 6.4.3 Heteroskedastic level-1 residuals over groups generate white = 1 - black - hisp eq het: white black hisp gllamm lwage black hisp union married exper yeart educt, i(nr) /// s(het) adapt *compare with xtmixed display .3270843^2 display exp([lns1]white) display exp([lns1]black) display exp([lns1]hisp) ******************* Chapter 7 *********************************** ** Data preparation: use http://www.stata-press.com/data/mlmus3/asian, clear label define g 1 "Boy" 2 "Girl" label values gender g *** 7.3 Models for nonlinear growth ** 7.3.1 Polynomial models * fitting the models recode gender 2=1 1=0, generate(girl) generate age2 = age^2 generate cons = 1 eq inter: cons eq slope: age gllamm weight girl age age2, i(id) nrf(2) eqs(inter slope) ip(m) nip(15) adapt estimates store RC eq slope2: age2 gllamm weight girl age age2, i(id) nrf(3) eqs(inter slope slope2) /// ip(m) nip(7) adapt estimates restore RC * compare with xtmixed: display .57233^2 display .5947313^2 display .5097091^2 * Predicting trajectories for individual children gllapred traj, linpred twoway (scatter weight age) (line traj age, sort) if girl==1, /// by(id, compact legend(off)) *** 7.5 Heteroskedasticity ** 7.5.1 Heteroskedasticity at level 1 generate boy = 1 - girl eq het: boy girl eq inter: cons eq slope: age gllamm weight girl age age2, i(id) nrf(2) eqs(inter slope) /// s(het) ip(m) nip(15) adapt *compare with xtmixed: display .6271286^2 display .4857393^2 display exp([lns1]boy) display exp([lns1]girl) lrtest RC . ** 7.5.2 Heteroskedasticity at level 2 generate age_boy = age*boy generate age_girl = age*girl eq intboy: boy eq slboy: age_boy eq intgirl: girl eq slgirl: age_girl gllamm weight girl age age2, i(id id) nrf(2 2) /// eqs(intboy slboy intgirl slgirl) ip(m) nip(5) adapt * compare with xtmixed: display .567544^2 display .5364599^2 display .6891223^2 display .6987662^2 display .218154^2 lrtest RC . * alternative model eq inter: boy girl eq slope: age_boy age_girl gllamm weight girl age age2, i(id) nrf(2) eqs(inter slope) /// ip(m) nip(15) adapt * intercept variance for girls display .26625045*1.4211837^2 * slope variance for girls display .45535006*.3644414^2 lrtest RC . ******************* Chapter 8 *********************************** ** Data preparation: use http://www.stata-press.com/data/mlmus3/pefr, clear reshape long wm wp, i(id) j(occasion) generate i = _n reshape long w, i(i) j(meth) string sort id meth occasion list id meth occasion w in 1/8, clean noobs encode meth, generate(method) recode method 2=0 *** 8.6 (g) Estimation using gllamm gllamm w, i(method id) adapt * compare with xtmixed display 17.75859^2 display 19.47623^2 display 108.6037^2 *** 8.7 Empirical Bayes prediction gllapred reff, u sort id method label define m 0 "Wright" 1 "Mini Wright" label values method m list id method reffm2 reffm1 if id<8 & occasion==1, noobs sepby(id) *** 8.8 Testing variance components estimates store THR quietly gllamm w, i(id) adapt lrtest THR . *** 8.9 Crossed versus nested random effects revisited gllamm w method, i(method id) adapt * compare with xtmixed display 17.75859^2 display 19.00386^2 display 108.6455^2 ** Data preparation use http://www.stata-press.com/data/mlmus3/kenya, clear egen pick_school = tag(schoolid) tabulate treatment, generate(treat) rename treat1 meat rename treat2 milk rename treat3 calorie quietly tabulate gender, generate(g) rename g1 boy *** 8.12 Three-level random-intercept model ** 8.12.3 (g) Estimation using gllamm generate meat_year = meat*relyear generate milk_year = milk*relyear generate calorie_year = calorie*relyear gllamm ravens meat milk calorie relyear meat_year milk_year calorie_year /// age_at_time0 boy, i(id schoolid) nip(5) adapt estimates store M1 * compare with xtmixed display 2.412956^2 display 1.532073^2 display .2209702^2 *xtmixed ravens ib4.treatment##c.relyear age_at_time0 boy || schoolid: || id:, mle *** 8.13 Three-level random-coefficient model ** 8.13.1 Random coefficient at the child level generate cons = 1 eq inter: cons eq slope: relyear gllamm ravens meat milk calorie relyear meat_year milk_year calorie_year /// age_at_time0 boy, i(id schoolid) nrf(2 1) eqs(inter slope inter) /// ip(m) nip(5) adapt estimates store M2 * compare with xtmixed display 2.315432^2 display 1.45822^2 display .864055^2 display .2546548^2 lrtest M1 M2 ** 8.13.2 Random coefficient at the child and school levels matrix a=e(b) gllamm ravens meat milk calorie relyear meat_year milk_year calorie_year /// age_at_time0 boy, i(id schoolid) nrf(2 2) eqs(inter slope inter slope) /// ip(m) nip(5) from(a) adapt * compare with xtmixed display 2.31526^2 display 1.452461^2 display .8583323^2 display .3172716^2 display .1110955^2 * (results closer with nip(7), but it takes longer) estimates restore M2 *** 8.14 Residual diagnostics and predictions predict eb, u predict RES, pearson replace RES = RES*exp([lns1]_cons) rename ebm1 RI2 rename ebm2 RC2 rename ebm3 RI3 replace RI3=. if pick_school!=1 replace RI2 =. if rn!=1 graph box RI3 RI2 RES, ascategory box(1, bstyle(outline)) /// yvaroptions(relabel(1 "School" 2 "Child" 3 "Occasion")) /// medline(lcolor(black)) scatter RC2 RI2 if rn==1, saving(yx, replace) /// xtitle("Random intercept") ytitle("Random slope") histogram RC2, freq horiz saving(hy, replace) /// yscale(alt) ytitle(" ") fxsize(35) normal histogram RI2, freq saving(hx, replace) /// xscale(alt) xtitle(" ") fysize(35) normal graph combine hx.gph yx.gph hy.gph, hole(2) imargin(0 0 0 0) gllapred ptraj, linpred sort schoolid id relyear twoway (line ptraj relyear, connect(ascending)), by(schoolid, compact) /// xtitle(Time in years) ytitle(Raven's score) summarize age_at_time0 if rn==1 replace age_at_time0 = r(mean) replace boy = 1 gllapred Means, xb twoway (line Means relyear if meat==1, sort lpatt(solid)) /// (line Means relyear if milk==1, sort lpatt(dash)) /// (line Means relyear if calorie==1, sort lpatt(shortdash)) /// (line Means relyear if treatment==4, sort lpatt(longdash_dot)), /// legend(order(1 "Meat" 2 "Milk" 3 "Calorie" 4 "Control")) /// xtitle(Time in years) ytitle(Predicted mean Raven's score)