R/brmsfit-methods.R

In mvuorre/vmisc: Vuorre Miscellaneous

#' Meta-analysis forest plots from brms output
#'
#' This function draws a forest plot from a random-effects meta-analysis
#' model fitted with brms.
#'
#' @param model A meta-analytic model estimated with brms.
#' @param level The "Confidence" level for the Credible Intervals.
#' Defaults to 0.95.
#' @param xlim Limits for the x-axis. The defaults are the smallest observed
#' ES - 4 Standard errors to the largest observed ES + 4 SEs.
#' @param show_data Logical; whether to show the observed effect size
#' and standard error below the meta-analytic estimates. Defaults to FALSE.
#' @param sort_estimates Logical; whether to sort the estimates in ascending
#' order of magnitude from bottom to top. Defaults to FALSE.
#' @param dens_fill String; color to fill the densities. Defaults to "grey60".
#' @param dens_col String; color for the outlines of the densities. Default: NA.
#'
#' @details
#'
#' \subsection{Requirements for brms model}{
#'
#' The meta-analytic model must be fitted with the following brms formula:
#' \code{yi | se(sei) ~ 1 + (1|study)}
#'
#' }
#'
#' @examples
#' \dontrun{
#' ## Use a data frame called d
#' fit <- brm(yi | se(sei) ~ 1 + (1|study), data = d)
#' forest(fit)
#' }
#' @method forest brmsfit
#' @export
forest <- function(model) {

    grouping <- model$ranef$group
    parameter <- model$ranef$coef

    # Obtain samples of each varying intercept
    samples_r <- coef(model, summary = FALSE)[[grouping]][,,parameter]
    samples_r <- tidyr::gather_(as.data.frame(samples_r),
                                key_col = grouping,
                                value_col = parameter,
                                gather_cols = model$data[[grouping]])

    # Add samples of population-level intercept
    samples_f <- as.data.frame(fixef(model, summary = FALSE))
    samples_f[[grouping]] <- "ME"

    # Combine
    samples <- rbind(samples_f, samples_r)

    # Summaries
    sum_r <- coef(model)[[grouping]][,,parameter]
    sum_r <- as.data.frame(sum_r)
    sum_r[[grouping]] <- row.names(sum_r)
    sum_f <- fixef(model)
    sum_f <- as.data.frame(sum_f)
    sum_f[[grouping]] <- "ME"
    samples_sum <- rbind(sum_r, sum_f)
    samples_sum[["Interval"]] <- paste0(round(samples_sum[["Estimate"]], 2),
                                        " [",
                                        round(samples_sum[["2.5%ile"]], 2),
                                        ", ",
                                        round(samples_sum[["97.5%ile"]], 2), "]")

    # Order effects
    samples_sum <- dplyr::arrange(samples_sum, Estimate)
    samples_sum[["order"]] <- row.names(samples_sum)
    # Put ME in bottom with an empty row
    samples_sum[["order"]] <- ifelse(samples_sum[[grouping]]=="ME", -1, samples_sum[["order"]])
    samples_sum[[grouping]] <- reorder(
        samples_sum[[grouping]],
        as.numeric(samples_sum[["order"]])
        )


    # Create graph
    g <- ggplot(samples_sum, aes_string(parameter, grouping))
    g <- g + geom_text(data = samples_sum,
                       aes(label = Interval, x = max(samples[[parameter]])),
                       hjust = "inward")
    g <- g + ggridges::geom_density_ridges(data=samples, rel_min_height = 0.01, scale = 0.95, col = NA)
    g <- g + geom_point(data = samples_sum, aes(Estimate))
    g <- g + geom_segment(data = samples_sum,
                          aes_string(yend = grouping,
                                     x = "`2.5%ile`",
                                     xend = "`97.5%ile`"))
    g
}