R/er-plot-layer.R

In erplots: Model-Agnostic Exposure-Response Plots

# layer_model ------------------------------------------------------------------

.layer_model <- function(object, model, stratify, conf_level, style,
                          predict_args = list(), dots = list()) {

  layer_model <- list()
  config <- list()

  # the fitted model, supplied by the caller (erplots never fits models
  # itself -- see `er_model_interface`)
  config$model <- model

  # confidence level
  config$conf_level <- conf_level

  # stashed so `.refresh_model_predictions()` (called from
  # `er_plot_build()`) can recompute `config$predictions` against
  # whatever `object$exposure$limits`/`object$strata` look like at build
  # time, rather than being stuck with the snapshot taken here -- see
  # that function's own comment for why this matters (issue #14: a later
  # `er_plot_theme(xlim = ...)` call must still be reflected in the
  # model curve/ribbon, not just the coordinate system)
  config$predict_args <- predict_args

  # model predictions, via the `er_predict()` generic. `predict_args`
  # (from `er_plot_add_model()`'s own argument of the same name, kept
  # separate from `dots`/`config$dots` below) is spliced into the
  # `er_predict()` call itself, for model-specific arguments beyond the
  # fixed `model`/`newdata`/`conf_level` contract -- see
  # `?er_model_interface` and `?er_plot_add_model`'s "Details". This is
  # an eager, add-time computation purely so a bad model/`predict_args`
  # combination fails immediately at the `er_plot_add_model()` call site
  # rather than silently, much later, inside `plot()`/`print()`; the
  # value computed here is unconditionally replaced by
  # `.refresh_model_predictions()` at build time (`er_plot_build()`), so
  # it never actually reaches a builder.
  config$predictions <- .get_model_predictions(
    config$model,
    config$conf_level,
    object$exposure,
    object$strata,
    stratify,
    object$data,
    predict_args = predict_args
  )

  # `style` is the escape hatch documented in `?er_style`: any function
  # matching the standard `er_style_*()` signature can be plugged in
  # without touching package internals. `er_plot_add_model()` has already
  # resolved a default when the caller didn't supply one, so this is
  # always a function here.
  config$style <- style

  # extra named arguments from `er_plot_add_model()`'s `...` -- see
  # `?er_style`'s "Passing extra arguments to a builder" section
  config$dots <- dots

  # store and return
  layer_model$stratify <- stratify
  layer_model$config <- config

  return(layer_model)
}


# layer_summary -----------------------------------------------------------------

.layer_summary <- function(object, model, stratify, style, conf_level = 0.95,
                            summary_args = list(), dots = list()) {

  layer_summary <- list()
  config <- list()

  # the fitted model, supplied by the caller. Unlike `.layer_model()`,
  # this is optional -- a summary builder doesn't have to be a *model*
  # summary (e.g. `er_style_summary_n()` is purely descriptive of the
  # data) -- so `model` may be `NULL` here. Use `[` (not `$`) so a `NULL`
  # model is retained as a named element rather than dropped from `config`.
  config["model"] <- list(model)

  # model summary, via the `er_summary()` generic, when a model was
  # supplied. Computed unconditionally (regardless of `stratify`) --
  # whether a given value makes sense to show when the layer is
  # stratified is a decision for the builder itself (see
  # `er_style_summary_pvalue()`), not this generic config-building step.
  # `config$summary` holds the full, raw `er_summary()` return value (see
  # `?er_model_interface` for its contract -- `p_value`/`coefficients`/
  # `glance`), so builders needing more than a bare p-value (e.g.
  # `er_style_summary_coefficients()`) can read it directly.
  # `config$p_value` is kept as a separate, extracted field alongside it
  # purely so `er_style_summary_pvalue()`'s existing, simpler read
  # (`config$p_value`) doesn't need to change.
  #
  # `conf_level` is now always forwarded (previously `er_summary(model)`
  # was called with no arguments at all, despite `conf_level` being part
  # of the documented `?er_model_interface` contract); `summary_args`
  # (from `er_plot_add_summary()`'s own argument of the same name, kept
  # separate from `dots`/`config$dots` below) is spliced in alongside it
  # for any further model-specific argument.
  config["summary"] <- list(NULL) # use `[` (not `$`) so the NULL is retained as a named element
  config["p_value"] <- list(NULL)
  if (!is.null(model)) {
    model_summary <- rlang::exec(er_summary, model = model, conf_level = conf_level, !!!summary_args)
    config$summary <- model_summary
    config$p_value <- model_summary$p_value
  }

  # visual distance from corners (used for placement of the summary
  # annotation), based on the raw observed data rather than any model's
  # fitted curve -- this works whether or not a model was supplied, and
  # avoids overlapping the raw points a data/quantile layer might also be
  # showing. See `.compute_corner_distance()` -- shared with
  # `.layer_quantile()`, which uses the same calculation to keep a
  # quantile boundary-vline label out of the summary annotation's corner
  # (see `?er_style_quantile`).
  config$corner_distance <- .compute_corner_distance(object$data, object$exposure, object$response)

  # `style` is the escape hatch documented in `?er_style`; `er_plot_add_summary()`
  # has already resolved a default when the caller didn't supply one.
  config$style <- style

  # extra named arguments from `er_plot_add_summary()`'s `...` -- see
  # `?er_style`'s "Passing extra arguments to a builder" section
  config$dots <- dots

  # store and return
  layer_summary$stratify <- stratify
  layer_summary$config <- config

  return(layer_summary)
}


