------------------------------------------------------------------------------- name: log: Q:\C-modelling\runmlwin\website\logfiles\2020-03-27\16\8_Running_a > _Simulation_Study_in_MLwiN.smcl log type: smcl opened on: 27 Mar 2020, 17:49:33 . **************************************************************************** . * 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 number of observations (_N) was 0, now 6 . . gen school = _n . . gen u = sqrt(30)*invnorm(uniform()) . . expand 18 (102 observations created) . . 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 MLwiN 3.05 multilevel model Number of obs = 108 Normal response model (hierarchical) Estimation algorithm: IGLS ----------------------------------------------------------- | No. of Observations per Group Level Variable | Groups Minimum Average Maximum ----------------+------------------------------------------ school | 6 18 18.0 18 ----------------------------------------------------------- Run time (seconds) = 0.72 Number of iterations = 2 Log likelihood = -371.45705 Deviance = 742.91411 ------------------------------------------------------------------------------ resp | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- cons | 9.127995 1.23475 7.39 0.000 6.707929 11.54806 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ Level 2: school | var(cons) | 6.179135 5.29773 -4.204225 16.56249 -----------------------------+------------------------------------------------ Level 1: pupil | var(cons) | 53.43324 7.482151 38.7685 68.09799 ------------------------------------------------------------------------------ . . . . * 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 { 2. // 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' 3. . quietly { 4. . drop _all 5. . set obs 6 6. . gen school = _n 7. . gen u = sqrt(30)*invnorm(uniform()) 8. // NOTE: For Stata 10 and above, you can use: . // gen u = rnormal(0,sqrt(30)) . . expand 18 9. . bysort school: gen pupil = _n 10. . gen cons = 1 11. . gen e = sqrt(40)*invnorm(uniform()) 12. // NOTE: For Stata 10 and above, you can use: . // gen e = rnormal(0,sqrt(40)) . . gen resp = 30*cons + u + e 13. . runmlwin resp cons, /// > level2(school: cons) /// > level1(pupil: cons) /// > nopause 14. . scalar b0_igls = [FP1]_b[cons] 15. . scalar s2u_igls = [RP2]_b[var(cons)] 16. . scalar s2e_igls = [RP1]_b[var(cons)] 17. . scalar cov_b0_igls = (inrange(30, b0_igls - 1.96*[FP1]_se[con > s], b0_igls + 1.96*[FP1]_se[cons])) 18. . scalar cov_s2u_igls = (inrange(10, s2u_igls - 1.96*[RP2]_se[v > ar(cons)], s2u_igls + 1.96*[RP2]_se[var(cons)])) 19. . scalar cov_s2e_igls = (inrange(40, s2e_igls - 1.96*[RP1]_se[v > ar(cons)], s2e_igls + 1.96*[RP1]_se[var(cons)])) 20. . runmlwin resp cons, /// > level2(school: cons) /// > level1(pupil: cons) /// > mcmc(burnin(500) chain(10001)) /// > initsprevious nopause 21. . mcmcsum, getchains 22. . rename FP1_cons b0 23. . rename RP2_var_cons_ s2u 24. . rename RP1_var_cons_ s2e 25. . foreach par in b0 s2u s2e { 26. . egen `par'_mcmc = mean(`par') 27. . egen `par'_lo_mcmc = pctile(`par'), p(2.5) 28. . egen `par'_hi_mcmc = pctile(`par'), p(97.5) 29. . } 30. . keep in 1 31. . gen cov_b0_mcmc = (inrange(30, b0_lo_mcmc, b0_hi_mcmc)) 32. . gen cov_s2u_mcmc = (inrange(10, s2u_lo_mcmc, s2u_hi_mcmc)) 33. . gen cov_s2e_mcmc = (inrange(40,s2e_lo_mcmc, s2e_hi_mcmc)) 34. . 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) 35. . } 36. . } Iteration: 1 Iteration: 2 Iteration: 3 Iteration: 4 Iteration: 5 Iteration: 6 Iteration: 7 Iteration: 8 Iteration: 9 Iteration: 10 . . 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) stats | b0_igls s2u_igls s2e_igls ---------+------------------------------ mean | 31.44 18.54 37.55 ---------------------------------------- . // 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) stats | b0_mcmc s2u_mcmc s2e_mcmc ---------+------------------------------ mean | 31.16 36.39 38.51 ---------------------------------------- . . tabstat s2u_mcmc, format(%3.2f) stat(median) variable | p50 -------------+---------- s2u_mcmc | 31.56 ------------------------ . . tabstat cov_b0_igls cov_s2u_igls cov_s2e_igls /// > cov_b0_mcmc cov_s2u_mcmc cov_s2e_mcmc, format(%3.2f) stats | cov_b0~s c~u_igls c~e_igls cov_b0~c c~u_mcmc c~e_mcmc ---------+------------------------------------------------------------ mean | 0.80 1.00 1.00 0.80 0.70 1.00 ---------------------------------------------------------------------- . . . . . * Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .110 . . . . . . **************************************************************************** . exit end of do-file