fix: improve calibration CIs, Gage R&R warnings, and COW dead code cleanup

- Calibration inverse prediction intervals now use proper ISO 11843 formula
  accounting for leverage and calibration center distance, not just sigma/slope
- Gage R&R now warns when negative variance components are clipped to zero,
  instead of silently forcing them
- Remove dead code line in COW alignment DP (prev_total was always 0)
- Update Gage R&R tests to suppressWarnings for random data scenarios

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: JamesHWade

R/align.R

@@ -1171,9 +1171,6 @@ tunable.step_measure_align_cow <- function(x, ...) {
     for (w in warp_options) {
       # Try all previous states
       for (prev_w in warp_options) {
-        prev_total <- (seg - 2) * 0 # Not tracking cumulative directly
-        # This approach: track segment-by-segment, each warp is local
-
         prev_state <- state_idx(prev_w)
         if (cost[seg - 1, prev_state] == Inf) {
           next

R/calibration-fit.R

@@ -537,29 +537,35 @@ measure_calibration_predict <- function(
   interval_type,
   level
 ) {
-  # Simplified interval calculation using delta method approximation
-  # For rigorous intervals, would need to implement Fieller's theorem
-
   fit <- object$model
   sigma <- summary(fit)$sigma
   n <- nrow(object$data)
-
-  # Get slope for delta method
-  coefs <- stats::coef(fit)
-  slope <- coefs[object$conc_col]
-
-  # Standard error of inverse prediction (approximate)
-  se_pred <- sigma / abs(slope)
-
-  # Additional uncertainty for prediction intervals
-  if (interval_type == "prediction") {
-    se_pred <- se_pred * sqrt(1 + 1 / n)
-  }
-
-  # Critical value
   df <- summary(fit)$df[2]
   t_crit <- stats::qt((1 + level) / 2, df)
 
+  if (object$model_type == "linear") {
+    # Classical inverse prediction intervals (ISO 11843 / Draper & Smith)
+    # Confidence: SE = (sigma/|b|) * sqrt(1/n + (x0 - xbar)^2 / Sxx)
+    # Prediction: SE = (sigma/|b|) * sqrt(1 + 1/n + (x0 - xbar)^2 / Sxx)
+    coefs <- stats::coef(fit)
+    slope <- coefs[object$conc_col]
+    cal_conc <- object$data[[object$conc_col]]
+    x_bar <- mean(cal_conc)
+    sxx <- sum((cal_conc - x_bar)^2)
+
+    extra <- if (interval_type == "prediction") 1 else 0
+    se_pred <- (sigma / abs(slope)) *
+      sqrt(extra + 1 / n + (pred_conc - x_bar)^2 / sxx)
+  } else {
+    # Quadratic models: use delta method approximation
+    coefs <- stats::coef(fit)
+    slope <- coefs[object$conc_col]
+    se_pred <- sigma / abs(slope)
+    if (interval_type == "prediction") {
+      se_pred <- se_pred * sqrt(1 + 1 / n)
+    }
+  }
+
   list(
     lower = pred_conc - t_crit * se_pred,
     upper = pred_conc + t_crit * se_pred

R/precision.R

@@ -409,13 +409,28 @@ measure_gage_rr <- function(
   ms_error <- anova_table["Residuals", "Mean Sq"]
 
   # Calculate variance components (EMS method)
+  # Negative estimates are set to zero (standard practice for ANOVA-based
+
+  # variance components), but we warn since this can indicate poor model fit
   var_repeatability <- ms_error
-  var_interaction <- max(0, (ms_interaction - ms_error) / n_replicates)
-  var_operator <- max(
-    0,
-    (ms_operator - ms_interaction) / (n_parts * n_replicates)
+  raw_interaction <- (ms_interaction - ms_error) / n_replicates
+  raw_operator <- (ms_operator - ms_interaction) / (n_parts * n_replicates)
+  raw_part <- (ms_part - ms_interaction) / (n_operators * n_replicates)
+
+  clipped <- c(
+    if (raw_interaction < 0) "interaction",
+    if (raw_operator < 0) "operator",
+    if (raw_part < 0) "part"
   )
-  var_part <- max(0, (ms_part - ms_interaction) / (n_operators * n_replicates))
+  if (length(clipped) > 0) {
+    cli::cli_warn(
+      "Negative variance component{?s} set to zero for {.val {clipped}}. This may indicate insufficient variation or poor model fit."
+    )
+  }
+
+  var_interaction <- max(0, raw_interaction)
+  var_operator <- max(0, raw_operator)
+  var_part <- max(0, raw_part)
 
   # Combine components
   var_reproducibility <- var_operator + var_interaction

tests/testthat/test-precision.R

@@ -216,12 +216,12 @@ test_that("measure_gage_rr calculates variance components", {
     ifelse(data$operator == "A", 0.5, ifelse(data$operator == "B", -0.3, 0)) +
     rnorm(nrow(data), 0, 0.5)
 
-  result <- measure_gage_rr(
+  result <- suppressWarnings(measure_gage_rr(
     data,
     response_col = "measurement",
     part_col = "part",
     operator_col = "operator"
-  )
+  ))
 
   expect_s3_class(result, "measure_gage_rr")
   expect_equal(nrow(result), 4)
@@ -240,12 +240,12 @@ test_that("measure_gage_rr calculates percentages correctly", {
   )
   data$measurement <- rnorm(nrow(data), 100, 2)
 
-  result <- measure_gage_rr(
+  result <- suppressWarnings(measure_gage_rr(
     data,
     response_col = "measurement",
     part_col = "part",
     operator_col = "operator"
-  )
+  ))
 
   # Percentages should sum appropriately
   expect_true(all(result$pct_contribution >= 0))
@@ -261,13 +261,13 @@ test_that("measure_gage_rr calculates tolerance percentage when provided", {
   )
   data$measurement <- 50 + (data$part - 3) * 5 + rnorm(nrow(data), 0, 1)
 
-  result <- measure_gage_rr(
+  result <- suppressWarnings(measure_gage_rr(
     data,
     response_col = "measurement",
     part_col = "part",
     operator_col = "operator",
     tolerance = 20
-  )
+  ))
 
   expect_false(anyNA(result$pct_tolerance))
   expect_true(all(result$pct_tolerance >= 0))
@@ -282,12 +282,12 @@ test_that("measure_gage_rr stores metadata", {
   )
   data$measurement <- rnorm(nrow(data))
 
-  result <- measure_gage_rr(
+  result <- suppressWarnings(measure_gage_rr(
     data,
     response_col = "measurement",
     part_col = "part",
     operator_col = "operator"
-  )
+  ))
 
   expect_equal(attr(result, "n_parts"), 10)
   expect_equal(attr(result, "n_operators"), 3)
@@ -318,13 +318,13 @@ test_that("measure_gage_rr print method works", {
   )
   data$measurement <- rnorm(nrow(data), 50, 2)
 
-  result <- measure_gage_rr(
+  result <- suppressWarnings(measure_gage_rr(
     data,
     response_col = "measurement",
     part_col = "part",
     operator_col = "operator",
     tolerance = 20
-  )
+  ))
 
   expect_output(print(result), "measure_gage_rr")
   expect_output(print(result), "Variance Components")
@@ -345,12 +345,12 @@ test_that("measure_gage_rr ndc calculation is reasonable", {
     (data$part - 5) * 10 + # Large part variation
     rnorm(nrow(data), 0, 0.5) # Small measurement error
 
-  result <- measure_gage_rr(
+  result <- suppressWarnings(measure_gage_rr(
     data,
     response_col = "measurement",
     part_col = "part",
     operator_col = "operator"
-  )
+  ))
 
   # With large part variation and small error, ndc should be high
   expect_true(attr(result, "ndc") >= 5)
@@ -378,12 +378,12 @@ test_that("tidy.measure_gage_rr returns tibble", {
   )
   data$measurement <- rnorm(nrow(data))
 
-  result <- measure_gage_rr(
+  result <- suppressWarnings(measure_gage_rr(
     data,
     response_col = "measurement",
     part_col = "part",
     operator_col = "operator"
-  )
+  ))
 
   tidy_result <- tidy(result)
   expect_s3_class(tidy_result, "tbl_df")