R/crosnma.summary.R

In TasnimHamza/crosnma:

#' Generic summary function for crosnma.run object
#'@param x A crosrun object generated by the crosnma.run function.
#'@param digits The number of significant digits printed. The default value is 3.
#'@param quantiles A numeric vector of probabilities to present posterior summaries. The default value is c(0.025, 0.5, 0.975) for the 95% credible interval and the median.
#'@param expo A logical value indicating whether to exponentiate the parameters of relative treatment effect and covariate effect. Default is TRUE.
#'@return \code{summary.crosnma} returns a matrix containing the following summary statistics (in columns) for each estimated parameter:
#'@return \code{Mean} the mean of the posterior distribution
#'@return \code{SD} the standard deviation of the posterior distribution
#'@return \code{2.5\%} (default) the 2.5\% quantile of the posterior distribution (the lower bound of the 95% credible interval)
#'@return \code{50\%} (default) the median of the posterior distribution
#'@return \code{97.5\%} (default) the 97.5\% quantile of the posterior distribution (the upper bound of the 95% credible interval)
#'@return \code{Rhat} Gelman-Rubin statistic. The further the value of Rhat from 1, the worse the mixing of chains and so the convergence.
#'@return \code{n.eff} an estimate of the effective sample size. The smaller the value of n.eff the greater the uncertainty associated with the corresponding parameter.
#'
#'@export
#'
#'
#'

summary.crosnma <- function(x, digits = 3,  quantiles = c(0.025, 0.5, 0.975), expo=TRUE, ...)
{

  # initial
  sum.fit <-summary((x$samples),digits=digits,quantiles = quantiles, ...)
  # mean and sd
  sum.fit.mat <- sum.fit$statistics[,c('Mean','SD')]
  # add quantiles
  sum.fit.mat <- cbind(sum.fit.mat, sum.fit$quantiles)
  if (expo) {
    sum.fit.mat[row.names(sum.fit.mat)[!startsWith(rownames(sum.fit.mat),"tau")],-2]<- exp(sum.fit.mat[row.names(sum.fit.mat)[!startsWith(rownames(sum.fit.mat),"tau")],-2])
    message("Mean and quantiles are exponentiated")
  }
  # add Rhat
  Rhat <- round(rhat(x$samples)$psrf[,'Point est.'],3)
  sum.fit.mat <- cbind(sum.fit.mat,Rhat)

  # add number of effective sample size
  n.eff <- round(effectiveSize(x$samples))
  sum.fit.mat <- cbind(sum.fit.mat,n.eff)


  # attach treatment names to d's
  row.names(sum.fit.mat)[startsWith(rownames(sum.fit.mat),"d")] <- paste("d",c(x$trt.key$trt.ini),sep=".")
  return(round(sum.fit.mat,digits=digits))
}