R/plotAverage.R

In Dani-Schulze/measurementInvariance: (Partial) measurement invariance analysis

### # Thu May 27 17:43:54 2021 ------------------------------

#' @title plotAverage
#'
#' @description Plotting the distribution of a parameter resulting from Bayesian
#' model averaging
#'
#' @param res_modelAveraging  Object generated by previous use
#' \code{\link{modelAveraging}}.
#' @param ... further arguments passed to \code{hist()}, except \code{ylim}.
#'
#' @return A histogram with 95% credibility interval lines.
#'
#' @usage plotAverage(res_modelAveraging = NULL)
#'
#' @export

plotAverage <- function(res_modelAveraging,
                        ...) {
  res <- res_modelAveraging

  range <- max(res$muAver) - min(res$muAver)
  suppressWarnings({
    h <- hist(res$muAver, plot = FALSE, ...) # to get height
  })
  # dealing with ellipsis arguments
  args <- list(main = "Model averaging",
               xlab = "Mean difference",
               ylab = "",
               xlim = c(min(res$muAver) - 0.2*range,
                        max(res$muAver) + 0.2*range),
               yaxt = "n")
  extraArgs <- list(...)
  if (any(names(extraArgs) == "ylim")) { # keeping ylim's defaults anytime
    extraArgs$ylim <- NULL
  }
  args[names(extraArgs)] <- extraArgs
  do.call(hist,
          c(list(x = res$muAver,
                 ylim = c(0, 1.3*h$counts[which.max(h$counts)])),
            args))

  abline(v = res$muAver[order(res$muAver)][0.025*length(res$muAver)],
         col = "blue",
         lwd = 2)
  abline(v = res$muAver[order(res$muAver)][0.975*length(res$muAver)],
         col = "blue",
         lwd = 2)
  abline(v = mean(res$muAver),
         col = "red",
         lwd = 2)

  rect(min(res$muAver),
       1.02*h$counts[which.max(h$counts)],
       max(res$muAver),
       1.4*h$counts[which.max(h$counts)],
       col = "white",
       border = NA)

  text(res$muAver[order(res$muAver)][0.025*length(res$muAver)],
       1.2*h$counts[which.max(h$counts)],
       round(res$muAver[order(res$muAver)][0.025*length(res$muAver)], 3))
  text(res$muAver[order(res$muAver)][0.975*length(res$muAver)],
       1.2*h$counts[which.max(h$counts)],
       round(res$muAver[order(res$muAver)][0.975*length(res$muAver)], 3))
  text(mean(res$muAver),
       1.2*h$counts[which.max(h$counts)],
       round(mean(res$muAver), 3))

}