diff --git a/NAMESPACE b/NAMESPACE
index 07530394a7..3a32db7308 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS.md b/NEWS.md
index 8a9f28d61a..a051af4f71 100644
--- a/NEWS.md
+++ b/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.
diff --git a/R/summarize_ancova.R b/R/summarize_ancova.R
index f740c43e46..cd131749fc 100644
--- a/R/summarize_ancova.R
+++ b/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")
     )
   }
diff --git a/R/survival_coxph_pairwise.R b/R/survival_coxph_pairwise.R
index 1132646768..96dc41e0d1 100644
--- a/R/survival_coxph_pairwise.R
+++ b/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 {
diff --git a/R/utils_default_stats_formats_labels.R b/R/utils_default_stats_formats_labels.R
index 937abc8892..63f8e15776 100644
--- a/R/utils_default_stats_formats_labels.R
+++ b/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",
diff --git a/man/summarize_ancova.Rd b/man/summarize_ancova.Rd
index 2b246a9b88..5484c8d293 100644
--- a/man/summarize_ancova.Rd
+++ b/man/summarize_ancova.Rd
@@ -91,7 +91,7 @@ times, to avoid warnings from \code{rtables}.}
 
 \item{.stats}{(\code{character})\cr statistics to select for the table.
 
-Options are: \verb{'n', 'lsmean', 'lsmean_diff', 'lsmean_diff_ci', 'pval'}}
+Options are: \verb{'n', 'lsmean', 'lsmean_se', 'lsmean_ci', 'lsmean_diff', 'lsmean_diff_ci', 'lsmean_diffci', 'pval'}}
 
 \item{.stat_names}{(\code{character})\cr names of the statistics that are passed directly to name single statistics
 (\code{.stats}). This option is visible when producing \code{\link[rtables:as_result_df]{rtables::as_result_df()}} with \code{make_ard = TRUE}.}
@@ -123,14 +123,18 @@ the statistics from \code{s_ancova()} to the table layout.
 }
 
 \itemize{
-\item \code{s_ancova()} returns a named list of 5 statistics:
+\item \code{s_ancova()} returns a named list of 8 statistics:
 \itemize{
 \item \code{n}: Count of complete sample size for the group.
 \item \code{lsmean}: Estimated marginal means in the group.
+\item \code{lsmean_se}: Adjusted mean with standard error as a 2-element vector \code{c(emmean, SE)}.
+\item \code{lsmean_ci}: Adjusted mean with confidence interval as a 3-element vector \code{c(emmean, lower.CL, upper.CL)}.
 \item \code{lsmean_diff}: Difference in estimated marginal means in comparison to the reference group.
 If working with the reference group, this will be empty.
 \item \code{lsmean_diff_ci}: Confidence level for difference in estimated marginal means in comparison
 to the reference group.
+\item \code{lsmean_diffci}: Difference in adjusted means with confidence interval as a 3-element vector
+\code{c(estimate, lower.CL, upper.CL)}. If working with the reference group, this will be \code{NA}.
 \item \code{pval}: p-value (not adjusted for multiple comparisons).
 }
 }
diff --git a/tests/testthat/_snaps/summarize_ancova.md b/tests/testthat/_snaps/summarize_ancova.md
index c9f1e45fd9..cb821b2390 100644
--- a/tests/testthat/_snaps/summarize_ancova.md
+++ b/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")
diff --git a/tests/testthat/test-summarize_ancova.R b/tests/testthat/test-summarize_ancova.R
index df6bc820f4..3d14ff85c1 100644
--- a/tests/testthat/test-summarize_ancova.R
+++ b/tests/testthat/test-summarize_ancova.R
@@ -231,3 +231,89 @@ testthat::test_that("summarize_ancova works with irregular arm levels", {
   res <- testthat::expect_silent(result3)
   testthat::expect_snapshot(res)
 })
+
+testthat::test_that("s_ancova returns lsmean_se, lsmean_ci, lsmean_diffci for non-ref column", {
+  df_col <- iris %>% dplyr::filter(Species == "versicolor")
+  df_ref <- iris %>% dplyr::filter(Species == "setosa")
+  result <- s_ancova(
+    df = df_col,
+    .var = "Sepal.Length",
+    .df_row = iris,
+    variables = list(arm = "Species", covariates = "Petal.Length"),
+    .ref_group = df_ref,
+    .in_ref_col = FALSE,
+    conf_level = 0.95
+  )
+
+  testthat::expect_length(result$lsmean_se, 2)
+  testthat::expect_length(result$lsmean_ci, 3)
+  testthat::expect_length(result$lsmean_diffci, 3)
+  testthat::expect_equal(as.numeric(result$lsmean_se[[1]]), as.numeric(result$lsmean))
+  testthat::expect_equal(as.numeric(result$lsmean_diffci[[1]]), as.numeric(result$lsmean_diff))
+  testthat::expect_equal(as.numeric(result$lsmean_diffci[2:3]), as.numeric(result$lsmean_diff_ci))
+})
+
+testthat::test_that("s_ancova returns lsmean_se and lsmean_ci for ref column", {
+  df_ref <- iris %>% dplyr::filter(Species == "setosa")
+  result <- s_ancova(
+    df = df_ref,
+    .var = "Sepal.Length",
+    .df_row = iris,
+    variables = list(arm = "Species", covariates = "Petal.Length"),
+    .ref_group = df_ref,
+    .in_ref_col = TRUE,
+    conf_level = 0.95
+  )
+
+  testthat::expect_length(result$lsmean_se, 2)
+  testthat::expect_length(result$lsmean_ci, 3)
+  testthat::expect_true(all(is.na(result$lsmean_diffci)))
+  testthat::expect_length(result$lsmean_diff, 0)
+})
+
+testthat::test_that("s_ancova returns NA diffs when .ref_group is NULL", {
+  df_col <- iris %>% dplyr::filter(Species == "versicolor")
+  result <- s_ancova(
+    df = df_col,
+    .var = "Sepal.Length",
+    .df_row = iris,
+    variables = list(arm = "Species", covariates = "Petal.Length"),
+    .ref_group = NULL,
+    .in_ref_col = FALSE,
+    conf_level = 0.95
+  )
+
+  testthat::expect_true(is.na(result$lsmean_diff))
+  testthat::expect_true(all(is.na(result$lsmean_diff_ci)))
+  testthat::expect_true(all(is.na(result$lsmean_diffci)))
+  testthat::expect_true(is.na(result$pval))
+})
+
+testthat::test_that("s_ancova respects weights_emmeans parameter", {
+  df_col <- iris %>% dplyr::filter(Species == "versicolor")
+  df_ref <- iris %>% dplyr::filter(Species == "setosa")
+
+  result_null <- s_ancova(
+    df = df_col,
+    .var = "Sepal.Length",
+    .df_row = iris,
+    variables = list(arm = "Species", covariates = NULL),
+    .ref_group = df_ref,
+    .in_ref_col = FALSE,
+    conf_level = 0.95,
+    weights_emmeans = NULL
+  )
+  result_cf <- s_ancova(
+    df = df_col,
+    .var = "Sepal.Length",
+    .df_row = iris,
+    variables = list(arm = "Species", covariates = NULL),
+    .ref_group = df_ref,
+    .in_ref_col = FALSE,
+    conf_level = 0.95,
+    weights_emmeans = "counterfactual"
+  )
+
+  # With no covariates, proportional and counterfactual weights give the same emmeans
+  testthat::expect_equal(result_null$lsmean, result_cf$lsmean, tolerance = 1e-6)
+})