fbseq: fbseq

#' @title Function \code{contrast_estimates}
#' @description Extracts estimated posterior means, standard deviations and credible intervals
#' of linear combinations of beta parameters.
#' @export
#' @return a data frame of marginal posterior parameter estimates: posterior means, poseterior standard deviations, and approximate credible intervals.
#' Parameters not updeted in the MCMC are excluded.
#' @param obj a \code{Chain} object or list of \code{Chain} objects generated by \code{fbseq()}.
#' @param level level of the credible intervals from 0 to 1
contrast_estimates = function(obj, level = 0.95){
  if(class(obj) != "list") obj = list(obj)
  n = sum(sapply(obj, function(ch) ch@iterations * ch@thin))
  Mean = rowSums(sapply(obj, function(ch) ch@contrastsPostMean * ch@iterations * ch@thin))/n
  MeanSq = rowSums(sapply(obj, function(ch) ch@contrastsPostMeanSquare * ch@iterations * ch@thin))/n
  Sd =  sqrt(n*(MeanSq - Mean^2)/(n - 1))
  p = 1 - (1 - level)/2
  lower = qnorm(1 - p, mean = Mean, sd = Sd)
  upper = qnorm(p, mean = Mean, sd = Sd)
  out = data.frame(Mean, Sd, lower, upper)
  colnames(out) = c("mean", "sd", paste0(c("lower", "upper"), 
        "_ci_", level))
  ns = names(Scenario(obj[[1]])@contrasts)
  gs = rownames(Scenario(obj[[1]])@counts)
  rownames(out) = paste0(rep(ns, each = obj[[1]]@G), "_", rep(gs, times = length(ns)))
  out
}