Add lsmean_se/lsmean_ci/lsmean_diffci stats, weights_emmeans param,

NAMESPACE

@@ -265,6 +265,7 @@ export(range_noinf)
 export(reapply_varlabels)
 export(ref_group_position)
 export(rtable2gg)
+export(s_ancova)
 export(s_bland_altman)
 export(s_compare)
 export(s_count_occurrences)

NEWS.md

@@ -2,11 +2,14 @@
 
 ### Miscellaneous
 * Updated `roxygen2` to 8.0.0 and added `@exportS3Method` tags for S3 methods in `decorate_grob.R` and `utils_grid.R`.
+* Converted  `s_ancova()`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 `lsmean_se`, `lsmean_ci`, and `lsmean_diffci` statistics to `s_ancova()`.
+* Added `weights_emmeans` parameter on `s_ancova()`.
 
 ### Bug Fixes
 * Fixed bug in `prop_diff_cmh()` which previously failed when strata combinations had 0 observations.

R/summarize_ancova.R

@@ -101,16 +101,20 @@ h_ancova <- function(.var,
 #'   of the investigated linear model.
 #'
 #' @return
-#' * `s_ancova()` returns a named list of 5 statistics:
+#' * `s_ancova()` returns a named list of 8 statistics:
 #'   * `n`: Count of complete sample size for the group.
 #'   * `lsmean`: Estimated marginal means in the group.
+#'   * `lsmean_se`: Adjusted mean with standard error as a 2-element vector `c(emmean, SE)`.
+#'   * `lsmean_ci`: Adjusted mean with confidence interval as a 3-element vector `c(emmean, lower.CL, upper.CL)`.
 #'   * `lsmean_diff`: Difference in estimated marginal means in comparison to the reference group.
 #'     If working with the reference group, this will be empty.
 #'   * `lsmean_diff_ci`: Confidence level for difference in estimated marginal means in comparison
 #'     to the reference group.
+#'   * `lsmean_diffci`: Difference in adjusted means with confidence interval as a 3-element vector
+#'     `c(estimate, lower.CL, upper.CL)`. If working with the reference group, this will be `NA`.
 #'   * `pval`: p-value (not adjusted for multiple comparisons).
 #'
-#' @keywords internal
+#' @export
 s_ancova <- function(df,
                      .var,
                      .df_row,
@@ -150,63 +154,98 @@ s_ancova <- function(df,
     # convert characters selected in interaction_y into the numeric order
     interaction_y <- which(sum_fit_level[[interaction_item]] == interaction_y)
     sum_fit_level <- sum_fit_level[interaction_y, ]
-    # if interaction is called, reset the index
-    ref_key <- seq(sum_fit[[arm]][unique(.ref_group[[arm]])])
-    ref_key <- tail(ref_key, n = 1)
-    ref_key <- (interaction_y - 1) * length(unique(.df_row[[arm]])) + ref_key
   } else {
     y <- df[[.var]]
-    # Get the index of the ref arm when interaction is not called
-    ref_key <- seq(sum_fit[[arm]][unique(.ref_group[[arm]])])
-    ref_key <- tail(ref_key, n = 1)
   }
 
   if (.in_ref_col) {
     list(
       n = length(y[!is.na(y)]),
       lsmean = formatters::with_label(sum_fit_level$emmean, "Adjusted Mean"),
+      lsmean_se = formatters::with_label(
+        c(sum_fit_level$emmean, sum_fit_level$SE),
+        "Adjusted Mean (SE)"
+      ),
+      lsmean_ci = formatters::with_label(
+        c(sum_fit_level$emmean, sum_fit_level$lower.CL, sum_fit_level$upper.CL),
+        f_conf_level(conf_level)
+      ),
       lsmean_diff = formatters::with_label(numeric(), "Difference in Adjusted Means"),
       lsmean_diff_ci = formatters::with_label(numeric(), f_conf_level(conf_level)),
+      lsmean_diffci = formatters::with_label(
+        c(NA_real_, NA_real_, NA_real_),
+        "Difference in Adjusted Means (95% CI)"
+      ),
       pval = formatters::with_label(numeric(), "p-value")
     )
   } else {
-    # Estimate the differences between the marginal means.
-    emmeans_contrasts <- emmeans::contrast(
-      emmeans_fit,
-      # Compare all arms versus the control arm.
-      method = "trt.vs.ctrl",
-      # Take the arm factor from .ref_group as the control arm.
-      ref = ref_key,
-      level = conf_level
-    )
-    sum_contrasts <- summary(
-      emmeans_contrasts,
-      # Derive confidence intervals, t-tests and p-values.
-      infer = TRUE,
-      # Do not adjust the p-values for multiplicity.
-      adjust = "none"
-    )
+    if (!is.null(.ref_group)) {
+      # Compute ref_key — index of the reference arm in the emmeans grid
+      if (interaction_y != FALSE) {
+        ref_key <- seq(sum_fit[[arm]][unique(.ref_group[[arm]])])
+        ref_key <- tail(ref_key, n = 1)
+        ref_key <- (interaction_y - 1) * length(unique(.df_row[[arm]])) + ref_key
+      } else {
+        ref_key <- seq(sum_fit[[arm]][unique(.ref_group[[arm]])])
+        ref_key <- tail(ref_key, n = 1)
+      }
+      # Estimate the differences between the marginal means.
+      emmeans_contrasts <- emmeans::contrast(
+        emmeans_fit,
+        # Compare all arms versus the control arm.
+        method = "trt.vs.ctrl",
+        # Take the arm factor from .ref_group as the control arm.
+        ref = ref_key,
+        level = conf_level
+      )
+      sum_contrasts <- summary(
+        emmeans_contrasts,
+        # Derive confidence intervals, t-tests and p-values.
+        infer = TRUE,
+        # Do not adjust the p-values for multiplicity.
+        adjust = "none"
+      )
 
-    contrast_lvls <- gsub(
-      "^\\(|\\)$", "", gsub(paste0(" - \\(*", .ref_group[[arm]][1], ".*"), "", sum_contrasts$contrast)
-    )
-    if (!is.null(interaction_item)) {
-      sum_contrasts_level <- sum_contrasts[grepl(sum_level, contrast_lvls, fixed = TRUE), ]
+      contrast_lvls <- gsub(
+        "^\\(|\\)$", "", gsub(paste0(" - \\(*", .ref_group[[arm]][1], ".*"), "", sum_contrasts$contrast)
+      )
+      if (!is.null(interaction_item)) {
+        sum_contrasts_level <- sum_contrasts[grepl(sum_level, contrast_lvls, fixed = TRUE), ]
+      } else {
+        sum_contrasts_level <- sum_contrasts[sum_level == contrast_lvls, ]
+      }
+      if (interaction_y != FALSE) {
+        sum_contrasts_level <- sum_contrasts_level[interaction_y, ]
+      }
     } else {
-      sum_contrasts_level <- sum_contrasts[sum_level == contrast_lvls, ]
-    }
-    if (interaction_y != FALSE) {
-      sum_contrasts_level <- sum_contrasts_level[interaction_y, ]
+      sum_contrasts_level <- list(
+        estimate = NA_real_,
+        lower.CL = NA_real_,
+        upper.CL = NA_real_,
+        p.value = NA_real_
+      )
     }
 
     list(
       n = length(y[!is.na(y)]),
       lsmean = formatters::with_label(sum_fit_level$emmean, "Adjusted Mean"),
+      lsmean_se = formatters::with_label(
+        c(sum_fit_level$emmean, sum_fit_level$SE),
+        "Adjusted Mean (SE)"
+      ),
+      lsmean_ci = formatters::with_label(
+        c(sum_fit_level$emmean, sum_fit_level$lower.CL, sum_fit_level$upper.CL),
+        f_conf_level(conf_level)
+      ),
       lsmean_diff = formatters::with_label(sum_contrasts_level$estimate, "Difference in Adjusted Means"),
       lsmean_diff_ci = formatters::with_label(
         c(sum_contrasts_level$lower.CL, sum_contrasts_level$upper.CL),
         f_conf_level(conf_level)
       ),
+      lsmean_diffci = formatters::with_label(
+        c(sum_contrasts_level$estimate, sum_contrasts_level$lower.CL, sum_contrasts_level$upper.CL),
+        "Difference in Adjusted Means (95% CI)"
+      ),
       pval = formatters::with_label(sum_contrasts_level$p.value, "p-value")
     )
   }

R/survival_coxph_pairwise.R

@@ -112,37 +112,37 @@ s_coxph_pairwise <- function(df,
 
   pval <- if (pval_method == "wald") {
     if (alternative == "two.sided") {
-        sum_cox$waldtest["pvalue"]
+      sum_cox$waldtest["pvalue"]
     } else {
-        # Need to calculate the signed Wald statistic.
-        beta_est <- unname(cox_fit$coefficients)
-        beta_se <- sqrt(cox_fit$var[1, 1])
-        signed_wald_stat <- beta_est / beta_se
-        if (alternative == "less") {
-          stats::pnorm(signed_wald_stat)
-        } else {
-          stats::pnorm(signed_wald_stat, lower.tail = FALSE)
-        }
-    }    
+      # Need to calculate the signed Wald statistic.
+      beta_est <- unname(cox_fit$coefficients)
+      beta_se <- sqrt(cox_fit$var[1, 1])
+      signed_wald_stat <- beta_est / beta_se
+      if (alternative == "less") {
+        stats::pnorm(signed_wald_stat)
+      } else {
+        stats::pnorm(signed_wald_stat, lower.tail = FALSE)
+      }
+    }
   } else if (pval_method == "log-rank") {
     if (alternative == "two.sided") {
-        original_survdiff$pvalue
+      original_survdiff$pvalue
     } else {
-        # Need to calculate the signed log-rank statistic, 
-        # which is not included in the original survdiff output.
-        otmp <- rowSums(original_survdiff$obs)
-        etmp <- rowSums(original_survdiff$exp)
-        signed_lr_stat <- (otmp[2] - etmp[2]) / sqrt(original_survdiff$var[2, 2])
+      # Need to calculate the signed log-rank statistic,
+      # which is not included in the original survdiff output.
+      otmp <- rowSums(original_survdiff$obs)
+      etmp <- rowSums(original_survdiff$exp)
+      signed_lr_stat <- (otmp[2] - etmp[2]) / sqrt(original_survdiff$var[2, 2])
 
-        if (alternative == "less") {
-          stats::pnorm(signed_lr_stat)
-        } else {
-          stats::pnorm(signed_lr_stat, lower.tail = FALSE)
-        }
+      if (alternative == "less") {
+        stats::pnorm(signed_lr_stat)
+      } else {
+        stats::pnorm(signed_lr_stat, lower.tail = FALSE)
+      }
     }
   } else if (pval_method == "likelihood") {
     if (alternative != "two.sided") {
-        stop("Likelihood ratio test does not support one-sided alternatives")
+      stop("Likelihood ratio test does not support one-sided alternatives")
     }
     sum_cox$logtest["pvalue"]
   } else {

R/utils_default_stats_formats_labels.R

@@ -580,7 +580,8 @@ tern_default_stats <- list(
   estimate_odds_ratio = c("or_ci", "n_tot"),
   estimate_proportion = c("n_prop", "prop_ci"),
   estimate_proportion_diff = c("diff", "diff_ci"),
-  summarize_ancova = c("n", "lsmean", "lsmean_diff", "lsmean_diff_ci", "pval"),
+  summarize_ancova = c("n", "lsmean", "lsmean_se", "lsmean_ci",
+                       "lsmean_diff", "lsmean_diff_ci", "lsmean_diffci", "pval"),
   summarize_coxreg = c("n", "hr", "ci", "pval", "pval_inter"),
   summarize_glm_count = c("n", "rate", "rate_ci", "rate_ratio", "rate_ratio_ci", "pval"),
   summarize_num_patients = c("unique", "nonunique", "unique_count"),
@@ -626,8 +627,11 @@ tern_default_formats <- c(
   lr_stat_df = "xx.xx (xx.)",
   iqr = "xx.x",
   lsmean = "xx.xx",
+  lsmean_ci = "xx.xx (xx.xx - xx.xx)",
   lsmean_diff = "xx.xx",
   lsmean_diff_ci = "(xx.xx, xx.xx)",
+  lsmean_diffci = "xx.xx (xx.xx - xx.xx)",
+  lsmean_se = "xx.xx (xx.xx)",
   mad = "xx.x",
   max = "xx.x",
   mean = "xx.x",

tests/testthat/_snaps/summarize_ancova.md

@@ -23,6 +23,16 @@
       attr(,"label")
       [1] "Adjusted Mean"
 
+      $lsmean_se
+      [1] 5.71740917 0.06680849
+      attr(,"label")
+      [1] "Adjusted Mean (SE)"
+
+      $lsmean_ci
+      [1] 5.717409 5.542996 5.891823
+      attr(,"label")
+      [1] "99% CI"
+
       $lsmean_diff
       [1] -0.4374138
       attr(,"label")
@@ -33,6 +43,11 @@
       attr(,"label")
       [1] "99% CI"
 
+      $lsmean_diffci
+      [1] -0.4374138 -1.4268150  0.5519873
+      attr(,"label")
+      [1] "Difference in Adjusted Means (95% CI)"
+
       $pval
       [1] 0.2503574
       attr(,"label")
@@ -52,6 +67,16 @@
       attr(,"label")
       [1] "Adjusted Mean"
 
+      $lsmean_se
+      [1] 4.36264203 0.08598191
+      attr(,"label")
+      [1] "Adjusted Mean (SE)"
+
+      $lsmean_ci
+      [1] 4.362642 4.192651 4.532633
+      attr(,"label")
+      [1] "95% CI"
+
       $lsmean_diff
       numeric(0)
       attr(,"label")
@@ -62,6 +87,11 @@
       attr(,"label")
       [1] "95% CI"
 
+      $lsmean_diffci
+      [1] NA NA NA
+      attr(,"label")
+      [1] "Difference in Adjusted Means (95% CI)"
+
       $pval
       numeric(0)
       attr(,"label")