# layer_quantile ---------------------------------------------------------------

.layer_quantile <- function(object, stratify, n_bins, conf_level, style, dots = list(),
                             ties = "upward", quantile_type = 7, labeller = NULL) {

  layer_quantile <- list()
  config <- list()

  config$n_quantiles <- n_bins
  config$conf_level <- conf_level

  binned <- object$data |>
    dplyr::mutate(
      response = .data[[object$response$name]],
      exposure_bins = cut_exposure_quantile(
        x = .data[[object$exposure$name]],
        n_bins = config$n_quantiles,
        ties = ties,
        quantile_type = quantile_type,
        labeller = labeller
      ),
      strata = .get_strata_values(.data, object$strata$name)
    )

  # quantile cutpoints (excluding placebo), for builders that draw
  # bin-boundary separators (e.g. `er_style_quantile_errorbar_vlines()`)
  # -- see `cut_exposure_quantile()`'s `"breaks"` attribute
  config$breaks <- attr(binned$exposure_bins, "breaks")
  # tie-break rule actually used, read back off the binned column's own
  # attribute (rather than re-storing the `ties` argument directly) so
  # `.check_exposure_binning_consistency()` compares what was actually
  # applied -- these two are always identical today, but keeping the
  # source-of-truth on the computed column avoids the two ever silently
  # drifting apart if that changes
  config$ties <- attr(binned$exposure_bins, "ties")

  # visual distance from corners, identical to `.layer_summary()`'s own
  # computation -- lets `er_style_quantile_errorbar_vlines()`/
  # `er_style_quantile_pointrange_vlines()`'s optional vline labels pick
  # the vertical half *opposite* wherever a summary annotation would
  # render (if one is present), without either layer needing to know
  # about the other. See `.compute_corner_distance()`.
  config$corner_distance <- .compute_corner_distance(object$data, object$exposure, object$response)

  # binary response: response *rate* per bin, via a Clopper-Pearson CI.
  # count response, when explicitly declared (`response_type = "count"`):
  # bin *mean*, via an exact Poisson interval. continuous (and, when not
  # explicitly declared "count", an approximation for count) response: bin
  # *mean*, via a t-interval. Label placement (y_lwr_lbl/y_upr_lbl/y_lbl) is
  # generalised across all branches below rather than duplicated, using
  # the response's own scale (`object$response$limits`) in place of the
  # binary-only [0, 1] assumption.
  if (object$response$type == "binary") {
    config$summary <- binned |>
      dplyr::summarise(
        n1 = sum(response == 1, na.rm = TRUE),
        n0 = sum(response == 0, na.rm = TRUE),
        x_mid = mean(.data[[object$exposure$name]], na.rm = TRUE),
        y_mid = n1 / (n0 + n1),
        y_mid_lbl = object$theme$format_percent(n1 / (n0 + n1)),
        ci_lower = ci_clopper_pearson(n1, n0 + n1, config$conf_level)["lower"],
        ci_upper = ci_clopper_pearson(n1, n0 + n1, config$conf_level)["upper"],
        .by = c("exposure_bins", "strata")
      )
  } else if (object$response$type == "count") {
    config$summary <- binned |>
      dplyr::summarise(
        n_units = sum(!is.na(response)),
        x_mid = mean(.data[[object$exposure$name]], na.rm = TRUE),
        y_mid = mean(response, na.rm = TRUE),
        y_mid_lbl = object$theme$format_number(mean(response, na.rm = TRUE)),
        ci_lower = ci_poisson(sum(response, na.rm = TRUE), n_units, config$conf_level)["lower"],
        ci_upper = ci_poisson(sum(response, na.rm = TRUE), n_units, config$conf_level)["upper"],
        .by = c("exposure_bins", "strata")
      ) |>
      dplyr::select(-n_units)
  } else {
    config$summary <- binned |>
      dplyr::summarise(
        x_mid = mean(.data[[object$exposure$name]], na.rm = TRUE),
        y_mid = mean(response, na.rm = TRUE),
        y_mid_lbl = object$theme$format_number(mean(response, na.rm = TRUE)),
        ci_lower = ci_t(response, config$conf_level)["lower"],
        ci_upper = ci_t(response, config$conf_level)["upper"],
        .by = c("exposure_bins", "strata")
      )
  }

  response_lo <- object$response$limits[1]
  response_hi <- object$response$limits[2]
  margin <- 0.05 * (response_hi - response_lo)

  config$summary <- config$summary |>
    dplyr::mutate(
      y_lwr_lbl = ci_lower - margin,
      y_upr_lbl = ci_upper + margin,
      y_lbl = dplyr::if_else(
        (y_lwr_lbl - response_lo) > (response_hi - y_upr_lbl),
        y_lwr_lbl,
        y_upr_lbl
      )
    )

  # see `?er_style` for the `style` escape hatch; `er_plot_add_quantiles()`
  # has already resolved a default when the caller didn't supply one
  config$style <- style

  # extra named arguments from `er_plot_add_quantiles()`'s `...` -- see
  # `?er_style`'s "Passing extra arguments to a builder" section
  config$dots <- dots

  # store and return
  layer_quantile$stratify <- stratify
  layer_quantile$config <- config

  return(layer_quantile)
}


