diff --git a/R/get_marginal_effect.R b/R/get_marginal_effect.R
index 89a785e..d395158 100644
--- a/R/get_marginal_effect.R
+++ b/R/get_marginal_effect.R
@@ -35,10 +35,11 @@
 #'
 #' @param contrast a string indicating choice of contrast. Defaults to 'diff' for a risk difference. See \link[beeca]{apply_contrast}.
 #'
-#' @param reference a string indicating which treatment group should be considered as
-#' the reference level. Accepted values are one of the levels in the treatment
-#' variable. Default to the first level used in the `glm` object. This parameter influences the calculation of treatment effects
-#' relative to the chosen reference group.
+#' @param reference a string or list of strings indicating which treatment
+#' group(s) to use as reference level for pairwise comparisons. Accepted values
+#' must be a subset of the levels in the treatment variable. Default to the
+#' first n-1 treatment levels used in the `glm` object. This parameter influences
+#' the calculation of treatment effects relative to the chosen reference group.
 #'
 #' @param mod for Ye's method, the implementation of open-source RobinCar package
 #' has an additional variance decomposition step when estimating the robust variance,
diff --git a/R/predict_counterfactuals.R b/R/predict_counterfactuals.R
index 64eca24..a12be17 100644
--- a/R/predict_counterfactuals.R
+++ b/R/predict_counterfactuals.R
@@ -17,19 +17,18 @@
 #' @return an updated `glm` object appended with an
 #' additional component `counterfactual.predictions`.
 #'
-#' This component contains a tibble with two columns: `cf_pred_0` and `cf_pred_1`,
-#' representing counterfactual predictions for each level of the
-#' treatment variable. A descriptive `label` attribute explains the counterfactual
-#' scenario associated with each column.
+#' This component contains a tibble with columns representing counterfactual
+#' predictions for each level of the treatment variable. A descriptive `label`
+#' attribute explains the counterfactual scenario associated with each column.
 #'
 #' @importFrom stats predict
 #' @importFrom dplyr as_tibble
 #' @export
 #'
 #' @details
-#' The function works by creating two new datasets from the original data used
-#' to fit the GLM model. In these datasets, the treatment variable
-#' is set to each of its levels across all records (e.g., patients).
+#' The function works by creating new datasets from the original data used
+#' to fit the GLM model. In these datasets, the treatment variable for all
+#' records (e.g., patients) is set to each possible treatment level.
 #'
 #' Predictions are then made for each dataset based on the fitted GLM model,
 #' simulating the response variable under each treatment condition.
