diff --git a/.Rbuildignore b/.Rbuildignore
index 79cee3a..801080a 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -17,4 +17,5 @@ copilot-instructions.md
 ^doc$
 ^Meta$
 ^icon$
+^examples$
 ^CRAN-SUBMISSION$
diff --git a/NAMESPACE b/NAMESPACE
index a6cb523..72528df 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -16,6 +16,7 @@ export(br_get_x)
 export(br_get_x2)
 export(br_get_y)
 export(br_pipeline)
+export(br_predict)
 export(br_rename_models)
 export(br_run)
 export(br_set_model)
@@ -29,6 +30,7 @@ export(br_show_forest_ggstats)
 export(br_show_forest_ggstatsplot)
 export(br_show_residuals)
 export(br_show_risk_network)
+export(br_show_survival_curves)
 export(br_show_table)
 export(br_show_table_gt)
 export(breg)
@@ -67,7 +69,13 @@ importFrom(rlang,is_string)
 importFrom(rlang,is_vector)
 importFrom(rlang,set_names)
 importFrom(stats,cor)
+importFrom(stats,pchisq)
+importFrom(stats,predict)
+importFrom(stats,quantile)
+importFrom(survival,Surv)
 importFrom(survival,coxph)
+importFrom(survival,survdiff)
+importFrom(survival,survfit)
 importFrom(utils,combn)
 importFrom(utils,packageDescription)
 importFrom(utils,packageVersion)
diff --git a/R/00-bregr-package.R b/R/00-bregr-package.R
index 5890802..898e3aa 100644
--- a/R/00-bregr-package.R
+++ b/R/00-bregr-package.R
@@ -5,11 +5,11 @@
 #' @import vctrs
 #' @import cli
 #' @importFrom mirai .progress
-#' @importFrom survival coxph
+#' @importFrom survival Surv survdiff survfit coxph
 #' @importFrom lifecycle deprecated
 #' @importFrom rlang .data .env
 #' @importFrom utils combn packageDescription packageVersion
-#' @importFrom stats cor
+#' @importFrom stats cor pchisq quantile predict
 #' @rawNamespace if (getRversion() < "4.3.0") importFrom("S7", "@")
 ## usethis namespace: end
 NULL
diff --git a/R/02-pipeline.R b/R/02-pipeline.R
index ee3b88c..3fb5702 100644
--- a/R/02-pipeline.R
+++ b/R/02-pipeline.R
@@ -99,18 +99,19 @@
 #' # 3. Customized model -----------
 #' dt <- data.frame(x = rnorm(100))
 #' dt$y <- rpois(100, exp(1 + dt$x))
+#' \donttest{
 #' m5 <- breg(dt) |>
 #'   br_set_y("y") |>
 #'   br_set_x("x") |>
 #'   br_set_model(method = 'quasi(variance = "mu", link = "log")') |>
 #'   br_run()
+#' }
 #'
 #' @testexamples
 #' assert_breg_obj(m)
 #' assert_breg_obj(m2)
 #' assert_breg_obj(m3)
 #' assert_breg_obj(m4)
-#' assert_breg_obj(m5)
 #' @seealso [accessors] for accessing `breg` object properties.
 NULL
 
