diff --git a/DESCRIPTION b/DESCRIPTION
index b8c49ea..8feb318 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: bregr
 Title: Easy and Efficient Batch Processing of Regression Models
-Version: 1.3.2
+Version: 1.4.0
 Authors@R: c(
     person("Shixiang", "Wang", , "w_shixiang@163.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-9855-7357")),
@@ -50,7 +50,7 @@ Suggests:
     lme4,
     merDeriv,
     parallel,
-    qs,
+    qs2,
     rmarkdown,
     testthat (>= 3.0.0),
     UCSCXenaShiny,
diff --git a/NAMESPACE b/NAMESPACE
index f636e8d..41212ba 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,9 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(base::print,br_diagnostics)
+S3method(print,breg_comparison)
 export(br_avail_method_config)
 export(br_avail_methods)
 export(br_avail_methods_use_exp)
+export(br_compare_models)
 export(br_diagnose)
 export(br_get_config)
 export(br_get_data)
@@ -30,6 +32,7 @@ export(br_show_fitted_line)
 export(br_show_fitted_line_2d)
 export(br_show_forest)
 export(br_show_forest_circle)
+export(br_show_forest_comparison)
 export(br_show_forest_ggstats)
 export(br_show_forest_ggstatsplot)
 export(br_show_nomogram)
diff --git a/NEWS.md b/NEWS.md
index 5bb2057..f42a344 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# bregr 1.4.0
+
+- Replaced the deprecated package `qs` with `qs2`.
+
 # bregr 1.3.2
 
 - Checked with `devtools::check(env_vars = c('_R_CHECK_DEPENDS_ONLY_' = "true"))` and fixed reported check issues.
