Sixway: save accuracy statistics / summary 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 >> http://www.bris.ac.uk/cmm/software/r2mlwin/
Post Reply
GerineLodder
Posts: 42
Joined: Thu Jul 28, 2016 11:33 am

Sixway: save accuracy statistics / summary diagnostics

Post by GerineLodder »

Is it possible to save the output from sixway function (accuracy statistics and summary diagnostics) without the figures in an output?

I have to run about a million analyses (slight exaggeration) and I would just like 1 output with all of the numbers. I already save the trajectories and sixway outputs to a different folder, so I also have the visual material, but I would like to have an output file with the numbers combined in 1 file as well.
ChrisCharlton
Posts: 1351
Joined: Mon Oct 19, 2009 10:34 am

Re: Sixway: save accuracy statistics / summary diagnostics

Post by ChrisCharlton »

There is no automatic way to do this, however if you look at the sixway function code (https://github.com/rforge/r2mlwin/blob/ ... R/sixway.R) it should be reasonably easy to see where all the numbers come from, so you could use this information to extract the information that you want to record.
GerineLodder
Posts: 42
Joined: Thu Jul 28, 2016 11:33 am

Re: Sixway: save accuracy statistics / summary diagnostics

Post by GerineLodder »

Thanks! Unfortunately, as I am not used to programming with R code, relatively easy is still not feasible for me (I can adopt a manual to my own data, but that's about it...).
Could you give an example of how to extract, for instance, Nhat and quartiles for a specific parameter are derived?
ChrisCharlton
Posts: 1351
Joined: Mon Oct 19, 2009 10:34 am

Re: Sixway: save accuracy statistics / summary diagnostics

Post by ChrisCharlton »

Note, that this code could be tidied up a bit, I have more of less copied it directly from the sixway function:

Code: Select all

# Convert chains object (corresponding to the chains argument of the sixway function) into a matrix (may not be necessary)
flatchain <- as.matrix(chain)

# Raftery-Lewis nhat for 0.025 quantile: (raftery.diag function is in the coda package)
RL1 <- raftery.diag(flatchain, q = 0.025, r = 0.005, s = 0.95, converge.eps = 0.001)
# nhat is:
RL1$resmatrix[1, "N"]

# Raftery-Lewis nhat for 0.975 quantile: (raftery.diag function is in the coda package)
RL2 <- raftery.diag(flatchain, q = 0.975, r = 0.005, s = 0.95, converge.eps = 0.001)
# nhat is:
RL2$resmatrix[1, "N"]

# Brooks-Draper nhat: (BD function is in the R2MLwiN package)
acf.maxlag <- 100
aa <- acf(flatchain, acf.maxlag, plot=FALSE)
rho <- aa$acf[2]
# nhat is
Ndb <- BD(mean(flatchain), var(flatchain), rho, k = 2, alpha = 0.05)

# Quantiles of the chain
quants <- quantile(flatchain, c(0.025, 0.05, 0.5, 0.95, 0.975))
Post Reply