@@ -120,10 +121,9 @@ NULL
 br_pipeline <- function(
     data, y, x, method, x2 = NULL,
     group_by = NULL,
-    n_workers = 1L, run_parallel = 1L,
+    n_workers = 1L, run_parallel = lifecycle::deprecated(),
     model_args = list(),
     run_args = list()) {
-
   if (lifecycle::is_present(run_parallel)) {
     lifecycle::deprecate_warn("1.1.0", "bregr::br_run(run_parallel = )", "bregr::br_run(n_workers = )")
     n_workers <- run_parallel
diff --git a/R/03-accessors.R b/R/03-accessors.R
index 54e0f09..6187fe3 100644
--- a/R/03-accessors.R
+++ b/R/03-accessors.R
@@ -21,6 +21,7 @@
 #' - `br_get_results()` returns modeling result `data.frame`.
 #'
 #' @name accessors
+#' @family accessors
 #' @seealso [pipeline] for building `breg` objects.
 #' @examples
 #' m <- br_pipeline(mtcars,
@@ -188,3 +189,83 @@ br_get_results <- function(obj, tidy = FALSE, ...) {
   }
   dplyr::filter(results, ...)
 }
+
+#' Predict method for breg objects
+#'
+#' @description
+#' `r lifecycle::badge('experimental')`
+#'
+#' Generate predictions from fitted models in a `breg` object.
+#' For Cox regression models, returns linear predictors (log relative hazard).
+#' For other models, returns predicted values.
+#'
+#' @param obj A `breg` object with fitted models.
+#' @param newdata Optional data frame for predictions. If NULL, uses original data.
+#' @param idx Model index, an integer or string.
+#' @param type Type of prediction. For Cox models: "lp" (linear predictor, default)
+#' or "risk" (relative risk). For other models: "response" (default) or "link".
+#' @returns Typically, a numeric vector of predictions.
+#' @export
+#' @family accessors
+#' @examples
+#' # Cox regression example
+#' if (requireNamespace("survival", quietly = TRUE)) {
+#'   lung <- survival::lung |> dplyr::filter(ph.ecog != 3)
+#'   mds <- br_pipeline(
+#'     lung,
+#'     y = c("time", "status"),
+#'     x = c("age", "ph.ecog"),
+#'     x2 = "sex",
+#'     method = "coxph"
+#'   )
+#'   scores <- br_predict(mds)
+#'   head(scores)
+#' }
+br_predict <- function(obj, newdata = NULL, idx = NULL, type = NULL) {
+  assert_breg_obj_with_results(obj)
+
+  # Get the model to use
+  if (is.null(idx)) {
+    cli::cli_inform("{.arg idx} not set, use the first model")
+    idx <- 1
+  } else {
+    if (length(idx) != 1) {
+      cli::cli_abort("please specify one model")
+    }
+  }
+  model <- br_get_models(obj, idx)
+
+  # Get data for prediction
+  if (is.null(newdata)) {
+    newdata <- br_get_data(obj)
+  }
+
+  # Determine prediction type
+  model_class <- class(model)[1]
+  if (is.null(type)) {
+    if (model_class == "coxph") {
+      # https://www.rdocumentation.org/packages/survival/versions/3.8-3/topics/predict.coxph
+      type <- "lp" # linear predictor (log relative hazard)
+    } else {
+      type <- "response" # predicted response
+    }
+    cli::cli_inform("{.arg type} is not specified, use {type} for the model")
+  }
+
+  # Generate predictions
+  tryCatch(
+    {
+      predictions <- predict(model, newdata = newdata, type = type)
+
+      # Handle missing values
+      if (any(is.na(predictions))) {
+        cli::cli_warn("Some predictions are NA, consider checking your data for missing values")
+      }
+
+      predictions
+    },
+    error = function(e) {
+      cli::cli_abort("Failed to generate predictions: {e$message}")
+    }
+  )
+}
diff --git a/R/04-show.R b/R/04-show.R
index f943b7c..1dfc86c 100644
--- a/R/04-show.R
+++ b/R/04-show.R
@@ -527,6 +527,174 @@ br_show_table_gt <- function(
   t
 }
 
+#' Show survival curves based on model scores
+#'
+#' @description
+#' `r lifecycle::badge('experimental')`
+#'
+#' Generate survival curves by grouping observations based on model prediction scores.
+#' This function is specifically designed for Cox regression models and creates
+#' survival curves comparing different risk groups.
+#'
+#' @param breg A `breg` object with fitted Cox regression models.
+#' @param idx Index or name of the model to use for prediction.
+#' If NULL, uses the first model.
+#' @param n_groups Number of groups to create based on score quantiles. Default is 3.
+#' @param group_labels Custom labels for the groups. If NULL, uses "Low", "Medium", "High"
+#' for 3 groups or "Q1", "Q2", etc. for other numbers.
+#' @param title Plot title. If NULL, generates automatic title.
+#' @param subtitle Plot subtitle.
+#' @returns A ggplot2 object showing survival curves.
+#' @export
+#' @family br_show
+#' @examples
+#' \donttest{
+#' # Cox regression example with survival curves
+#' if (requireNamespace("survival", quietly = TRUE)) {
+#'   lung <- survival::lung |> dplyr::filter(ph.ecog != 3)
+#'   mds <- br_pipeline(
+#'     lung,
+#'     y = c("time", "status"),
+#'     x = c("age", "ph.ecog"),
+#'     x2 = "sex",
+#'     method = "coxph"
+#'   )
+#'   p <- br_show_survival_curves(mds)
+#'   print(p)
+#' }
+#' }
+br_show_survival_curves <- function(breg,
+                                    idx = NULL,
+                                    n_groups = 3,
+                                    group_labels = NULL,
+                                    title = NULL,
+                                    subtitle = NULL) {
+  assert_breg_obj_with_results(breg)
+
+  # Get the model to use
+  if (is.null(idx)) {
+    cli::cli_inform("{.arg idx} not set, use the first model")
+    idx <- 1
+  } else {
+    if (length(idx) != 1) {
+      cli::cli_abort("please specify one model")
+    }
+  }
+  model <- br_get_models(breg, idx)
+
+  # Check if we have survival models
+  if (!inherits(model, "coxph")) {
+    cli::cli_abort("survival curves are only supported for Cox regression models (coxph)")
+  }
+
+  # Get predictions (risk scores)
+  scores <- br_predict(breg, idx = idx, type = "lp")
+  data <- br_get_data(breg)
+
+  # Get survival variables
+  y_vars <- br_get_y(breg)
+  if (length(y_vars) != 2) {
+    cli::cli_abort("Expected 2 survival variables (time, status), got {length(y_vars)}")
+  }
+
+  time_var <- y_vars[1]
+  status_var <- y_vars[2]
+
+  # Create score groups
+  score_quantiles <- quantile(scores, probs = seq(0, 1, length.out = n_groups + 1), na.rm = TRUE)
+
+  # Handle case where scores have little variation
+  if (length(unique(score_quantiles)) <= 2) {
+    cli::cli_warn("Score values have limited variation, grouping may not be meaningful")
+  }
+
+  score_groups <- cut(scores,
+    breaks = score_quantiles,
+    include.lowest = TRUE,
+    labels = FALSE
+  )
+
+  # Create group labels
+  if (is.null(group_labels)) {
+    if (n_groups == 3) {
+      group_labels <- c("Low Risk", "Medium Risk", "High Risk")
+    } else if (n_groups == 2) {
+      group_labels <- c("Low Risk", "High Risk")
+    } else {
+      group_labels <- paste0("Q", 1:n_groups)
+    }
+  } else if (length(group_labels) != n_groups) {
+    cli::cli_abort("Length of group_labels ({length(group_labels)}) must match n_groups ({n_groups})")
+  }
+
+  # Prepare data for survival analysis
+  surv_data <- data |>
+    dplyr::mutate(
+      .score = scores,
+      .score_group = factor(score_groups, labels = group_labels),
+      .time = .data[[time_var]],
+      .status = .data[[status_var]]
+    ) |>
+    dplyr::filter(!is.na(.data$.score_group))
+
+  # Create survival object
+  surv_obj <- survival::Surv(surv_data$.time, surv_data$.status)
+
+  # Fit Kaplan-Meier curves for each group
+  km_fit <- survival::survfit(surv_obj ~ .score_group, data = surv_data)
+
+  # Create survival curve data for plotting
+  surv_summary <- summary(km_fit)
+
+  plot_data <- data.frame(
+    time = surv_summary$time,
+    surv = surv_summary$surv,
+    group = surv_summary$strata,
+    upper = surv_summary$upper,
+    lower = surv_summary$lower
+  )
+
+  # Clean group names
+  plot_data$group <- gsub(".*=", "", plot_data$group)
+
+  # Create the plot
+  if (is.null(title)) {
+    focal_var <- names(br_get_models(breg))[idx %||% 1]
+    title <- paste("Survival Curves by", focal_var, "Model Score")
+  }
+
+  p <- ggplot2::ggplot(plot_data, ggplot2::aes(x = .data$time, y = .data$surv, color = .data$group)) +
+    ggplot2::geom_step(linewidth = 1) +
+    ggplot2::geom_ribbon(ggplot2::aes(ymin = .data$lower, ymax = .data$upper, fill = .data$group),
+      alpha = 0.3, color = NA
+    ) +
+    ggplot2::scale_y_continuous(limits = c(0, 1), labels = function(x) paste0(x * 100, "%")) +
+    ggplot2::labs(
+      title = title,
+      subtitle = subtitle,
+      x = "Time",
+      y = "Survival Probability",
+      color = "Risk Group",
+      fill = "Risk Group"
+    ) +
+    ggplot2::theme_minimal() +
+    ggplot2::theme(
+      legend.position = "bottom",
+      plot.title = ggplot2::element_text(hjust = 0.5),
+      plot.subtitle = ggplot2::element_text(hjust = 0.5)
+    )
+
+  # Add log-rank test p-value if requested
+  if (n_groups > 1) {
+    logrank_test <- survival::survdiff(surv_obj ~ .score_group, data = surv_data)
+    p_value <- pchisq(logrank_test$chisq, df = length(logrank_test$n) - 1, lower.tail = FALSE)
+
+    p <- p + ggplot2::labs(caption = paste("Log-rank test p-value:", format.pval(p_value, digits = 3)))
+  }
+
+  return(p)
+}
+
 #' Show residuals vs fitted plot for regression models
 #'
 #' @description
diff --git a/README.Rmd b/README.Rmd
index 2bae063..f262d4d 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -113,7 +113,7 @@ mds <- br_pipeline(
 
 Run in parallel:
 
-```{r}
+```{r warning=FALSE}
 mds_p <- br_pipeline(
   lung,
   y = c("time", "status"),
@@ -183,6 +183,25 @@ mds2 <- br_pipeline(
 br_show_risk_network(mds2)
 ```
 
+#### Model Score Prediction and Survival Curves
+
+For Cox-PH models, you can generate model predictions (risk scores) and create survival curves grouped by these scores:
+
+```{r}
+# Generate model predictions
+scores <- br_predict(mds2, idx = "ph.ecog")
+head(scores)
+```
+
+```{r dpi=150, fig.height=6}
+# Create survival curves based on model scores
+br_show_survival_curves(
+  mds2,
+  idx = "ph.ecog",
+  n_groups = 3,
+  title = "Survival Curves by 'ph.ecog' Model Risk Score"
+)
+```
 
 ### Table
 
diff --git a/README.md b/README.md
index 1dd85b4..2292e49 100644
--- a/README.md
+++ b/README.md
@@ -136,8 +136,10 @@ mds_p <- br_pipeline(
   method = "coxph",
   n_workers = 3
 )
-#> Warning: running in parallel is typically not recommended for small number (<100) of focal terms
 #> exponentiate estimates of model(s) constructed from coxph method at default
+#> ■■■■■■■                           20% | ETA:  6s
+#> ■■■■■■■■■■■■■■■■■■■               60% | ETA:  7s
+#>                                                  
 ```
 
 ``` r
@@ -328,6 +330,36 @@ br_show_risk_network(mds2)
 
 <img src="man/figures/README-unnamed-chunk-10-1.png" width="100%" />
 
+#### Model Score Prediction and Survival Curves
+
+For Cox-PH models, you can generate model predictions (risk scores) and
+create survival curves grouped by these scores:
+
+``` r
+# Generate model predictions
+scores <- br_predict(mds2, idx = "ph.ecog")
+#> `type` is not specified, use lp for the model
+#> Warning: Some predictions are NA, consider checking your data for missing
+#> values
+head(scores)
+#>          1          2          3          4          5          6 
+#>  0.3692998 -0.1608293 -0.2936304  0.1811648 -0.2493634  0.3692998
+```
+
+``` r
+# Create survival curves based on model scores
+br_show_survival_curves(
+  mds2,
+  idx = "ph.ecog",
+  n_groups = 3,
+  title = "Survival Curves by 'ph.ecog' Model Risk Score"
+)
+#> Warning: Some predictions are NA, consider checking your data for missing
+#> values
+```
+
+<img src="man/figures/README-unnamed-chunk-12-1.png" width="100%" />
+
 ### Table
 
 Show tidy table result as pretty table:
@@ -394,12 +426,12 @@ site](https://wanglabcsu.github.io/bregr/).
 
 ``` r
 covr::package_coverage()
-#> bregr Coverage: 77.37%
-#> R/98-utils.R: 37.18%
-#> R/02-pipeline.R: 78.01%
-#> R/04-show.R: 78.45%
+#> bregr Coverage: 76.26%
+#> R/98-utils.R: 51.14%
+#> R/03-accessors.R: 72.50%
+#> R/04-show.R: 75.31%
+#> R/02-pipeline.R: 76.32%
 #> R/06-avail.R: 78.57%
-#> R/03-accessors.R: 80.00%
 #> R/01-class.R: 90.70%
 #> R/99-zzz.R: 90.91%
 #> R/05-polar.R: 95.19%
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 24079df..bd3282f 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -32,6 +32,8 @@ reference:
   - title: Model Diagnostics and Use
     desc: Inspect and visualize models for diagnostic purposes, and enable effective model utilization.
     contents:
+    - br_predict
+    - br_show_survival_curves
     - br_show_fitted_line
     - br_show_fitted_line_2d
     - br_show_residuals
diff --git a/man/accessors.Rd b/man/accessors.Rd
index 9a34230..78e018b 100644
--- a/man/accessors.Rd
+++ b/man/accessors.Rd
@@ -104,4 +104,8 @@ br_get_results(m, tidy = TRUE, term == "cyl")
 }
 \seealso{
 \link{pipeline} for building \code{breg} objects.
+
+Other accessors: 
+\code{\link{br_predict}()}
 }
+\concept{accessors}
diff --git a/man/br_predict.Rd b/man/br_predict.Rd
new file mode 100644
index 0000000..6442f3d
--- /dev/null
+++ b/man/br_predict.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/03-accessors.R
+\name{br_predict}
+\alias{br_predict}
+\title{Predict method for breg objects}
+\usage{
+br_predict(obj, newdata = NULL, idx = NULL, type = NULL)
+}
+\arguments{
+\item{obj}{A \code{breg} object with fitted models.}
+
+\item{newdata}{Optional data frame for predictions. If NULL, uses original data.}
+
+\item{idx}{Model index, an integer or string.}
+
+\item{type}{Type of prediction. For Cox models: "lp" (linear predictor, default)
+or "risk" (relative risk). For other models: "response" (default) or "link".}
+}
+\value{
+Typically, a numeric vector of predictions.
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+Generate predictions from fitted models in a \code{breg} object.
+For Cox regression models, returns linear predictors (log relative hazard).
+For other models, returns predicted values.
+}
+\examples{
+# Cox regression example
+if (requireNamespace("survival", quietly = TRUE)) {
+  lung <- survival::lung |> dplyr::filter(ph.ecog != 3)
+  mds <- br_pipeline(
+    lung,
+    y = c("time", "status"),
+    x = c("age", "ph.ecog"),
+    x2 = "sex",
+    method = "coxph"
+  )
+  scores <- br_predict(mds)
+  head(scores)
+}
+}
+\seealso{
+Other accessors: 
+\code{\link{accessors}}
+}
+\concept{accessors}
diff --git a/man/br_show_fitted_line.Rd b/man/br_show_fitted_line.Rd
index c7fae81..f7a9abc 100644
--- a/man/br_show_fitted_line.Rd
+++ b/man/br_show_fitted_line.Rd
@@ -47,6 +47,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_residuals}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()},
 \code{\link{br_show_table_gt}()}
 }
