R/variance_reliability.R

In wavecov: Overlap-Aware Workforce Wave Comparisons

#' Reliability-based variance for change in means
#'
#' @keywords internal
variance_reliability <- function(N0, N1, nS,
                                 V0_ind = NULL,
                                 V1_ind = NULL,
                                 sigma2 = NULL,
                                 reliability) {
  if (is.null(reliability) || !is.numeric(reliability) || length(reliability) != 1L) {
    stop("Reliability must be a single numeric value.")
  }
  if (reliability < 0 || reliability > 1) {
    stop("Reliability must lie in [0, 1].")
  }
  
  if (!is.null(V0_ind) && !is.null(V1_ind)) {
    if (is.null(sigma2)) {
      sigma2 <- ((N0 - 1) * V0_ind + (N1 - 1) * V1_ind) / (N0 + N1 - 2)
    } else {
      if (!is.numeric(sigma2) || length(sigma2) != 1L || sigma2 <= 0) {
        stop("sigma2 must be a positive numeric scalar.")
      }
      ratio <- max(V0_ind, V1_ind) / min(V0_ind, V1_ind)
      if (ratio > 2) {
        warning("Wave-specific individual-level variances differ by more than a factor of 2; ",
                "the common-variance reliability approximation may be fragile.")
      }
    }
  } else if (is.null(sigma2)) {
    stop("Reliability method requires either V0_ind and V1_ind or a common sigma2.")
  } else if (!is.numeric(sigma2) || length(sigma2) != 1L || sigma2 <= 0) {
    stop("sigma2 must be a positive numeric scalar.")
  }
  
  V0_bar <- sigma2 / N0
  V1_bar <- sigma2 / N1
  
  C01_ind <- reliability * sigma2
  C01_bar <- nS * C01_ind / (N0 * N1)
  
  var_Delta <- V0_bar + V1_bar - 2 * C01_bar
  if (var_Delta < 0) {
    warning("Computed variance under the reliability method is negative; setting to NA.")
    var_Delta <- NA_real_
  }
  
  list(
    sigma2    = sigma2,
    V0_bar    = V0_bar,
    V1_bar    = V1_bar,
    C01_ind   = C01_ind,
    C01_bar   = C01_bar,
    var_Delta = var_Delta,
    reliability = reliability
  )
}