R/tabmeta-plot.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ============================================================================
# R4VN tabmeta plotting methods
# ============================================================================

.r4vn_meta_analysis_ref <- function(x, ref = NULL) {
  if (is.null(ref)) {
    if (x$config$scale == "log") return(0)
    if (x$config$scale == "zcor") return(0)
    return(x$config$ref)
  }
  ref <- as.numeric(ref)[1L]
  if (!is.finite(ref)) stop("`ref` must be finite.", call. = FALSE)
  if (x$config$scale == "log") {
    if (ref <= 0) stop("`ref` must be > 0 for ratio/rate measures.",
                       call. = FALSE)
    return(log(ref))
  }
  if (x$config$scale == "zcor") {
    if (abs(ref) >= 1) stop("Correlation `ref` must be between -1 and 1.",
                            call. = FALSE)
    return(atanh(ref))
  }
  ref
}

.r4vn_meta_natural_ref <- function(x, ref = NULL) {
  if (!is.null(ref)) return(as.numeric(ref)[1L])
  x$config$ref
}

.r4vn_meta_merge_plot_args <- function(defaults, engine_args,
                                       protected = character()) {
  if (is.null(engine_args)) engine_args <- list()
  if (!is.list(engine_args) ||
      (length(engine_args) && is.null(names(engine_args)))) {
    stop("`engine_args` must be a named list.", call. = FALSE)
  }
  engine_args[intersect(names(engine_args), protected)] <- NULL
  utils::modifyList(defaults, engine_args)
}

.r4vn_meta_plot_file <- function(
    x, type, file, width = 1800, height = 1400, res = 180,
    moderator = NULL, contour = FALSE, title = NULL, subtitle = NULL,
    caption = NULL, xlab = NULL,
    ref = NULL, xlim = NULL, ticks = NULL, color = NULL,
    font_family = "sans", text_size = 0.82, axis_size = 0.9,
    title_size = 1.05, title_color = "black",
    subtitle_size = 0.9, subtitle_color = "gray30",
    caption_size = 0.75, caption_color = "gray40",
    margins = NULL, background = "white",
    point_color = NULL, point_bg = "white", ci_color = NULL,
    summary_color = NULL, summary_border = NULL,
    point_shape = NULL, point_size = NULL,
    line_type = 1, line_width = 1,
    ref_color = "gray40", ref_type = 2, ref_width = 1,
    show_weights = TRUE, show_prediction = FALSE,
    weight_title = "Weight", estimate_title = NULL,
    show_abcd = FALSE, abcd_titles = c("a", "b", "c", "d"),
    prediction_style = "bar", row_shade = "zebra",
    shade_color = "gray95", digits = NULL, study_order = NULL,
    header = NULL, annotate = TRUE, yaxis = "sei", ylab = NULL,
    contour_levels = c(90, 95, 99),
    contour_colors = c("#FEE2E2", "#FEF3C7", "#E5E7EB"),
    funnel_label = FALSE, funnel_legend = FALSE,
    bubble_ci = TRUE, bubble_pi = FALSE,
    bubble_shade = c("#DBEAFE", "#E5E7EB"), grid = FALSE,
    engine_args = list()) {

  grDevices::png(
    filename = file, width = width, height = height, res = res,
    bg = background
  )
  on.exit(grDevices::dev.off(), add = TRUE)
  .r4vn_meta_draw(
    x = x, type = type, moderator = moderator, contour = contour,
    title = title, subtitle = subtitle, caption = caption,
    xlab = xlab, ref = ref, xlim = xlim, ticks = ticks,
    color = color, font_family = font_family, text_size = text_size,
    axis_size = axis_size, title_size = title_size,
    title_color = title_color, subtitle_size = subtitle_size,
    subtitle_color = subtitle_color, caption_size = caption_size,
    caption_color = caption_color, margins = margins,
    background = background,
    point_color = point_color, point_bg = point_bg,
    ci_color = ci_color, summary_color = summary_color,
    summary_border = summary_border, point_shape = point_shape,
    point_size = point_size, line_type = line_type,
    line_width = line_width, ref_color = ref_color,
    ref_type = ref_type, ref_width = ref_width,
    show_weights = show_weights, show_prediction = show_prediction,
    weight_title = weight_title, estimate_title = estimate_title,
    show_abcd = show_abcd, abcd_titles = abcd_titles,
    prediction_style = prediction_style, row_shade = row_shade,
    shade_color = shade_color, digits = digits,
    study_order = study_order, header = header, annotate = annotate,
    yaxis = yaxis, ylab = ylab, contour_levels = contour_levels,
    contour_colors = contour_colors, funnel_label = funnel_label,
    funnel_legend = funnel_legend, bubble_ci = bubble_ci,
    bubble_pi = bubble_pi, bubble_shade = bubble_shade, grid = grid,
    engine_args = engine_args
  )
  invisible(file)
}

