Support for multiarm ancova

NAMESPACE

@@ -46,6 +46,7 @@ export(Stack)
 export(add_class)
 export(analyse)
 export(ancova)
+export(ancova_m_groups)
 export(as_class)
 export(as_vcov)
 export(control_bayes)

NEWS.md

@@ -4,6 +4,8 @@
 
 * All covariance structures are now also supported for Bayesian multiple imputation: `method_bayes()` gained additional `covariance` and `prior_cov` arguments to allow users to specify the covariance structure and prior for the Bayesian imputation model. Please see the updated statistical specifications vignette for details. (#501, #518)
 * New function `mcse()` to calculate the Monte Carlo standard error for pooled estimates from (approximate) Bayesian imputation. (#493)
+* Added support for multi-group ANCOVA analysis with `ancova_m_groups()` function. This enables analysis of covariance with more than two treatment groups, where each non-reference group is compared against the reference group. (#520)
+
 
 # rbmi 1.4.1
 

R/analyse.R

@@ -58,6 +58,11 @@
 #' a linear transformation of the outcomes.
 #' Thus care is required when applying alternative analysis functions in this setting.
 #'
+#' Note that as of verion 1.5.0 rbmi provides the `ancova_m_groups()` function; this is an extension to
+#' the `ancova()` function which supports the analysis of more than 2 treatment arms.
+#' The `ancova()` function is left as the default for backwards compatibility as `ancova_m_groups()`
+#' uses a different naming structure.
+#'
 #' The `delta` argument can be used to specify offsets to be applied
 #' to the outcome variable in the imputed datasets prior to the analysis.
 #' This is typically used for sensitivity or tipping point analyses. The

R/ancova.R

@@ -51,6 +51,7 @@
 #'
 #' @inheritSection lsmeans Weighting
 #'
+#' @seealso [ancova_m_groups()] for multi-group analysis
 #' @seealso [analyse()]
 #' @seealso [stats::lm()]
 #' @seealso [set_vars()]
@@ -60,6 +61,120 @@ ancova <- function(
     vars,
     visits = NULL,
     weights = c("counterfactual", "equal", "proportional_em", "proportional")
+) {
+    ancova_core(
+        data = data,
+        vars = vars,
+        visits = visits,
+        weights = weights,
+        ancova_fun = ancova_single
+    )
+}
+
+
+#' Multi-Group Analysis of Covariance
+#'
+#' Performs an analysis of covariance with support for multiple treatment groups,
+#' returning treatment effects (contrasts between each non-reference group and the
+#' reference group) and least square means estimates for all groups.
+#'
+#' @inheritParams ancova
+#'
+#' @details
+#' This function extends [ancova()] to support more than two treatment groups.
+#' The function works as follows:
+#'
+#' 1. Select the first value from `visits`.
+#' 2. Subset the data to only the observations that occurred on this visit.
+#' 3. Fit a linear model as `vars$outcome ~ vars$group + vars$covariates`.
+#' 4. Extract treatment effects (each non-reference group vs reference) and
+#'    least square means for all treatment groups.
+#' 5. Repeat for all other values in `visits`.
+#'
+#' The reference group is determined by the first factor level of `vars$group`.
+#' Only pairwise comparisons between each non-reference group and the reference
+#' group are computed (not all pairwise comparisons between groups).
+#'
+#' @section Output Structure:
+#' Results are returned as a named list with elements suffixed by visit name:
+#' ```
+#' list(
+#'     trt_L1_L0_visit_1 = list(est = ...), # Group L1 vs L0 (reference)
+#'     trt_L2_L0_visit_1 = list(est = ...), # Group L2 vs L0 (reference)
+#'     lsm_L0_visit_1 = list(est = ...),    # LSMean for reference group L0
+#'     lsm_L1_visit_1 = list(est = ...),    # LSMean for group L1
+#'     lsm_L2_visit_1 = list(est = ...),    # LSMean for group L2
+#'     ...
+#' )
+#' ```
+#' Group levels are standardized to "L0", "L1", "L2", etc., where "L0" represents
+#' the reference group (first factor level).
+#'
+#' @inheritSection lsmeans Weighting
+#'
+#' @examples
+#' \dontrun{
+#' # 3-group analysis
+#' data$treatment <- factor(data$treatment, levels = c("Control", "Low", "High"))
+#' vars <- set_vars(
+#'     outcome = "response",
+#'     group = "treatment",
+#'     visit = "visit",
+#'     covariates = c("age", "sex")
+#' )
+#' result <- ancova_m_groups(data, vars)
+#'
+#' # With group interactions
+#' vars_int <- set_vars(
+#'     outcome = "response",
+#'     group = "treatment",
+#'     visit = "visit",
+#'     covariates = c("age", "sex", "treatment*age")
+#' )
+#' result_int <- ancova_m_groups(data, vars_int)
+#' }
+#'
+#' @seealso [ancova()] for two-group analysis
+#' @seealso [analyse()]
+#' @seealso [stats::lm()]
+#' @seealso [set_vars()]
+#' @export
+ancova_m_groups <- function(
+    data,
+    vars,
+    visits = NULL,
+    weights = c("counterfactual", "equal", "proportional_em", "proportional")
+) {
+    ancova_core(
+        data = data,
+        vars = vars,
+        visits = visits,
+        weights = weights,
+        ancova_fun = ancova_single_m_group
+    )
+}
+
+
+
+#' Core ANCOVA Implementation
+#'
+#' Internal function that provides common functionality for both two-group and
+#' multi-group ANCOVA analyses. This function handles input validation, visit
+#' processing, and delegates the actual analysis to the specified ancova function.
+#'
+#' @inheritParams ancova
+#' @param ancova_fun Function to perform the single-visit ANCOVA analysis.
+#'   Should be either [ancova_single()] for two-group analysis or
+#'   [ancova_single_m_group()] for multi-group analysis.
+#'
+#'
+#' @keywords internal
+ancova_core <- function(
+    data,
+    vars,
+    visits = NULL,
+    weights = c("counterfactual", "equal", "proportional_em", "proportional"),
+    ancova_fun
 ) {
     outcome <- vars[["outcome"]]
     group <- vars[["group"]]
@@ -130,7 +245,7 @@ ancova <- function(
         visits,
         function(x) {
             data2 <- data[data[[visit]] == x, ]
-            res <- ancova_single(data2, outcome, group, covariates, weights)
+            res <- ancova_fun(data2, outcome, group, covariates, weights)
             names(res) <- paste0(names(res), "_", x)
             return(res)
         }
@@ -164,6 +279,7 @@ ancova <- function(
 #' }
 #' @seealso [ancova()]
 #' @importFrom stats lm coef vcov df.residual
+#' @keywords internal
 ancova_single <- function(
     data,
     outcome,
@@ -206,3 +322,110 @@ ancova_single <- function(
     )
     return(x)
 }
+
+
+
+#' Multi-Group ANCOVA Implementation for Single Visit
+#'
+#' Internal function that implements analysis of covariance for multiple treatment
+#' groups within a single visit. This function extends [ancova_single()] to handle
+#' more than two groups.
+#'
+#' @inheritParams ancova_single
+#'
+#' @details
+#' This function:
+#' - Accepts factor variables with 2 or more levels for the group parameter
+#' - Standardizes group level names to "L0", "L1", "L2", etc. for consistent output
+#' - Calculates treatment effects as contrasts between each non-reference group and the reference
+#' - Computes least square means for all groups using the specified weighting strategy
+#'
+#' The reference group is always the first factor level (standardized to "L0").
+#' Treatment effects compare each subsequent group ("L1", "L2", ...) against "L0".
+#'
+#' @section Reserved Variable Names:
+#' This function uses `"rbmiGroup"` as an internal variable name. If your dataset
+#' contains a variable with this name, the function will error.
+#'
+#' @return A named list containing:
+#' \itemize{
+#'   \item `trt_L1_L0`, `trt_L2_L0`, etc.: Treatment effects (each group vs reference)
+#'   \item `lsm_L0`, `lsm_L1`, `lsm_L2`, etc.: Least square means for all groups
+#' }
+#' Each element contains `est` (estimate), `se` (standard error), and `df` (degrees of freedom).
+#'
+#' @examples
+#' \dontrun{
+#' # Internal function - typically called via ancova_m_groups()
+#' data$grp <- factor(c("Control", "Treatment1", "Treatment2"))
+#' result <- ancova_single_m_group(
+#'     data = data,
+#'     outcome = "response",
+#'     group = "grp",
+#'     covariates = c("age", "sex"),
+#'     weights = "counterfactual"
+#' )
+#' }
+#'
+#' @seealso [ancova_single()] for two-group analysis
+#' @seealso [ancova_m_groups()] for the user-facing multi-group function
+#' @importFrom stats lm coef vcov df.residual
+#' @keywords internal
+ancova_single_m_group <- function(
+    data,
+    outcome,
+    group,
+    covariates,
+    weights = c("counterfactual", "equal", "proportional_em", "proportional")
+) {
+    weights <- match.arg(weights)
+    assert_that(
+        is.factor(data[[group]]),
+        length(levels(data[[group]])) >= 2,
+        length(unique(data[[group]])) == length(levels(data[[group]])),
+        msg = "`data[[vars$group]]` must be a factor variable with >=2 levels"
+    )
+
+    data2 <- data[, c(extract_covariates(covariates), outcome, group)]
+
+    # Standardise level names and group variable name to make extraction reliable
+    levels(data2[[group]]) <- paste0("L", seq_along(levels(data2[[group]])) - 1)
+    assert_that(
+        !"rbmiGroup" %in% names(data),
+        msg = paste(
+            "`rbmiGroup` is a reserved variable name for internal use, please rename",
+            "your variable to avoid conflicts"
+        )
+    )
+    data2[["rbmiGroup"]] <- data2[[group]]
+    data2 <- data2[!names(data2) %in% group]
+
+    frm <- as_simple_formula(outcome, c(covariates, group))
+    frm_std <- frm_find_and_replace(frm, as.name(group), as.name("rbmiGroup"))
+
+    mod <- lm(formula = frm_std, data = data2)
+    args <- list(model = mod, .weights = weights)
+
+    lsm <- lapply(
+        levels(data2[["rbmiGroup"]]),
+        function(x) {
+            args[["rbmiGroup"]] <- x
+            do.call(lsmeans, args)
+        }
+    )
+    names(lsm) <- sprintf("lsm_%s", levels(data2[["rbmiGroup"]]))
+
+    trt <- lapply(
+        levels(data2[["rbmiGroup"]])[-1],
+        function(x) {
+            term <- sprintf("rbmiGroup%s", x)
+            list(
+                est = coef(mod)[[term]],
+                se = sqrt(vcov(mod)[term, term]),
+                df = df.residual(mod)
+            )
+        }
+    )
+    names(trt) <- sprintf("trt_%s_L0", levels(data2[["rbmiGroup"]])[-1])
+    append(trt, lsm)
+}

R/utilities.R

@@ -912,3 +912,33 @@ set_options <- function() {
         }
     }
 }
+
+
+#' Recursively find and replace symbols in a language object
+#'
+#' This function traverses a language object (such as an expression or call)
+#' and recursively replaces all occurrences of a specified symbol with another symbol.
+#'
+#' @param frm A language object (e.g., call, expression, or list of calls) to search and modify.
+#' @param find_sym A symbol (as a name) to find within \code{frm}.
+#' @param replace_sym A symbol (as a name) to replace \code{find_sym} with.
+#'
+#' @return The modified language object with all instances of \code{find_sym} replaced by \code{replace_sym}.
+#'
+#' @examples
+#' \dontrun{
+#' expr <- quote(a + b * c)
+#' frm_find_and_replace(expr, as.name("b"), as.name("x"))
+#' # Returns: a + x * c
+#' }
+#' @keywords internal
+frm_find_and_replace <- function(frm, find_sym, replace_sym) {
+    for (i in seq_along(frm)) {
+        if (is.call(frm[[i]])) {
+            frm[[i]] <- frm_find_and_replace(frm[[i]], find_sym, replace_sym)
+        } else if (is.name(frm[[i]])) {
+            frm[[i]] <- ifelse(frm[[i]] == find_sym, replace_sym, frm[[i]])
+        }
+    }
+    frm
+}