1482: Add CMH test with Sato variance estimator

.lintr

@@ -1,6 +1,5 @@
 linters: linters_with_defaults(
   line_length_linter = line_length_linter(120),
-  cyclocomp_linter = NULL,
   object_usage_linter = NULL,
   object_length_linter = NULL,
   pipe_consistency_linter = NULL

NEWS.md

@@ -1,19 +1,21 @@
 # tern 0.9.10.9012
 
-### Miscellaneous
-* Updated `roxygen2` to 8.0.0 and added `@exportS3Method` tags for S3 methods in `decorate_grob.R` and `utils_grid.R`.
-* Converted  `s_surv_time()`to exported functions.
-
 ### Enhancements
 * Added `alternative` argument to `s_coxph_pairwise()` to allow one-sided hypothesis testing. 
 * Added `lr_stat_df` to the parameters return list of `s_coxph_pairwise()`.
 * Added `uncond_exact_diff` method to `estimate_proportion_diff()` for the unconditional exact confidence interval for the difference in proportions by inverting one-sided tail tests over a nuisance parameter.
+* Added `cmh_sato` method to `test_proportion_diff()` for CMH testing with the Sato variance estimator. (#1482)
 * Added `range_with_cens_info` statistic to `s_surv_time()`.
 * Added `lsmean_se`, `lsmean_ci`, and `lsmean_diffci` statistics to `s_ancova()`.
 * Added `s_ancova()` to exported functions.
 
 ### Bug Fixes
 * Fixed bug in `prop_diff_cmh()` which previously failed when strata combinations had 0 observations.
+* Fixed one-sided p-values in `prop_cmh()` with Wilson-Hilferty transformation — the sign of the effect was lost, producing incorrect p-values for `alternative = "less"` and `alternative = "greater"`.
+
+### Miscellaneous
+* Updated `roxygen2` to 8.0.0 and added `@exportS3Method` tags for S3 methods in `decorate_grob.R` and `utils_grid.R`.
+* Converted `s_surv_time()` to exported functions.
 
 # tern 0.9.10
 

R/coxph.R

@@ -25,8 +25,7 @@ univariate <- function(x) {
 # Get the right-hand-term of a formula
 rht <- function(x) {
   checkmate::assert_formula(x)
-  y <- as.character(rev(x)[[1]])
-  return(y)
+  as.character(rev(x)[[1]])
 }
 
 #' Hazard ratio estimation in interactions
@@ -109,7 +108,7 @@ estimate_coef <- function(variable, given,
   design_mat <- apply(
     X = design_mat, MARGIN = 1, FUN = function(x) {
       mmat[names(mmat) %in% x[-which(names(x) == "given")]] <- 1
-      return(mmat)
+      mmat
     }
   )
   colnames(design_mat) <- interaction_names
@@ -124,7 +123,7 @@ estimate_coef <- function(variable, given,
     y <- vcov[vcov_el, vcov_el]
     y <- sum(y)
     y <- sqrt(y)
-    return(y)
+    y
   })
 
   q_norm <- stats::qnorm((1 + conf_level) / 2)
@@ -135,7 +134,7 @@ estimate_coef <- function(variable, given,
     x["lcl"] <- exp(x["coef"] - q_norm * x["se(coef)"])
     x["ucl"] <- exp(x["coef"] + q_norm * x["se(coef)"])
 
-    return(x)
+    x
   })
 
   y <- t(y)
@@ -147,7 +146,7 @@ estimate_coef <- function(variable, given,
     " hazard ratio given the level of ", given, " compared to ",
     variable, " level ", lvl_var[1], "."
   )
-  return(y)
+  y
 }
 
 #' `tryCatch` around `car::Anova`
@@ -195,11 +194,10 @@ try_car_anova <- function(mod,
         invokeRestart("muffleWarning")
       }
     ),
-    finally = {
-    }
+    finally = {}
   )
 
-  return(y)
+  y
 }
 
 #' Fit a Cox regression model and ANOVA
@@ -241,7 +239,7 @@ fit_n_aov <- function(formula,
   y <- list(mod = mod, msum = msum, aov = aov)
   attr(y, "message") <- warn_attr
 
-  return(y)
+  y
 }
 
 # argument_checks
@@ -266,7 +264,7 @@ name_covariate_names <- function(covariates) {
   no_names <- is.null(names(covariates))
   if (any(miss_names)) names(covariates)[miss_names] <- vapply(covariates[miss_names], FUN = rht, FUN.VALUE = "name")
   if (no_names) names(covariates) <- vapply(covariates, FUN = rht, FUN.VALUE = "name")
-  return(covariates)
+  covariates
 }
 
 check_increments <- function(increments, covariates) {

R/h_step.R

@@ -171,8 +171,7 @@ h_step_survival_est <- function(formula,
         invokeRestart("muffleWarning")
       }
     ),
