Feature eair: Extra combined statistic required to be able to present "n (EAIR per 100 SY)"

.github/workflows/check_CRAN.yaml

@@ -53,6 +53,7 @@ jobs:
           needs: check
 
       - uses: r-lib/actions/check-r-package@v2
+        timeout-minutes: 10
         with:
           upload-snapshots: true
           build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")'

DESCRIPTION

@@ -72,12 +72,15 @@ Suggests:
     pharmaverseadamjnj (>= 0.0.2),
     multcomp,
     lifecycle,
-    dplyr
+    dplyr,
+    ratesci (>= 1.1.0)
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 Remotes:
     insightsengineering/formatters@main,
     insightsengineering/tern@main,
+    petelaud/ratesci@master,
+    johnsonandjohnson/pharmaverseadamjnj@dev,
     insightsengineering/rtables@main,
     insightsengineering/rtables.officer@main
 Config/roxygen2/version: 8.0.0

NAMESPACE

@@ -7,6 +7,7 @@ export(a_cmhrms_j_with_exclude)
 export(a_coxph_hr)
 export(a_diff_means)
 export(a_eair100_j)
+export(a_eair_strat_j)
 export(a_event_free)
 export(a_freq_combos_j)
 export(a_freq_j)
@@ -160,6 +161,7 @@ importFrom(formatters,wrap_string_ttype)
 importFrom(generics,tidy)
 importFrom(rtables,in_rows)
 importFrom(rtables,rcell)
+importFrom(stats,ave)
 importFrom(stats,complete.cases)
 importFrom(stats,setNames)
 importFrom(stats,t.test)

NEWS.md