# layer_data -------------------------------------------------------------------

.layer_data <- function(object, stratify, panel, style, dots = list()) {

  layer_data <- list()

  config <- list()
  config$layout <- "panel"
  config$panel <- panel
  # `NULL` unless the caller passes `seed = <value>` through
  # `er_plot_add_data()`'s own `...` -- no seed management happens
  # otherwise, so the jitter draws from the ambient RNG stream like any
  # other jittered geom and differs across repeated `plot()` calls on
  # the same object. This used to be a hard-coded literal (`1234L`),
  # which CRAN policy forbids (a package must not fix a seed without
  # the user's consent); rather than swap it for an auto-generated
  # random seed (which would silently reintroduce a form of
  # always-on seed management the user never asked for), seeding is
  # opt-in only, matching the CRAN Cookbook's recommended pattern.
  # Assigning via `[<-`/`list()` (rather than `$<-`) keeps the `seed`
  # name present in `config` even when the value is `NULL` -- `$<-`
  # would silently drop the element instead of storing `NULL`.
  config["seed"] <- list(dots$seed)
  # `er_plot_add_data()` has already resolved `style` (and confirmed
  # its layout is "panel") before calling here -- see `?er_style` for
  # the `style`/`er_style_tag()` escape hatch
  config$style <- style
  # extra named arguments from `er_plot_add_data()`'s `...` -- see
  # `?er_style`'s "Passing extra arguments to a builder" section
  config$dots <- dots

  # `panels` is a named list of panels to build, keyed by panel name, in
  # build order; `panel_position` records where each named panel sits
  # relative to the base plot ("above"/"below"), which the composition
  # helpers (R/er-plot-compose.R) use instead of hardcoding "upper"/"lower".
  # `color_role` tags what the layer's `colour` aesthetic means --
  # "strata" (the usual case, dispatched to via the shared strata legend)
  # or "response" (the continuous/count variant's colour-encoded response
  # value, which needs its own label/legend and isn't deduplicated across
  # stratum panels) -- consumed by `.polish_labels()`/`.polish_legends()`
  # in R/er-plot-compose.R.
  if (object$response$type == "binary") {
    config$color_role <- "strata"

    panels <- character(0)
    if (panel %in% c("upper", "both")) panels <- c(panels, "upper")
    if (panel %in% c("lower", "both")) panels <- c(panels, "lower")
    config$panels <- panels
    config$panel_position <- c(upper = "above", lower = "below")[panels]

  } else {
    # continuous/count response: a single panel, points coloured
    # continuously by the response value, in place of the binary
    # upper/lower partition -- `er_plot_add_data()` guards `panel` to
    # "both" for this response type, since there's no upper/lower
    # partition to select from. When stratified, the colour channel is
    # already spoken for by the response, so stratification becomes one
    # panel per stratum level instead (all placed "below" the base plot).
    config$color_role <- "response"

    if (stratify) {
      panels <- as.character(object$strata$limits)
    } else {
      panels <- "data"
    }
    config$panels <- panels
    config$panel_position <- stats::setNames(rep("below", length(panels)), panels)
  }

  layer_data$stratify <- stratify
  layer_data$config <- config

  return(layer_data)
}


