***************************************************************************** * Stata do-file to replicate results given in the following presentation: * * TITLE: Running MLwiN from within Stata: The runmlwin command * AUTHORS: George Leckie and Chris Charlton * MEETING: MUGS, Department of Social Medicine * LOCATION: Bristol * DATE: 14-06-2011 * * George Leckie and Chris Charlton * Centre for Multilevel Modelling, 2011 ***************************************************************************** * The examples in this do-file are adapted from: * * Rabe-Hesketh, S. and Skrondal, A. (2008) Multilevel and Longitudinal * Modeling Using Stata, 2nd Edition, Stata Press. * ***************************************************************************** * To run this do-file in Stata you must be working with: * * (1) Stata 9.2 or higher * * (2) The latest version of MLwiN * * (3) The latest version of the runmlwin command * * runmlwin will inform you if you are working with out of date versions of * Stata or the MLwiN software package. * * You must also tell runmlwin where MLwiN is installed. One way you can do * this is by specifying the MLwiN path by defining a global macro called * MLwiN_path. For example: * * . global MLwiN_path C:\Program Files\MLwiN v2.25\mlwin.exe * * See the "Remarks on installation instructions" section of the runmlwin help * file for further information: * * . help runmlwin * ***************************************************************************** * Change this global macro to point to your MLwiN path ///global MLwiN_path C:\Program Files (x86)\MLwiN v2.25\mlwin.exe ***************************************************************************** * (1) TWO-LEVEL CONTINUOUS RESPONSE MODELS ***************************************************************************** * Open the data use "http://www.stata-press.com/data/mlmus2/asian.dta", clear * Plot child weight against age graph twoway (connect weight age, connect(ascending)), /// ytitle("Weight in Kg") xtitle("Age in years") * Generate a variable cons to act as the constant in the models generate cons = 1 * Generate age squared generate age2 = age^2 * Fit a two-level random intercepts model to child weight runmlwin weight cons age age2, /// level2(id: cons) /// level1(occ: cons) * Calculate the predicted weight of a 2 year-old lincom 1*cons + 2*age + 4*age2 * Fit a random intercepts and random slopes model to child weight runmlwin weight cons age age2, /// level2(id: cons age) /// level1(occ: cons) nopause * Display the variance parameters as standard deviations runmlwin, sd * Display the correlation between the random intercepts and random slopes runmlwin, correlation * Plot the child-level varaince function twoway (function [RP2]var(cons) + 2*[RP2]cov(cons\age)*x /// + [RP2]var(age)*x^2, range(age)), /// ytitle("Between-child variance") xtitle("Age in years") * Re-fit model, but request the posterior estimates for the random effects runmlwin weight cons age age2, /// level2(id: cons age, residuals(u)) /// level1(occ: cons) nopause * Tag each child once egen pickone = tag(id) * Plot the slope residuals against the intercept residuals scatter u1 u0 if pickone==1, yline(0) xline(0) /// aspectratio(1) ylabel(-2(1)2) xlabel(-2(1)2) * Generate the predicted child weights generate prediction = _b[cons]*cons + _b[age]*age + _b[age2]*age2 /// + u0 + u1*age * Plot the predicted child weights against age line prediction age, connect(a) /// ytitle("Predicted weight in Kg") xtitle("Age in years") * Drop redundant variables drop u0 u1 prediction * Define a label for gender label define genderlabel 1 "Boy" 2 "Girl" * Attach the gender label to the gender variable label values gender genderlabel * Plot child weight against age, by gender graph twoway (connect weight age, connect(ascending)), by(gender) /// ytitle("Weight in Kg") xtitle("Age in years") * Generate a dummy variable for girls generate girl = (gender==2) * Add a gender main effect runmlwin weight cons age age2 girl, /// level2(id: cons age) /// level1(occ: cons) nopause * Generate an interaction term between girl and age generate girlXage = girl*age * Add the interaction term to the model runmlwin weight cons age age2 girl girlXage, /// level2(id: cons age) /// level1(occ: cons) nopause * Store the model results estimates store m1 * Generate a dummy variable for boys generate boy = 1 - girl * Allow girls and boys to have different residual variances runmlwin weight cons age age2 girl girlXage, /// level2(id: cons age) /// level1(occ: girl boy, diagonal) nopause * Store the model results estimates store m2 * Wald test of whether girls and boys have different residual variances test [RP1]var(girl) = [RP1]var(boy) * LR test of whether girls and boys have different residual variances lrtest m1 m2 * Generate an interaction term between boy and age generate boyXage = boy*age * Specify an elements on/off vector for the following model matrix a = (1,1,1,0,0,1,0,0,1,1) * Allow boys and girls to have separate intercept and slope (cov)ariances runmlwin weight boy boyXage girl girlXage age2, /// level2(id: boy boyXage girl girlXage, elements(a)) /// level1(occ: boy girl, diagonal) nopause ***************************************************************************** * (2) MULTILEVEL BINARY RESPONSE MODELS ***************************************************************************** use "http://www.stata-press.com/data/mlmus2/guatemala.dta", clear generate cons = 1 runmlwin immun cons kid2p rural pcInd81, /// level3(cluster: cons) /// level2(mom: cons) /// level1(kid:) /// discrete(distribution(binomial) link(logit) denominator(cons)) /// nopause estimates store mql1 runmlwin immun cons kid2p rural pcInd81, /// level3(cluster: cons) /// level2(mom: cons) /// level1(kid:) /// discrete(distribution(binomial) link(logit) denominator(cons) pql2) /// initsprevious maxiterations(500) nopause estimates store pql2 estimates table mql1 pql2 ***************************************************************************** * (3) SIMULATION STUDIES ARE EASY ***************************************************************************** ***************************************************************************** * (4) MCMC ESTIMATION ***************************************************************************** runmlwin immun cons kid2p rural pcInd81, /// level3(cluster: cons) /// level2(mom: cons) /// level1(kid:) /// discrete(distribution(binomial) link(logit) denominator(cons)) /// mcmc(on) initsprevious nopause estimates store mcmc estimates table mql1 pql2 mcmc mcmcsum, trajectories mcmcsum [RP2]var(cons), fiveway mcmcsum [RP2]var(cons) runmlwin immun cons kid2p rural pcInd81, /// level3(cluster: cons) /// level2(mom: cons) /// level1(kid:) /// discrete(distribution(binomial) link(logit) denominator(cons)) /// mcmc(burnin(500) chain(5000) thinning(1) hc(3)) /// initsprevious nopause ***************************************************************************** * (5) WORK EFFICIENTLY ***************************************************************************** ***************************************************************************** * (6) RESOURCES TO HELP YOU LEARN RUNMLWIN ***************************************************************************** ***************************************************************************** exit