@@ -23,6 +23,7 @@
 - Reduce sampling of `rbmi` test to make tests shorter (#323)
 - Optimize shared tables in `test-tt_to_tblfile` (#323)
 - Changed forked `h_ancova` for the `tern` one
+- Add extra statistics to `a_eair100_j` and introduce scaling factor `num_p_year` (default = 100) (#361)
 - Removed ellipsis argument from `a_freq_resp_var_j` (#236) 
 - Removed `dplyr` from `junco` and replaced by base R #201
 - Replaced `assertthat` by `checkmate` for consistency #201
@@ -36,6 +37,7 @@
 - Use new exported splv_extra and value_expr accessors (insightsengineering#1098)
 - Update new exported calls from rtables.officer
 - update documentation to `roxygen2` 8.0.0 
+- Add extra statistics to `a_eair100_j` and introduce scaling factor `num_p_year` (default = 100) (#361)
 
 ### Added
 - Added `categorize_pval()` for assigning p-values to validated, user-defined categories.

R/a_eair_diff.R

@@ -0,0 +1,142 @@
+.h_s_eair_diff_scalar <- function(x1,
+                                  n1,
+                                  x2,
+                                  n2,
+                                  conf_type = c("wald", "waldcc", "mn", "scas"),
+                                  conf_level = 0.95) {
+  conf_type <- match.arg(conf_type)
+  diff_est <- switch(conf_type,
+    waldcc = h_s_eair_diff_wald(x1, n1, x2, n2, conf_level, cc = TRUE),
+    wald   = h_s_eair_diff_wald(x1, n1, x2, n2, conf_level, cc = FALSE),
+    mn     = h_s_eair_diff_mn(x1, n1, x2, n2, conf_level, skew = FALSE),
+    scas   = h_s_eair_diff_mn(x1, n1, x2, n2, conf_level, skew = TRUE)
+  )
+  list(diff_est = diff_est$est, diff_lcl = diff_est$ci[1], diff_ucl = diff_est$ci[2])
+}
+
+h_s_eair_diff <- function(x1,
+                          n1,
+                          x2,
+                          n2,
+                          conf_type,
+                          conf_level = 0.95) {
+  results <- mapply(
+    .h_s_eair_diff_scalar,
+    x1, n1, x2, n2,
+    MoreArgs = list(conf_type = conf_type, conf_level = conf_level),
+    SIMPLIFY = FALSE
+  )
+  data.frame(
+    diff_est = sapply(results, `[[`, "diff_est"),
+    diff_lcl = sapply(results, `[[`, "diff_lcl"),
+    diff_ucl = sapply(results, `[[`, "diff_ucl")
+  )
+}
+
+h_s_eair_diff_wald <- function(x1, n1, x2, n2, conf_level, cc = TRUE) {
+  alpha <- 1 - conf_level
+  est <- x1 / n1 - x2 / n2
+  se <- sqrt(x1 / n1^2 + x2 / n2^2)
+
+  coeff <- stats::qnorm(1 - alpha / 2)
+  term2 <- coeff * se
+  if (cc) {
+    term2 <- term2 + 0.5 * (1 / n1 + 1 / n2)
+  }
+  lcl <- est - term2
+  ucl <- est + term2
+
+  eair_diff <- list()
+  eair_diff$est <- est
+  eair_diff$ci <- c(lcl, ucl)
+  eair_diff
+}
+
+h_s_eair_diff_mn <- function(x1, n1, x2, n2, conf_level, skew = FALSE) {
+  alpha <- 1 - conf_level
+  z <- stats::qnorm(1 - alpha / 2)
+
+  est <- x1 / n1 - x2 / n2
+
+  lcl <- .conf_poi_rd(x1, n1, x2, n2, z, lower = TRUE, skew = skew, level = conf_level)
+  ucl <- .conf_poi_rd(x1, n1, x2, n2, z, lower = FALSE, skew = skew, level = conf_level)
+
+  eair_diff <- list()
+  eair_diff$est <- est
+  eair_diff$ci <- c(lcl, ucl)
+  eair_diff
+}
+
+.score_poi_rd <- function(x1, n1, x2, n2, theta, skew = FALSE, level = 0.95) {
+  alpha <- 1 - level
+  p1hat <- x1 / n1
+  p2hat <- x2 / n2
+  Stheta <- (p1hat - p2hat) - theta
+  if (isTRUE(abs(Stheta) == 0)) {
+    return(0)
+  }
+  x <- x1 + x2
+  N <- n1 + n2
+  A <- N
+  B <- N * theta - x
+  C_ <- -x2 * theta
+  num <- (-B + sqrt(B^2 - 4 * A * C_))
+  p2d <- ifelse(isTRUE(num == 0), 0, num / (2 * A))
+  p1d <- p2d + theta
+  V <- p1d / n1 + p2d / n2
+  score1 <- Stheta / sqrt(V)
+
+  if (!skew) {
+    return(score1)
+  }
+
+  # Skewness correction (SCAS)
+  mu3 <- p1d / (n1^2) - p2d / (n2^2)
+  scterm <- mu3 / (6 * V^(3 / 2))
+  if (isTRUE(scterm == 0)) {
+    return(score1)
+  }
+
+  A_q <- scterm
+  B_q <- 1
+  C_q <- -(score1 + scterm)
+  dsct <- B_q^2 - 4 * A_q * C_q
+  if (dsct < 0) {
+    # Fallback to simplified skewness correction
+    qtnorm <- stats::qnorm(1 - alpha / 2)
+    score1 <- score1 - (qtnorm^2 - 1) * scterm
+  } else {
+    score1 <- (-B_q + sqrt(dsct)) / (2 * A_q)
+  }
+
+  return(score1)
+}
+
+.conf_poi_rd <- function(x1, n1, x2, n2, z, lower = FALSE,
+                         skew = FALSE, level = 0.95) {
+  p1 <- x1 / n1
+  p2 <- x2 / n2
+  p_hat <- p1 - p2
+  dp <- 1 + ifelse(lower, 1, -1) * p_hat
+  i <- 1
+  while (i <= 50) {
+    dp <- 0.5 * dp
+    y <- p_hat + ifelse(lower, -1, 1) * dp
+    score <- .score_poi_rd(x1, n1, x2, n2, y, skew = skew, level = level)
+    # For lower bound: score decreases as y moves left, find where score = +z
+    #   -> move y further left if score > z (still too close to estimate)
+    # For upper bound: score decreases as y moves right, find where score = -z
+    #   -> move y further right if score > -z (still too close to estimate)
+    if (lower) {
+      if (score < z) p_hat <- y
+    } else {
+      if (score > -z) p_hat <- y
+    }
+    if ((dp < 1e-07) || (abs(score - ifelse(lower, z, -z)) < 1e-06)) {
+      break
+    } else {
+      i <- i + 1
+    }
+  }
+  return(y)
+}