Question about hierarchical centering and MCMC diagnostics

Welcome to the forum for R2MLwiN users. Feel free to post your question about R2MLwiN here. The Centre for Multilevel Modelling take no responsibility for the accuracy of these posts, we are unable to monitor them closely. Do go ahead and post your question and thank you in advance if you find the time to post any answers!

Go to R2MLwiN: Running MLwiN from within R >>
Post Reply
Posts: 13
Joined: Tue Sep 26, 2017 2:35 pm

Question about hierarchical centering and MCMC diagnostics

Post by divasong »


I have a problem with MCMC in R2MLwiN. I want to build up three-level model with hierarchical centering (HC), but now MLwiN only allows to select one level (level-2 or level-3) to incorporate HC. The MCMC-web mentioned that it's okay to apply HC with three-level 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 Brooks-Draper 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.


Posts: 1165
Joined: Mon Oct 19, 2009 10:34 am

Re: Question about hierarchical centering and MCMC diagnostics

Post by ChrisCharlton »

The R replication code for the examples in the MCMC manual can be found at ... CGuide25.R. If you prefer to specify the hierarchical centring at level-3 instead of level-2 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:

Code: Select all

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)
Taking the formula from BD.R the last line could be replaced with:

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))

Post Reply