Add MOVER (Zou 2008) confidence-interval option for RERI

addint() and addintlist() gain `ci_method = c("delta", "mover")`. The default
("delta") keeps the existing symmetric Wald interval for backward compatibility;
"mover" uses the Method Of Variance Estimates Recovery interval of Zou (2008)
<doi:10.1093/aje/kwn104>, which recovers the variance of RERI from the component
risk ratios' confidence limits and generally has better coverage for this
bounded, skewed quantity. AP and S continue to use the delta method.

The implementation builds the (log RR11, log RR10, log RR01) contrast from the
model's coefficients and 3x3 vcov (handling both the joint and interaction
parameterisations) and applies Zou's MOVER-R formula with the component
correlations. Validated to match the interactionR package exactly on a shared
glm fit.

Adds a MOVER regression test (85 tests pass). R CMD check stays clean
(0 errors, 0 warnings, 1 expected NOTE).

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

Author: ehsanx

NAMESPACE

@@ -83,6 +83,7 @@ importFrom(scales,comma)
 importFrom(stats,as.formula)
 importFrom(stats,coef)
 importFrom(stats,confint)
+importFrom(stats,cov2cor)
 importFrom(stats,fitted)
 importFrom(stats,formula)
 importFrom(stats,model.frame)

R/addint.R

@@ -21,6 +21,9 @@
 #' S = (RR11 - 1) / ((RR10 - 1) + (RR01 - 1))
 #'
 #' Confidence intervals for S are calculated on the log scale and then exponentiated.
+#' By default the RERI interval uses the symmetric delta-method (Wald) interval;
+#' set \code{ci_method = "mover"} for the MOVER interval of Zou (2008)
+#' <doi:10.1093/aje/kwn104>, which generally has better coverage.
 #'
 #' \strong{Scale caveat:} these measures are defined on the risk-ratio scale.
 #' When the model supplies odds ratios (logistic) or hazard ratios (Cox), RERI,
@@ -44,13 +47,19 @@
 #' @param measures A character vector specifying which measures to calculate.
 #'   Options: `"RERI"`, `"AP"`, `"S"`, or `"all"`. Default is `"all"`.
 #' @param conf.level Confidence level for the interval (default 0.95).
+#' @param ci_method Character string giving the confidence-interval method for
+#'   RERI: `"delta"` (default) for the symmetric delta-method Wald interval, or
+#'   `"mover"` for the MOVER interval of Zou (2008)
+#'   <doi:10.1093/aje/kwn104>, which recovers the variance from the component
+#'   risk ratios' confidence limits and generally has better coverage. The AP and
+#'   S intervals always use the delta method.
 #'
 #' @return A list containing named vectors for each requested measure.
 #'   Each vector includes Estimate, SE (SE for RERI and AP, SE_log for S),
 #'   LowerCI, UpperCI. Returns `NULL` or a partial list with `NA` values if
 #'   calculation fails (e.g., missing coefficients, delta method error).
 #'
-#' @importFrom stats coef vcov qnorm pnorm na.omit
+#' @importFrom stats coef vcov qnorm pnorm na.omit cov2cor
 #' @importFrom msm deltamethod
 #'
 #' @export
