diff --git a/NAMESPACE b/NAMESPACE index 38e6e6e1..41acd84e 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index 142d8c05..cb1814f4 100644 --- a/NEWS.md +++ b/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 diff --git a/R/analyse.R b/R/analyse.R index b765fb4d..c347cda7 100644 --- a/R/analyse.R +++ b/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 diff --git a/R/ancova.R b/R/ancova.R index 0306f5b0..8a520138 100644 --- a/R/ancova.R +++ b/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) +} diff --git a/R/utilities.R b/R/utilities.R index 36759285..34217441 100644 --- a/R/utilities.R +++ b/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 +} diff --git a/man/analyse.Rd b/man/analyse.Rd index b69a5440..5eac6d21 100644 --- a/man/analyse.Rd +++ b/man/analyse.Rd @@ -89,6 +89,11 @@ method (\code{method = method_condmean()} in \code{\link[=draws]{draws()}}) reli 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 \code{ancova_m_groups()} function; this is an extension to +the \code{ancova()} function which supports the analysis of more than 2 treatment arms. +The \code{ancova()} function is left as the default for backwards compatibility as \code{ancova_m_groups()} +uses a different naming structure. + The \code{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 diff --git a/man/ancova.Rd b/man/ancova.Rd index dd5926e2..2ef6d4e9 100644 --- a/man/ancova.Rd +++ b/man/ancova.Rd @@ -120,6 +120,8 @@ for \code{weights = "counterfactual"}. } \seealso{ +\code{\link[=ancova_m_groups]{ancova_m_groups()}} for multi-group analysis + \code{\link[=analyse]{analyse()}} \code{\link[stats:lm]{stats::lm()}} diff --git a/man/ancova_core.Rd b/man/ancova_core.Rd new file mode 100644 index 00000000..227d7b3a --- /dev/null +++ b/man/ancova_core.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ancova.R +\name{ancova_core} +\alias{ancova_core} +\title{Core ANCOVA Implementation} +\usage{ +ancova_core( + data, + vars, + visits = NULL, + weights = c("counterfactual", "equal", "proportional_em", "proportional"), + ancova_fun +) +} +\arguments{ +\item{data}{A \code{data.frame} containing the data to be used in the model.} + +\item{vars}{A \code{vars} object as generated by \code{\link[=set_vars]{set_vars()}}. Only the \code{group}, +\code{visit}, \code{outcome} and \code{covariates} elements are required. See details.} + +\item{visits}{An optional character vector specifying which visits to +fit the ancova model at. If \code{NULL}, a separate ancova model will be fit to the +outcomes for each visit (as determined by \code{unique(data[[vars$visit]])}). +See details.} + +\item{weights}{Character, either \code{"counterfactual"} (default), \code{"equal"}, +\code{"proportional_em"} or \code{"proportional"}. +Specifies the weighting strategy to be used when calculating the lsmeans. +See the weighting section for more details.} + +\item{ancova_fun}{Function to perform the single-visit ANCOVA analysis. +Should be either \code{\link[=ancova_single]{ancova_single()}} for two-group analysis or +\code{\link[=ancova_single_m_group]{ancova_single_m_group()}} for multi-group analysis.} +} +\description{ +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. +} +\keyword{internal} diff --git a/man/ancova_m_groups.Rd b/man/ancova_m_groups.Rd new file mode 100644 index 00000000..d6006d94 --- /dev/null +++ b/man/ancova_m_groups.Rd @@ -0,0 +1,152 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ancova.R +\name{ancova_m_groups} +\alias{ancova_m_groups} +\title{Multi-Group Analysis of Covariance} +\usage{ +ancova_m_groups( + data, + vars, + visits = NULL, + weights = c("counterfactual", "equal", "proportional_em", "proportional") +) +} +\arguments{ +\item{data}{A \code{data.frame} containing the data to be used in the model.} + +\item{vars}{A \code{vars} object as generated by \code{\link[=set_vars]{set_vars()}}. Only the \code{group}, +\code{visit}, \code{outcome} and \code{covariates} elements are required. See details.} + +\item{visits}{An optional character vector specifying which visits to +fit the ancova model at. If \code{NULL}, a separate ancova model will be fit to the +outcomes for each visit (as determined by \code{unique(data[[vars$visit]])}). +See details.} + +\item{weights}{Character, either \code{"counterfactual"} (default), \code{"equal"}, +\code{"proportional_em"} or \code{"proportional"}. +Specifies the weighting strategy to be used when calculating the lsmeans. +See the weighting section for more details.} +} +\description{ +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. +} +\details{ +This function extends \code{\link[=ancova]{ancova()}} to support more than two treatment groups. +The function works as follows: +\enumerate{ +\item Select the first value from \code{visits}. +\item Subset the data to only the observations that occurred on this visit. +\item Fit a linear model as \code{vars$outcome ~ vars$group + vars$covariates}. +\item Extract treatment effects (each non-reference group vs reference) and +least square means for all treatment groups. +\item Repeat for all other values in \code{visits}. +} + +The reference group is determined by the first factor level of \code{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: + +\if{html}{\out{
}}\preformatted{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 + ... +) +}\if{html}{\out{
}} + +Group levels are standardized to "L0", "L1", "L2", etc., where "L0" represents +the reference group (first factor level). +} + +\section{Weighting}{ + +\subsection{Counterfactual}{ + +For \code{weights = "counterfactual"} (the default) the lsmeans are obtained by +taking the average of the predicted values for each patient after assigning all patients +to each arm in turn. +This approach is equivalent to standardization or g-computation. +In comparison to \code{emmeans} this approach is equivalent to: + +\if{html}{\out{
}}\preformatted{emmeans::emmeans(model, specs = "", counterfactual = "") +}\if{html}{\out{
}} + +Note that to ensure backwards compatibility with previous versions of \code{rbmi} +\code{weights = "proportional"} is an alias for \code{weights = "counterfactual"}. +To get results consistent with \code{emmeans}'s \code{weights = "proportional"} +please use \code{weights = "proportional_em"}. +} + +\subsection{Equal}{ + +For \code{weights = "equal"} the lsmeans are obtained by taking the model fitted +value of a hypothetical patient whose covariates are defined as follows: +\itemize{ +\item Continuous covariates are set to \code{mean(X)} +\item Dummy categorical variables are set to \code{1/N} where \code{N} is the number of levels +\item Continuous * continuous interactions are set to \code{mean(X) * mean(Y)} +\item Continuous * categorical interactions are set to \code{mean(X) * 1/N} +\item Dummy categorical * categorical interactions are set to \code{1/N * 1/M} +} + +In comparison to \code{emmeans} this approach is equivalent to: + +\if{html}{\out{
}}\preformatted{emmeans::emmeans(model, specs = "", weights = "equal") +}\if{html}{\out{
}} +} + +\subsection{Proportional}{ + +For \code{weights = "proportional_em"} the lsmeans are obtained as per \code{weights = "equal"} +except instead of weighting each observation equally they are weighted by the proportion +in which the given combination of categorical values occurred in the data. +In comparison to \code{emmeans} this approach is equivalent to: + +\if{html}{\out{
}}\preformatted{emmeans::emmeans(model, specs = "", weights = "proportional") +}\if{html}{\out{
}} + +Note that this is not to be confused with \code{weights = "proportional"} which is an alias +for \code{weights = "counterfactual"}. +} +} + +\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{ +\code{\link[=ancova]{ancova()}} for two-group analysis + +\code{\link[=analyse]{analyse()}} + +\code{\link[stats:lm]{stats::lm()}} + +\code{\link[=set_vars]{set_vars()}} +} diff --git a/man/ancova_single.Rd b/man/ancova_single.Rd index 61cecb29..377b37d5 100644 --- a/man/ancova_single.Rd +++ b/man/ancova_single.Rd @@ -98,3 +98,4 @@ ancova_single(iris2, "Sepal.Length", "Species", c("Petal.Length * Petal.Width")) \seealso{ \code{\link[=ancova]{ancova()}} } +\keyword{internal} diff --git a/man/ancova_single_m_group.Rd b/man/ancova_single_m_group.Rd new file mode 100644 index 00000000..98789b59 --- /dev/null +++ b/man/ancova_single_m_group.Rd @@ -0,0 +1,80 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ancova.R +\name{ancova_single_m_group} +\alias{ancova_single_m_group} +\title{Multi-Group ANCOVA Implementation for Single Visit} +\usage{ +ancova_single_m_group( + data, + outcome, + group, + covariates, + weights = c("counterfactual", "equal", "proportional_em", "proportional") +) +} +\arguments{ +\item{data}{A \code{data.frame} containing the data to be used in the model.} + +\item{outcome}{Character, the name of the outcome variable in \code{data}.} + +\item{group}{Character, the name of the group variable in \code{data}.} + +\item{covariates}{Character vector containing the name of any additional covariates +to be included in the model as well as any interaction terms.} + +\item{weights}{Character, either \code{"counterfactual"} (default), \code{"equal"}, +\code{"proportional_em"} or \code{"proportional"}. +Specifies the weighting strategy to be used when calculating the lsmeans. +See the weighting section for more details.} +} +\value{ +A named list containing: +\itemize{ +\item \code{trt_L1_L0}, \code{trt_L2_L0}, etc.: Treatment effects (each group vs reference) +\item \code{lsm_L0}, \code{lsm_L1}, \code{lsm_L2}, etc.: Least square means for all groups +} +Each element contains \code{est} (estimate), \code{se} (standard error), and \code{df} (degrees of freedom). +} +\description{ +Internal function that implements analysis of covariance for multiple treatment +groups within a single visit. This function extends \code{\link[=ancova_single]{ancova_single()}} to handle +more than two groups. +} +\details{ +This function: +\itemize{ +\item Accepts factor variables with 2 or more levels for the group parameter +\item Standardizes group level names to "L0", "L1", "L2", etc. for consistent output +\item Calculates treatment effects as contrasts between each non-reference group and the reference +\item 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 \code{"rbmiGroup"} as an internal variable name. If your dataset +contains a variable with this name, the function will error. +} + +\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{ +\code{\link[=ancova_single]{ancova_single()}} for two-group analysis + +\code{\link[=ancova_m_groups]{ancova_m_groups()}} for the user-facing multi-group function +} +\keyword{internal} diff --git a/man/frm_find_and_replace.Rd b/man/frm_find_and_replace.Rd new file mode 100644 index 00000000..29df03c3 --- /dev/null +++ b/man/frm_find_and_replace.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utilities.R +\name{frm_find_and_replace} +\alias{frm_find_and_replace} +\title{Recursively find and replace symbols in a language object} +\usage{ +frm_find_and_replace(frm, find_sym, replace_sym) +} +\arguments{ +\item{frm}{A language object (e.g., call, expression, or list of calls) to search and modify.} + +\item{find_sym}{A symbol (as a name) to find within \code{frm}.} + +\item{replace_sym}{A symbol (as a name) to replace \code{find_sym} with.} +} +\value{ +The modified language object with all instances of \code{find_sym} replaced by \code{replace_sym}. +} +\description{ +This function traverses a language object (such as an expression or call) +and recursively replaces all occurrences of a specified symbol with another symbol. +} +\examples{ +\dontrun{ +expr <- quote(a + b * c) +frm_find_and_replace(expr, as.name("b"), as.name("x")) +# Returns: a + x * c +} +} +\keyword{internal} diff --git a/tests/testthat/test-ancova.R b/tests/testthat/test-ancova.R index 5226131f..8e9340be 100644 --- a/tests/testthat/test-ancova.R +++ b/tests/testthat/test-ancova.R @@ -236,3 +236,385 @@ test_that("ancova", { regex = "`k`" ) }) + + +test_that("ancova_m_groups - basic 3-group functionality", { + set.seed(201) + + n <- 900 + dat <- tibble( + visit = "vis1", + age = rnorm(n), + sex = factor(sample(c("M", "F"), size = n, replace = TRUE)), + grp = factor(sample(c("Control", "Trt1", "Trt2"), size = n, replace = TRUE)), + out = rnorm(n, mean = 50 + 2 * age + 3 * f2n(sex) + + 4 * (grp == "Trt1") + 6 * (grp == "Trt2"), sd = 10) + ) + + # Manual model for comparison + mod <- lm(out ~ age + sex + grp, data = dat) + + # Test ancova_m_groups + result_actual <- ancova_m_groups( + dat, + list( + outcome = "out", + group = "grp", + covariates = c("age", "sex"), + visit = "visit" + ) + ) + + # Check output structure + expected_names <- c( + "trt_L1_L0_vis1", "trt_L2_L0_vis1", # Treatment effects vs reference + "lsm_L0_vis1", "lsm_L1_vis1", "lsm_L2_vis1" # LSMeans for all groups + ) + expect_equal(sort(names(result_actual)), sort(expected_names)) + + + expected_trt_L1_L0_vis1 <- list( + est = coef(mod)[["grpTrt1"]], + se = sqrt(vcov(mod)["grpTrt1", "grpTrt1"]), + df = df.residual(mod) + ) + expect_equal(result_actual$trt_L1_L0_vis1, expected_trt_L1_L0_vis1) + + + expected_trt_L2_L0_vis1 <- list( + est = coef(mod)[["grpTrt2"]], + se = sqrt(vcov(mod)["grpTrt2", "grpTrt2"]), + df = df.residual(mod) + ) + expect_equal(result_actual$trt_L2_L0_vis1, expected_trt_L2_L0_vis1) +}) + + +test_that("ancova_m_groups - 4-group functionality", { + set.seed(301) + + n <- 800 + grp_levels <- c("Placebo", "Low", "Medium", "High") + dat <- tibble( + visit = "baseline", + baseline_score = rnorm(n, mean = 30, sd = 8), + grp = factor(sample(grp_levels, size = n, replace = TRUE), levels = grp_levels), + out = rnorm(n, mean = 20 + 0.7 * baseline_score + + 2 * (grp == "Low") + + 5 * (grp == "Medium") + + 8 * (grp == "High"), sd = 6) + ) + + # Manual model for comparison + mod <- lm(out ~ baseline_score + grp, data = dat) + + result_actual <- ancova_m_groups( + dat, + list( + outcome = "out", + group = "grp", + covariates = "baseline_score", + visit = "visit" + ) + ) + + # Check we have all expected components + expected_names <- c( + "trt_L1_L0_baseline", "trt_L2_L0_baseline", "trt_L3_L0_baseline", # 3 treatment effects + "lsm_L0_baseline", "lsm_L1_baseline", "lsm_L2_baseline", "lsm_L3_baseline" # 4 LSMeans + ) + expect_equal(sort(names(result_actual)), sort(expected_names)) + + # Verify treatment effects + expect_equal(result_actual$trt_L1_L0_baseline$est, coef(mod)[["grpLow"]], tolerance = 1e-10) + expect_equal(result_actual$trt_L2_L0_baseline$est, coef(mod)[["grpMedium"]], tolerance = 1e-10) + expect_equal(result_actual$trt_L3_L0_baseline$est, coef(mod)[["grpHigh"]], tolerance = 1e-10) +}) + + +test_that("ancova_m_groups - LSMeans with group interactions", { + set.seed(401) + + grp_levels <- c("Control", "Treatment1", "Treatment2") + sex_levels <- c("Male", "Female") + n <- 600 + dat <- tibble( + visit = "week12", + age = rnorm(n, mean = 45, sd = 12), + sex = factor(sample(sex_levels, size = n, replace = TRUE), levels = sex_levels), + grp = factor(sample(grp_levels, size = n, replace = TRUE), levels = grp_levels), + # Create interaction effect between group and sex + out = rnorm(n, mean = 40 + 0.5 * age + + 3 * (sex == "Female") + + 5 * (grp == "Treatment1") + + 7 * (grp == "Treatment2") + + # Group-sex interactions + 2 * (grp == "Treatment1" & sex == "Female") + + -1 * (grp == "Treatment2" & sex == "Female"), + sd = 8) + ) + + # Test with group*sex interaction in covariates + result_actual <- ancova_m_groups( + dat, + list( + outcome = "out", + group = "grp", + covariates = c("age", "sex", "grp*sex"), # Note: using rbmiGroup for interaction + visit = "visit" + ) + ) + + mod <- lm(out ~ age + sex + grp + grp:sex, data = dat) + + + # Check structure + expected_names <- c( + "trt_L1_L0_week12", "trt_L2_L0_week12", + "lsm_L0_week12", "lsm_L1_week12", "lsm_L2_week12" + ) + expect_equal(sort(names(result_actual)), sort(expected_names)) + + + suppressMessages({ + emmean_counter <- as.data.frame( + emmeans::emmeans(mod, "grp", counterfactual = "grp") + ) + }) + + expected <- as.list(emmean_counter[emmean_counter$grp == "Control", c("emmean", "SE", "df")]) + names(expected) <- c("est", "se", "df") + expect_equal(result_actual$lsm_L0_week12, expected) + + expected <- as.list(emmean_counter[emmean_counter$grp == "Treatment1", c("emmean", "SE", "df")]) + names(expected) <- c("est", "se", "df") + expect_equal(result_actual$lsm_L1_week12, expected) + + expected <- as.list(emmean_counter[emmean_counter$grp == "Treatment2", c("emmean", "SE", "df")]) + names(expected) <- c("est", "se", "df") + expect_equal(result_actual$lsm_L2_week12, expected) + + + # Same again but with a different weighting scheme + result_actual <- ancova_m_groups( + dat, + list( + outcome = "out", + group = "grp", + covariates = c("age", "sex", "grp*sex"), # Note: using rbmiGroup for interaction + visit = "visit" + ), + weights = "equal" + ) + suppressMessages({ + emmean_counter <- as.data.frame( + emmeans::emmeans(mod, "grp", "grp") + ) + }) + expected <- as.list(emmean_counter[emmean_counter$grp == "Control", c("emmean", "SE", "df")]) + names(expected) <- c("est", "se", "df") + expect_equal(result_actual$lsm_L0_week12, expected) + + expected <- as.list(emmean_counter[emmean_counter$grp == "Treatment1", c("emmean", "SE", "df")]) + names(expected) <- c("est", "se", "df") + expect_equal(result_actual$lsm_L1_week12, expected) + + expected <- as.list(emmean_counter[emmean_counter$grp == "Treatment2", c("emmean", "SE", "df")]) + names(expected) <- c("est", "se", "df") + expect_equal(result_actual$lsm_L2_week12, expected) + +}) + + +test_that("ancova_m_groups - multiple visits", { + set.seed(601) + + n_per_visit <- 200 + grp_levels <- c("A", "B", "C") + dat <- bind_rows( + tibble( + visit = "visit1", + age = rnorm(n_per_visit), + grp = factor(sample(grp_levels, size = n_per_visit, replace = TRUE), levels = grp_levels), + out = rnorm(n_per_visit, mean = 20 + 2 * age + + 3 * (grp == "B") + 5 * (grp == "C"), sd = 4) + ), + tibble( + visit = "visit2", + age = rnorm(n_per_visit), + grp = factor(sample(c("A", "B", "C"), size = n_per_visit, replace = TRUE)), + out = rnorm(n_per_visit, mean = 25 + 2 * age + + 4 * (grp == "B") + 7 * (grp == "C"), sd = 4) + ) + ) + + result_actual <- ancova_m_groups( + dat, + list( + outcome = "out", + group = "grp", + covariates = "age", + visit = "visit" + ) + ) + + # Check we have results for both visits + visit1_names <- grep("_visit1$", names(result_actual), value = TRUE) + visit2_names <- grep("_visit2$", names(result_actual), value = TRUE) + + expect_length(visit1_names, 5) # 2 treatment effects + 3 LSMeans + expect_length(visit2_names, 5) # 2 treatment effects + 3 LSMeans + + expected_names <- c( + "trt_L1_L0_visit1", "trt_L2_L0_visit1", "lsm_L0_visit1", "lsm_L1_visit1", "lsm_L2_visit1", + "trt_L1_L0_visit2", "trt_L2_L0_visit2", "lsm_L0_visit2", "lsm_L1_visit2", "lsm_L2_visit2" + ) + expect_true(all(expected_names %in% names(result_actual))) +}) + + +test_that("ancova_m_groups - no covariates", { + set.seed(701) + + n <- 300 + grp_levels <- c("Control", "Low", "High") + dat <- tibble( + visit = "final", + grp = factor(sample(grp_levels, size = n, replace = TRUE), levels = grp_levels), + out = rnorm(n, mean = 40 + 3 * (grp == "Low") + 8 * (grp == "High"), sd = 6) + ) + + result_actual <- ancova_m_groups( + dat, + list( + outcome = "out", + group = "grp", + covariates = NULL, + visit = "visit" + ) + ) + + # Manual model for comparison + mod <- lm(out ~ grp, data = dat) + + # Check treatment effects + expect_equal( + result_actual$trt_L1_L0_final$est, + coef(mod)[["grpLow"]], + tolerance = 1e-10 + ) + expect_equal( + result_actual$trt_L2_L0_final$est, + coef(mod)[["grpHigh"]], + tolerance = 1e-10 + ) + + # Check we have all expected components + expected_names <- c( + "trt_L1_L0_final", "trt_L2_L0_final", + "lsm_L0_final", "lsm_L1_final", "lsm_L2_final" + ) + expect_equal(sort(names(result_actual)), sort(expected_names)) +}) + + +test_that("ancova_m_groups - error handling", { + set.seed(801) + + n <- 100 + dat <- tibble( + visit = "v1", + age = rnorm(n), + grp = factor(sample(c("A", "B"), size = n, replace = TRUE)), # Only 2 groups + out = rnorm(n) + ) + + # Should work with 2 groups (minimum requirement) + expect_no_error({ + result <- ancova_m_groups( + dat, + list(outcome = "out", group = "grp", covariates = "age", visit = "visit") + ) + }) + + # Test with single group - should fail + dat_single_group <- dat + dat_single_group$grp <- factor(rep("A", n)) + + expect_error( + ancova_m_groups( + dat_single_group, + list(outcome = "out", group = "grp", covariates = "age", visit = "visit") + ), + regexp = "must be a factor variable with >=2 levels" + ) + + # Test reserved variable name conflict + dat_conflict <- dat + dat_conflict$rbmiGroup <- 1 + + expect_error( + ancova_m_groups( + dat_conflict, + list(outcome = "out", group = "grp", covariates = c("age", "rbmiGroup"), visit = "visit") + ), + regexp = "rbmiGroup.*reserved variable name" + ) +}) + + +test_that("ancova_m_groups - backward compatibility with 2-group case", { + set.seed(901) + + n <- 400 + grp_levels <- c("Control", "Treatment") + sex_levels <- c("M", "F") + dat <- tibble( + visit = "week8", + age = rnorm(n), + sex = factor(sample(sex_levels, size = n, replace = TRUE), levels = sex_levels), + grp = factor(sample(grp_levels, size = n, replace = TRUE), levels = grp_levels), + out = rnorm(n, mean = 30 + 2 * age + 3 * f2n(sex) + 5 * f2n(grp), sd = 7) + ) + + # Compare ancova_m_groups with traditional ancova for 2-group case + result_traditional <- ancova( + dat, + list( + outcome = "out", + group = "grp", + covariates = c("age", "sex"), + visit = "visit" + ) + ) + + result_multigroup <- ancova_m_groups( + dat, + list( + outcome = "out", + group = "grp", + covariates = c("age", "sex"), + visit = "visit" + ) + ) + + expect_equal( + result_traditional$trt_week8, + result_multigroup$trt_L1_L0_week8, + tolerance = 1e-10 + ) + + expect_equal( + result_traditional$lsm_ref_week8, + result_multigroup$lsm_L0_week8, + tolerance = 1e-10 + ) + + expect_equal( + result_traditional$lsm_alt_week8, + result_multigroup$lsm_L1_week8, + tolerance = 1e-10 + ) + + expect_true(length(result_traditional) == length(result_multigroup)) +}) diff --git a/tests/testthat/test-utilities.R b/tests/testthat/test-utilities.R index b9a932fd..098e36ae 100644 --- a/tests/testthat/test-utilities.R +++ b/tests/testthat/test-utilities.R @@ -342,3 +342,33 @@ test_that("get_stan_model works as expected depending on covariance and prior on model <- expect_silent(get_stan_model("us", "lkj")) expect_snapshot(model) }) + + + + +test_that("frm_find_and_replace works as expected", { + # Things being tested here include: + # - Can change values on both sides of the formula + # - Interactions e.g. `*` and `:` + # - No arg functions k() + # - 1-arg functions h(z) + # - 2-arg functions f(z, z) + # - Constants e.g. + 1 + # - Data modification functions I(z^2) + # - Name subset e.g. zz doesn't get renamed (where we are searching for z) + frm <- x + z ~ 1 + z + a + b + zz + z:a + a * z + a * b + h(z) + f(z, z) + k() + I(z^2) + actual <- frm_find_and_replace(frm, as.name("z"), as.name("P")) + expected <- x + P ~ 1 + P + a + b + zz + P:a + a * P + a * b + h(P) + f(P, P) + k() + I(P^2) + environment(actual) <- globalenv() + environment(expected) <- globalenv() + expect_equal(actual, expected) + + + # Special names / special characters + frm <- ~ ` .. !abc & `:x - 1 * x + actual <- frm_find_and_replace(frm, as.name(" .. !abc & "), as.name("bob")) + expected <- ~ bob:x - 1 * x + environment(actual) <- globalenv() + environment(expected) <- globalenv() + expect_equal(actual, expected) +}) diff --git a/vignettes/FAQ.Rmd b/vignettes/FAQ.Rmd index df8681e7..7a460599 100644 --- a/vignettes/FAQ.Rmd +++ b/vignettes/FAQ.Rmd @@ -123,3 +123,24 @@ anaObj <- rbmi::analyse( ```
+ + + +## Is ANCOVA on more than 2 treatment groups supported? + +This can be implemented by using the `ancova_m_groups()` function instead of the default `ancova()` function e.g. + +```r +anaObj <- rbmi::analyse( + imputeObj, + ancova_m_groups, + vars = vars + ) +``` + +Please note that the object created by this function uses a different naming convention to +that of the `ancova()` function which means any manual extraction code looking for specific +names such as `"lsm_alt"` will need to be updated. Please see the documentation for `ancova_m_groups()` +for more details. + +
diff --git a/vignettes/FAQ.html b/vignettes/FAQ.html index 82dec069..70de0b49 100644 --- a/vignettes/FAQ.html +++ b/vignettes/FAQ.html @@ -362,6 +362,7 @@

Alessandro Noci, Craig Gower-Page and Marcel Wolbers

  • 1.4 How to handle missing data in baseline covariates in rbmi?
  • 1.5 Why does rbmi by default use an ANCOVA analysis model and not an MMRM analysis model?
  • 1.6 How can I analyse the change-from-baseline in the analysis model when imputation was done on the original outcomes?
    • +
  • 1.7 Is ANCOVA on more than 2 treatment groups supported?
    • @@ -447,6 +448,20 @@

    1.6 How can I analyse the change- vars = vars )


    + +
    +

    1.7 Is ANCOVA on more than 2 treatment groups supported?

    +

    This can be implemented by using the ancova_m_groups() function instead of the default ancova() function e.g.

    +
    anaObj <- rbmi::analyse(
+    imputeObj,
+    ancova_m_groups,
+    vars = vars
+ )
    +

    Please note that the object created by this function uses a different naming convention to +that of the ancova() function which means any manual extraction code looking for specific +names such as "lsm_alt" will need to be updated. Please see the documentation for ancova_m_groups() +for more details.

    +


    White, Ian R, and Simon G Thompson. 2005. “Adjusting for Partially Missing Baseline Measurements in Randomized Trials.” Statistics in Medicine 24 (7): 993–1007. diff --git a/vignettes/quickstart.Rmd b/vignettes/quickstart.Rmd index 20e52a4e..8cf4f989 100644 --- a/vignettes/quickstart.Rmd +++ b/vignettes/quickstart.Rmd @@ -229,6 +229,17 @@ function which determines the names of the key variables within the data and the Please also note that the names of the analysis estimates contain "ref" and "alt" to refer to the two treatment arms. In particular "ref" refers to the first factor level of `vars$group` which does not necessarily coincide with the control arm. In this example, since `levels(dat[[vars$group]]) = c("DRUG", PLACEBO`), the results associated with "ref" correspond to the intervention arm, while those associated with "alt" correspond to the control arm. +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. That is, instead of `"ref"` and "`alt`" the values `"L0"`, `"L1"`, `"L2"`, ... +are used to represent the different arms where `"L0"` is the reference level and `"L1"` is the first +offset group (e.g. the second factor level). That is the resulting list will contain entries such as: +`lsm_L1_day40` which contains the least square mean for the second factor level at the `"day40"` +visit. For treatment effects it will contain entries such as `trt_L1_L0_day40` which contains +the treatment effect for the second factor level compared to the reference level at the `"day40"` +visit. See the documentation for `ancova_m_groups()` for more details. + Additionally, we can use the `delta` argument of `analyse()` to perform a delta adjustments of the imputed datasets prior to the analysis. In brief, this is implemented by specifying a data.frame that contains the amount of adjustment to be added to each longitudinal outcome for each subject and visit, i.e. diff --git a/vignettes/quickstart.html b/vignettes/quickstart.html index b9b89f9d..d1174679 100644 --- a/vignettes/quickstart.html +++ b/vignettes/quickstart.html @@ -360,17 +360,9 @@

    2 The Data

    We use a publicly available example dataset from an antidepressant clinical trial of an active drug versus placebo. The relevant endpoint is the Hamilton 17-item depression rating scale (HAMD17) which was assessed at baseline and at weeks 1, 2, 4, and 6. Study drug discontinuation occurred in 24% of subjects from the active drug and 26% of subjects from placebo. All data after study drug discontinuation are missing and there is a single additional intermittent missing observation.

    library(rbmi)
 library(dplyr)
-#> 
-#> Attaching package: 'dplyr'
-#> The following objects are masked from 'package:stats':
-#> 
-#>     filter, lag
-#> The following objects are masked from 'package:base':
-#> 
-#>     intersect, setdiff, setequal, union
-
-data("antidepressant_data")
-dat <- antidepressant_data
    + +data("antidepressant_data") +dat <- antidepressant_data

    We consider an imputation model with the mean change from baseline in the HAMD17 score as the outcome (variable CHANGE in the dataset). The following covariates are included in the imputation model: the treatment group (THERAPY), the (categorical) visit (VISIT), treatment-by-visit interactions, the baseline HAMD17 score (BASVAL), and baseline HAMD17 score-by-visit interactions. A common unstructured covariance matrix structure is assumed for both groups. The analysis model is an ANCOVA model with the treatment group as the primary factor and adjustment for the baseline HAMD17 score.

    rbmi expects its input dataset to be complete; that is, there must be one row per subject for each visit. Missing outcome values should be coded as NA, while missing covariate values are not allowed. If the dataset is incomplete, then the expand_locf() helper function can be used to add any missing rows, using LOCF imputation to carry forward the observed baseline covariate values to visits with missing outcomes. Rows corresponding to missing outcomes are not present in the antidepressant trial dataset. To address this we will therefore use the expand_locf() function as follows:

    
@@ -469,14 +461,16 @@ 
    3 Draws
    
 #> Imputation Type: random
 #> Method:
 #>     name: Bayes
-#>     same_cov: TRUE
-#>     n_samples: 150
-#> Controls:
-#>     warmup: 200
-#>     thin: 5
-#>     chains: 1
-#>     init: mmrm
-#>     seed: 214647346
    +#> covariance: us +#> same_cov: TRUE +#> n_samples: 150 +#> prior_cov: default +#> Controls: +#> warmup: 200 +#> thin: 5 +#> chains: 1 +#> init: mmrm +#> seed: 1074096309

    Note the use of set_vars() which specifies the names of the key variables within the dataset and the imputation model. Additionally, note that whilst vars$group and vars$visit are added as terms to the imputation model by default, their interaction is not, @@ -603,6 +597,16 @@

    5 Analyse

    (in addition to the treatment group) for which the analysis model will be adjusted.

    Please also note that the names of the analysis estimates contain “ref” and “alt” to refer to the two treatment arms. In particular “ref” refers to the first factor level of vars$group which does not necessarily coincide with the control arm. In this example, since levels(dat[[vars$group]]) = c("DRUG", PLACEBO), the results associated with “ref” correspond to the intervention arm, while those associated with “alt” correspond to the control arm.

    +

    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. That is, instead of "ref" and “alt” the values "L0", "L1", "L2", … +are used to represent the different arms where "L0" is the reference level and "L1" is the first +offset group (e.g. the second factor level). That is the resulting list will contain entries such as: +lsm_L1_day40 which contains the least square mean for the second factor level at the "day40" +visit. For treatment effects it will contain entries such as trt_L1_L0_day40 which contains +the treatment effect for the second factor level compared to the reference level at the "day40" +visit. See the documentation for ancova_m_groups() for more details.

    Additionally, we can use the delta argument of analyse() to perform a delta adjustments of the imputed datasets prior to the analysis. In brief, this is implemented by specifying a data.frame that contains the amount of adjustment to be added to each longitudinal outcome for each subject and visit, i.e.  @@ -701,34 +705,34 @@

    6 Pool

    #> trt_4 -0.092 0.683 -1.439 1.256 0.893 #> lsm_ref_4 -1.616 0.486 -2.576 -0.656 0.001 #> lsm_alt_4 -1.708 0.475 -2.645 -0.77 <0.001 -#> trt_5 1.335 0.923 -0.488 3.158 0.15 -#> lsm_ref_5 -4.154 0.659 -5.455 -2.853 <0.001 -#> lsm_alt_5 -2.819 0.646 -4.095 -1.542 <0.001 -#> trt_6 1.954 0.998 -0.016 3.925 0.052 -#> lsm_ref_6 -6.087 0.717 -7.503 -4.671 <0.001 -#> lsm_alt_6 -4.133 0.703 -5.521 -2.745 <0.001 -#> trt_7 2.178 1.129 -0.053 4.41 0.056 -#> lsm_ref_7 -6.974 0.821 -8.597 -5.351 <0.001 -#> lsm_alt_7 -4.796 0.784 -6.345 -3.246 <0.001 +#> trt_5 1.332 0.922 -0.489 3.154 0.151 +#> lsm_ref_5 -4.146 0.658 -5.446 -2.845 <0.001 +#> lsm_alt_5 -2.813 0.644 -4.085 -1.542 <0.001 +#> trt_6 1.946 0.995 -0.02 3.912 0.052 +#> lsm_ref_6 -6.097 0.713 -7.505 -4.689 <0.001 +#> lsm_alt_6 -4.151 0.7 -5.534 -2.768 <0.001 +#> trt_7 2.172 1.118 -0.038 4.382 0.054 +#> lsm_ref_7 -6.987 0.817 -8.603 -5.372 <0.001 +#> lsm_alt_7 -4.816 0.782 -6.362 -3.269 <0.001 #> --------------------------------------------------

    The table of values shown in the print message for poolObj can also be extracted using the as.data.frame() function:

    as.data.frame(poolObj)
-#>    parameter         est        se        lci        uci         pval
-#> 1      trt_4 -0.09180645 0.6826279 -1.4394968  1.2558839 8.931772e-01
-#> 2  lsm_ref_4 -1.61581996 0.4862316 -2.5757714 -0.6558685 1.093708e-03
-#> 3  lsm_alt_4 -1.70762640 0.4749573 -2.6453193 -0.7699335 4.262148e-04
-#> 4      trt_5  1.33531976 0.9230204 -0.4876244  3.1582639 1.499512e-01
-#> 5  lsm_ref_5 -4.15415338 0.6588023 -5.4553174 -2.8529894 2.730742e-09
-#> 6  lsm_alt_5 -2.81883362 0.6463331 -4.0954650 -1.5422022 2.331688e-05
-#> 7      trt_6  1.95436277 0.9976289 -0.0164307  3.9251562 5.191672e-02
-#> 8  lsm_ref_6 -6.08732723 0.7166468 -7.5032662 -4.6713882 1.747812e-14
-#> 9  lsm_alt_6 -4.13296446 0.7025598 -5.5211647 -2.7447642 2.534885e-08
-#> 10     trt_7  2.17814646 1.1288440 -0.0532519  4.4095448 5.564659e-02
-#> 11 lsm_ref_7 -6.97364970 0.8207107 -8.5966654 -5.3506340 3.049429e-14
-#> 12 lsm_alt_7 -4.79550324 0.7838188 -6.3448217 -3.2461848 8.537927e-09
    +#> parameter est se lci uci pval +#> 1 trt_4 -0.09180645 0.6826279 -1.43949684 1.2558839 8.931772e-01 +#> 2 lsm_ref_4 -1.61581996 0.4862316 -2.57577141 -0.6558685 1.093708e-03 +#> 3 lsm_alt_4 -1.70762640 0.4749573 -2.64531931 -0.7699335 4.262148e-04 +#> 4 trt_5 1.33244564 0.9223725 -0.48916744 3.1540587 1.505325e-01 +#> 5 lsm_ref_5 -4.14553184 0.6584403 -5.44594687 -2.8451168 2.851010e-09 +#> 6 lsm_alt_5 -2.81308620 0.6437591 -4.08452555 -1.5416468 2.238611e-05 +#> 7 trt_6 1.94579706 0.9951403 -0.02003552 3.9116296 5.235157e-02 +#> 8 lsm_ref_6 -6.09703116 0.7126425 -7.50494478 -4.6891175 1.158175e-14 +#> 9 lsm_alt_6 -4.15123411 0.6998827 -5.53407714 -2.7683911 1.976048e-08 +#> 10 trt_7 2.17178546 1.1181330 -0.03812541 4.3816963 5.403151e-02 +#> 11 lsm_ref_7 -6.98740689 0.8171118 -8.60323818 -5.3715756 2.188701e-14 +#> 12 lsm_alt_7 -4.81562143 0.7823681 -6.36209275 -3.2691501 7.119562e-09

    These outputs gives an estimated difference of -2.178 (95% CI -0.053 to 4.410) -between the two groups at the last visit with an associated p-value of 0.056.

    +2.172 (95% CI -0.038 to 4.382) +between the two groups at the last visit with an associated p-value of 0.054.

    7 Code