R/plot_meta.R

In rcalinjageman/esci: Estimation Statistics with Confidence Intervals

#' Generates a forest plot displaying results of a meta-analysis
#'
#' @description `plot_meta' returns a ggplot2 object visualizing the results of
#'   a meta-analysis, showing each study effect size and CI, the overall effect
#'   size and CI as a diamond, effect sizes estimated at each moderator level
#'   (if defined), and (optionally) prediction intervals for subsequent studies.
#'   This function requires as input an esci_estimate object generated by an
#'   esci meta-analysis function: [esci::meta_any()], [esci::meta_d1()],
#'   [esci::meta_d2()], [esci::meta_mdiff_two()], [esci::meta_mean()],
#'   [esci::meta_pdiff_two()], [esci::meta_proportion()], and [esci::meta_r()].
#'
#'
#' @inherit plot_describe details
#'
#'
#' @param estimate - an esci_estimate object generated by an esci meta_ function
#'
#' @param mark_zero - Boolean; defaults to TRUE to include a dotted line
#'   indicated no effect (effect_size = 0)
#' @param include_PIs - Boolean; defaults to FALSE; set to TRUE to include
#'   prediction intervals for the overall effect and each moderator level (if
#'   defined)
#' @param report_CIs - Boolean; defaults to FALSE; set to TRUE to include
#'   printed representation of each study effect size and CI along the right-
#'   hand of the figure
#' @param explain_DR - Boolean; defaults to FALSE; set to TRUE if no moderator
#'   is defined to show both the RE and FE effect sizes to represent how the
#'   diamond ration measure of effect-size heterogeneity is calculated
#' @param meta_diamond_height - Optional real number > 0 to indicate height that
#'   each meta-analytic diamond should be drawn; defaults to 0.35
#' @param ggtheme - Optional ggplot2 theme object to control overall styling;
#'   defaults to [ggplot2::theme_classic()]
#'
#'
#' @inherit plot_describe return
#'
#'
#' @inherit meta_mdiff_two examples
#'
#' @export
plot_meta <- function(
    estimate,
    mark_zero = TRUE,
    include_PIs = FALSE,
    report_CIs = FALSE,
    explain_DR = FALSE,
    meta_diamond_height = 0.35,
    ggtheme = ggplot2::theme_classic()
) {

  UL <- LL <- effect_size <- weight <- yendv <- PI_LL <- PI_UL <- pline <- NULL

  if (include_PIs & is.na(estimate$es_meta$PI_LL[[1]])) {
    include_PIs <- FALSE
    warning("No prediction intervals present, ignoring include_PIs = TRUE")
  }
  if (explain_DR & include_PIs) {
    explain_DR <- FALSE
    warning("Prediction intervals requested; ignoring explain_DR = TRUE")
  }
  if (explain_DR & !is.null(estimate$raw_data$moderator)) {
    warning("Meta-analysis has a moderator; ignoring explain_DR = TRUE")
    explain_DR <- FALSE
  }

  # ------------------------------------
  # Basic plot
  myplot <- ggplot2::ggplot() + ggtheme

  # Mark 0 ?
  if (mark_zero) {
    myplot <- myplot + ggplot2::geom_vline(
      xintercept = 0,
      color = "black",
      linetype = "dashed"
    )
    myplot <- esci_plot_layers(myplot, "ref_zero")
  }

  # --------------------------------------
  # Raw data prep
  rdata <- estimate$raw_data

  # Match moderator levels to contrast levels
  if (is.null(rdata$moderator)) {
    rdata$type <- "Reference"
  } else {
    contrast <- estimate$properties$contrast
    reference_levels <- names(contrast[contrast < 0])
    comparison_levels <- names(contrast[contrast > 0])
    unused_levels <- names(contrast[contrast == 0])

    rdata[rdata$moderator %in% reference_levels, "type"] <- "Reference"
    rdata[rdata$moderator %in% comparison_levels, "type"] <- "Comparison"
    if (length(unused_levels) > 0) {
      rdata[rdata$moderator %in% unused_levels, "type"] <- "Unused"
    }
  }

  rows_raw <- nrow(rdata)
  rdata$line <- -seq(1:rows_raw)
  rdata$CI_label <- paste(
    "[",
    format(rdata$LL, digits = 2),
    ", ",
    format(rdata$UL, digits = 2),
    "]",
    sep = ""
  )

  # Raw data plot
  plot_levels <- c(
    "#008DF9" = "Reference",
    "#009F81" = "Comparison",
    "gray65" = "Unused"
  )
  for (x in 1:length(plot_levels)) {
    next_up <- plot_levels[x]
    this_data <- rdata[rdata$type == next_up, ]
    if (nrow(this_data) > 0) {
      myplot <- myplot + ggplot2::geom_errorbarh(
        data = this_data,
        ggplot2::aes(
          xmin = LL,
          xmax = UL,
          y = line
        ),
        height = 0
      )
      myplot <- esci_plot_layers(myplot, paste("raw_", next_up, "_error", sep = ""))

      myplot <- myplot + ggplot2::geom_point(
        data = this_data,
        ggplot2::aes(
          x = effect_size,
          y = line,
          size = weight
        ),
        shape = 22,
        colour = "black",
        fill = names(plot_levels[x])
      )
      myplot <- esci_plot_layers(myplot, paste("raw_", next_up, "_point", sep = ""))
    }
  }



  # -------------------------------------------
  # Group data
  gdata <- estimate$es_meta
  if (is.null(gdata$moderator_variable_level)) {
    gdata$type <- "Overall"
    gdata$moderator_variable_level <- "Overall"

    if (explain_DR) {
      gdata <- rbind(gdata, gdata)
      gdata$moderator_variable_level <- c(
        "RE: Overall",
        "FE: Overall"
      )
      gdata$effect_size[[1]] <- gdata$RE_effect_size[[1]]
      gdata$effect_size[[2]] <- gdata$FE_effect_size[[1]]
      gdata$LL[[1]] <- gdata$RE_effect_size[[1]] - gdata$RE_CI_width[[1]]/2
      gdata$UL[[1]] <- gdata$RE_effect_size[[1]] + gdata$RE_CI_width[[1]]/2
      gdata$LL[[2]] <- gdata$FE_effect_size[[1]] - gdata$FE_CI_width[[1]]/2
      gdata$UL[[2]] <- gdata$FE_effect_size[[1]] + gdata$FE_CI_width[[1]]/2
    }

  } else {
    gdata$type <- "Overall"
    gdata[gdata$moderator_variable_level %in% reference_levels, "type"] <- "Reference"
    gdata[gdata$moderator_variable_level %in% comparison_levels, "type"] <- "Comparison"
    if (length(unused_levels) > 0) {
      gdata[gdata$moderator_variable_level %in% unused_levels, "type"] <- "Unused"
    }
  }

  # Set lines of the meta-analysis for group data
  #   Ordering overall first, unused next, reference, and then comparison
  rows_start <- rows_raw+1
  rows_total <- rows_raw
  line_groups <- c("Overall", "Unused", "Reference", "Comparison")
  for (up_next in line_groups) {
    if (nrow(gdata[gdata$type == up_next, ]) > 0) {
      increment <- if(include_PIs) 2 else 1
      rows_total <- rows_total + nrow(gdata[gdata$type == up_next, ]) * increment
      gdata[gdata$type == up_next, "line"] <- -seq(
        from = rows_start,
        to = rows_total,
        by = increment
      )
      rows_start <- rows_total + 1
    }
  }

  gdata$CI_label <- paste(
    "[",
    format(gdata$LL, digits = 2),
    ", ",
    format(gdata$UL, digits = 2),
    "]",
    sep = ""
  )

  if (explain_DR) {
    gdata$CI_label <- paste(
      "Length: ",
      format(gdata$UL - gdata$LL, digits = 2)
    )
  }

  gdata$PI_label <- paste(
    "(",
    format(gdata$PI_LL, digits = 2),
    ", ",
    format(gdata$PI_UL, digits = 2),
    ")",
    sep = ""
  )

  if (include_PIs) {
    gdata$pline <- gdata$line - 1
  }

  # Connecting lines
  # Connecting lines for differences
  # Prep connecting lines
  if (!is.null(estimate$es_meta_difference)) {
    plot_levels <- c(
      "dashed" = "Reference",
      "dotted" = "Comparison"
    )
    gdata$yendv <- -Inf
    gdata[gdata$type == "Comparison", "yendv"] <- -(rows_total + 2)

    for (x in 1:length(plot_levels)) {
      next_up <- plot_levels[x]
      dline <- gdata[gdata$type == next_up, ]
      if (nrow(dline) > 0) {
        myplot <- myplot + ggplot2::geom_segment(
          data = dline,
          ggplot2::aes(
            x = effect_size,
            xend = effect_size,
            y = line,
            yend = yendv
          ),
          colour = "black",
          linetype = names(plot_levels[x])
        )
        myplot <- esci_plot_layers(myplot, paste("ref_", next_up, "_lines", sep = ""))
      }
    }
  }

  # Plot group data
  plot_levels <- c(
    "Black" = "Overall",
    "#008DF9" = "Reference",
    "#009F81" = "Comparison",
    "gray65" = "Unused"
  )
  for (x in 1:length(plot_levels)) {
    next_up <- plot_levels[x]
    this_data <- gdata[gdata$type == next_up, ]
    if (nrow(this_data) > 0) {

      myplot <- myplot + geom_meta_diamond_h(
        data = this_data,
        ggplot2::aes(
          x = effect_size,
          xmin = LL,
          xmax = UL,
          y = line
        ),
        height = meta_diamond_height,
        color = "black",
        fill = names(plot_levels[x]),
      )
      myplot <- esci_plot_layers(myplot, paste("group_", next_up, "_diamond", sep = ""))

      if(include_PIs) {
        myplot <- myplot + geom_segment(
          data = this_data,
          ggplot2::aes(
            x = PI_LL,
            xend = PI_UL,
            y = pline,
            yend = pline
          ),
          color = names(plot_levels[x]),
          linewidth = 2,
          alpha = 0.5
        )
        myplot <- esci_plot_layers(myplot, paste("group_", next_up, "_PI", sep = ""))
      }

    }
  }


  # Difference data
  dlabel <- NULL
  difference_CI_label <- NULL
  expand_bottom <- 0.05
  if (!is.null(estimate$es_meta_difference)) {
    expand_bottom <- 0

    # Prep difference data
    ddata <- estimate$es_meta_difference["Difference", ]
    difference_CI_label <- paste(
      "[",
      format(ddata$LL, digits = 2),
      ", ",
      format(ddata$UL, digits = 2),
      "]",
      sep = ""
    )
    difference_CI_label <- c("", difference_CI_label)

    reference <- estimate$es_meta_difference["Reference", "effect_size"]
    ddata$effect_size <- ddata$effect_size + reference
    ddata$LL <- ddata$LL + reference
    ddata$UL <- ddata$UL + reference

    rows_total <- rows_total + 2
    ddata$line <- -rows_total

    # Plot reference data


    myplot <- myplot + ggplot2::geom_errorbarh(
      data = ddata,
      ggplot2::aes(
        xmin = LL,
        xmax = UL,
        y = line
      ),
      height = 0
    )
    myplot <- esci_plot_layers(myplot, "group_Difference_line")

    myplot <- myplot + ggplot2::geom_point(
      data = ddata,
      ggplot2::aes(
        x = effect_size,
        y = line,
      ),
      size = 6,
      shape = 24,
      colour = "black",
      fill = "black"
    )

    #
    #
    # myplot <- myplot + geom_meta_diamond_h(
    #   data = ddata,
    #   ggplot2::aes(
    #     x = effect_size,
    #     xmin = LL,
    #     xmax = UL,
    #     y = line
    #   ),
    #   height = meta_diamond_height,
    #   colour = "Black",
    #   fill = "White"
    # )
    myplot <- esci_plot_layers(myplot, "group_Difference_diamond")

    dlabel <- c("", ddata$moderator_level)
  }

  # ----------------------------------------------
  # Axis
  # X axis labels and setup
  xlab <- paste(
    estimate$properties$effect_size_name_ggplot,
    ": ",
    estimate$es_meta$effect_label[[1]],
    sep = ""
  )
  myplot <- myplot + ggplot2::scale_x_continuous(
    name = xlab,
    position = "top"
  )
  myplot$esci_xlab <- xlab

  # If needed, difference axis
  if (!is.null(estimate$es_meta_difference)) {
    myplot <- esci_plot_difference_axis_x(
      myplot,
      estimate$es_meta_difference, d_lab = "Difference axis"
    )
  }

  # Y axis labels and setup
  gdata_labels <- gdata[order(-gdata$line), ]$moderator_variable_level

  gdata_labels[which(gdata_labels == "Overall")] <- if (estimate$properties$model == "fixed_effects") "FE overall" else "RE overall"

  cils <- paste(
    gdata_labels,
    " ",
    estimate$properties$conf_level * 100,
    "% CI",
    sep = ""
  )

  if (include_PIs) {

    pils <- paste(
      gdata_labels,
      " ",
      estimate$properties$conf_level * 100,
      "% PI",
      sep = ""
    )
    gdata_labels <- c(rbind(cils, pils))
  } else {
    gdata_labels <- cils
  }

  all_labels <- c(
    rdata$label,
    gdata_labels,
    dlabel
  )

  if (is.null(estimate$es_meta_difference)) {
    expand_bottom <- 0.05
  } else {
    expand_bottom <- 0
  }

  if (report_CIs) {
    glabel <- gdata$CI_label
    if (include_PIs) {
      glabel <- c(rbind(gdata$CI_label, gdata$PI_label))
    }

    all_CI_labels <- c(
      rdata$CI_label,
      glabel,
      difference_CI_label
    )

    sec_axis_CIs <- ggplot2::sec_axis(
      name = NULL,
      breaks = -seq(1:rows_total),
      trans = ~.-0,
      labels = all_CI_labels
    )

    myplot$esci_sec_axis_CIs <- sec_axis_CIs

    myplot <- myplot + ggplot2::scale_y_continuous(
      name = NULL,
      breaks = -seq(1:rows_total),
      labels = all_labels,
      sec.axis = sec_axis_CIs #,
      #expand = ggplot2::expansion(mult = c(expand_bottom, 0.05), add = c(meta_diamond_height, 0))
    )

  } else {
    myplot <- myplot + ggplot2::scale_y_continuous(
      name = NULL,
      breaks = -seq(1:rows_total),
      labels = all_labels #,
      #expand = ggplot2::expansion(mult = c(expand_bottom, 0.05), add = c(meta_diamond_height, 0))
    )
  }

  # -------------------------------
  # Theme tweaks
  myplot <- myplot + ggplot2::theme(
    legend.position = "none",
    axis.line.y.right = element_blank(),
    axis.ticks.y.right = element_blank(),
    axis.title.x.top = ggtext::element_markdown()
  )

  return(myplot)
}

