****************************************************************************
*     MLwiN MCMC Manual 
*
* 8   Running a Simulation Study in MLwiN . . . . . . . . . . . . . . . . 97
*
*     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/
****************************************************************************

* 8.1 JSP dataset simulation study . . . . . . . . . . . . . . . . . . . .97

* 8.2 Setting up the structure of the dataset . . . . . . . . . . . . . . 98

clear

set seed 12345

set obs 6

gen school = _n

gen u = sqrt(30)*invnorm(uniform())

expand 18

bysort school: gen pupil = _n

gen cons = 1

gen e = sqrt(40)*invnorm(uniform())

gen resp = 10*cons + u + e

runmlwin resp cons, ///
	level2(school: cons) ///
	level1(pupil: cons) ///
	nopause



* 8.3 Generating simulated datasets based on true values . . . . . . . . 102

postutil clear

postfile simstudy ///
	b0_igls s2u_igls s2e_igls ///
	cov_b0_igls cov_s2u_igls cov_s2e_igls ///
	b0_mcmc s2u_mcmc s2e_mcmc ///
	cov_b0_mcmc cov_s2u_mcmc cov_s2e_mcmc ///
	using "simstudy.dta", replace
	
forvalues iteration = 1/10 {
// Note: To replicate the analysis presented in the manual as cloesly as
// possible, increase the maximum number of iterations from 10 to 100.

	display "Iteration: " `iteration'
	
	quietly {

		drop _all

		set obs 6

		gen school = _n

		gen u = sqrt(30)*invnorm(uniform())
		// NOTE: For Stata 10 and above, you can use: 
		//       gen u = rnormal(0,sqrt(30))

		expand 18

		bysort school: gen pupil = _n

		gen cons = 1

		gen e = sqrt(40)*invnorm(uniform())
		// NOTE: For Stata 10 and above, you can use: 
		//       gen e = rnormal(0,sqrt(40))

		gen resp = 30*cons + u + e

		runmlwin resp cons, ///
			level2(school: cons) ///
			level1(pupil: cons) ///
			nopause

		scalar b0_igls = [FP1]_b[cons]
		
		scalar s2u_igls = [RP2]_b[var(cons)]

		scalar s2e_igls = [RP1]_b[var(cons)]

		scalar cov_b0_igls = (inrange(30, b0_igls - 1.96*[FP1]_se[cons], b0_igls + 1.96*[FP1]_se[cons]))

		scalar cov_s2u_igls = (inrange(10, s2u_igls - 1.96*[RP2]_se[var(cons)], s2u_igls + 1.96*[RP2]_se[var(cons)]))

		scalar cov_s2e_igls = (inrange(40, s2e_igls - 1.96*[RP1]_se[var(cons)], s2e_igls + 1.96*[RP1]_se[var(cons)]))

		runmlwin resp cons, ///
			level2(school: cons) ///
			level1(pupil: cons) ///
			mcmc(burnin(500) chain(10001)) ///
			initsprevious nopause

		mcmcsum, getchains

		rename FP1_cons b0

		rename RP2_var_cons_ s2u

		rename RP1_var_cons_ s2e

		foreach par in b0 s2u s2e {

			egen `par'_mcmc = mean(`par')

			egen `par'_lo_mcmc = pctile(`par'), p(2.5)

			egen `par'_hi_mcmc = pctile(`par'), p(97.5)

		}
		
		keep in 1
		
		gen cov_b0_mcmc = (inrange(30, b0_lo_mcmc, b0_hi_mcmc))

		gen cov_s2u_mcmc = (inrange(10, s2u_lo_mcmc, s2u_hi_mcmc))

		gen cov_s2e_mcmc = (inrange(40,s2e_lo_mcmc, s2e_hi_mcmc))

		post simstudy ///
			(b0_igls) (s2u_igls) (s2e_igls) ///
			(cov_b0_igls) (cov_s2u_igls) (cov_s2e_igls) ///
			(b0_mcmc) (s2u_mcmc) (s2e_mcmc) ///
			(cov_b0_mcmc) (cov_s2u_mcmc) (cov_s2e_mcmc)

	}

}

postclose simstudy



* 8.4 Fitting the model to the simulated datasets . . . . . . . . . . . .106
// Note: The analysis that appears in this section of the manual is carried
// out under Section 8.3 of this do-file



* 8.5 Analysing the simulation results . . . . . . . . . . . . . . . . . 109

use "simstudy.dta", clear

tabstat b0_igls s2u_igls s2e_igls , format(%3.2f)
// Note: These results are smililar, but not identical to those reported in
// the MCMC manual. The reason for the discrepancies will be due to the
// different seed used in Stata to generate the data from the seed used in
// MLwiN to generate the data. This comment also applies to the results
// presented below.

tabstat	b0_mcmc s2u_mcmc s2e_mcmc, format(%3.2f)

tabstat s2u_mcmc, format(%3.2f) stat(median)

tabstat cov_b0_igls cov_s2u_igls cov_s2e_igls ///
	cov_b0_mcmc cov_s2u_mcmc cov_s2e_mcmc, format(%3.2f)




* Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .110





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