R/compare_waves.R

In wavecov: Overlap-Aware Workforce Wave Comparisons

#' Compare two workforce waves for a single metric
#'
#' Compute overlap-aware inference for the change in a repeated organizational
#' metric between two workforce waves, using linked stayers, a reliability
#' assumption, or Cauchy-based bounds.
#'
#' @param data Data frame with at least id, wave, and metric columns.
#' @param wave_0 Value in `wave_col` for the first wave.
#' @param wave_1 Value in `wave_col` for the second wave.
#' @param id_col Column name for the employee identifier.
#' @param wave_col Column name for the wave identifier.
#' @param metric_col Column name for the metric.
#' @param method Method to use: "auto", "linked", "reliability", or "bounds".
#'   In "auto" mode, `wavecov` prefers linked, then reliability, then bounds.
#' @param reliability Same-person cross-wave reliability for the relevant
#'   stayer population (used only for the reliability method).
#' @param sigma2 Optional common individual-level variance. When supplied with
#'   the reliability method, overrides the pooled variance from the two waves.
#' @param sigma2S Optional pooled individual-level stayer variance. When
#'   supplied with the reliability method, is treated as a candidate `sigma2`.
#'   When supplied with the bounds method and `overlap_mode = "overlap_aware"`,
#'   is used as a common stayer variance to tighten the Cauchy bounds.
#' @param cov_sign Covariance sign restriction for bounds: "unrestricted" or
#'   "nonnegative".
#' @param overlap_mode Whether to use generic Cauchy bounds ("generic") or
#'   overlap-aware tightening ("overlap_aware").
#' @param level Confidence level for intervals (point-estimate methods) or
#'   sensitivity envelopes (bounded results).
#' @param ... Reserved for future extensions.
#'
#' @return A "wavecov_result" object (linked or reliability) or a
#'   "wavecov_bounds_result" object (bounds). Use `print()`, `tidy()`,
#'   and `confint()` for summaries.
#'
#' @export
compare_waves <- function(data,
                          wave_0,
                          wave_1,
                          id_col,
                          wave_col,
                          metric_col,
                          method = c("auto", "linked", "reliability", "bounds"),
                          reliability = NULL,
                          sigma2 = NULL,
                          sigma2S = NULL,
                          cov_sign = c("unrestricted", "nonnegative"),
                          overlap_mode = c("generic", "overlap_aware"),
                          level = 0.95,
                          ...) {
  method       <- base::match.arg(method)
  cov_sign     <- base::match.arg(cov_sign)
  overlap_mode <- base::match.arg(overlap_mode)
  
  if (missing(data) || is.null(data)) {
    stop("data must be supplied and cannot be NULL.")
  }
  
  summ <- summarize_panel(
    data       = data,
    wave_0     = wave_0,
    wave_1     = wave_1,
    id_col     = id_col,
    wave_col   = wave_col,
    metric_col = metric_col
  )
  
  # Effective sigma2 for reliability: explicit sigma2 wins, else sigma2S,
  # else pooled from V0_ind/V1_ind inside variance_reliability().
  sigma2_eff <- if (!is.null(sigma2)) sigma2 else sigma2S
  
  chosen_method <- method
  if (identical(method, "auto")) {
    if (!is.na(summ$C01_ind) && summ$nS > 1L) {
      chosen_method <- "linked"
    } else if (!is.null(reliability)) {
      chosen_method <- "reliability"
    } else {
      chosen_method <- "bounds"
    }
  }
  
  if (identical(chosen_method, "linked")) {
    
    var_out <- variance_linked(
      N0      = summ$N0,
      N1      = summ$N1,
      nS      = summ$nS,
      V0_ind  = summ$V0_ind,
      V1_ind  = summ$V1_ind,
      C01_ind = summ$C01_ind
    )
    
    res <- new_wavecov_result(
      summary = summ,
      var_out = var_out,
      method  = "linked",
      level   = level
    )
    
    # Add normal-approximation confidence interval fields
    alpha <- 1 - level
    z     <- stats::qnorm(1 - alpha / 2)
    res$ci_lower <- res$estimate - z * res$se
    res$ci_upper <- res$estimate + z * res$se
    
  } else if (identical(chosen_method, "reliability")) {
    
    var_out <- variance_reliability(
      N0         = summ$N0,
      N1         = summ$N1,
      nS         = summ$nS,
      V0_ind     = summ$V0_ind,
      V1_ind     = summ$V1_ind,
      sigma2     = sigma2_eff,
      reliability = reliability
    )
    
    res <- new_wavecov_result(
      summary = summ,
      var_out = var_out,
      method  = "reliability",
      level   = level
    )
    
    # Add normal-approximation confidence interval fields
    alpha <- 1 - level
    z     <- stats::qnorm(1 - alpha / 2)
    res$ci_lower <- res$estimate - z * res$se
    res$ci_upper <- res$estimate + z * res$se
    
  } else if (identical(chosen_method, "bounds")) {
    
    V0_bar <- summ$V0_ind / summ$N0
    V1_bar <- summ$V1_ind / summ$N1
    
    bounds <- variance_bounds(
      N0          = summ$N0,
      N1          = summ$N1,
      V0_bar      = V0_bar,
      V1_bar      = V1_bar,
      nS          = if (identical(overlap_mode, "overlap_aware")) summ$nS else NULL,
      V0S_ind     = if (identical(overlap_mode, "overlap_aware")) summ$V0S_ind else NULL,
      V1S_ind     = if (identical(overlap_mode, "overlap_aware")) summ$V1S_ind else NULL,
      sigma2S     = if (identical(overlap_mode, "overlap_aware")) sigma2S else NULL,
      cov_sign    = cov_sign,
      overlap_mode = overlap_mode
    )
    
    res <- new_wavecov_bounds_result(
      summary = summ,
      bounds  = bounds,
      level   = level
    )
    
    # For bounds objects, provide placeholders so $ci_lower / $ci_upper always exist.
    # Actual uncertainty is represented by `bounds`, not simple symmetric CIs.
    res$ci_lower <- NA_real_
    res$ci_upper <- NA_real_
    
  } else {
    stop("Unsupported method: ", chosen_method)
  }
  
  res
}