****************************************************************************
*     MLwiN MCMC Manual 
*
* 4   Other Features of Variance Components Models . . . . . . . . . . . .45
*
*     Browne, W. J. (2009). MCMC Estimation in MLwiN, v2.26. Centre for
*     Multilevel Modelling, University of Bristol.
****************************************************************************
*     Stata do-file to replicate all analyses using runmlwin
*
*     George Leckie and Chris Charlton, 
*     Centre for Multilevel Modelling, 2012
*     http://www.bristol.ac.uk/cmm/software/runmlwin/
****************************************************************************

* 4.1 Metropolis Hastings (MH) sampling for the variance components model 46

use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear

runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons) ///
	nopause

estimates store igls

matrix b = e(b)

matrix V = e(V)

quietly runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons)  ///
	mcmc(on) initsb(b) initsv(V) nopause 

estimates store gibbs



* 4.2 Metropolis-Hastings settings . . . . . . . . . . . . . . . . . . . .47



* 4.3 Running the variance components with Metropolis Hastings . . . . . .48

quietly runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons)  ///
	mcmc(femethod(univariatemh) remethod(univariatemh)) ///
	initsb(b) initsv(V) nopause

estimates store mh_adaptive

mcmcsum [FP1]standlrt, detail

mcmcsum [FP1]standlrt, fiveway

quietly runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons)  ///
	mcmc(femethod(univariatemh) remethod(univariatemh) noadapt) ///
	initsb(b) initsv(V) nopause

estimates store mh_scale	

estimates table igls gibbs mh_scale mh_adaptive, b(%4.3f)




* 4.4 MH cycles per Gibbs iteration . . . . . . . . . . . . . . . . . . . 49



* 4.5 Block updating MH sampling . . . . . . . . . . . . . . . . . . . . .49

runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons)  ///
	mcmc(femethod(multivariatemh) remethod(univariatemh) acceptance(0.4)) ///
	initsb(b) initsv(V) nopause

mcmcsum in 4951/5000, trajectories



* 4.6 Residuals in MCMC . . . . . . . . . . . . . . . . . . . . . . . . . 51

quietly runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons) ///
	nopause

runmlwin normexam cons standlrt, ///
	level2(school: cons, residuals(u, savechains("u.dta", replace))) ///
	level1(student: cons) ///
	mcmc(chain(5001)) ///
	initsprevious nopause

use "u.dta", clear

keep iteration school value

rename value u0_

reshape wide u0, i(iteration) j(school)

mcmcsum u0_1, variables detail

mcmcsum u0_1, variables fiveway



* 4.7 Comparing two schools . . . . . . . . . . . . . . . . . . . . . . . 54

gen diff = u0_1 - u0_2

mcmcsum diff, variables detail

mcmcsum diff, variables fiveway

gen prop = (diff>0)

sum prop



* 4.8 Calculating ranks of schools . . . . . . . . . . . . . . . . . . . .55

use "u.dta", clear

keep iteration school value

rename value u0

bysort iteration: egen u0rank = rank(u0), unique

drop u0

bysort school: egen u0rankmn = mean(u0rank)

bysort school: egen u0ranklo = pctile(u0rank), p(2.5)

bysort school: egen u0rankmd = pctile(u0rank), p(50)

bysort school: egen u0rankhi = pctile(u0rank), p(97.5)

keep school u0rankmn u0ranklo u0rankmd u0rankhi

duplicates drop

twoway (scatter u0rankmd school) (rcap u0rankhi u0ranklo school), ///
	ytitle("Rank") xtitle("School") legend(off)

egen rank = rank(u0rankmn), unique

twoway (scatter u0rankmn rank) (rcap u0rankhi u0ranklo rank), ///
	ytitle("Rank") xtitle("School") legend(off)



* 4.9 Estimating a function of parameters . . . . . . . . . . . . . . . . 58

use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear

quietly runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons) ///
	nopause

quietly runmlwin normexam cons standlrt, ///
	level2(school: cons) ///
	level1(student: cons) ///
	mcmc(chain(5001)) ///
	initsprevious nopause

mcmcsum, getchains

gen isc = RP2_var_cons_/( RP2_var_cons_ +  RP1_var_cons_)

mcmcsum isc, variables detail

mcmcsum isc, variables fiveway



* Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . . 60





****************************************************************************
exit