Sixway: save accuracy statistics / summary diagnostics

[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]

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]

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]

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