***************************************************************************** * 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 * CONFERENCE: CCSR/ISC Seminar series * LOCATION: Manchester * DATE: 25-09-2012 * * George Leckie and Chris Charlton * Centre for Multilevel Modelling, 2012 ***************************************************************************** * To run this do-file in Stata you must be working with: * * (1) Stata 12.1 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.26\mlwin.exe **************************************************************************** * 1. TWO-LEVEL MULTILEVEL MODELS **************************************************************************** * Open the tutorial data set use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear * Fit a two-level (students within schools) variance components model to * a continuous educational response variable, normexam. Note, you will need * to click the "Resume Macro" button twice in MLwiN to return the model * results to the Stata output window. runmlwin normexam cons, /// level2(school: cons) /// level1(student: cons) * Store the model estimates estimates store model1 * Generate a boy dummy variable generate boy = 1 - girl * Extend the previous model to include fixed part covariates, a random school * level slope and separate level 1 residuals for boys and girls. The runmlwin * command also requests that runmlwin extracts the predicted values for the * school level residuals from MLwiN and returns them to Stata. The nopause * option prevents MLwiN from pausing before and after model estimation and so * returns the model results automatically to Stata. runmlwin normexam cons standlrt girl, /// level2(school: cons standlrt, residuals(u)) /// level1(student: girl boy, diagonal) nopause * Store the model estimates estimates store model2 * Perform a likelihood ratio test to compare the boy and girl residual * variances lrtest model1 model2 * Perform a Wald test to compare the boy and girl residual variances test [RP1]var(girl) = [RP1]var(boy) * Calculate the VPC for the average boy nlcom (Boy_VPC_xis0: [RP2]var(cons)/([RP2]var(cons) + [RP1]var(boy))) * Generate predicted school lines for boys generate yhat = [FP1]cons + [FP1]standlrt*standlrt + u0 + u1*standlrt sort school standlrt line yhat standlrt, connect(ascending) * Keep the first student in each school bysort school: keep if _n==1 * Generate the ranks of the school level intercept residuals egen u0rank = rank(u0) * Produce a caterpillar plot for school level intercept residuals serrbar u0 u0se u0rank, scale(1.96) yline(0) **************************************************************************** * 2. MULTILEVEL MODELS FOR BINARY RESPONSES **************************************************************************** * Open the tutorial data set use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear * Generate a binary pass/fail version of the continuous response gen passexam = (normexam>0) * Fit a two-level random intercepts random slopes model to the binary * response with the same fixed part and school level random part as * before runmlwin passexam cons standlrt girl, /// level2(school: cons standlrt) /// level1(student:) /// discrete(distribution(binomial) link(logit) denominator(cons)) /// nopause estimates store mql1 * Fit the previous model using PQL2 (as opposed to the default of MQL1) runmlwin passexam cons standlrt girl, /// level2(school: cons standlrt) /// level1(student:) /// discrete(dist(binomial) link(logit) denom(cons) pql2) /// initsprevious nopause estimates store pql2 estimates table mql1 pql2, /// stats(ll N) b(%4.3f) stfmt(%4.0f) varwidth(18) newpanel **************************************************************************** * 4 MCMC ESTIMATION **************************************************************************** * Fit the previous model using MCMC estimation with a burnin of 500 * iterations and a chain of 5000 iterations runmlwin passexam cons standlrt girl, /// level2(school: cons standlrt) /// level1(student:) /// discrete(d(binomial) l(logit) de(cons)) /// mcmc(burnin(500) chain(5000)) /// initsprevious nopause * Produce the MCMC trajectories plot for each parameter mcmcsum, trajectories * Produce the MCMC kernel densities plot for each parameter mcmcsum, densities * Produce the MCMC fiveway diagnostic plot for school level slope variance * parameter mcmcsum [RP2]var(standlrt), fiveway * Produce MCMC summary diagnostics for the school level slope variance * parameter mcmcsum [RP2]var(standlrt) * Redisplay the model results. Present the fixed part parameters as odds * ratios. Present the variance parameters as standard deviations and the * covariance parameters as correlations. Present the modes of each chain * rather than the means, runmlwin, nogroup mode or sd correlation **************************************************************************** * 5 EXPORT MODELS TO WINBUGS **************************************************************************** * Generate the WinBUGS model, initial values and data files for this model runmlwin passexam cons standlrt girl, /// level2(school: cons standlrt) /// level1(student:) /// discrete(d(binomial) l(logit) de(cons)) /// mcmc(b(500) c(5000) savewinbugs(model(m.txt) inits(i.txt) data(d.txt) nofit)) /// initsprevious nogroup nopause * View the WinBUGS model file view m.txt ***************************************************************************** exit