@@ -89,9 +98,11 @@ addint <- function(model,
                    type = c("joint", "interaction"),
                    coef_names,
                    measures = "all",
-                   conf.level = 0.95) {
+                   conf.level = 0.95,
+                   ci_method = c("delta", "mover")) {
 
   type <- match.arg(type)
+  ci_method <- match.arg(ci_method)
 
   if (!requireNamespace("msm", quietly = TRUE)) {
     stop("Package 'msm' needed for this function to work. Please install it.", call. = FALSE)
@@ -224,6 +235,22 @@ addint <- function(model,
         lower_ci <- est_reri - z * se_reri
         upper_ci <- est_reri + z * se_reri
       }
+      # Optionally replace the symmetric Wald interval with the MOVER interval
+      # (Zou 2008), which recovers the variance from the component risk ratios'
+      # confidence limits and generally has better coverage for RERI.
+      if (identical(ci_method, "mover")) {
+        mover <- tryCatch(
+          .reri_mover_ci(selected_beta, selected_V, type, z, est_reri),
+          error = function(e) {
+            warning("MOVER interval failed for RERI; keeping the Wald interval: ",
+                    e$message)
+            NULL
+          })
+        if (!is.null(mover) && all(is.finite(mover))) {
+          lower_ci <- mover[1]
+          upper_ci <- mover[2]
+        }
+      }
     }
     results_list$RERI <- c(Estimate = est_reri, SE = se_reri, LowerCI = lower_ci, UpperCI = upper_ci)
     names(results_list$RERI) <- c("RERI_Estimate", "RERI_SE", paste0("RERI_CI", floor(conf.level*100), "_low"), paste0("RERI_CI", floor(conf.level*100), "_upp"))
@@ -289,3 +316,45 @@ addint <- function(model,
 
   return(results_list)
 }
+
+# Internal: MOVER (Method Of Variance Estimates Recovery) confidence interval for
+# RERI = RR11 - RR10 - RR01 + 1, following Zou (2008) <doi:10.1093/aje/kwn104>.
+# `selected_beta`/`selected_V` are the three model coefficients and their 3x3
+# variance-covariance matrix, in the order implied by `type`. Returns
+# c(lower, upper). Not exported.
+.reri_mover_ci <- function(selected_beta, selected_V, type, z, est_reri) {
+  # Contrast matrix mapping the coefficients to (log RR11, log RR10, log RR01).
+  if (identical(type, "interaction")) {
+    # selected_beta = (bA, bB, bAB): RR10 = exp(bA), RR01 = exp(bB),
+    # RR11 = exp(bA + bB + bAB).
+    C <- rbind(c(1, 1, 1), c(1, 0, 0), c(0, 1, 0))
+  } else {
+    # joint: selected_beta = (b10, b01, b11): RR10 = exp(b10), RR01 = exp(b01),
+    # RR11 = exp(b11).
+    C <- rbind(c(0, 0, 1), c(1, 0, 0), c(0, 1, 0))
+  }
+
+  log_rr <- as.numeric(C %*% as.numeric(selected_beta))
+  Sigma  <- C %*% selected_V %*% t(C)
+  se     <- sqrt(diag(Sigma))
+  if (any(!is.finite(se)) || any(se <= 0)) return(NULL)
+
+  R <- stats::cov2cor(Sigma)
+  r12 <- R[1, 2]; r13 <- R[1, 3]; r23 <- R[2, 3]   # (RR11,RR10),(RR11,RR01),(RR10,RR01)
+
+  theta <- exp(log_rr)                              # RR11, RR10, RR01
+  l <- exp(log_rr - z * se)
+  u <- exp(log_rr + z * se)
+
+  # RERI = theta1 - theta2 - theta3 + 1 (a1 = +1, a2 = a3 = -1).
+  var_lower <- (theta[1] - l[1])^2 + (u[2] - theta[2])^2 + (u[3] - theta[3])^2 -
+    2 * r12 * (theta[1] - l[1]) * (u[2] - theta[2]) -
+    2 * r13 * (theta[1] - l[1]) * (u[3] - theta[3]) +
+    2 * r23 * (u[2] - theta[2]) * (u[3] - theta[3])
+  var_upper <- (u[1] - theta[1])^2 + (theta[2] - l[2])^2 + (theta[3] - l[3])^2 -
+    2 * r12 * (u[1] - theta[1]) * (theta[2] - l[2]) -
+    2 * r13 * (u[1] - theta[1]) * (theta[3] - l[3]) +
+    2 * r23 * (theta[2] - l[2]) * (theta[3] - l[3])
+
+  c(est_reri - sqrt(max(0, var_lower)), est_reri + sqrt(max(0, var_upper)))
+}

R/addintlist.R

@@ -27,6 +27,8 @@
 #'   CIs in the final table. Default is 3. (Set to NULL for no rounding).
 #' @param verbose Logical. If `TRUE`, prints a progress message for each
 #'   factor-level combination as it is calculated. Defaults to `FALSE`.
+#' @param ci_method Confidence-interval method for RERI passed to
+#'   \code{\link{addint}}: `"delta"` (default) or `"mover"` (Zou 2008).
 #'
 #' @return A `tibble` (data frame) summarizing the additive interaction results.
 #'   Columns include:
@@ -96,7 +98,10 @@ addintlist <- function(model,
                        measures = "all",
                        conf.level = 0.95,
                        digits = 3,
-                       verbose = FALSE) {
+                       verbose = FALSE,
+                       ci_method = c("delta", "mover")) {
+
+  ci_method <- match.arg(ci_method)
 
   # Check if the core calculation function 'addint' exists
   if (!exists("addint") || !is.function(addint)) {
@@ -227,7 +232,8 @@ addintlist <- function(model,
         type = "interaction",
         coef_names = current_coef_names_list,
         measures = measures_to_calc, # Use the validated list
-        conf.level = conf.level
+        conf.level = conf.level,
+        ci_method = ci_method
       )
 
       # --- Process results into a flat tibble for this combination ---

man/addint.Rd

@@ -64,6 +64,36 @@ test_that("addint warns and returns NULL when a required coefficient is missing"
   expect_null(res)
 })
 
+test_that("addint MOVER interval for RERI is valid and asymmetric (Zou 2008)", {
+  beta <- c(A = 0.4, B = 0.3, `A:B` = 0.5)
+  V <- matrix(c(0.02, 0.005, 0.004,
+                0.005, 0.03, 0.006,
+                0.004, 0.006, 0.05), 3, 3)
+  dimnames(V) <- list(names(beta), names(beta))
+  m <- make_mock(beta, V)
+  cn <- list(exp1_coef = "A", exp2_coef = "B", inter_coef = "A:B")
+
+  delta <- addint(m, type = "interaction", coef_names = cn, measures = "RERI",
+                  ci_method = "delta")$RERI
+  mover <- addint(m, type = "interaction", coef_names = cn, measures = "RERI",
+                  ci_method = "mover")$RERI
+
+  # Same point estimate; MOVER interval differs from the symmetric Wald one.
+  expect_equal(unname(mover[["RERI_Estimate"]]),
+               unname(delta[["RERI_Estimate"]]), tolerance = 1e-10)
+  expect_lt(mover[["RERI_CI95_low"]], mover[["RERI_Estimate"]])
+  expect_gt(mover[["RERI_CI95_upp"]], mover[["RERI_Estimate"]])
+
+  # Asymmetric around the estimate (unlike the delta interval).
+  upper_half <- mover[["RERI_CI95_upp"]] - mover[["RERI_Estimate"]]
+  lower_half <- mover[["RERI_Estimate"]] - mover[["RERI_CI95_low"]]
+  expect_false(isTRUE(all.equal(upper_half, lower_half)))
+
+  # Pinned values, cross-checked against the interactionR package.
+  expect_equal(unname(mover[["RERI_CI95_low"]]), 0.3191, tolerance = 1e-3)
+  expect_equal(unname(mover[["RERI_CI95_upp"]]), 4.4499, tolerance = 1e-3)
+})
+
 test_that("addint validates the measures argument", {
   beta <- c(A = 0.4, B = 0.3, `A:B` = 0.5)
   V <- diag(3)