diff --git a/man/br_show_fitted_line_2d.Rd b/man/br_show_fitted_line_2d.Rd
index 9c5d3f6..d9bbbc5 100644
--- a/man/br_show_fitted_line_2d.Rd
+++ b/man/br_show_fitted_line_2d.Rd
@@ -44,6 +44,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_residuals}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()},
 \code{\link{br_show_table_gt}()}
 }
diff --git a/man/br_show_forest.Rd b/man/br_show_forest.Rd
index acf8dae..e50a184 100644
--- a/man/br_show_forest.Rd
+++ b/man/br_show_forest.Rd
@@ -73,6 +73,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_residuals}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()},
 \code{\link{br_show_table_gt}()}
 }
diff --git a/man/br_show_forest_ggstats.Rd b/man/br_show_forest_ggstats.Rd
index 1fe9296..c8014d5 100644
--- a/man/br_show_forest_ggstats.Rd
+++ b/man/br_show_forest_ggstats.Rd
@@ -41,6 +41,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_residuals}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()},
 \code{\link{br_show_table_gt}()}
 }
diff --git a/man/br_show_forest_ggstatsplot.Rd b/man/br_show_forest_ggstatsplot.Rd
index ab56eba..66d139b 100644
--- a/man/br_show_forest_ggstatsplot.Rd
+++ b/man/br_show_forest_ggstatsplot.Rd
@@ -43,6 +43,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstats}()},
 \code{\link{br_show_residuals}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()},
 \code{\link{br_show_table_gt}()}
 }