diff --git a/README.md b/README.md
index b2e2ee3..735b995 100644
--- a/README.md
+++ b/README.md
@@ -24,7 +24,7 @@ Motivated by the recent [FDA guidance (2023)](https://www.fda.gov/regulatory-inf
 
 ## Scope
 
-The package is designed to estimate marginal (unconditional) estimands in a binary endpoint setting with covariate adjustment. It is suited for 2-arm clinical trials with or without covariate adaptive (stratified permuted block or biased coin) randomization where the summary measure of the marginal estimand is one of (risk difference, odds ratio, risk ratio, log odds ratio, log risk ratio). For practical considerations on the implications of covariate adjustment in superiority vs non-inferiority trials, please see [Nicholas et al. (2015)](https://doi.org/10.1002/sim.6447) and [Morris et al. (2022)](https://doi.org/10.1186/s13063-022-06097-z).
+The package is designed to estimate marginal (unconditional) estimands in a binary endpoint setting with covariate adjustment. It is suited for clinical trials with or without covariate adaptive (stratified permuted block or biased coin) randomization where the summary measure of the marginal estimand is one of (risk difference, odds ratio, risk ratio, log odds ratio, log risk ratio). For practical considerations on the implications of covariate adjustment in superiority vs non-inferiority trials, please see [Nicholas et al. (2015)](https://doi.org/10.1002/sim.6447) and [Morris et al. (2022)](https://doi.org/10.1186/s13063-022-06097-z).
 
 ## Example
 
diff --git a/man/get_marginal_effect.Rd b/man/get_marginal_effect.Rd
index 481eb6a..5770627 100644
--- a/man/get_marginal_effect.Rd
+++ b/man/get_marginal_effect.Rd
@@ -40,10 +40,11 @@ model-based, HC3, HC, HC1, HC2, HC4, HC4m, and HC5. See \link[sandwich]{vcovHC}
 
 \item{contrast}{a string indicating choice of contrast. Defaults to 'diff' for a risk difference. See \link[beeca]{apply_contrast}.}
 
-\item{reference}{a string indicating which treatment group should be considered as
-the reference level. Accepted values are one of the levels in the treatment
-variable. Default to the first level used in the \code{glm} object. This parameter influences the calculation of treatment effects
-relative to the chosen reference group.}
+\item{reference}{a string or list of strings indicating which treatment
+group(s) to use as reference level for pairwise comparisons. Accepted values
+must be a subset of the levels in the treatment variable. Default to the
+first n-1 treatment levels used in the \code{glm} object. This parameter influences
+the calculation of treatment effects relative to the chosen reference group.}
 
 \item{mod}{for Ye's method, the implementation of open-source RobinCar package
 has an additional variance decomposition step when estimating the robust variance,
diff --git a/man/predict_counterfactuals.Rd b/man/predict_counterfactuals.Rd
index 9853084..41d050c 100644
--- a/man/predict_counterfactuals.Rd
+++ b/man/predict_counterfactuals.Rd
@@ -18,10 +18,9 @@ in fitting the \code{object}.}
 an updated \code{glm} object appended with an
 additional component \code{counterfactual.predictions}.
 
-This component contains a tibble with two columns: \code{cf_pred_0} and \code{cf_pred_1},
-representing counterfactual predictions for each level of the
-treatment variable. A descriptive \code{label} attribute explains the counterfactual
-scenario associated with each column.
+This component contains a tibble with columns representing counterfactual
+predictions for each level of the treatment variable. A descriptive \code{label}
+attribute explains the counterfactual scenario associated with each column.
 }
 \description{
 This function calculates counterfactual predictions for each level of a
@@ -30,9 +29,9 @@ It is designed to aid in the assessment of treatment effects by predicting
 outcomes under different treatments under causal inference framework.
 }
 \details{
-The function works by creating two new datasets from the original data used
-to fit the GLM model. In these datasets, the treatment variable
-is set to each of its levels across all records (e.g., patients).
+The function works by creating new datasets from the original data used
+to fit the GLM model. In these datasets, the treatment variable for all
+records (e.g., patients) is set to each possible treatment level.
 
 Predictions are then made for each dataset based on the fitted GLM model,
 simulating the response variable under each treatment condition.
diff --git a/vignettes/estimand_and_implementations.Rmd b/vignettes/estimand_and_implementations.Rmd
index d22fb1e..a86fd95 100644
--- a/vignettes/estimand_and_implementations.Rmd
+++ b/vignettes/estimand_and_implementations.Rmd
@@ -56,11 +56,10 @@ See [Magirr et al. (2024)](https://osf.io/9mp58/) for discussion of the estimand
 
 ## Analysis example
 
-In the following example, we will use an example CDISC clinical trial dataset in ADaM format (`?trial02_cdisc`).
-We will only include the two treatment arms of interest for the analysis: "Xanomeline High Dose" and "Placebo".
-Below we show how to describe and perform an adjusted analysis targeting marginal risk difference between the two arms.  
+In the following example, we will use an example CDISC clinical trial dataset in ADaM format (`?trial02_cdisc`) with three treatment arms.
+Below we show how to describe and perform an adjusted analysis targeting marginal risk difference comparing the two active arms against placebo.  
 
-A logistic regression model will be used to estimate the average treatment effect between Xanomeline High Dose and Placebo. The model will adjust for baseline values of patient sex, race and age. Difference in marginal response proportions with p-value and corresponding 95% confidence interval will be estimated from the logistic regression model using the methodology described in Ge et al. 2011 with sandwich variance estimator (Liu et al. 2023). The sandwich estimator provides robustness to model-misspecification.
+A logistic regression model will be used to estimate the average treatment effect for Xanomeline High Dose vs Placebo and Xanomeline Low Dose vs Placebo. The model will adjust for baseline values of patient sex, race and age. Difference in marginal response proportions with p-value and corresponding 95% confidence interval will be estimated from the logistic regression model using the methodology described in Ge et al. 2011 with sandwich variance estimator (Liu et al. 2023). The sandwich estimator provides robustness to model-misspecification.
 
 ```{r, message=FALSE, warning=FALSE}
 library(beeca)