Page 1 of 1

Question about hierarchical centering and MCMC diagnostics

Posted: Tue Oct 24, 2017 1:44 pm
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.


Re: Question about hierarchical centering and MCMC diagnostics

Posted: Tue Oct 24, 2017 3:08 pm
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))