diff --git a/man/br_show_residuals.Rd b/man/br_show_residuals.Rd
index 7ee5699..fe54d0e 100644
--- a/man/br_show_residuals.Rd
+++ b/man/br_show_residuals.Rd
@@ -56,6 +56,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstats}()},
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()},
 \code{\link{br_show_table_gt}()}
 }
diff --git a/man/br_show_risk_network.Rd b/man/br_show_risk_network.Rd
index ee98828..f1af6df 100644
--- a/man/br_show_risk_network.Rd
+++ b/man/br_show_risk_network.Rd
@@ -38,6 +38,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstats}()},
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_residuals}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()},
 \code{\link{br_show_table_gt}()}
 
diff --git a/man/br_show_survival_curves.Rd b/man/br_show_survival_curves.Rd
new file mode 100644
index 0000000..4eaed7c
--- /dev/null
+++ b/man/br_show_survival_curves.Rd
@@ -0,0 +1,70 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/04-show.R
+\name{br_show_survival_curves}
+\alias{br_show_survival_curves}
+\title{Show survival curves based on model scores}
+\usage{
+br_show_survival_curves(
+  breg,
+  idx = NULL,
+  n_groups = 3,
+  group_labels = NULL,
+  title = NULL,
+  subtitle = NULL
+)
+}
+\arguments{
+\item{breg}{A \code{breg} object with fitted Cox regression models.}
+
+\item{idx}{Index or name of the model to use for prediction.
+If NULL, uses the first model.}
+
+\item{n_groups}{Number of groups to create based on score quantiles. Default is 3.}
+
+\item{group_labels}{Custom labels for the groups. If NULL, uses "Low", "Medium", "High"
+for 3 groups or "Q1", "Q2", etc. for other numbers.}
+
+\item{title}{Plot title. If NULL, generates automatic title.}
+
+\item{subtitle}{Plot subtitle.}
+}
+\value{
+A ggplot2 object showing survival curves.
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+Generate survival curves by grouping observations based on model prediction scores.
+This function is specifically designed for Cox regression models and creates
+survival curves comparing different risk groups.
+}
+\examples{
+\donttest{
+# Cox regression example with survival curves
+if (requireNamespace("survival", quietly = TRUE)) {
+  lung <- survival::lung |> dplyr::filter(ph.ecog != 3)
+  mds <- br_pipeline(
+    lung,
+    y = c("time", "status"),
+    x = c("age", "ph.ecog"),
+    x2 = "sex",
+    method = "coxph"
+  )
+  p <- br_show_survival_curves(mds)
+  print(p)
+}
+}
+}
+\seealso{
+Other br_show: 
+\code{\link{br_show_fitted_line}()},
+\code{\link{br_show_fitted_line_2d}()},
+\code{\link{br_show_forest}()},
+\code{\link{br_show_forest_ggstats}()},
+\code{\link{br_show_forest_ggstatsplot}()},
+\code{\link{br_show_residuals}()},
+\code{\link{br_show_risk_network}()},
+\code{\link{br_show_table}()},
+\code{\link{br_show_table_gt}()}
+}
+\concept{br_show}
diff --git a/man/br_show_table.Rd b/man/br_show_table.Rd
index 44ae9cf..c989d3d 100644
--- a/man/br_show_table.Rd
+++ b/man/br_show_table.Rd
@@ -52,6 +52,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_residuals}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table_gt}()}
 }
 \concept{br_show}
