Dears,
I have a problem with MCMC in R2MLwiN. I want to build up threelevel model with hierarchical centering (HC), but now MLwiN only allows to select one level (level2 or level3) to incorporate HC. The MCMCweb mentioned that it's okay to apply HC with threelevel model, but I don't know how to incorporate it in R2MLwiN. Is there anyone could give some suggestions or provide some sample code?
Another problem is about the MCMC diagnostic result. When I use sixway(), it will automatically plot the diagnostic results in Graph Device. However, I have many models to check and want to extract some diagnostic indices within the graph (especially BrooksDraper diagnostic statistics, BD). As BD() is an internal function within sixway(), just wander how to generate this statistics or is there any method to extract the BD index from the graph and then output it into a .csv file. I know that there is a mcmcsum() function in Stata which could generate BD index, is there any similar function in R?
Thanks and regards.
Diva
Question about hierarchical centering and MCMC diagnostics

 Posts: 1071
 Joined: Mon Oct 19, 2009 10:34 am
Re: Question about hierarchical centering and MCMC diagnostics
The R replication code for the examples in the MCMC manual can be found at http://www.bristol.ac.uk/cmm/media/r2ml ... CGuide25.R. If you prefer to specify the hierarchical centring at level3 instead of level2 you would simply set the option hcen to equal 3 instead of 2.
Extracting the code to calculate BD from the sixway function you could do something like the following:
Taking the formula from BD.R the last line could be replaced with:
Extracting the code to calculate BD from the sixway function you could do something like the following:
Code: Select all
library(R2MLwiN)
data(tutorial, package = "R2MLwiN")
(mymodel < runMLwiN(normexam ~ 1 + standlrt + (1  student), estoptions = list(EstM = 1), data = tutorial))
paramchain < mymodel@chains[, "FP_Intercept"]
rho < acf(paramchain, 100, plot=FALSE)$acf[2]
R2MLwiN:::BD(mean(paramchain), var(paramchain), rho)
Code: Select all
ceiling(4 * qnorm(1  0.05/2)^2 * (sqrt(var(paramchain))/(10^(floor(log10(abs(mean(paramchain))))  2 + 1)))^2 * (1 + rho)/(1  rho))