Skip to content

Instantly share code, notes, and snippets.

@dwinter
Created May 24, 2013 03:57
Show Gist options
  • Save dwinter/5641202 to your computer and use it in GitHub Desktop.
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
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