diff --git a/man/br_show_table_gt.Rd b/man/br_show_table_gt.Rd
index 7be4542..d9fe4de 100644
--- a/man/br_show_table_gt.Rd
+++ b/man/br_show_table_gt.Rd
@@ -50,6 +50,7 @@ Other br_show:
 \code{\link{br_show_forest_ggstatsplot}()},
 \code{\link{br_show_residuals}()},
 \code{\link{br_show_risk_network}()},
+\code{\link{br_show_survival_curves}()},
 \code{\link{br_show_table}()}
 }
 \concept{br_show}
diff --git a/man/figures/README-unnamed-chunk-12-1.png b/man/figures/README-unnamed-chunk-12-1.png
new file mode 100644
index 0000000..c05d21f
Binary files /dev/null and b/man/figures/README-unnamed-chunk-12-1.png differ
diff --git a/man/pipeline.Rd b/man/pipeline.Rd
index 318b32a..e06fa0b 100644
--- a/man/pipeline.Rd
+++ b/man/pipeline.Rd
@@ -18,6 +18,7 @@ br_pipeline(
   x2 = NULL,
   group_by = NULL,
   n_workers = 1L,
+  run_parallel = lifecycle::deprecated(),
   model_args = list(),
   run_args = list()
 )
@@ -167,11 +168,13 @@ m4 <- br_pipeline(mtcars,
 # 3. Customized model -----------
 dt <- data.frame(x = rnorm(100))
 dt$y <- rpois(100, exp(1 + dt$x))
+\donttest{
 m5 <- breg(dt) |>
   br_set_y("y") |>
   br_set_x("x") |>
   br_set_model(method = 'quasi(variance = "mu", link = "log")') |>
   br_run()
+}
 
 }
 \seealso{
diff --git a/tests/testthat/test-predict-survival.R b/tests/testthat/test-predict-survival.R
new file mode 100644
index 0000000..9b37c0b
--- /dev/null
+++ b/tests/testthat/test-predict-survival.R
@@ -0,0 +1,99 @@
+test_that("br_predict works for Cox regression", {
+  skip_if_not_installed("survival")
+
+  # Create test data similar to survival::lung
+  set.seed(123)
+  n <- 50
+  test_data <- data.frame(
+    time = runif(n, 1, 1000),
+    status = rbinom(n, 1, 0.6),
+    age = rnorm(n, 60, 10),
+    sex = factor(sample(c("M", "F"), n, replace = TRUE))
+  )
+
+  # Create breg object with Cox regression
+  breg_obj <- br_pipeline(
+    test_data,
+    y = c("time", "status"),
+    x = "age",
+    x2 = "sex",
+    method = "coxph"
+  )
+
+  # Test predictions
+  predictions <- br_predict(breg_obj)
+
+  expect_true(is.numeric(predictions))
+  expect_equal(length(predictions), n)
+  expect_true(!any(is.na(predictions)))
+})
+
+test_that("br_show_survival_curves works for Cox regression", {
+  skip_if_not_installed("survival")
+  skip_if_not_installed("ggplot2")
+
+  # Create test data
+  set.seed(123)
+  n <- 100
+  test_data <- data.frame(
+    time = rexp(n, 0.01),
+    status = rbinom(n, 1, 0.7),
+    age = rnorm(n, 60, 10),
+    sex = factor(sample(c("M", "F"), n, replace = TRUE))
+  )
+
+  # Create breg object with Cox regression
+  breg_obj <- br_pipeline(
+    test_data,
+    y = c("time", "status"),
+    x = "age",
+    x2 = "sex",
+    method = "coxph"
+  )
+
+  # Test survival curves
+  p <- br_show_survival_curves(breg_obj, n_groups = 2)
+
+  expect_s3_class(p, "ggplot")
+})
+
+test_that("br_predict fails gracefully for non-coxph models", {
+  # Create test data for linear regression
+  test_data <- data.frame(
+    y = rnorm(50),
+    x = rnorm(50)
+  )
+
+  # Create breg object with linear regression
+  breg_obj <- br_pipeline(
+    test_data,
+    y = "y",
+    x = "x",
+    method = "gaussian"
+  )
+
+  # Test predictions work for non-coxph
+  predictions <- br_predict(breg_obj)
+  expect_true(is.numeric(predictions))
+})
+
+test_that("br_show_survival_curves fails for non-coxph models", {
+  # Create test data for linear regression
+  test_data <- data.frame(
+    y = rnorm(50),
+    x = rnorm(50)
+  )
+
+  # Create breg object with linear regression
+  breg_obj <- br_pipeline(
+    test_data,
+    y = "y",
+    x = "x",
+    method = "gaussian"
+  )
+
+  # Test that survival curves fail for non-coxph
+  expect_error(
+    br_show_survival_curves(breg_obj)
+  )
+})
diff --git a/tests/testthat/test-roxytest-testexamples-02-pipeline.R b/tests/testthat/test-roxytest-testexamples-02-pipeline.R
index 737be6c..9d37e84 100644
--- a/tests/testthat/test-roxytest-testexamples-02-pipeline.R
+++ b/tests/testthat/test-roxytest-testexamples-02-pipeline.R
@@ -2,7 +2,7 @@
 
 # File R/"02-pipeline.R": @testexamples
 
-test_that("[unknown alias] @ L115", {
+test_that("[unknown alias] @ L116", {
   
   library(bregr)
   # 1. Pipeline -------------------------
@@ -47,16 +47,17 @@ test_that("[unknown alias] @ L115", {
   # 3. Customized model -----------
   dt <- data.frame(x = rnorm(100))
   dt$y <- rpois(100, exp(1 + dt$x))
+  
   m5 <- breg(dt) |>
     br_set_y("y") |>
     br_set_x("x") |>
     br_set_model(method = 'quasi(variance = "mu", link = "log")') |>
     br_run()
   
+  
   assert_breg_obj(m)
   assert_breg_obj(m2)
   assert_breg_obj(m3)
   assert_breg_obj(m4)
-  assert_breg_obj(m5)
 })
 
diff --git a/tests/testthat/test-roxytest-testexamples-03-accessors.R b/tests/testthat/test-roxytest-testexamples-03-accessors.R
index f189f23..ff50903 100644
--- a/tests/testthat/test-roxytest-testexamples-03-accessors.R
+++ b/tests/testthat/test-roxytest-testexamples-03-accessors.R
@@ -2,7 +2,7 @@
 
 # File R/"03-accessors.R": @testexamples
 
-test_that("[unknown alias] @ L49", {
+test_that("[unknown alias] @ L50", {
   
   m <- br_pipeline(mtcars,
     y = "mpg",
diff --git a/tests/testthat/test-roxytest-testexamples-04-show.R b/tests/testthat/test-roxytest-testexamples-04-show.R
index 1a5e6ff..04728c6 100644
--- a/tests/testthat/test-roxytest-testexamples-04-show.R
+++ b/tests/testthat/test-roxytest-testexamples-04-show.R
@@ -120,7 +120,7 @@ test_that("Function br_show_table_gt() @ L500", {
 })
 
 
-test_that("Function br_show_residuals() @ L569", {
+test_that("Function br_show_residuals() @ L737", {
   
   m <- br_pipeline(mtcars,
     y = "mpg",