.r4vn_meta_draw <- function(
    x, type = "forest", moderator = NULL, contour = FALSE,
    title = NULL, subtitle = NULL, caption = NULL, xlab = NULL,
    ref = NULL, xlim = NULL, ticks = NULL,
    color = NULL, font_family = "sans", text_size = 0.82,
    axis_size = 0.9, title_size = 1.05, title_color = "black",
    subtitle_size = 0.9, subtitle_color = "gray30",
    caption_size = 0.75, caption_color = "gray40",
    margins = NULL, background = "white",
    point_color = NULL, point_bg = "white",
    ci_color = NULL, summary_color = NULL, summary_border = NULL,
    point_shape = NULL, point_size = NULL, line_type = 1,
    line_width = 1, ref_color = "gray40", ref_type = 2,
    ref_width = 1, show_weights = TRUE, show_prediction = FALSE,
    weight_title = "Weight", estimate_title = NULL,
    show_abcd = FALSE, abcd_titles = c("a", "b", "c", "d"),
    prediction_style = "bar", row_shade = "zebra",
    shade_color = "gray95", digits = NULL, study_order = NULL,
    header = NULL, annotate = TRUE, yaxis = "sei", ylab = NULL,
    contour_levels = c(90, 95, 99),
    contour_colors = c("#FEE2E2", "#FEF3C7", "#E5E7EB"),
    funnel_label = FALSE, funnel_legend = FALSE,
    bubble_ci = TRUE, bubble_pi = FALSE,
    bubble_shade = c("#DBEAFE", "#E5E7EB"), grid = FALSE,
    engine_args = list()) {

  .r4vn_meta_need()
  type <- tolower(as.character(type)[1L])
  allowed <- c(
    "forest", "funnel", "baujat", "influence", "leaveout",
    "cumulative", "radial", "labbe", "bubble", "trimfill"
  )
  if (!type %in% allowed) {
    stop(
      "`type` must be one of: ", paste(allowed, collapse = ", "), ".",
      call. = FALSE
    )
  }

  if (is.null(title)) {
    title <- switch(
      type,
      forest = x$title,
      funnel = "Funnel plot",
      trimfill = "Trim-and-fill funnel plot",
      baujat = "Baujat plot",
      influence = "Influence diagnostics",
      leaveout = "Leave-one-out sensitivity analysis",
      cumulative = "Cumulative meta-analysis",
      radial = "Radial plot",
      labbe = "L'Abbe plot",
      bubble = if (!is.null(moderator))
        paste0("Meta-regression: ", moderator) else "Meta-regression"
    )
  }
  if (is.null(xlab)) xlab <- x$config$label
  col_main <- if (is.null(color)) "black" else as.character(color)[1L]
  if (is.null(point_color)) point_color <- col_main
  if (is.null(ci_color)) ci_color <- point_color
  if (is.null(summary_color)) summary_color <- col_main
  if (is.null(summary_border)) summary_border <- summary_color
  if (is.null(digits)) digits <- x$digits
  if (!is.null(show_prediction)) {
    if (!is.logical(show_prediction) || length(show_prediction) != 1L ||
        is.na(show_prediction)) {
      stop("`show_prediction` must be TRUE, FALSE, or NULL.", call. = FALSE)
    }
  }
  if (!is.logical(show_abcd) || length(show_abcd) != 1L ||
      is.na(show_abcd)) {
    stop("`show_abcd` must be TRUE or FALSE.", call. = FALSE)
  }
  if (isTRUE(show_abcd) &&
      !all(c(".a", ".b", ".c", ".d") %in% names(x$analysis_data))) {
    warning(
      "`show_abcd=TRUE` is available when `tabmeta()` receives a, b, c, d input.",
      call. = FALSE
    )
    show_abcd <- FALSE
  }
  if (length(abcd_titles) != 4L) {
    stop("`abcd_titles` must contain four column headings.", call. = FALSE)
  }
  old_par <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(old_par), add = TRUE)
  graphics::par(
    family = as.character(font_family)[1L], bg = background,
    lwd = as.numeric(line_width)[1L],
    cex.main = title_size, col.main = title_color,
    cex.lab = axis_size, cex.axis = axis_size
  )
  if (!is.null(margins)) {
    margins <- as.numeric(margins)
    if (length(margins) != 4L || any(!is.finite(margins))) {
      stop("`margins` must contain four finite values.", call. = FALSE)
    }
    graphics::par(mar = margins)
  }
  back <- .r4vn_meta_back_fun(x)
  finish <- function() {
    if (!is.null(subtitle) && length(subtitle) &&
        nzchar(as.character(subtitle)[1L])) {
      graphics::mtext(
        as.character(subtitle)[1L], side = 3, line = 0.5,
        cex = subtitle_size, col = subtitle_color
      )
    }
    if (!is.null(caption) && length(caption) &&
        nzchar(as.character(caption)[1L])) {
      graphics::mtext(
        as.character(caption)[1L], side = 1, line = 3.2,
        cex = caption_size, col = caption_color
      )
    }
    invisible(x)
  }

  if (type == "forest") {
    prediction_on <- isTRUE(show_prediction) && isTRUE(x$random)
    forest_header <- if (is.null(header)) {
      c(
        "Study",
        if (is.null(estimate_title)) {
          paste0(x$config$label, " (", round(x$ci * 100), "% CI)")
        } else as.character(estimate_title)[1L]
      )
    } else header
    forest_method <- tryCatch(
      utils::getS3method(
        "forest", "rma", optional = TRUE,
        envir = asNamespace("metafor")
      ),
      error = function(e) NULL
    )
    supports_colci <- !is.null(forest_method) &&
      "colci" %in% names(formals(forest_method))
    supports_ilab_lab <- !is.null(forest_method) &&
      "ilab.lab" %in% names(formals(forest_method))
    ci_lty <- line_type[1L]
    prediction_lty <- if (length(line_type) >= 2L) {
      line_type[2L]
    } else {
      ref_type[1L]
    }
    heading_lty <- if (length(line_type) >= 3L) {
      line_type[3L]
    } else {
      ci_lty
    }
    args <- list(
      x = x$primary_model,
      slab = x$analysis_data$.study,
      atransf = back,
      xlab = xlab,
      refline = NA_real_,
      header = forest_header,
      annotate = isTRUE(annotate),
      showweights = FALSE,
      addpred = prediction_on,
      shade = row_shade,
      colshade = shade_color,
      pch = if (is.null(point_shape)) 15 else point_shape,
      colout = point_color,
      col = summary_color,
      border = summary_border,
      lty = c(ci_lty, prediction_lty, heading_lty),
      fonts = font_family,
      cex = text_size,
      cex.lab = axis_size,
      cex.axis = axis_size,
      digits = digits,
      annosym = c(" (", ", ", ")")
    )
    if (prediction_on) args$predstyle <- prediction_style

    zcrit <- stats::qnorm(1 - (1 - x$ci) / 2)
    observed_limits <- range(
      c(
        x$analysis_data$yi - zcrit * sqrt(x$analysis_data$vi),
        x$analysis_data$yi + zcrit * sqrt(x$analysis_data$vi),
        as.numeric(x$primary_model$ci.lb),
        as.numeric(x$primary_model$ci.ub)
      ),
      finite = TRUE
    )
    if (length(observed_limits) != 2L ||
        any(!is.finite(observed_limits)) ||
        observed_limits[1L] >= observed_limits[2L]) {
      observed_limits <- c(-1, 1)
    }

    effect_limits <- NULL
    if (!is.null(xlim)) {
      effect_limits <- as.numeric(xlim)
      if (length(effect_limits) != 2L ||
          any(!is.finite(effect_limits)) ||
          effect_limits[1L] >= effect_limits[2L]) {
        stop("`xlim` must contain two increasing finite values.",
             call. = FALSE)
      }
      if (x$config$scale == "log") {
        if (any(effect_limits <= 0)) {
          stop("Ratio/rate `xlim` values must be greater than 0.",
               call. = FALSE)
        }
        effect_limits <- log(effect_limits)
      }
      if (x$config$scale == "zcor") {
        if (any(abs(effect_limits) >= 1)) {
          stop("Correlation `xlim` must lie between -1 and 1.",
               call. = FALSE)
        }
        effect_limits <- atanh(effect_limits)
      }
    } else {
      effect_limits <- range(pretty(observed_limits, n = 5L))
    }
    args$alim <- effect_limits

    effect_span <- diff(effect_limits)
    if (!is.finite(effect_span) || effect_span <= 0) effect_span <- 1
    ilab_columns <- list()
    ilab_labels <- character()
    ilab_positions <- numeric()

    if (isTRUE(show_abcd)) {
      abcd <- x$analysis_data[, c(".a", ".b", ".c", ".d"), drop = FALSE]
      abcd[] <- lapply(abcd, .r4vn_meta_fmt, digits = 0)
      ilab_columns <- c(ilab_columns, unname(as.list(abcd)))
      ilab_labels <- c(ilab_labels, as.character(abcd_titles))
      ilab_positions <- c(
        ilab_positions,
        effect_limits[1L] - c(1.45, 1.15, 0.85, 0.55) * effect_span
      )
    }

    if (isTRUE(show_weights)) {
      weight_values <- .r4vn_meta_weights(x$primary_model)
      if (length(weight_values) != nrow(x$analysis_data)) {
        weight_values <- rep(NA_real_, nrow(x$analysis_data))
      }
      weight_text <- ifelse(
        is.finite(weight_values),
        paste0(.r4vn_meta_fmt(weight_values, 1), "%"), ""
      )
      ilab_columns <- c(ilab_columns, list(weight_text))
      ilab_labels <- c(ilab_labels, as.character(weight_title)[1L])
      # Keep Weight to the right of metafor's exact estimate annotation.
      ilab_positions <- c(
        ilab_positions, effect_limits[2L] + 1.45 * effect_span
      )
    }

    if (length(ilab_columns)) {
      args$ilab <- do.call(cbind, ilab_columns)
      if (supports_ilab_lab) args$ilab.lab <- ilab_labels
      args$ilab.xpos <- ilab_positions
      args$ilab.pos <- rep(2L, length(ilab_positions))
      estimate_position <- effect_limits[2L] + 0.75 * effect_span
      args$xlim <- c(
        effect_limits[1L] - if (isTRUE(show_abcd)) 2.10 * effect_span else
          0.75 * effect_span,
        effect_limits[2L] + if (isTRUE(show_weights)) 1.95 * effect_span else
          1.15 * effect_span
      )
      args$textpos <- c(args$xlim[1L], estimate_position)
    }
    if (supports_colci) {
      args$colci <- ci_color
    } else if (!identical(ci_color, point_color)) {
      # metafor versions before `colci` use `colout` for both study symbols
      # and their confidence intervals. Prefer the requested CI color and
      # avoid forwarding an unsupported graphical argument through `...`.
      args$colout <- ci_color
    }
    if (!is.null(point_size)) args$psize <- point_size

    if (!is.null(study_order)) {
      args$order <- study_order
    } else if (!is.null(x$by) && x$by %in% names(x$analysis_data)) {
      args$order <- x$analysis_data[[x$by]]
    }

    args <- .r4vn_meta_merge_plot_args(
      args, engine_args, protected = c("x", "atransf", "refline")
    )

    if (!is.null(ticks)) {
      tick_values <- as.numeric(ticks)
      if (x$config$scale == "log") {
        if (any(tick_values <= 0)) {
          stop("Ratio/rate tick values must be greater than 0.",
               call. = FALSE)
        }
        tick_values <- log(tick_values)
      }
      if (x$config$scale == "zcor") {
        if (any(abs(tick_values) >= 1)) {
          stop("Correlation ticks must lie between -1 and 1.",
               call. = FALSE)
        }
        tick_values <- atanh(tick_values)
      }
      args$at <- tick_values
    }

    forest_info <- do.call(metafor::forest, args)
    if (length(ilab_columns) && !supports_ilab_lab &&
        !is.null(args$ilab.xpos)) {
      header_y <- x$primary_model$k + 2
      if (is.list(forest_info) && !is.null(forest_info$rows) &&
          any(is.finite(forest_info$rows))) {
        header_y <- max(forest_info$rows, na.rm = TRUE) + 2
      }
      graphics::text(
        x = args$ilab.xpos, y = rep(header_y, length(args$ilab.xpos)),
        labels = ilab_labels,
        pos = 2, font = 2, cex = text_size
      )
    }
    graphics::abline(
      v = .r4vn_meta_analysis_ref(x, ref), col = ref_color,
      lty = ref_type, lwd = ref_width
    )
    graphics::title(
      main = title, line = 2.2, cex.main = title_size,
      col.main = title_color
    )
    return(finish())
  }

  if (type == "funnel") {
    funnel_ref <- if (is.null(ref)) {
      as.numeric(x$primary_model$beta[1L])
    } else .r4vn_meta_analysis_ref(x, ref)
    args <- list(
      x = x$primary_model,
      yaxis = yaxis,
      atransf = back,
      xlab = xlab,
      ylab = ylab,
      digits = digits,
      pch = if (is.null(point_shape)) 19 else point_shape,
      col = point_color,
      bg = point_bg,
      refline = funnel_ref,
      lty = ref_type,
      label = funnel_label,
      legend = funnel_legend
    )
    if (!is.null(point_size)) args$cex <- point_size
    if (isTRUE(contour)) {
      if (length(contour_colors) != length(contour_levels)) {
        stop("`contour_colors` must match `contour_levels` in length.",
             call. = FALSE)
      }
      args$level <- contour_levels
      args$shade <- contour_colors
    }
    args <- .r4vn_meta_merge_plot_args(
      args, engine_args, protected = c("x", "atransf", "refline")
    )
    do.call(metafor::funnel, args)
    graphics::abline(
      v = funnel_ref, col = ref_color, lty = ref_type, lwd = ref_width
    )
    graphics::title(main = title, cex.main = title_size,
                    col.main = title_color)
    return(finish())
  }

  if (type == "trimfill") {
    if (is.null(x$bias) || is.null(x$bias$trimfill)) {
      stop("Trim-and-fill results are unavailable. Request `bias=TRUE` and ",
           "include `trimfill` in `bias_methods`.", call. = FALSE)
    }
    funnel_ref <- if (is.null(ref)) {
      as.numeric(x$bias$trimfill$beta[1L])
    } else .r4vn_meta_analysis_ref(x, ref)
    args <- list(
      x = x$bias$trimfill,
      yaxis = yaxis,
      atransf = back,
      xlab = xlab,
      ylab = ylab,
      digits = digits,
      pch = if (is.null(point_shape)) 19 else point_shape,
      pch.fill = if (is.null(point_shape)) 21 else point_shape,
      col = point_color,
      bg = point_bg,
      refline = funnel_ref,
      lty = ref_type,
      label = funnel_label,
      legend = funnel_legend
    )
    if (!is.null(point_size)) args$cex <- point_size
    if (isTRUE(contour)) {
      if (length(contour_colors) != length(contour_levels)) {
        stop("`contour_colors` must match `contour_levels` in length.",
             call. = FALSE)
      }
      args$level <- contour_levels
      args$shade <- contour_colors
    }
    args <- .r4vn_meta_merge_plot_args(
      args, engine_args, protected = c("x", "atransf", "refline")
    )
    do.call(metafor::funnel, args)
    graphics::abline(
      v = funnel_ref, col = ref_color, lty = ref_type, lwd = ref_width
    )
    graphics::title(main = title, cex.main = title_size,
                    col.main = title_color)
    return(finish())
  }

  if (type == "baujat") {
    args <- .r4vn_meta_merge_plot_args(
      list(
        x = x$primary_model, main = title,
        pch = if (is.null(point_shape)) 21 else point_shape,
        col = point_color, bg = point_bg
      ),
      engine_args, protected = "x"
    )
    do.call(metafor::baujat, args)
    return(finish())
  }

  if (type == "radial") {
    args <- .r4vn_meta_merge_plot_args(
      list(
        x = x$primary_model, main = title,
        pch = if (is.null(point_shape)) 21 else point_shape,
        col = point_color, bg = point_bg
      ),
      engine_args, protected = "x"
    )
    do.call(metafor::radial, args)
    return(finish())
  }

  if (type == "labbe") {
    required <- c(".event1", ".n1", ".event0", ".n0")
    if (!all(required %in% names(x$analysis_data))) {
      stop(
        "L'Abbe plots require raw binary event and sample-size input.",
        call. = FALSE
      )
    }
    d <- x$analysis_data
    pc <- d$.event0 / d$.n0
    pt <- d$.event1 / d$.n1
    w <- .r4vn_meta_weights(x$primary_model)
    if (!length(w) || all(!is.finite(w))) w <- rep(1, nrow(d))
    w[!is.finite(w)] <- min(w[is.finite(w)], na.rm = TRUE)
    cex <- 0.7 + 1.8 * sqrt(w / max(w, na.rm = TRUE))
    lim <- range(c(pc, pt, 0, 1), finite = TRUE)
    graphics::plot(
      pc, pt,
      xlim = lim, ylim = lim,
      xlab = "Control-group risk",
      ylab = "Treatment-group risk",
      main = title,
      pch = if (is.null(point_shape)) 21 else point_shape,
      bg = point_bg,
      col = point_color,
      cex = if (is.null(point_size)) cex else point_size,
      cex.main = title_size, col.main = title_color,
      cex.lab = axis_size, cex.axis = axis_size
    )
    graphics::abline(
      a = 0, b = 1, col = ref_color, lty = ref_type, lwd = ref_width
    )
    graphics::text(
      pc, pt, labels = d$.study, pos = 3,
      cex = max(0.45, text_size * 0.7), col = point_color
    )
    return(finish())
  }

  if (type == "influence") {
    if (is.null(x$influence) || is.null(x$influence$cook)) {
      stop(
        "Influence diagnostics are unavailable. Re-run with ",
        "`influence=TRUE` or `full=TRUE`.",
        call. = FALSE
      )
    }
    v <- x$influence$cook
    graphics::barplot(
      v,
      names.arg = x$analysis_data$.study,
      las = 2,
      ylab = "Cook's distance",
      main = title,
      col = shade_color,
      border = point_color,
      cex.names = max(0.5, text_size * 0.82),
      cex.main = title_size, col.main = title_color,
      cex.lab = axis_size, cex.axis = axis_size
    )
    graphics::abline(
      h = 4 / length(v), col = ref_color,
      lty = ref_type, lwd = ref_width
    )
    return(finish())
  }

  if (type == "leaveout") {
    z <- x$leaveout
    if (is.null(z) || is.null(z$table)) {
      stop(
        "Leave-one-out results are unavailable. Re-run with ",
        "`leaveout=TRUE` or `full=TRUE`.",
        call. = FALSE
      )
    }
    y <- rev(seq_along(z$estimate))
    xr <- range(c(z$lower, z$upper), finite = TRUE)
    if (!length(xr) || any(!is.finite(xr))) xr <- c(0, 1)
    graphics::plot(
      NA,
      xlim = if (is.null(xlim)) xr else xlim,
      ylim = c(0.5, length(y) + 0.5),
      yaxt = "n",
      xlab = xlab,
      ylab = "",
      main = title,
      cex.main = title_size, col.main = title_color,
      cex.lab = axis_size, cex.axis = axis_size
    )
    graphics::axis(
      2, at = y, labels = rev(z$table[[1L]]),
      las = 2, cex.axis = 0.66
    )
    graphics::segments(
      z$lower, y, z$upper, y, col = ci_color,
      lty = line_type, lwd = line_width
    )
    graphics::points(
      z$estimate, y, pch = if (is.null(point_shape)) 15 else point_shape,
      col = point_color, bg = point_bg,
      cex = if (is.null(point_size)) 1 else point_size
    )
    graphics::abline(
      v = x$overall$estimate, col = ref_color,
      lty = ref_type, lwd = ref_width
    )
    return(finish())
  }

  if (type == "cumulative") {
    z <- x$cumulative
    if (is.null(z) || is.null(z$plot_data)) {
      stop(
        "Cumulative results are unavailable. Supply `cumulative=` in ",
        "`tabmeta()`.",
        call. = FALSE
      )
    }
    p <- z$plot_data
    y <- rev(seq_len(nrow(p)))
    xr <- range(c(p$Lower, p$Upper), finite = TRUE)
    if (!length(xr) || any(!is.finite(xr))) xr <- c(0, 1)
    graphics::plot(
      NA,
      xlim = if (is.null(xlim)) xr else xlim,
      ylim = c(0.5, nrow(p) + 0.5),
      yaxt = "n",
      xlab = xlab,
      ylab = "",
      main = title,
      cex.main = title_size, col.main = title_color,
      cex.lab = axis_size, cex.axis = axis_size
    )
    graphics::axis(
      2, at = y, labels = rev(p[[1L]]),
      las = 2, cex.axis = 0.72
    )
    graphics::segments(
      p$Lower, y, p$Upper, y, col = ci_color,
      lty = line_type, lwd = line_width
    )
    graphics::points(
      p$Estimate, y, pch = if (is.null(point_shape)) 15 else point_shape,
      col = point_color, bg = point_bg,
      cex = if (is.null(point_size)) 1 else point_size
    )
    graphics::abline(
      v = .r4vn_meta_natural_ref(x, ref), col = ref_color,
      lty = ref_type, lwd = ref_width
    )
    return(finish())
  }

  if (type == "bubble") {
    if (is.null(x$regression) || is.null(x$regression$model)) {
      stop(
        "Meta-regression results are unavailable. Supply `reg=` in `tabmeta()`.",
        call. = FALSE
      )
    }
    if (is.null(moderator)) {
      if (length(x$reg) == 1L) {
        moderator <- x$reg[1L]
      } else {
        stop("Specify one moderator with `moderator=`.", call. = FALSE)
      }
    }
    moderator <- as.character(moderator)[1L]
    if (!moderator %in% x$reg) {
      stop(
        "Moderator `", moderator, "` was not included in `reg`.",
        call. = FALSE
      )
    }
    if (!is.numeric(x$analysis_data[[moderator]])) {
      stop(
        "Bubble plots currently require a numeric moderator.",
        call. = FALSE
      )
    }

    # With a single numeric moderator the first moderator coefficient is mod=1.
    args <- list(
      x = x$regression$model, mod = 1, atransf = back,
      xlab = moderator, ylab = xlab,
      refline = .r4vn_meta_analysis_ref(x, ref), main = title,
      ci = isTRUE(bubble_ci), pi = isTRUE(bubble_pi),
      shade = bubble_shade, grid = grid,
      pch = if (is.null(point_shape)) 21 else point_shape,
      col = point_color, bg = point_bg,
      lcol = c(summary_color, ci_color, ci_color),
      lwd = line_width, lty = line_type
    )
    if (!is.null(point_size)) args$psize <- point_size
    args <- .r4vn_meta_merge_plot_args(
      args, engine_args, protected = c("x", "atransf", "mod")
    )
    do.call(metafor::regplot, args)
    return(finish())
  }

  invisible(x)
}

