R/JAGSrun.R

In bayesmix: Bayesian Mixture Models with JAGS

JAGSrun <-  function(y, prefix = yname, model = BMMmodel(k = 2),
                     control = JAGScontrol(variables = c("mu", "tau", "eta")),
                     tmp = TRUE, cleanup = TRUE, ...) {
  yname <- deparse(substitute(y))
  if (!is.null(dim(y))) {
    if (dim(y)[1] == 1) y <- y[1,]
    else if (dim(y)[2] == 1) y <- y[,1]
    else stop("Only univariate response allowed")
  }
  y <- as.numeric(y)
  cl <- match.call()
  if (tmp) {
    dir <- getwd()
    tmpdir <- tempdir()
    if (!file.exists(tmpdir)){
      if (!dir.create(tmpdir)) stop("Error creating tmp directory")
    }
    setwd(tmpdir)
    on.exit(setwd(dir))
  }
  results <- JAGScall(model, y, prefix, control, ...)
  if (cleanup) {
    unlink(paste(prefix, ".bug", sep = ""))
  }
  z <- list(call = cl, results = results$results, model = results$model,
            variables = results$variables, data = y)
  class(z) <- "JAGSrun"
  z
}

# Print method for JAGSrun objects adapted from print.mcmc and
# summary.mcmc in package coda written by Martyn Plummer, Nicky Best,
# Kate Cowles, Karen Vines

print.JAGSrun <- function(x, ...) {
  cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
  if (inherits(x$model, "BMMmodel")) {
    if (is.null(x$results)) cat("No results!\n")
    else {
      cat("Markov Chain Monte Carlo (MCMC) output:\nStart =", stats::start(x$results), 
          "\nEnd =", stats::end(x$results), "\nThinning interval =", coda::thin(x$results), 
          "\n")
      for (i in x$variables) {
        y <- x$results[,grep(paste("^", i, sep = ""), colnames(x$results)), drop = FALSE]
        if(dim(y)[2] <=  x$model$data$k) {
          yout <- summaryShort_mcmc(y)
          class(yout) <- "summaryShort_mcmc"
          cat(paste("\n Empirical mean, standard deviation and 95% CI for", i, "\n"))
          print(yout, ...)
        }
      }
    }
  }
}

summaryShort_mcmc <- function (object, quantiles = c(0.025, 0.975), 
                               ...) {
  x <- methods::as(object, "mcmc")
  statnames <- c("Mean", "SD")
  varstats <- matrix(nrow = coda::nvar(x), ncol = length(statnames), 
                     dimnames = list(coda::varnames(x), statnames))
  if (is.matrix(x)) {
    xmean <- apply(x, 2, mean)
    xvar <- apply(x, 2, stats::var)
    varquant <- t(apply(x, 2, stats::quantile, quantiles))
  }
  else {
    xmean <- mean(x, na.rm = TRUE)
    xvar <- stats::var(x, na.rm = TRUE)
    varquant <- stats::quantile(x, quantiles)
  }
  varstats[, 1] <- xmean
  varstats[, 2] <- sqrt(xvar)
  varstats <- drop(varstats)
  varquant <- drop(varquant)
  out <- list(statistics = varstats, quantiles = varquant,
              start = stats::start(x), end = stats::end(x),
              thin = coda::thin(x), nchain = 1)
  class(out) <- "summaryShort_mcmc"
  return(out)
}

print.summaryShort_mcmc <- function(x, digits = max(3, .Options$digits - 3), ...) {
  if (is.matrix(x$statistics)) {
    print(cbind(x$statistics, x$quantiles), digits = digits, ...)
  }
  else print(c(x$statistics, x$quantiles), digits = digits, ...)
}