# layer_overlay ------------------------------------------------------------------

.layer_overlay <- function(object, stratify, style, dots = list()) {

  layer_overlay <- list()

  config <- list()
  # see the matching comment in `.layer_data()` above: seeding is
  # opt-in only, via a user-supplied `seed`; `NULL` otherwise.
  config["seed"] <- list(dots$seed)

  # unlike `.layer_data()`, there's a single builder regardless of
  # response type -- `er_style_data_overlay()` only needs to know the
  # response type to decide how much vertical jitter to apply (binary
  # responses get a small nudge so 0/1 points don't overplot into two
  # solid lines; continuous/count responses get none). `er_plot_add_data()`
  # has already resolved `style` (and confirmed its layout is "overlay")
  # before calling here -- see `?er_style` for the `style`/`er_style_tag()`
  # escape hatch.
  config$response_type <- object$response$type
  config$style <- style
  # extra named arguments from `er_plot_add_data()`'s `...` -- see
  # `?er_style`'s "Passing extra arguments to a builder" section
  config$dots <- dots

  layer_overlay$stratify <- stratify
  layer_overlay$config <- config

  return(layer_overlay)
}


# layer_group ------------------------------------------------------------------

.layer_group <- function(object, group_cols, stratify, n_bins, style, dots = list(),
                          ties = "upward", quantile_type = 7, labeller = NULL) {

  # grouping by the plot's own stratification variable while also
  # keeping strata (`stratify == TRUE`) bakes the same column name into
  # `config$groupings` twice (`c(g, object$strata$name)`), which makes
  # the `dplyr::left_join()` below fail with "Join columns in `x` must
  # be unique" -- catch it here with a message that names the actual
  # problem, rather than letting the join error surface uninformatively
  if (stratify && !is.null(object$strata$name) && any(group_cols %in% object$strata$name)) {
    offenders <- group_cols[group_cols %in% object$strata$name]
    rlang::abort(c(
      "`group_by` cannot include the plot's own stratification variable when `keep_strata = TRUE`.",
      "x" = paste0(
        "`", paste(offenders, collapse = "`, `"), "` is already used to stratify this plot ",
        "(see `stratify_by` in `er_plot()`)."
      ),
      "i" = "Set `keep_strata = FALSE` in `er_plot_add_groups()`, or group by a different variable."
    ))
  }

  layer_group <- list()
  layer_group$stratify <- stratify
  layer_group$config <- list()

  for(g in group_cols) {

    config <- list()
    # see `?er_style` for the `style` escape hatch; `er_plot_add_groups()`
    # has already resolved a default when the caller didn't supply one
    config$style <- style
    # extra named arguments from `er_plot_add_groups()`'s `...`, shared
    # identically across every grouping variable added by this call --
    # see `?er_style`'s "Passing extra arguments to a builder" section
    config$dots <- dots

    # data
    dat <- object$data

    # create factor from continuous grouping variables
    if (is.numeric(dat[[g]])) {
      new_g <- paste0(".", g, "_quantile")
      new_g_sym <- dplyr::sym(new_g)
      if (g == object$exposure$name) {
        dat <- dat |>
          dplyr::mutate(
            {{new_g_sym}} := .data[[g]] |>
              cut_exposure_quantile(n_bins = n_bins %||% 4, ties = ties, quantile_type = quantile_type, labeller = labeller) |>
              .set_label(.get_label(dat[[g]]) %||% g)
          )

      } else {
        dat <- dat |>
          dplyr::mutate(
            {{new_g_sym}} := .data[[g]] |>
              cut_quantile(n_bins = n_bins %||% 4, ties = ties, quantile_type = quantile_type, labeller = labeller) |>
              .set_label(.get_label(dat[[g]]) %||% g)
          )
      }
      # tie-break rule actually used and quantile cutpoints (exposure
      # variable only -- `cut_quantile()` has no `"breaks"` attribute),
      # read back off the binned column's own attributes for
      # `.check_exposure_binning_consistency()` to compare against the
      # quantile layer's own `config$breaks`/`config$ties`
      config$ties <- attr(dat[[new_g]], "ties")
      config$breaks <- attr(dat[[new_g]], "breaks")
      g <- new_g
    }

    # store the variable names used for grouping
    if (stratify)  config$groupings <- c(g, object$strata$name)
    if (!stratify) config$groupings <- g

    # `er_plot_add_groups()` is additive -- each call may pass a different
    # `keep_strata`, so `stratify` is baked into each group's own config
    # (rather than only the shared `layer_group$stratify` used pre-additivity)
    # and `.build_group_plot()` reads it from here per group
    config$stratify <- stratify

    # store information about the y-axis variable
    config$y <- .plot_variable(
      name = g,
      label = .get_label(dat[[g]]) %||% g,
      role = paste("group", g, sep = "_")
    )

    # store sample size information (for merge into plot labels)
    config$counts <- dat |>
      dplyr::summarise(
        n   = sum(!is.na(.data[[object$exposure$name]])),
        lbl = paste0("N=", n),
        .by = config$groupings
      ) |>
      dplyr::mutate(lvl = paste0(.data[[g]], " (", lbl, ")")) |>
      dplyr::arrange(.data[[g]])

    # store the number of groups plotted on the y-axis
    config$n_groups <- nrow(config$counts)

    # store a modified data set to use for plotting
    config$data <- dat |>
      dplyr::select(dplyr::all_of(c(config$groupings, object$exposure$name))) |>
      dplyr::left_join(config$counts, by = config$groupings)

    layer_group$config[[g]] <- config
  }

  return(layer_group)
}