#' Plot an R4VN meta-analysis
#'
#' Draw forest, funnel, trim-and-fill, influence, cumulative, diagnostic, and
#' meta-regression figures. Every figure can be drawn in the R/RStudio Plot
#' pane or saved independently in a publication format.
#'
#' @param x An object created by `tabmeta()`.
#' @param type Plot type: `"forest"`, `"funnel"`, `"trimfill"`, `"baujat"`,
#'   `"influence"`, `"leaveout"`, `"cumulative"`, `"radial"`, `"labbe"`,
#'   or `"bubble"`.
#' @param moderator Moderator for a bubble/meta-regression plot.
#' @param contour Logical; create a contour-enhanced funnel plot.
#' @param title,subtitle,caption Main title, subtitle, and figure caption.
#' @param xlab Plot x-axis label.
#' @param ref Reference line on the natural effect scale.
#' @param xlim,ticks Optional effect-scale x limits and tick positions. A study
#'   estimate or confidence interval outside `xlim` is indicated by an arrow;
#'   its exact, unclipped value remains in the numerical estimate column.
#' @param color Backward-compatible overall plotting color.
#' @param font_family,text_size,axis_size,title_size,title_color Font family,
#'   study-label size, axis size, title size, and title color.
#' @param subtitle_size,subtitle_color,caption_size,caption_color Subtitle and
#'   caption sizes and colors.
#' @param margins Optional four-value base-graphics margin vector.
#' @param background Figure background color.
#' @param point_color,point_bg,ci_color,summary_color,summary_border Colors for
#'   study points, point fill, confidence intervals, pooled diamond, and its
#'   border.
#' @param point_shape,point_size Point symbol and optional fixed point size.
#' @param line_type,line_width Confidence-interval line type and width.
#' @param ref_color,ref_type,ref_width Reference-line color, type, and width.
#' @param show_weights Show study weights in a forest plot.
#' @param show_prediction Show the prediction interval in the forest plot.
#'   Default `FALSE`; the prediction interval remains available in the
#'   numerical results when requested in `tabmeta()`. Set `TRUE` to draw it
#'   for a random-effects model.
#' @param weight_title,estimate_title Headings for the separate weight and
#'   numerical estimate columns. `estimate_title=NULL` uses the effect measure
#'   and confidence level. Weight is placed to the right of the numerical
#'   estimate column.
#' @param show_abcd Show the four original binary cells in separate forest
#'   columns. This is available when `tabmeta()` was called with `a`, `b`, `c`,
#'   and `d`.
#' @param abcd_titles Four headings used when `show_abcd=TRUE`.
#' @param prediction_style Prediction display: `"line"`, `"polygon"`,
#'   `"bar"`, `"shade"`, or `"dist"`.
#' @param row_shade,shade_color Forest-row shading style and color.
#' @param digits Number of displayed decimals.
#' @param study_order Optional forest ordering vector or metafor order keyword.
#' @param header Forest headings; `NULL` uses R4VN headings.
#' @param annotate Show effect and confidence-interval annotations.
#' @param yaxis,ylab Funnel-plot y-axis definition and label.
#' @param contour_levels,contour_colors Funnel contour levels and colors.
#' @param funnel_label Label funnel points (`FALSE`, `TRUE`, `"all"`, `"out"`,
#'   or a number of extreme points).
#' @param funnel_legend Funnel legend control, including a position such as
#'   `"topright"`.
#' @param bubble_ci,bubble_pi,bubble_shade,grid Meta-regression confidence and
#'   prediction bands, their shading colors, and grid display.
#' @param engine_args Named list passed to the underlying metafor plotting
#'   method. This provides access to advanced engine-specific controls.
#' @param file Optional output file. Supported extensions are PNG, JPEG, TIFF,
#'   PDF, and SVG.
#' @param width,height Device width and height. Raster units are pixels.
#' @param res Raster resolution.
#' @param ... Additional named arguments merged into `engine_args`.
#'
#' @return The meta-analysis object invisibly.
#' @export
#' @examples
#' \donttest{
#' if (requireNamespace("metafor", quietly = TRUE)) {
#'   dat <- read.csv(system.file("extdata", "meta_example.csv", package = "R4VN"))
#'   m <- tabmeta(
#'     dat, study, effect = OR, lower = LCI, upper = UCI, or = TRUE,
#'     show = FALSE
#'   )
#'
#'   # Journal-style forest plot in the Plot pane.
#'   plot(
#'     m, type = "forest", font_family = "sans",
#'     subtitle = "Random-effects model with 95% confidence intervals",
#'     caption = "Square size reflects study weight; diamond is pooled effect.",
#'     margins = c(5.5, 4.2, 5.0, 2.0),
#'     point_color = "#1F4E79", ci_color = "#5B9BD5",
#'     summary_color = "#C00000", summary_border = "#7F0000",
#'     row_shade = "zebra", shade_color = "#F5F7FA",
#'     show_weights = TRUE, weight_title = "Weight",
#'     estimate_title = "OR (95% CI)", show_prediction = FALSE,
#'     xlim = c(0.2, 2), ticks = c(0.25, 0.5, 1, 1.5, 2)
#'   )
#'
#'   # Contour-enhanced funnel and trim-and-fill plots.
#'   plot(
#'     m, type = "funnel", contour = TRUE,
#'     point_shape = 21, point_color = "#1F4E79", point_bg = "#D9EAF7",
#'     contour_levels = c(90, 95, 99),
#'     contour_colors = c("#FFF2CC", "#FCE4D6", "#E2F0D9"),
#'     funnel_label = "out", funnel_legend = "topright"
#'   )
#'   if (!is.null(m$models$trimfill)) {
#'     plot(m, type = "trimfill", contour = TRUE)
#'   }
#'
#'   # Save a 300-dpi TIFF independently.
#'   forest_file <- tempfile(fileext = ".tiff")
#'   plot(
#'     m, type = "forest", file = forest_file,
#'     width = 2400, height = 1800, res = 300,
#'     font_family = "sans", point_color = "#1F4E79",
#'     ci_color = "#5B9BD5", summary_color = "#C00000"
#'   )
#'   unlink(forest_file)
#'
#'   # Advanced metafor controls.
#'   plot(m, type = "forest",
#'        engine_args = list(efac = c(1, 1.2), plim = c(0.6, 1.8)))
#' }
#' }
plot.r4vn_meta <- function(
    x, type = "forest", moderator = NULL, contour = FALSE,
    title = NULL, subtitle = NULL, caption = NULL, xlab = NULL,
    ref = NULL, xlim = NULL, ticks = NULL,
    color = NULL, font_family = "sans", text_size = 0.82,
    axis_size = 0.9, title_size = 1.05, title_color = "black",
    subtitle_size = 0.9, subtitle_color = "gray30",
    caption_size = 0.75, caption_color = "gray40",
    margins = NULL, background = "white",
    point_color = NULL, point_bg = "white",
    ci_color = NULL, summary_color = NULL, summary_border = NULL,
    point_shape = NULL, point_size = NULL, line_type = 1,
    line_width = 1, ref_color = "gray40", ref_type = 2,
    ref_width = 1, show_weights = TRUE, show_prediction = FALSE,
    weight_title = "Weight", estimate_title = NULL,
    show_abcd = FALSE, abcd_titles = c("a", "b", "c", "d"),
    prediction_style = "bar", row_shade = "zebra",
    shade_color = "gray95", digits = NULL, study_order = NULL,
    header = NULL, annotate = TRUE, yaxis = "sei", ylab = NULL,
    contour_levels = c(90, 95, 99),
    contour_colors = c("#FEE2E2", "#FEF3C7", "#E5E7EB"),
    funnel_label = FALSE, funnel_legend = FALSE,
    bubble_ci = TRUE, bubble_pi = FALSE,
    bubble_shade = c("#DBEAFE", "#E5E7EB"), grid = FALSE,
    engine_args = list(), file = NULL, width = 1800, height = 1400,
    res = 180, ...) {

  if (!inherits(x, "r4vn_meta")) {
    stop("`x` must be created by `tabmeta()`.", call. = FALSE)
  }

  dots <- list(...)
  if (length(dots) && is.null(names(dots))) {
    stop("Additional plot arguments must be named.", call. = FALSE)
  }
  engine_args <- .r4vn_meta_merge_plot_args(engine_args, dots)
  draw_args <- list(
    x = x, type = type, moderator = moderator, contour = contour,
    title = title, subtitle = subtitle, caption = caption,
    xlab = xlab, ref = ref, xlim = xlim, ticks = ticks,
    color = color, font_family = font_family, text_size = text_size,
    axis_size = axis_size, title_size = title_size,
    title_color = title_color, subtitle_size = subtitle_size,
    subtitle_color = subtitle_color, caption_size = caption_size,
    caption_color = caption_color, margins = margins,
    background = background,
    point_color = point_color, point_bg = point_bg,
    ci_color = ci_color, summary_color = summary_color,
    summary_border = summary_border, point_shape = point_shape,
    point_size = point_size, line_type = line_type,
    line_width = line_width, ref_color = ref_color,
    ref_type = ref_type, ref_width = ref_width,
    show_weights = show_weights, show_prediction = show_prediction,
    weight_title = weight_title, estimate_title = estimate_title,
    show_abcd = show_abcd, abcd_titles = abcd_titles,
    prediction_style = prediction_style, row_shade = row_shade,
    shade_color = shade_color, digits = digits,
    study_order = study_order, header = header, annotate = annotate,
    yaxis = yaxis, ylab = ylab, contour_levels = contour_levels,
    contour_colors = contour_colors, funnel_label = funnel_label,
    funnel_legend = funnel_legend, bubble_ci = bubble_ci,
    bubble_pi = bubble_pi, bubble_shade = bubble_shade, grid = grid,
    engine_args = engine_args
  )

  if (is.null(file)) return(do.call(.r4vn_meta_draw, draw_args))

  ext <- tolower(tools::file_ext(file))
  allowed <- c("png", "jpg", "jpeg", "tif", "tiff", "pdf", "svg")
  if (!ext %in% allowed) {
    stop(
      "Unsupported plot file extension. Use: ",
      paste(allowed, collapse = ", "), ".",
      call. = FALSE
    )
  }

  if (ext == "png") {
    grDevices::png(file, width = width, height = height, res = res,
                   bg = background)
  } else if (ext %in% c("jpg", "jpeg")) {
    grDevices::jpeg(file, width = width, height = height, res = res,
                    bg = background)
  } else if (ext %in% c("tif", "tiff")) {
    grDevices::tiff(file, width = width, height = height, res = res,
                    bg = background)
  } else if (ext == "pdf") {
    grDevices::pdf(file, width = width / res, height = height / res,
                   bg = background)
  } else if (ext == "svg") {
    grDevices::svg(file, width = width / res, height = height / res,
                   bg = background)
  }
  on.exit(grDevices::dev.off(), add = TRUE)

  do.call(.r4vn_meta_draw, draw_args)
  invisible(x)
}