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.
Sixway: save accuracy statistics / summary diagnostics

 Posts: 42
 Joined: Thu Jul 28, 2016 11:33 am

 Posts: 1111
 Joined: Mon Oct 19, 2009 10:34 am
Re: Sixway: save accuracy statistics / summary diagnostics
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.

 Posts: 42
 Joined: Thu Jul 28, 2016 11:33 am
Re: Sixway: save accuracy statistics / summary diagnostics
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?
Could you give an example of how to extract, for instance, Nhat and quartiles for a specific parameter are derived?

 Posts: 1111
 Joined: Mon Oct 19, 2009 10:34 am
Re: Sixway: save accuracy statistics / summary diagnostics
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)
# RafteryLewis 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"]
# RafteryLewis 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"]
# BrooksDraper 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))