# Warns (at `er_plot_build()` time, so it catches either add-order) when
# `er_plot_add_quantiles()` and an `er_plot_add_groups()` call that groups
# by *the exposure variable itself* end up binning it differently.
# Deliberately narrow: two `er_plot_add_groups()` calls for two different
# (non-exposure) covariates are never compared against each other -- there's
# no reason they'd need to agree, so nothing there is checked. Compares
# `breaks`/`ties` (read back off each layer's own binned column, stored on
# `config$breaks`/`config$ties` by `.layer_quantile()`/`.layer_group()`
# above) rather than the raw `n_bins`/`ties`/`quantile_type` arguments, so a
# `quantile_type` difference that happens to produce identical breaks
# doesn't spuriously warn.
#' @noRd
.check_exposure_binning_consistency <- function(object) {
  quantile_config <- object$layer$quantile$config
  if (is.null(quantile_config)) return(invisible(NULL))

  exposure_key <- paste0(".", object$exposure$name, "_quantile")
  group_config <- object$layer$group$config[[exposure_key]]
  if (is.null(group_config)) return(invisible(NULL))

  breaks_differ <- !isTRUE(all.equal(
    unname(quantile_config$breaks), unname(group_config$breaks)
  ))
  ties_differ <- !identical(quantile_config$ties, group_config$ties)

  if (breaks_differ || ties_differ) {
    rlang::warn(c(
      sprintf(
        "`er_plot_add_quantiles()` and `er_plot_add_groups()` bin `%s` differently.",
        object$exposure$name
      ),
      "i" = sprintf(
        "Quantile layer: %d bin%s, ties = \"%s\". Group layer: %d bin%s, ties = \"%s\".",
        length(quantile_config$breaks) - 1, if (length(quantile_config$breaks) == 2) "" else "s",
        quantile_config$ties,
        length(group_config$breaks) - 1, if (length(group_config$breaks) == 2) "" else "s",
        group_config$ties
      ),
      "i" = "Pass matching `n_bins`/`ties`/`quantile_type` to both calls if you want the two panels' bins to line up, or ignore this warning if the difference is intentional."
    ))
  }
  invisible(NULL)
}


# miscellaneous helpers -------------------------------------------------------

.plot_variable <- function(name = NULL, label = NULL, limits = NULL, role = NULL, type = NULL) {
  list(name = name, label = label, limits = limits, role = role, type = type)
}