-    finally = {
-    }
+    finally = {}
   )
   if (!is.null(coxph_warnings)) {
     warning(paste(
@@ -285,8 +284,7 @@ h_step_rsp_est <- function(formula,
         invokeRestart("muffleWarning")
       }
     ),
-    finally = {
-    }
+    finally = {}
   )
   if (!is.null(fit_warnings)) {
     warning(paste(

R/prop_diff.R

@@ -596,6 +596,81 @@ prop_diff_nc <- function(rsp,
   )
 }
 
+#' @describeIn h_prop_diff Helper function to calculate the CMH weighted
+#' difference in proportions.
+#'
+#' @param tbl (`array`)\cr 3-dimensional array with dimensions corresponding to
+#'   group, response, and strata. The second dimension (response) should have names
+#'   "TRUE" and "FALSE".
+#'
+#' @keywords internal
+h_diff_cmh <- function(tbl) {
+  checkmate::assert_array(tbl)
+  checkmate::assert_integer(c(ncol(tbl), nrow(tbl)), lower = 2, upper = 2)
+  checkmate::assert_integer(length(dim(tbl)), lower = 3, upper = 3)
+  checkmate::assert_true(identical(dimnames(tbl)[[2]], c("TRUE", "FALSE")))
+
+  n1 <- colSums(tbl[1, 1:2, ])
+  n2 <- colSums(tbl[2, 1:2, ])
+  x1 <- tbl[1, 1, ]
+  p1 <- x1 / n1
+  x2 <- tbl[2, 1, ]
+  p2 <- x2 / n2
+  # CMH weights
+  use_stratum <- (n1 > 0) & (n2 > 0)
+  n1 <- n1[use_stratum]
+  n2 <- n2[use_stratum]
+  p1 <- p1[use_stratum]
+  p2 <- p2[use_stratum]
+  wt <- (n1 * n2 / (n1 + n2))
+  wt_normalized <- wt / sum(wt)
+  est1 <- sum(wt_normalized * p1)
+  est2 <- sum(wt_normalized * p2)
+  diff_est <- est2 - est1
+
+  list(
+    x1 = x1,
+    n1 = n1,
+    p1 = p1,
+    x2 = x2,
+    n2 = n2,
+    p2 = p2,
+    wt = wt,
+    wt_normalized = wt_normalized,
+    est1 = est1,
+    est2 = est2,
+    diff_est = diff_est
+  )
+}
+
+#' @describeIn h_prop_diff Helper function to calculate the standard error for the
+#'   CMH weighted difference in proportions.
+#'
+#' @param cmh_results (`list`)\cr output of [h_diff_cmh()].
+#'
+#' @keywords internal
+h_diff_cmh_se <- function(cmh_results, diff_se = c("standard", "sato")) {
+  checkmate::assert_list(cmh_results, types = "numeric", any.missing = FALSE, names = "unique")
+  diff_se <- match.arg(diff_se)
+  l <- cmh_results # For easier readability of the formulas below.
+  if (diff_se == "standard") {
+    terms1 <- l$p1 * (1 - l$p1) / l$n1
+    terms2 <- l$p2 * (1 - l$p2) / l$n2
+    sqrt(sum((terms1 + terms2) * l$wt_normalized^2))
+  } else {
+    # Sato variance estimator.
+    p_terms_num <- l$n2^2 * l$x1 -
+      l$n1^2 * l$x2 +
+      l$n1 * l$n2 * (l$n1 - l$n2) / 2
+    p_terms <- p_terms_num / (l$n1 + l$n2)^2
+    q_terms_num <- l$x1 * (l$n2 - l$x2) + l$x2 * (l$n1 - l$x1)
+    q_terms <- q_terms_num / (2 * (l$n1 + l$n2))
+    num <- l$diff_est * sum(p_terms) + sum(q_terms)
+    denom <- sum(l$wt)^2
+    sqrt(num / denom)
+  }
+}
+
 #' @describeIn h_prop_diff Calculates the weighted difference. This is defined as the difference in
 #'   response rates between the experimental treatment group and the control treatment group, adjusted
 #'   for stratification factors by applying Cochran-Mantel-Haenszel (CMH) weights. For the CMH chi-squared
@@ -645,94 +720,71 @@ prop_diff_cmh <- function(rsp,
     warning("Less than 5 observations in some strata.")
   }
 
-  # first dimension: FALSE, TRUE
-  # 2nd dimension: CONTROL, TX
+  # first dimension: CONTROL, TX
+  # 2nd dimension: TRUE, FALSE
   # 3rd dimension: levels of strata
-  # rsp as factor rsp to handle edge case of no FALSE (or TRUE) rsp records
+  # Note: rsp needs to be a factor to handle edge case of
+  #       no FALSE (or TRUE) rsp records.
   t_tbl <- table(
-    factor(rsp, levels = c("FALSE", "TRUE")),
     grp,
+    factor(rsp, levels = c("TRUE", "FALSE")),
     strata
   )
-  n1 <- colSums(t_tbl[1:2, 1, ])
-  n2 <- colSums(t_tbl[1:2, 2, ])
-  x1 <- t_tbl[2, 1, ]
-  p1 <- x1 / n1
-  x2 <- t_tbl[2, 2, ]
-  p2 <- x2 / n2
-  # CMH weights
-  use_stratum <- (n1 > 0) & (n2 > 0)
-  n1 <- n1[use_stratum]
-  n2 <- n2[use_stratum]
-  p1 <- p1[use_stratum]
-  p2 <- p2[use_stratum]
-  wt <- (n1 * n2 / (n1 + n2))
-  wt_normalized <- wt / sum(wt)
-  est1 <- sum(wt_normalized * p1)
-  est2 <- sum(wt_normalized * p2)
-  estimate <- c(est1, est2)
+  cmh <- h_diff_cmh(t_tbl)
+
+  estimate <- c(cmh$est1, cmh$est2)
   names(estimate) <- levels(grp)
-  se1 <- sqrt(sum(wt_normalized^2 * p1 * (1 - p1) / n1))
-  se2 <- sqrt(sum(wt_normalized^2 * p2 * (1 - p2) / n2))
+  se1 <- sqrt(sum(cmh$wt_normalized^2 * cmh$p1 * (1 - cmh$p1) / cmh$n1))
+  se2 <- sqrt(sum(cmh$wt_normalized^2 * cmh$p2 * (1 - cmh$p2) / cmh$n2))
   z <- stats::qnorm((1 + conf_level) / 2)
   err1 <- z * se1
   err2 <- z * se2
-  ci1 <- c((est1 - err1), (est1 + err1))
-  ci2 <- c((est2 - err2), (est2 + err2))
+  ci1 <- c((cmh$est1 - err1), (cmh$est1 + err1))
+  ci2 <- c((cmh$est2 - err2), (cmh$est2 + err2))
   estimate_ci <- list(ci1, ci2)
   names(estimate_ci) <- levels(grp)
-  diff_est <- est2 - est1
 
   if (diff_se %in% c("standard", "sato")) {
-    se_diff <- if (diff_se == "standard") {
-      sqrt(sum(((p1 * (1 - p1) / n1) + (p2 * (1 - p2) / n2)) * wt_normalized^2))
-    } else {
-      # Sato variance estimator.
-      p_terms <- (n2^2 * x1 - n1^2 * x2 + n1 * n2 * (n1 - n2) / 2) / (n1 + n2)^2
-      q_terms <- (x1 * (n2 - x2) + x2 * (n1 - x1)) / (2 * (n1 + n2))
-      num <- diff_est * sum(p_terms) + sum(q_terms)
-      denom <- sum(wt)^2
-      sqrt(num / denom)
-    }
-    diff_ci <- c(diff_est - z * se_diff, diff_est + z * se_diff)
+    se_diff <- h_diff_cmh_se(cmh, diff_se = diff_se)
+    diff_ci <- c(cmh$diff_est - z * se_diff, cmh$diff_est + z * se_diff)
   } else {
     # Miettinen and Nurminen method is used.
     z_stat_fun <- function(delta) {
       var_est <- h_miettinen_nurminen_var_est(
-        n1 = n1, n2 = n2,
-        x1 = x1, x2 = x2,
+        n1 = cmh$n1, n2 = cmh$n2,
+        x1 = cmh$x1, x2 = cmh$x2,
         diff_par = delta
       )$var_est
-      num <- sum(wt * (p2 - p1 - delta))
-      denom <- sqrt(sum(wt^2 * var_est))
+      num <- sum(cmh$wt * (cmh$p2 - cmh$p1 - delta))
+      denom <- sqrt(sum(cmh$wt^2 * var_est))
       num / denom
     }
     # Find upper and lower confidence limits by root finding such that
     # z_stat_fun(limit) = +/- z quantile:
     root_lower <- function(delta) z_stat_fun(delta) - z
     root_upper <- function(delta) z_stat_fun(delta) + z
     diff_ci <- c(
-      stats::uniroot(root_lower, interval = c(-0.99, diff_est))$root,
-      stats::uniroot(root_upper, interval = c(diff_est, 0.99))$root
+      stats::uniroot(root_lower, interval = c(-0.99, cmh$diff_est))$root,
+      stats::uniroot(root_upper, interval = c(cmh$diff_est, 0.99))$root
     )
     # Calculate the standard error separately.
     var_est <- h_miettinen_nurminen_var_est(
-      n1 = n1, n2 = n2,
-      x1 = x1, x2 = x2,
-      diff_par = diff_est
+      n1 = cmh$n1, n2 = cmh$n2,
+      x1 = cmh$x1, x2 = cmh$x2,
+      diff_par = cmh$diff_est
     )$var_est
-    se_diff <- sqrt(sum(wt_normalized^2 * var_est))
+    se_diff <- sqrt(sum(cmh$wt_normalized^2 * var_est))
   }
 
   list(
     prop = estimate,
     prop_ci = estimate_ci,
-    diff = diff_est,
+    diff = cmh$diff_est,
     diff_ci = diff_ci,
     se_diff = se_diff,
-    weights = wt_normalized,
-    n1 = n1,
-    n2 = n2
+    weights = cmh$wt_normalized,
+    n1 = cmh$n1,
+    n2 = cmh$n2
   )
 }