diff --git a/R/02-pipeline.R b/R/02-pipeline.R
index 0fa2c4d..ce67109 100644
--- a/R/02-pipeline.R
+++ b/R/02-pipeline.R
@@ -532,7 +532,7 @@ runner_ <- function(m, data, dots, opts = NULL) {
   }
 
   if (isTRUE(as.logical(getOption("bregr.save_model", default = FALSE)))) {
-    rlang::check_installed(c("ids", "qs"))
+    rlang::check_installed(c("ids", "qs2"))
     md_path <- getOption("bregr.path", default = "")
     if (md_path == "") {
       if (requireNamespace("fs", quietly = TRUE)) {
@@ -560,7 +560,7 @@ runner_ <- function(m, data, dots, opts = NULL) {
     } else {
       md_file <- file.path(md_path, paste0(dg, ".qs"))
     }
-    qs::qsave(model, file = md_file)
+    qs2::qs_save(model, file = md_file)
     model <- md_file
   }
 
diff --git a/R/03-accessors.R b/R/03-accessors.R
index 72da8ae..e7d54e7 100644
--- a/R/03-accessors.R
+++ b/R/03-accessors.R
@@ -124,8 +124,8 @@ br_get_models <- function(obj, idx = NULL, auto_drop = TRUE) {
       if (len > 1000) {
         cli::cli_inform("directly retrieve >1000 models may resource-consuming, subsetting with {.arg idx} is more recommended")
       }
-      rlang::check_installed("qs")
-      mds <- map(mds, qs::qread)
+      rlang::check_installed("qs2")
+      mds <- map(mds, qs2::qs_read)
     }
   } else {
     if (is.numeric(idx)) {
@@ -142,8 +142,8 @@ br_get_models <- function(obj, idx = NULL, auto_drop = TRUE) {
 
     mds <- mds[idx]
     if (!insight::is_model(mds[[1]]) && requireNamespace("fs", quietly = TRUE) && fs::is_file(mds[[1]])) {
-      rlang::check_installed("qs")
-      mds <- map(mds, qs::qread)
+      rlang::check_installed("qs2")
+      mds <- map(mds, qs2::qs_read)
     }
     if (length(idx) == 1 && auto_drop) mds <- mds[[1]]
   }
diff --git a/R/08-compare.R b/R/08-compare.R
new file mode 100644
index 0000000..526abe7
--- /dev/null
+++ b/R/08-compare.R
@@ -0,0 +1,401 @@
+# Comparison utilities for univariate vs multivariate models
+#
+# Provides functions to compare univariate and multivariate regression models
+# side by side, similar to the functionality shown in autoReg package.
+# =====================
+
+
+#' Compare univariate and multivariate models
+#'
+#' @description
+#' `r lifecycle::badge('experimental')`
+#'
+#' This function builds both univariate models (each predictor separately) and
+#' a multivariate model (all predictors together), then combines the results
+#' for comparison. This is useful for understanding how predictor effects change
+#' when accounting for other variables.
+#'
+#' @param data A `data.frame` containing all necessary variables for analysis.
+#' @param y Character vector specifying dependent variables (response variables).
+#' For GLM models, this is typically a single character (e.g., `"outcome"`).
+#' For Cox-PH models, it should be a length-2 vector in the format `c("time", "status")`.
+#' @param x Character vector specifying focal independent terms (predictors).
+#' These will be modeled both individually (univariate) and together (multivariate).
+#' @param x2 Character vector specifying control independent terms (predictors, optional).
+#' These are included in all models (both univariate and multivariate).
+#' @param method Method for model construction. See [br_set_model()] for details.
+#' @param ... Additional arguments passed to [br_run()].
+#' @param n_workers Integer, indicating number of workers for parallel processing.
+#' @param model_args A list of arguments passed to `br_set_model()`.
+#' @param run_args A list of arguments passed to `br_run()`.
+#'
+#' @returns A list with class `breg_comparison` containing:
+#' - `univariate`: breg object with univariate model results
+#' - `multivariate`: breg object with multivariate model results
+#' - `combined_results`: Combined results data frame with a `mode` column
+#' - `combined_results_tidy`: Combined tidy results with a `mode` column
+#'
+#' @export
+#' @family br_compare
+#' @examples
+#' # Compare univariate vs multivariate for Cox models
+#' lung <- survival::lung |>
+#'   dplyr::filter(ph.ecog != 3)
+#' lung$ph.ecog <- factor(lung$ph.ecog)
+#'
+#' comparison <- br_compare_models(
+#'   lung,
+#'   y = c("time", "status"),
+#'   x = c("ph.ecog", "ph.karno", "pat.karno"),
+#'   x2 = c("age", "sex"),
+#'   method = "coxph"
+#' )
+#'
+#' # View combined results
+#' comparison$combined_results_tidy
+#'
+#' # Create forest plot comparison
+#' br_show_forest_comparison(comparison)
+#' @testexamples
+#' expect_s3_class(comparison, "breg_comparison")
+#' expect_true("mode" %in% colnames(comparison$combined_results_tidy))
+br_compare_models <- function(
+  data,
+  y,
+  x,
+  x2 = NULL,
+  method,
+  ...,
+  n_workers = 1L,
+  model_args = list(),
+  run_args = list()
+) {
+  # Validate inputs
+  assert_character(y, allow_na = FALSE)
+  assert_character(x, allow_na = FALSE)
+  assert_character(x2, allow_na = FALSE, allow_null = TRUE)
+
+  if (length(x) < 2) {
+    cli::cli_abort("{.arg x} must contain at least 2 variables for comparison")
+  }
+
+  # Build univariate models (each x separately with x2 as controls)
+  cli::cli_inform("Building univariate models...")
+  univariate <- br_pipeline(
+    data = data,
+    y = y,
+    x = x,
+    x2 = x2,
+    method = method,
+    n_workers = n_workers,
+    model_args = model_args,
+    run_args = run_args,
+    ...
+  )
+
+  # Build multivariate model (all x together with x2 as controls)
+  cli::cli_inform("Building multivariate model...")
+  # For multivariate, we treat all x as x2 (controls) with a dummy focal variable
+  # Actually, we need to create a single model with all x variables
+  # We can do this by setting x to just the first variable name and x2 to include all others
+
+  # Create a multivariate model by combining all x into the model
+  # We'll use a different approach: build using the same pipeline but with all x combined
+  x_combined <- paste(x, collapse = " + ")
+
+  # Build the multivariate model with all focal variables
+  multivariate_br <- breg(data) |>
+    br_set_y(y) |>
+    br_set_x(x[1]) # Use first as placeholder
+
+  # Set x2 to include all x except first, plus original x2
+  if (is.null(x2)) {
+    multivariate_br <- multivariate_br |> br_set_x2(x[-1])
+  } else {
+    multivariate_br <- multivariate_br |> br_set_x2(c(x[-1], x2))
+  }
+
+  # Now we need to hack this to create a single model with all x variables
+  # Let's take a different approach - manually construct the multivariate model
+
+  # Actually, let's use a cleaner approach:
+  # Build a single-model breg object with all x as focal (using a dummy approach)
+  multivariate <- build_multivariate_model(data, y, x, x2, method, run_args)
+
+  # Add mode column to distinguish results
+  univariate_results <- br_get_results(univariate)
+  univariate_results$mode <- "univariate"
+
+  multivariate_results <- br_get_results(multivariate)
+  multivariate_results$mode <- "multivariate"
+
+  univariate_results_tidy <- br_get_results(univariate, tidy = TRUE)
+  univariate_results_tidy$mode <- "univariate"
+
+  multivariate_results_tidy <- br_get_results(multivariate, tidy = TRUE)
+  multivariate_results_tidy$mode <- "multivariate"
+
+  # Filter to only focal variables for comparison
+  univariate_results_focal <- univariate_results |>
+    dplyr::filter(.data$Focal_variable == .data$variable)
+
+  multivariate_results_focal <- multivariate_results |>
+    dplyr::filter(.data$variable %in% x)
+
+  univariate_results_tidy_focal <- univariate_results_tidy |>
+    dplyr::filter(.data$Focal_variable == .data$term)
+
+  multivariate_results_tidy_focal <- multivariate_results_tidy |>
+    dplyr::filter(.data$term %in% x)
+
+  # Combine results
+  combined_results <- dplyr::bind_rows(
+    univariate_results_focal,
+    multivariate_results_focal
+  )
+
+  combined_results_tidy <- dplyr::bind_rows(
+    univariate_results_tidy_focal,
+    multivariate_results_tidy_focal
+  )
+
+  # Create comparison object
+  comparison <- list(
+    univariate = univariate,
+    multivariate = multivariate,
+    combined_results = combined_results,
+    combined_results_tidy = combined_results_tidy
+  )
+
+  class(comparison) <- c("breg_comparison", "list")
+  attr(comparison, "exponentiate") <- attr(univariate, "exponentiate")
+
+  comparison
+}
+
+
+#' Build a multivariate model with all focal variables
+#' @keywords internal
+#' @noRd
+build_multivariate_model <- function(data, y, x, x2, method, run_args) {
+  # Create a breg object with all x in a single model
+  # We'll use the first x as focal and rest as x2
+  br <- breg(data) |>
+    br_set_y(y) |>
+    br_set_x(x[1])
+
+  # Combine remaining x with x2
+  if (length(x) > 1) {
+    combined_x2 <- c(x[-1], x2)
+  } else {
+    combined_x2 <- x2
+  }
+
+  if (!is.null(combined_x2) && length(combined_x2) > 0) {
+    br <- br |> br_set_x2(combined_x2)
+  }
+
+  br <- br |> br_set_model(method)
+
+  # Run with specified arguments
+  if (length(run_args) > 0) {
+    br <- do.call(br_run, c(list(obj = br), run_args))
+  } else {
+    br <- br_run(br)
+  }
+
+  # Modify the results to treat all x variables as focal
+  # We need to extract results for all x variables and mark them as focal
+  results <- br@results
+  results_tidy <- br@results_tidy
+
+  # Update Focal_variable for all x variables
+  results <- results |>
+    dplyr::mutate(
+      Focal_variable = dplyr::if_else(
+        .data$variable %in% x,
+        .data$variable,
+        .data$Focal_variable
+      )
+    )
+
+  results_tidy <- results_tidy |>
+    dplyr::mutate(
+      Focal_variable = dplyr::if_else(
+        .data$term %in% x,
+        .data$term,
+        .data$Focal_variable
+      )
+    )
+
+  br@results <- results
+  br@results_tidy <- results_tidy
+
+  br
+}
+
+
+#' Show forest plot for model comparison
+#'
+#' @description
+#' `r lifecycle::badge('experimental')`
+#'
+#' Creates a forest plot comparing univariate and multivariate model results
+#' side by side. Each variable shows estimates from both modeling approaches.
+#'
+#' @param comparison A `breg_comparison` object from [br_compare_models()].
+#' @param ... Additional arguments passed to [forestploter::forest()].
+#' @param xlim Numeric vector of length 2 specifying x-axis limits.
+#' @param rm_controls If `TRUE`, show only focal variables (default).
+#'
+#' @returns A forest plot object.
+#' @export
+#' @family br_compare
+#' @examples
+#' lung <- survival::lung |>
+#'   dplyr::filter(ph.ecog != 3)
+#' lung$ph.ecog <- factor(lung$ph.ecog)
+#'
+#' comparison <- br_compare_models(
+#'   lung,
+#'   y = c("time", "status"),
+#'   x = c("ph.ecog", "ph.karno", "pat.karno"),
+#'   x2 = c("age", "sex"),
+#'   method = "coxph"
+#' )
+#'
+#' br_show_forest_comparison(comparison)
+#' @testexamples
+#' expect_s3_class(br_show_forest_comparison(comparison), "forestplot")
+br_show_forest_comparison <- function(
+  comparison,
+  ...,
+  xlim = NULL,
+  rm_controls = TRUE
+) {
+  if (!inherits(comparison, "breg_comparison")) {
+    cli::cli_abort("{.arg comparison} must be a {.cls breg_comparison} object from {.fn br_compare_models}")
+  }
+
+  dots <- rlang::list2(...)
+  exponentiate <- attr(comparison, "exponentiate")
+
+  # Get combined results
+  dt <- comparison$combined_results
+
+  if (rm_controls) {
+    # Only show focal variables
+    dt <- dt |>
+      dplyr::filter(.data$Focal_variable == .data$variable)
+  }
+
+  # Format the data for forest plot
+  dt <- dt |>
+    dplyr::arrange(.data$variable, .data$mode) |>
+    dplyr::mutate(
+      ` ` = paste(rep(" ", 20), collapse = " "),
+      `Estimate (95% CI)` = dplyr::case_when(
+        dt$reference_row ~ "Reference",
+        is.na(dt$std.error) ~ "",
+        TRUE ~
+          sprintf(
+            "%.2f (%.2f to %.2f)",
+            .data$estimate,
+            .data$conf.low,
+            .data$conf.high
+          )
+      ),
+      P = if_else(
+        is.na(.data$p.value),
+        "",
+        format.pval(.data$p.value, digits = 2, eps = 0.001)
+      ),
+      conf.low = if_else(is.na(.data$conf.low), .data$estimate, .data$conf.low),
+      conf.high = if_else(
+        is.na(.data$conf.high),
+        .data$estimate,
+        .data$conf.high
+      ),
+      mode = tools::toTitleCase(.data$mode)
+    )
+
+  # Calculate xlim if not provided
+  if (is.null(xlim)) {
+    xlim <- c(
+      floor(min(dt$conf.low, na.rm = TRUE)),
+      ceiling(max(dt$conf.high, na.rm = TRUE))
+    )
+    if (is.infinite(xlim[1])) {
+      cli_warn("infinite CI detected, set a minimal value -100")
+      xlim[1] <- -100
+    }
+    if (is.infinite(xlim[2])) {
+      cli_warn("infinite CI detected, set a maximal value 100")
+      xlim[2] <- 100
+    }
+  }
+
+  # Group by variable to show them together
+  dt <- dt |>
+    dplyr::group_by(.data$variable) |>
+    dplyr::mutate(
+      variable = if_else(
+        dplyr::row_number() == 1,
+        .data$variable,
+        ""
+      )
+    ) |>
+    dplyr::ungroup()
+
+  # Select and rename columns for display
+  dt <- dt |>
+    dplyr::select(
+      Variable = "variable",
+      Level = "label",
+      Mode = "mode",
+      N = "n_obs",
+      " ",
+      "Estimate (95% CI)",
+      c("P", "estimate", "conf.low", "conf.high")
+    )
+
+  # Set reference line
+  if (exponentiate && !("ref_line" %in% names(dots))) {
+    dots[["ref_line"]] <- 1L
+  }
+
+  # Create forest plot
+  p <- do.call(
+    forestploter::forest,
+    c(
+      list(
+        data = dt[, 1:7],
+        est = dt$estimate,
+        lower = dt$conf.low,
+        upper = dt$conf.high,
+        ci_column = 5,
+        xlim = xlim
+      ),
+      dots
+    )
+  )
+
+  p
+}
+
+
+#' Print method for breg_comparison object
+#' @param x A `breg_comparison` object.
+#' @param ... Additional arguments (not used).
+#' @export
+#' @method print breg_comparison
+print.breg_comparison <- function(x, ...) {
+  cli::cli_text("A {.cls breg_comparison} object")
+  cli::cli_text("")
+  cli::cli_text("Univariate models: {nrow(x$univariate@results_tidy)} terms from {length(x$univariate@models)} models")
+  cli::cli_text("Multivariate model: {nrow(x$multivariate@results_tidy)} terms from 1 model")
+  cli::cli_text("")
+  cli::cli_text("Use {.fn br_show_forest_comparison} to visualize the comparison")
+  cli::cli_text("Access results via {.code $combined_results} or {.code $combined_results_tidy}")
+  invisible(x)
+}
diff --git a/README.Rmd b/README.Rmd
index 03a9287..e3aca95 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -167,6 +167,26 @@ br_show_forest(
 We also provide some interfaces from other packages for plotting constructed model(s), e.g., `br_show_forest_ggstats()`, `br_show_forest_ggstatsplot()`, `br_show_fitted_line()`,
 and `br_show_fitted_line_2d()`.
 
+#### Comparing Univariate vs Multivariate Models
+
+A common analysis task is to compare how predictor effects change when modeled individually (univariate) versus together (multivariate). The `br_compare_models()` function builds both types of models and displays them side-by-side:
+
+```{r dpi=150, fig.height=5}
+# Compare univariate and multivariate models
+comparison <- br_compare_models(
+  lung,
+  y = c("time", "status"),
+  x = c("ph.ecog", "ph.karno", "pat.karno"),
+  x2 = c("age", "sex"),
+  method = "coxph"
+)
+
+# Show forest plot with both models
+br_show_forest_comparison(comparison)
+```
+
+This allows you to see how adjusting for other predictors affects the estimates for each variable.
+
 For Cox-PH modeling results (focal variables must be continuous type), we provide a risk network plotting function.
 
 ```{r}
@@ -229,7 +249,7 @@ All functions are documented in the [package reference](https://wanglabcsu.githu
 
 ## Coverage
 
-```{r}
+```{r eval=FALSE}
 covr::package_coverage()
 ```
 
diff --git a/README.md b/README.md
index d41b25c..f3bc9a3 100644
--- a/README.md
+++ b/README.md
@@ -44,7 +44,6 @@ A simplified overview of batch regression modeling is given below for
 illustration:
 
 <p align="center">
-
 <img src="man/figures/breg_formula.svg">
 </p>
 
@@ -78,7 +77,7 @@ Load package(s):
 library(bregr)
 #> Welcome to 'bregr' package!
 #> =======================================================================
-#> You are using bregr version 1.3.0
+#> You are using bregr version 1.3.2
 #> 
 #> Project home : https://github.com/WangLabCSU/bregr
 #> Documentation: https://wanglabcsu.github.io/bregr/
@@ -140,8 +139,9 @@ mds_p <- br_pipeline(
   n_workers = 3
 )
 #> exponentiate estimates of model(s) constructed from coxph method at default
-#> ■■■■■■■■■■■■■                     40% | ETA: 16s
-#>                                                  
+#> ■■■■■■■                           20% | ETA: 13s
+#> 
+#> ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
 ```
 
 ``` r
@@ -287,7 +287,7 @@ modeling results.
 br_show_forest(mds)
 ```
 
-<img src="man/figures/README-unnamed-chunk-7-1.png" width="100%" />
+<img src="man/figures/README-unnamed-chunk-7-1.png" alt="" width="100%" />
 
 We can tune the plot to only keep focal variables and adjust the limits
 of x axis.
@@ -304,13 +304,43 @@ br_show_forest(
 )
 ```
 
-<img src="man/figures/README-unnamed-chunk-8-1.png" width="100%" />
+<img src="man/figures/README-unnamed-chunk-8-1.png" alt="" width="100%" />
 
 We also provide some interfaces from other packages for plotting
 constructed model(s), e.g., `br_show_forest_ggstats()`,
 `br_show_forest_ggstatsplot()`, `br_show_fitted_line()`, and
 `br_show_fitted_line_2d()`.
 
+#### Comparing Univariate vs Multivariate Models
+
+A common analysis task is to compare how predictor effects change when
+modeled individually (univariate) versus together (multivariate). The
+`br_compare_models()` function builds both types of models and displays
+them side-by-side:
+
+``` r
+# Compare univariate and multivariate models
+comparison <- br_compare_models(
+  lung,
+  y = c("time", "status"),
+  x = c("ph.ecog", "ph.karno", "pat.karno"),
+  x2 = c("age", "sex"),
+  method = "coxph"
+)
+#> Building univariate models...
+#> exponentiate estimates of model(s) constructed from coxph method at default
+#> Building multivariate model...
+#> exponentiate estimates of model(s) constructed from coxph method at default
+
+# Show forest plot with both models
+br_show_forest_comparison(comparison)
+```
+
+<img src="man/figures/README-unnamed-chunk-9-1.png" alt="" width="100%" />
+
+This allows you to see how adjusting for other predictors affects the
+estimates for each variable.
+
 For Cox-PH modeling results (focal variables must be continuous type),
 we provide a risk network plotting function.
 
@@ -330,7 +360,7 @@ br_show_risk_network(mds2)
 #> please note only continuous focal terms analyzed and visualized
 ```
 
-<img src="man/figures/README-unnamed-chunk-10-1.png" width="100%" />
+<img src="man/figures/README-unnamed-chunk-11-1.png" alt="" width="100%" />
 
 #### Model Score Prediction and Survival Curves
 
@@ -360,7 +390,7 @@ br_show_survival_curves(
 #> values
 ```
 
-<img src="man/figures/README-unnamed-chunk-12-1.png" width="100%" />
+<img src="man/figures/README-unnamed-chunk-13-1.png" alt="" width="100%" />
 
 ### Table
 
@@ -428,17 +458,6 @@ site](https://wanglabcsu.github.io/bregr/).
 
 ``` r
 covr::package_coverage()
-#> bregr Coverage: 67.03%
-#> R/98-utils.R: 58.17%
-#> R/04-show-nomogram-helpers.R: 60.00%
-#> R/01-class.R: 61.19%
-#> R/07-diagnostics.R: 63.41%
-#> R/04-show.R: 66.52%
-#> R/03-accessors.R: 75.31%
-#> R/02-pipeline.R: 75.74%
-#> R/06-avail.R: 78.57%
-#> R/99-zzz.R: 92.31%
-#> R/05-polar.R: 92.37%
 ```
 
 ## Citation
diff --git a/_pkgdown.yml b/_pkgdown.yml
index fc39695..2393ecc 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -14,10 +14,12 @@ reference:
     contents:
     - breg
     - print.breg
+    - print.breg_comparison
   - title: Workflow
     desc: Core workflow and pipeline wrapper.
     contents:
     - br_pipeline
+    - br_compare_models
   - title: Accessor
     desc: Attributes and data accessors from [breg] objects.
     contents:
@@ -30,6 +32,7 @@ reference:
     - br_show_risk_network
     - br_show_forest_ggstats
     - br_show_forest_ggstatsplot
+    - br_show_forest_comparison
   - title: Model Diagnostics and Use
     desc: Inspect and visualize models for diagnostic purposes, and enable effective model utilization.
     contents:
diff --git a/copilot-instructions.md b/copilot-instructions.md
index ba2a893..d85202d 100644
--- a/copilot-instructions.md
+++ b/copilot-instructions.md
@@ -169,6 +169,7 @@ br_show_table(result)
 - Use `pak` instead of `install.packages()` for faster package management
 - Skip vignette building during development iterations
 - Use `_R_CHECK_FORCE_SUGGESTS_=false` for faster checking
+- Use `_R_CHECK_DEPENDS_ONLY_=true` for depends checking
 - Consider parallel testing for large test suites
 
 ## Common Issues and Solutions
@@ -234,4 +235,4 @@ rig default 4.5.1
 R --slave -e "sessionInfo()"
 ```
 
-This setup provides a modern, efficient development environment for the bregr R package using best practices and contemporary tools.
\ No newline at end of file
+This setup provides a modern, efficient development environment for the bregr R package using best practices and contemporary tools.
diff --git a/man/br_compare_models.Rd b/man/br_compare_models.Rd
new file mode 100644
index 0000000..031e424
--- /dev/null
+++ b/man/br_compare_models.Rd
@@ -0,0 +1,83 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/08-compare.R
+\name{br_compare_models}
+\alias{br_compare_models}
+\title{Compare univariate and multivariate models}
+\usage{
+br_compare_models(
+  data,
+  y,
+  x,
+  x2 = NULL,
+  method,
+  ...,
+  n_workers = 1L,
+  model_args = list(),
+  run_args = list()
+)
+}
+\arguments{
+\item{data}{A \code{data.frame} containing all necessary variables for analysis.}
+
+\item{y}{Character vector specifying dependent variables (response variables).
+For GLM models, this is typically a single character (e.g., \code{"outcome"}).
+For Cox-PH models, it should be a length-2 vector in the format \code{c("time", "status")}.}
+
+\item{x}{Character vector specifying focal independent terms (predictors).
+These will be modeled both individually (univariate) and together (multivariate).}
+
+\item{x2}{Character vector specifying control independent terms (predictors, optional).
+These are included in all models (both univariate and multivariate).}
+
+\item{method}{Method for model construction. See \code{\link[=br_set_model]{br_set_model()}} for details.}
+
+\item{...}{Additional arguments passed to \code{\link[=br_run]{br_run()}}.}
+
+\item{n_workers}{Integer, indicating number of workers for parallel processing.}
+
+\item{model_args}{A list of arguments passed to \code{br_set_model()}.}
+
+\item{run_args}{A list of arguments passed to \code{br_run()}.}
+}
+\value{
+A list with class \code{breg_comparison} containing:
+\itemize{
+\item \code{univariate}: breg object with univariate model results
+\item \code{multivariate}: breg object with multivariate model results
+\item \code{combined_results}: Combined results data frame with a \code{mode} column
+\item \code{combined_results_tidy}: Combined tidy results with a \code{mode} column
+}
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+This function builds both univariate models (each predictor separately) and
+a multivariate model (all predictors together), then combines the results
+for comparison. This is useful for understanding how predictor effects change
+when accounting for other variables.
+}
+\examples{
+# Compare univariate vs multivariate for Cox models
+lung <- survival::lung |>
+  dplyr::filter(ph.ecog != 3)
+lung$ph.ecog <- factor(lung$ph.ecog)
+
+comparison <- br_compare_models(
+  lung,
+  y = c("time", "status"),
+  x = c("ph.ecog", "ph.karno", "pat.karno"),
+  x2 = c("age", "sex"),
+  method = "coxph"
+)
+
+# View combined results
+comparison$combined_results_tidy
+
+# Create forest plot comparison
+br_show_forest_comparison(comparison)
+}
+\seealso{
+Other br_compare: 
+\code{\link{br_show_forest_comparison}()}
+}
+\concept{br_compare}
diff --git a/man/br_show_forest_comparison.Rd b/man/br_show_forest_comparison.Rd
new file mode 100644
index 0000000..b1fac4c
--- /dev/null
+++ b/man/br_show_forest_comparison.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/08-compare.R
+\name{br_show_forest_comparison}
+\alias{br_show_forest_comparison}
+\title{Show forest plot for model comparison}
+\usage{
+br_show_forest_comparison(comparison, ..., xlim = NULL, rm_controls = TRUE)
+}
+\arguments{
+\item{comparison}{A \code{breg_comparison} object from \code{\link[=br_compare_models]{br_compare_models()}}.}
+
+\item{...}{Additional arguments passed to \code{\link[forestploter:forest]{forestploter::forest()}}.}
+
+\item{xlim}{Numeric vector of length 2 specifying x-axis limits.}
+
+\item{rm_controls}{If \code{TRUE}, show only focal variables (default).}
+}
+\value{
+A forest plot object.
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+Creates a forest plot comparing univariate and multivariate model results
+side by side. Each variable shows estimates from both modeling approaches.
+}
+\examples{
+lung <- survival::lung |>
+  dplyr::filter(ph.ecog != 3)
+lung$ph.ecog <- factor(lung$ph.ecog)
+
+comparison <- br_compare_models(
+  lung,
+  y = c("time", "status"),
+  x = c("ph.ecog", "ph.karno", "pat.karno"),
+  x2 = c("age", "sex"),
+  method = "coxph"
+)
+
+br_show_forest_comparison(comparison)
+}
+\seealso{
+Other br_compare: 
+\code{\link{br_compare_models}()}
+}
+\concept{br_compare}
diff --git a/man/figures/README-unnamed-chunk-11-1.png b/man/figures/README-unnamed-chunk-11-1.png
index dfac099..601d8fc 100644
Binary files a/man/figures/README-unnamed-chunk-11-1.png and b/man/figures/README-unnamed-chunk-11-1.png differ
diff --git a/man/figures/README-unnamed-chunk-13-1.png b/man/figures/README-unnamed-chunk-13-1.png
index 8dd62a6..3aa19bd 100644
Binary files a/man/figures/README-unnamed-chunk-13-1.png and b/man/figures/README-unnamed-chunk-13-1.png differ
diff --git a/man/figures/README-unnamed-chunk-7-1.png b/man/figures/README-unnamed-chunk-7-1.png
index ed91162..91f3c2a 100644
Binary files a/man/figures/README-unnamed-chunk-7-1.png and b/man/figures/README-unnamed-chunk-7-1.png differ
diff --git a/man/figures/README-unnamed-chunk-8-1.png b/man/figures/README-unnamed-chunk-8-1.png
index e445fd4..829d0b5 100644
Binary files a/man/figures/README-unnamed-chunk-8-1.png and b/man/figures/README-unnamed-chunk-8-1.png differ
diff --git a/man/print.breg_comparison.Rd b/man/print.breg_comparison.Rd
new file mode 100644
index 0000000..2bb9aae
--- /dev/null
+++ b/man/print.breg_comparison.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/08-compare.R
+\name{print.breg_comparison}
+\alias{print.breg_comparison}
+\title{Print method for breg_comparison object}
+\usage{
+\method{print}{breg_comparison}(x, ...)
+}
+\arguments{
+\item{x}{A \code{breg_comparison} object.}
+
+\item{...}{Additional arguments (not used).}
+}
+\description{
+Print method for breg_comparison object
+}
diff --git a/tests/testthat/test-compare.R b/tests/testthat/test-compare.R
new file mode 100644
index 0000000..52c3233
--- /dev/null
+++ b/tests/testthat/test-compare.R
@@ -0,0 +1,118 @@
+test_that("br_compare_models works with continuous variables", {
+  library(bregr)
+
+  lung <- survival::lung |>
+    dplyr::filter(ph.ecog != 3)
+
+  comparison <- br_compare_models(
+    lung,
+    y = c("time", "status"),
+    x = c("ph.karno", "pat.karno"),
+    x2 = c("age", "sex"),
+    method = "coxph"
+  )
+
+  # Check object class
+  expect_s3_class(comparison, "breg_comparison")
+
+  # Check that mode column exists
+  expect_true("mode" %in% colnames(comparison$combined_results))
+  expect_true("mode" %in% colnames(comparison$combined_results_tidy))
+
+  # Check that we have both univariate and multivariate results
+  expect_true("univariate" %in% comparison$combined_results$mode)
+  expect_true("multivariate" %in% comparison$combined_results$mode)
+
+  # Check that we have results for all focal variables
+  expect_true(all(c("ph.karno", "pat.karno") %in% comparison$combined_results$variable))
+})
+
+test_that("br_compare_models works with categorical variables", {
+  library(bregr)
+
+  lung <- survival::lung |>
+    dplyr::filter(ph.ecog != 3)
+  lung$ph.ecog <- factor(lung$ph.ecog)
+
+  comparison <- br_compare_models(
+    lung,
+    y = c("time", "status"),
+    x = c("ph.ecog", "ph.karno"),
+    x2 = c("age", "sex"),
+    method = "coxph"
+  )
+
+  expect_s3_class(comparison, "breg_comparison")
+  expect_true("mode" %in% colnames(comparison$combined_results_tidy))
+})
+
+test_that("br_compare_models requires at least 2 variables", {
+  library(bregr)
+
+  lung <- survival::lung |>
+    dplyr::filter(ph.ecog != 3)
+
+  expect_error(
+    br_compare_models(
+      lung,
+      y = c("time", "status"),
+      x = "ph.karno", # Only one variable
+      x2 = c("age", "sex"),
+      method = "coxph"
+    ),
+    "must contain at least 2 variables"
+  )
+})
+
+test_that("br_show_forest_comparison creates a forest plot", {
+  library(bregr)
+
+  lung <- survival::lung |>
+    dplyr::filter(ph.ecog != 3)
+
+  comparison <- br_compare_models(
+    lung,
+    y = c("time", "status"),
+    x = c("ph.karno", "pat.karno"),
+    x2 = c("age", "sex"),
+    method = "coxph"
+  )
+
+  p <- br_show_forest_comparison(comparison)
+
+  expect_s3_class(p, "forestplot")
+})
+
+test_that("br_show_forest_comparison validates input", {
+  library(bregr)
+
+  # Create a regular breg object (not a comparison)
+  m <- br_pipeline(
+    mtcars,
+    y = "mpg",
+    x = c("cyl", "disp"),
+    method = "gaussian"
+  )
+
+  expect_error(
+    br_show_forest_comparison(m),
+    "breg_comparison"
+  )
+})
+
+test_that("print method for breg_comparison works", {
+  library(bregr)
+
+  lung <- survival::lung |>
+    dplyr::filter(ph.ecog != 3)
+
+  comparison <- br_compare_models(
+    lung,
+    y = c("time", "status"),
+    x = c("ph.karno", "pat.karno"),
+    x2 = c("age", "sex"),
+    method = "coxph"
+  )
+
+  expect_s3_class(comparison, "breg_comparison")
+})
diff --git a/tests/testthat/test-roxytest-testexamples-08-compare.R b/tests/testthat/test-roxytest-testexamples-08-compare.R
new file mode 100644
index 0000000..130a8d0
--- /dev/null
+++ b/tests/testthat/test-roxytest-testexamples-08-compare.R
@@ -0,0 +1,47 @@
+# Generated by roxytest: do not edit by hand!
+
+# File R/"08-compare.R": @testexamples
+
+test_that("Function br_compare_models() @ L62", {
+  
+  # Compare univariate vs multivariate for Cox models
+  lung <- survival::lung |>
+    dplyr::filter(ph.ecog != 3)
+  lung$ph.ecog <- factor(lung$ph.ecog)
+  
+  comparison <- br_compare_models(
+    lung,
+    y = c("time", "status"),
+    x = c("ph.ecog", "ph.karno", "pat.karno"),
+    x2 = c("age", "sex"),
+    method = "coxph"
+  )
+  
+  # View combined results
+  comparison$combined_results_tidy
+  
+  # Create forest plot comparison
+  br_show_forest_comparison(comparison)
+  expect_s3_class(comparison, "breg_comparison")
+  expect_true("mode" %in% colnames(comparison$combined_results_tidy))
+})
+
+
+test_that("Function br_show_forest_comparison() @ L270", {
+  
+  lung <- survival::lung |>
+    dplyr::filter(ph.ecog != 3)
+  lung$ph.ecog <- factor(lung$ph.ecog)
+  
+  comparison <- br_compare_models(
+    lung,
+    y = c("time", "status"),
+    x = c("ph.ecog", "ph.karno", "pat.karno"),
+    x2 = c("age", "sex"),
+    method = "coxph"
+  )
+  
+  br_show_forest_comparison(comparison)
+  expect_s3_class(br_show_forest_comparison(comparison), "forestplot")
+})
+