#' Add a difference axis to the x axis of an esci forest plot
#'
#'
#' `esci_plot_difference_axis_x` can be used to redraw the
#' difference axis from a forest plot created with [esci::plot_meta].
#' You must pass the plot returned from plot_meta and the effect size
#' table containing the estimated difference.
#'
#'
#' @param  myplot required ggplot2 plot returned from a plot_meta
#' function
#' @param difference_table required data frame from an esci-estimate that
#' has a difference-based effect size
#' @param dlim Optional 2-item vector to provide the lower and upper
#' boundaries of the difference axis.  Defaults to c(NA, NA) which is to
#' auto-set both boundaries.
#' @param d_n.breaks Optional numeric > 2 to suggest number of breaks for the
#' difference axis; defaults to NULL in which case number of breaks
#' is handled automatically by ggplot
#' @param d_lab Optional character to serve as the label for the difference
#' axis; defaults to NULL
#' @export
esci_plot_difference_axis_x <- function(
    myplot,
    difference_table,
    dlim = c(NA, NA),
    d_n.breaks = NULL,
    d_lab = NULL
) {

  # From debugging
  # myplot <- plot_meta(estimate)
  # difference_table <- estimate$es_meta_difference

  # If previous difference axis exists, remove it
  if (!is.null(myplot$layers$ref_difference_axis)) {
    myplot$layers$ref_difference_axis <- NULL
  }

  # Store y axis info
  y_labels <- ggplot2::ggplot_build(myplot)$layout$panel_params[[1]]$y$get_labels()
  y_breaks <- ggplot2::ggplot_build(myplot)$layout$panel_params[[1]]$y$get_breaks()
  y_min <- min(ggplot2::layer_scales(myplot)$y$range$range)
  y_max <- max(ggplot2::layer_scales(myplot)$y$range$range)
  xlim <- ggplot2::ggplot_build(myplot)$layout$panel_params[[1]]$x$limits
  xbreaks <- ggplot2::ggplot_build(myplot)$layout$panel_params[[1]]$x$breaks

  # Get reference value and center LL, UL, and effect size on it
  reference_value <- difference_table["Reference", "effect_size"]
  adj_r <- c("Comparison", "Reference")
  adj_c <- c("effect_size", "LL", "UL")
  difference_table[adj_r, adj_c] <- difference_table[adj_r, adj_c] - reference_value

  # Get initial lower and upper, either specified or based on difference table
  if (is.na(dlim[[1]])) {
    lower <- min(0, difference_table$LL, na.rm = TRUE)
  } else {
    lower <- dlim[[1]]
  }

  if (is.na(dlim[[2]])) {
    upper <- max(0, difference_table$UL, na.rm = TRUE)
  } else {
    upper <- dlim[[2]]
  }

  # Get initial distance between breaks, either for set number of based on SE
  if (is.null(d_n.breaks)) {
    bdist <- difference_table["Difference", "SE"]
    d_n.breaks <- 4
  } else {
    bdist <- (upper - lower) / d_n.breaks
  }

  # Adjust upper and lower to ensure 0 is included
  new_upper <- ceiling(upper/bdist) * bdist
  new_lower <- floor(lower/bdist) * bdist

  # Make difference axis breaks and labels
  d_breaks <- pretty(c(new_lower, new_upper), d_n.breaks)
  new_upper <- d_breaks[[1]]
  new_lower <- d_breaks[length(d_breaks)]
  # d_breaks <- seq(from = new_lower, to = new_upper, by = bdist)
  d_labels <- format(d_breaks, digits = 2)

  # Assemble and apply difference axis
  my_sec_axis <- ggplot2::sec_axis(
    name = d_lab,
    trans = ~.-reference_value,
    breaks = d_breaks,
    labels = d_labels
  )


  myplot <- myplot + ggplot2::scale_x_continuous(
    limits = xlim,
    breaks = xbreaks,
    position = "top",
    name = myplot$esci_xlab,
    sec.axis = my_sec_axis
  )

  # Just plot a truncated axis
  myplot <- myplot + ggplot2::guides(
    x.sec = ggh4x::guide_axis_truncated(
      trunc_lower = new_lower + reference_value,
      trunc_upper = new_upper + reference_value
    )
  )

  if (!is.null(d_lab)) {
    nx_min <- min(ggplot2::layer_scales(myplot)$x$range$range)
    nx_max <- max(ggplot2::layer_scales(myplot)$x$range$range)

    x_range <- nx_max - nx_min
    daxis_top <- (new_upper + reference_value) - nx_min
    daxis_top_location <- daxis_top / x_range

    daxis_half_length <- (new_upper - new_lower)/2
    daxis_half_location <- daxis_half_length / x_range

    myplot <- myplot + ggplot2::theme(
      axis.title.x.bottom = element_text(hjust = daxis_top_location - daxis_half_location)
    )

  }


  if (is.null(myplot$esci_sec_axis_CIs)) {
    myplot <- myplot + ggplot2::scale_y_continuous(
      name = NULL,
      breaks = y_breaks,
      labels = y_labels #,
      #expand = ggplot2::expansion(mult = c(0, 0.05), add = c(0, 0))
    )
  } else {
    myplot <- myplot + ggplot2::scale_y_continuous(
      name = NULL,
      breaks = y_breaks,
      labels = y_labels,
      sec.axis = myplot$esci_sec_axis_CIs #,
      #expand = ggplot2::expansion(mult = c(0, 0.05), add = c(0, 0))
    )
  }


  return(myplot)

}