R/utils.R

In FPScausal: Functional Propensity Score for Causal Inference

# Internal utility functions for FPScausal

#' @importFrom utils globalVariables
# Declare ggplot2 aesthetics to avoid R CMD check NOTEs
utils::globalVariables(c(
  "variable", "value", "weighted", "FPC", "confounder",
  "correlation", "type", "time", "subject", "Component",
  "Eigenvalue", "CumulativePVE", "PVE", "weight",
  "s", "t", "effect", "lower", "upper", "significant",
  "beta_type", "domain_val",
  "label", "group", "region", "sig_label",
  "start", "end", "effect_uw",
  "curve_type", "effect_true"
))

#' Detect the type of an outcome variable
#'
#' @param outcome A numeric vector or matrix.
#' @return One of `"scalar"`, `"binary"`, or `"functional"`.
#' @keywords internal
.detect_outcome_type <- function(outcome) {
  if (is.matrix(outcome) || inherits(outcome, "fd")) {
    return("functional")
  }
  outcome_vec <- as.numeric(outcome)
  unique_vals <- unique(outcome_vec[!is.na(outcome_vec)])
  if (length(unique_vals) <= 2 && all(unique_vals %in% c(0, 1))) {
    return("binary")
  }
  "scalar"
}

#' Auto-select number of B-spline basis functions
#'
#' @param t_grid Numeric vector of observation time points.
#' @return Integer: suggested nbasis.
#' @keywords internal
.auto_nbasis <- function(t_grid) {
  max(10L, round(length(t_grid) * 0.6))
}

#' Check for domain overlap between treatment and functional outcome
#'
#' Emits a warning when the treatment and outcome share the same domain name
#' and the domains numerically overlap, because the historical constraint
#' (non-anticipativity) may be needed but is not enforced.
#'
#' @param domain_t Numeric c(a, b) treatment domain.
#' @param name_t   Character treatment domain name.
#' @param domain_o Numeric c(a, b) outcome domain.
#' @param name_o   Character outcome domain name.
#' @keywords internal
.check_domains_overlap <- function(domain_t, name_t, domain_o, name_o) {
  if (identical(name_t, name_o)) {
    overlap_start <- max(domain_t[1], domain_o[1])
    overlap_end   <- min(domain_t[2], domain_o[2])
    if (overlap_start < overlap_end) {
      warning(
        "Treatment domain (", name_t, ") and outcome domain (", name_o, ") ",
        "share the same name and overlap on [", overlap_start, ", ", overlap_end, "]. ",
        "If causal ordering requires treatment to precede outcome (non-anticipativity), ",
        "consider using non-overlapping domains."
      )
    }
  }
}

#' Extract treatment matrix from fps_weighting or fps_effect_estimation object
#'
#' @param fps_obj An `fps_weighting` or `fps_effect_estimation` object.
#' @return The `fpca_treatment` list element.
#' @keywords internal
.get_fpca_treatment <- function(fps_obj) {
  if (inherits(fps_obj, "fps_weighting")) {
    return(fps_obj$fpca_treatment)
  }
  if (inherits(fps_obj, "fps_effect_estimation")) {
    return(fps_obj$fpca_treatment)
  }
  stop("Object must be of class 'fps_weighting' or 'fps_effect_estimation'.")
}

#' Check that an object is of class fps_weighting
#' @keywords internal
.check_fps_weighting <- function(fps_object) {
  if (!inherits(fps_object, "fps_weighting")) {
    stop("'fps_object' must be an object returned by fps_weighting().")
  }
}

#' Reflected-percentile bootstrap confidence interval
#'
#' @param obs_val Observed statistic (scalar).
#' @param boot_vals Numeric vector of bootstrap replicates.
#' @param alpha Significance level.
#' @return Named numeric vector c(lwr, obs, upr).
#' @keywords internal
.reflected_ci <- function(obs_val, boot_vals, alpha = 0.05) {
  q_hi <- stats::quantile(boot_vals, 1 - alpha / 2, names = FALSE, na.rm = TRUE)
  q_lo <- stats::quantile(boot_vals, alpha / 2,     names = FALSE, na.rm = TRUE)
  c(
    lwr = obs_val - (q_hi - obs_val),
    obs = obs_val,
    upr = obs_val - (q_lo - obs_val)
  )
}

#' Compute MISE, AISE, and ISB between estimated and true beta
#'
#' @param beta_hat Numeric vector (or matrix for functional case).
#' @param beta_true Numeric vector (or matrix) of the same shape.
#' @return Named numeric vector with MISE, AISE, ISB (for single estimate, MISE=AISE=ISB=ISE).
#' @keywords internal
.compute_error_metrics <- function(beta_hat, beta_true) {
  if (is.matrix(beta_hat)) {
    ise <- mean((beta_hat - beta_true)^2)
    bias2 <- mean((beta_hat - beta_true))^2
  } else {
    ise <- mean((beta_hat - beta_true)^2)
    bias2 <- mean((beta_hat - beta_true))^2
  }
  c(ISE = ise, ISB = bias2)
}

#' Find contiguous intervals where a logical mask is TRUE
#'
#' Returns a data frame with `start` and `end` columns (values from t_grid)
#' for each run of TRUE values in `mask`.
#' @param t_grid Numeric vector.
#' @param mask   Logical vector of same length as t_grid.
#' @return data.frame with columns start, end (or 0-row frame if no TRUE).
#' @keywords internal
.contiguous_intervals <- function(t_grid, mask) {
  if (!any(mask, na.rm = TRUE)) {
    return(data.frame(start = numeric(0), end = numeric(0)))
  }
  r      <- rle(as.logical(mask))
  ends   <- cumsum(r$lengths)
  starts <- c(1L, ends[-length(ends)] + 1L)
  idx    <- which(r$values)
  data.frame(start = t_grid[starts[idx]], end = t_grid[ends[idx]])
}