Created
May 24, 2013 03:57
-
-
Save dwinter/5641202 to your computer and use it in GitHub Desktop.
Quckly run a bunch of diagnostic tests on a list of MCMC runs
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
model_diagnostics <- function(modlist, plot=TRUE){ | |
get_param_ranges <- function(x) { | |
res <- apply((sapply(x, posterior.mode)),1,range) | |
rownames(res) <- c("min", "max") | |
return(res) | |
} | |
n <- length(modlist) | |
on.exit(devAskNewPage(devAskNewPage())) #ie send it back the way it was | |
random_mcmc <- mcmc.list(lapply(modlist, "[[", "VCV")) | |
fixed_mcmc <- mcmc.list(lapply(modlist, "[[", "Sol")) | |
#gelman diagnostics | |
gelman_random <- gelman.diag(random_mcmc) | |
gelman_fixed <- gelman.diag(fixed_mcmc) | |
#paramater estimates | |
fixed_ests <- get_param_ranges(fixed_mcmc) | |
random_ests <- get_param_ranges(random_mcmc) | |
#autocorrelations | |
autocorr_random <- autocorr(random_mcmc) | |
autocorr_fixed <- autocorr(fixed_mcmc) | |
res <- list(gelman_random=gelman_random, gelman_fixed=gelman_fixed, | |
autocorr_random=autocorr_random, autocorr_fixed=autocorr_fixed, | |
fixed_ests = fixed_ests, random_ests=random_ests, n=n) | |
class(res) <- "MCMCglmm_diagnostic" | |
if (plot){ | |
devAskNewPage(TRUE) | |
gelman.plot(random_mcmc) | |
gelman.plot(fixed_mcmc) | |
autocorr.plot(random_mcmc) | |
autocorr.plot(fixed_mcmc) | |
} | |
return(res) | |
} | |
print.MCMCglmm_diagnostic <- function(x, ...){ | |
on.exit(cat("\n")) | |
cat("MCMC diagnostics from", x$n, "model runs\n\n") | |
cat("Gelman Diagnostics for fixed effects:\n") | |
print(x$gelman_fixed) | |
cat("\nGelman Diagnostics for random effects\n") | |
print(x$gelman_fixed) | |
cat("Autocorrelation stats for fixed effects:\n") | |
print(x$autocorr_fixed) | |
cat("\nAutocorrelation stats for random effects\n") | |
print(x$autocorr_fixed) | |
cat("Ranges of values for fixed-effect parameter estimates\n") | |
print(x$random_ests) | |
cat("\nRanges of values for random-effect parameter estimates\n") | |
print(x$fixed_ests) | |
} | |
#example | |
x <- rnorm(250) | |
y <- x * 0.5 + rnorm(250,15, 1) | |
df <- data.frame(x,y,boxes=rep(letters[1:5]), each=50) | |
mods <- replicate(4, MCMCglmm(y~x, random=~boxes, data=df), simplify=FALSE) | |
print(model_diagnostics(mods, plot=FALSE)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment