Page 1 of 1
Sixway: save accuracy statistics / summary diagnostics
Posted: Mon Jan 22, 2018 3:02 pm
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.
Re: Sixway: save accuracy statistics / summary diagnostics
Posted: Tue Jan 23, 2018 9:44 am
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.
Re: Sixway: save accuracy statistics / summary diagnostics
Posted: Tue Jan 23, 2018 1:10 pm
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?
Re: Sixway: save accuracy statistics / summary diagnostics
Posted: Tue Jan 23, 2018 5:57 pm
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))