R/summary.metadiag.R

In bamdit: Bayesian Meta-Analysis of Diagnostic Test Data

#' Generic summary function for metadiag object in bamdit
#'
#' @param object The object generated by the metadiag function.
#'
#' @param digits The number of significant digits printed. The default value is 3.
#'
#' @param ... \dots
#'
#' @export
#'

summary.metadiag = function(object, digits = 3, ...) {
  bugs.output = object$BUGSoutput
  bugs.summary = bugs.output$summary

  summary.m = list()

  model.spec = list()
  model.spec$model.type = object$re.model
  model.spec$link = object$link
  model.spec$re = object$re
  if (model.spec$re == "sm") {
    model.spec$split.w = object$split.w
    model.spec$df = object$df
  }
  summary.m$model.specification = model.spec

  if (model.spec$model.type == "DS"){
    summary.m$summary.fe = rbind(bugs.summary["mu.D",],
                                 bugs.summary["mu.S",])
    row.names(summary.m$summary.fe) = c("mu.D", "mu.S")
    summary.m$summary.re = rbind(bugs.summary["sigma.D",],
                                 bugs.summary["sigma.S",],
                                 bugs.summary["rho",])
    row.names(summary.m$summary.re) = c("sigma.D", "sigma.S", "rho")
    summary.m$summary.cor = cor(bugs.output$sims.list$mu.D,
                                bugs.output$sims.list$mu.S)
  } else {
    summary.m$summary.fe = rbind(bugs.summary["mu.Se",],
                                 bugs.summary["mu.Sp",])
    row.names(summary.m$summary.fe) = c("mu.Se", "mu.Sp")
    summary.m$summary.re = rbind(bugs.summary["sigma.Se",],
                                 bugs.summary["sigma.Sp",],
                                 bugs.summary["rho",])
    row.names(summary.m$summary.re) = c("sigma.Se", "sigma.Sp", "rho")
    summary.m$summary.cor = cor(bugs.output$sims.list$mu.Se,
                                bugs.output$sims.list$mu.Sp)
  }

  summary.m$summary.expected.accuracy = rbind(bugs.summary["se.pool",],
                                              bugs.summary["sp.pool",])
  row.names(summary.m$summary.expected.accuracy) = c("se.pool", "sp.pool")

  summary.m$summary.pred.expected.accuracy = rbind(bugs.summary["se.new",],
                                                   bugs.summary["sp.new",])
  row.names(summary.m$summary.pred.expected.accuracy) = c("se.new", "sp.new")

  summary.m$DIC = bugs.output$DIC
  summary.m$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.m$mcmc.setup = mcmc.setup

  class(summary.m) = "summary.metadiag"

  print(summary.m, digits=digits)
  #return(summary.m)
}

print.summary.metadiag = function(x, digits, ...) {
  cat('Model specification:\n')
  model.spec = x$model.specification
  cat(paste('  Parametrization: ', model.spec$model.type, sep = ''))
  cat('\n')
  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('  Degrees of freedom: ', model.spec$df, sep = ''))
    cat('\n')
  }

  cat('\n')
  cat('Fixed effects: \n')
  print(round(x$summary.fe, digits))
  if (model.spec$model.type == "DS") {
    cat(paste('Correlation between mu.D and mu.S is ', round(x$summary.cor, digits), '.\n', sep = ''))
  } else {
    cat(paste('Correlation between mu.Se and mu.Sp is ', round(x$summary.cor, digits), '.\n', sep = ''))
  }

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

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

  cat('Expected accuracy:\n')
  print(round(x$summary.expected.accuracy, digits))
  cat('\n')
  cat('Predictive accuracy:\n')
  print(round(x$summary.pred.expected.accuracy, 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')

  if (model.spec$model.type == "DS") {
    cat('The current output can NOT be used for ReviewManager (RevMan) software. Please, run the model with re.model = "SeSp"')
  } else {

    cat('The current output can be used for ReviewManager (RevMan) software for writing Cochrane Reviews.')
  }
}