R/summary.hmr.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

#' Generic summary function for hmr object in jarbes
#' @param object The object generated by the hmr function.
#'
#' @param digits The number of significant digits printed. The default value is 3.
#' @param ... \dots
#'
#' @export
#'
summary.hmr = function(object, digits = 3, ...) {
  bugs.output = object$BUGSoutput
  bugs.summary = bugs.output$summary

  summary.h = list()

  model.spec = list()
  model.spec$model.type = object$re.model
  model.spec$link = object$link
  model.spec$re = object$re
  model.spec$split.w = object$split.w
  model.spec$df.estimate = object$df.estimate
  summary.h$model.specification = model.spec

  # FE's
  summary.h$summary.fe = bugs.summary[c("mu.1", "mu.2", "beta.0", "beta.1",
                                        "Odds.pool", "P_control.pool"),]
  row.names(summary.h$summary.fe)[c(3, 4)] = c("intercept", "slope")

  # RE's
  rows = c("sigma.1", "sigma.2", "rho")
  if (model.spec$re == "sm" & model.spec$df.estimate == TRUE)
    rows = c(rows, "df")
  summary.h$summary.re = bugs.summary[rows,]

  # bias parameter
  summary.h$summary.bias.par = matrix(bugs.summary["mu.phi",], nrow = 1,
                                      dimnames = list(c("mu.phi"),
                                                      colnames(bugs.summary)))

  # predictive effects
  summary.h$summary.predictive.effects = bugs.summary[c("Odds.new",
                                                        "P_control.new"),]

  # IPD Predictors
  ind.names = grep("beta.IPD", row.names(bugs.summary))
  summary.h$summary.IPD.predictors = rbind(bugs.summary["sigma.beta",],
                                           bugs.summary[ind.names,])
  row.names(summary.h$summary.IPD.predictors) = c("sigma.beta (regularization parameter)",
                                                  object$beta.names)

  # MCMC info
  summary.h$DIC = bugs.output$DIC
  summary.h$pD = bugs.output$pD

  mcmc.setup = list()
  mcmc.setup$n.chains = bugs.output$n.chains
  mcmc.setup$n.iter = bugs.output$n.iter
  mcmc.setup$n.burnin = bugs.output$n.burnin
  summary.h$mcmc.setup = mcmc.setup

  class(summary.h) = "summary.hmr"

  print(summary.h, digits = digits)
}

print.summary.hmr = function(x, digits, ...) {
  cat('Model specification:\n')
  model.spec = x$model.specification
  cat(paste('  Random effects: ', model.spec$re, sep = ''))
  cat('\n')
  cat(paste('  Link function: ', model.spec$link, sep = ''))
  cat('\n')
  if (model.spec$re == "sm") {
    cat(paste('  Split weights: ', model.spec$split.w, sep = ''))
    cat('\n')
    cat(paste('  Estimate degrees of freedom: ', model.spec$df, sep = ''))
    cat('\n')
  }

  cat('\n')
  cat('Fixed effects: \n')
  print(round(x$summary.fe, digits))

  cat('\n')
  cat('Random effects: \n')
  print(round(x$summary.re, digits))

  cat('\n')
  cat('Bias parameters: \n')
  print(round(x$summary.bias.par, digits))

  cat('\n-------------------\n')

  cat('Predictive effects:\n')
  print(round(x$summary.predictive.effects, digits))

  cat('\n-------------------\n')

  cat('Idividual Participant Data Predictors:\n')
  print(round(x$summary.IPD.predictors, digits))

  cat('\n-------------------\n')

  mcmc = x$mcmc.setup
  cat(paste('MCMC setup (fit using jags): ', mcmc$n.chains, ' chains, each with ', mcmc$n.iter, ' iterations (first ', mcmc$n.burnin, ' discarded)', sep = ''))
  cat('\n')
  cat(paste('DIC: ', round(x$DIC, digits), sep = ''))
  cat('\n')
  cat(paste('pD: ', round(x$pD, digits), sep = ''))
  cat('\n')
}