.get_strata_values <- function(data, name) {
  if (is.null(name)) return(NA)
  data[[name]]
}

# Recomputes the model layer's `config$predictions` from scratch, using
# `object`'s *current* `exposure`/`strata` rather than whatever they were
# when `er_plot_add_model()` was first called. Called unconditionally from
# `er_plot_build()` (before `.build_base_plot()`) whenever a model layer is
# present.
#
# Why this exists (issue #14): `.layer_model()` computes `config$predictions`
# eagerly, at add-layer time, over a `seq()` grid spanning
# `object$exposure$limits` as it stood *then*. A later
# `er_plot_theme(xlim = ...)` call mutates `object$exposure$limits` for the
# coordinate system, but never revisits an already-added model layer's
# cached predictions -- so the model curve/ribbon kept spanning the old
# (typically wider, data-derived) range while the panel's `coord_cartesian()`
# window narrowed to the new `xlim`. Because `.build_base_plot()` uses
# `clip = "off"` (needed so a point sitting exactly on a supplied limit
# isn't clipped at the panel edge -- see `AGENTS.md`), that stale excess
# wasn't cropped at the panel border either: it was drawn straight through
# it. Recomputing at build time -- right before the geoms are actually
# built -- makes the model layer's grid track `exposure$limits`/`strata`
# however they end up, regardless of what order `er_plot_add_model()` and
# `er_plot_theme(xlim = ...)` were called in (and, symmetrically, stretches
# to fill a *widened* `xlim` rather than leaving a short line that stops
# short of the new panel edge).
#
# `er_predict()` is a prediction call against an already-fitted model, not
# a refit -- re-running it here is cheap and doesn't touch model internals,
# consistent with "erplots never calls a model-fitting function".
.refresh_model_predictions <- function(object) {
  layer_model <- object$layer$model
  if (is.null(layer_model)) return(object)

  layer_model$config$predictions <- .get_model_predictions(
    layer_model$config$model,
    layer_model$config$conf_level,
    object$exposure,
    object$strata,
    layer_model$stratify,
    object$data,
    predict_args = layer_model$config$predict_args
  )

  object$layer$model <- layer_model
  object
}

.get_model_predictions <- function(model, conf_level, exposure, strata, stratify, data,
                                    predict_args = list()) {

  pred_dat <- seq(exposure$limits[1], exposure$limits[2], length.out = 300L) |>
    data.frame() |> .set_names(exposure$name)

  if (stratify) pred_dat <- pred_dat |>
    dplyr::cross_join(data.frame(strata$limits) |> .set_names(strata$name))

  # a fitted model's formula may reference covariates beyond the exposure
  # (and, when `stratify`, strata) variable -- erplots never inspects a
  # model's formula (it's deliberately model-agnostic), so instead of
  # guessing which covariates matter, fill in *every* other column of the
  # original fitting data at a single reference value. An `er_predict()`
  # method ignores extra columns it doesn't need, so this is harmless when
  # the model has no such covariates, and avoids the "object not found"
  # crash inside `predict()` when it does
  pred_dat <- .fill_reference_covariates(pred_dat, data)

  model_predictions <- rlang::exec(
    er_predict, model = model, newdata = pred_dat, conf_level = conf_level, !!!predict_args
  )
  return(model_predictions)
}

# fills every column of `data` not already present in `pred_dat` with a
# single reference value (recycled to `pred_dat`'s row count by ordinary
# vector recycling on assignment)
.fill_reference_covariates <- function(pred_dat, data) {
  extra_cols <- setdiff(names(data), names(pred_dat))
  for (col in extra_cols) {
    pred_dat[[col]] <- .reference_value(data[[col]])
  }
  pred_dat
}

# a single "typical" value for a covariate: first level for a factor
# (preserving its levels, so a model's contrasts still resolve correctly),
# first value (alphabetically) for a character vector, the mean for a
# numeric vector, `FALSE` for a logical vector, and the first observed
# value as a fallback for anything else (e.g. Date/POSIXct)
.reference_value <- function(x) {
  if (is.factor(x)) return(factor(levels(x)[1], levels = levels(x)))
  if (is.character(x)) return(sort(unique(x))[1])
  if (is.logical(x)) return(FALSE)
  if (is.numeric(x)) return(mean(x, na.rm = TRUE))
  x[1]
}