diff --git a/DESCRIPTION b/DESCRIPTION
index 59410f74..01ed895d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: sccomp
 Type: Package
 Title: Differential Composition and Variability Analysis for Single-Cell Data
-Version: 2.1.33
+Version: 2.1.34
 Date: 2026-05-11
 Authors@R: c(person("Stefano", "Mangiola", email = "stefano.mangiola@unimelb.edu.au", role = c("aut", "cre")), person("Alexandra J.", "Roth-Schulze", role = "aut"), person("Marie", "Trussart", role = "aut"), person("Enrique", "Zozaya-Valdés", role = "aut"), person("Mengyao", "Ma", role = "aut"), person("Zijie", "Gao", role = "aut"), person("Alan F.", "Rubin", role = "aut"), person("Terence P.", "Speed", role = "aut"), person("Heejung", "Shim", role = "aut"), person("Anthony T.", "Papenfuss", role = "aut"))
 Description: Comprehensive R package for differential composition and variability analysis in single-cell RNA sequencing, CyTOF, and microbiome data. Provides robust Bayesian modeling with outlier detection, random effects, and advanced statistical methods for cell type proportion analysis. Features include probabilistic outlier identification, mixed-effect modeling, differential variability testing, and comprehensive visualization tools. Perfect for cancer research, immunology, developmental biology, and single-cell genomics applications.
diff --git a/NAMESPACE b/NAMESPACE
index 7354520f..c5866835 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -29,6 +29,7 @@ export(sccomp_remove_outliers)
 export(sccomp_remove_unwanted_effects)
 export(sccomp_remove_unwanted_variation)
 export(sccomp_replicate)
+export(sccomp_scatterplot)
 export(sccomp_stan_models_cache_dir)
 export(sccomp_test)
 export(sccomp_theme)
diff --git a/R/plot.R b/R/plot.R
index 8abc72ef..f8238fdd 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -58,11 +58,6 @@ plot.sccomp_tbl <- function(
   significance_statistic <- match.arg(significance_statistic)
   sort_by <- match.arg(sort_by)
 
-  # Quosures from estimate (same names as legacy `.cell_group` args; stored on tbl attributes)
-  .sample <- attr(x, ".sample")
-  .cell_group <- attr(x, ".cell_group")
-  .count <- attr(x, ".count")
-
   # Define the variables as NULL to avoid CRAN NOTES
   v_effect <- NULL
 
@@ -72,30 +67,6 @@ plot.sccomp_tbl <- function(
   if(x |> select(ends_with("FDR")) |> ncol() |> equals(0))
     stop("sccomp says: to produce plots, you need to run the function sccomp_test() on your estimates.")
 
-  data_proportion =
-    x |>
-    
-    # Otherwise does not work
-    select(-`factor`) |>
-    
-    pivot_wider(names_from = parameter, values_from = c(contains("c_"), contains("v_"))) |>
-    left_join(
-      attr(x, "count_data"),
-      by = join_by(!!.cell_group)
-    ) |> 
-    with_groups(!!.sample, ~ mutate(.x, proportion = (!!.count)/sum(!!.count)) )
-  
-  
-  # If I don't have outliers add them
-  if(x |> attr("outliers") |> is.null() |> not())
-    data_proportion = 
-    data_proportion |> 
-    left_join(x |> attr("outliers"), by = join_by(!!.sample, !!.cell_group))
- else 
-  data_proportion = data_proportion |> mutate(outlier = FALSE) 
-  
- 
-
   # Select the factors to plot  
   my_factors =
     x |>
@@ -114,14 +85,12 @@ plot.sccomp_tbl <- function(
         ~ {
           
           # If variable is continuous
-          if(data_proportion |> select(all_of(.x)) |> pull(1) |> is("numeric"))
+          if(attr(x, "count_data") |> select(all_of(.x)) |> pull(1) |> is("numeric"))
             my_plot =
-              plot_scatterplot(
+              sccomp_scatterplot(
                 .data = x,
-                data_proportion = data_proportion,
-                factor_of_interest = .x,
-                significance_threshold = significance_threshold,
-                my_theme = sccomp_theme()
+                factor = .x,
+                significance_threshold = significance_threshold
               )
           
           # If discrete
@@ -135,8 +104,7 @@ plot.sccomp_tbl <- function(
               ) 
           
           # Return
-          my_plot +
-            ggtitle(sprintf("Grouped by %s (for multi-factor models, associations could be hardly observable with unidimensional data stratification)", .x))
+          my_plot
         } )
     
   }
diff --git a/R/plot_scatterplot.R b/R/plot_scatterplot.R
index d09b73e1..680b09c7 100644
--- a/R/plot_scatterplot.R
+++ b/R/plot_scatterplot.R
@@ -1,3 +1,123 @@
+#' sccomp_scatterplot
+#'
+#' @description
+#' Creates a scatterplot of cell-group proportions against a continuous covariate,
+#' optionally overlaying the posterior predictive trend from the fitted model.
+#' Use this visualisation when the factor of interest is numeric; for discrete
+#' factors prefer [sccomp_boxplot()].
+#'
+#' @import ggplot2
+#' @importFrom tidyr pivot_longer pivot_wider unite unnest
+#' @importFrom dplyr select with_groups mutate left_join rename
+#'
+#' @param .data A tibble containing the results from `sccomp_estimate` and
+#'   `sccomp_test`, including the columns: cell_group name, sample name, read
+#'   counts, factor(s), p-values, and significance indicators.
+#' @param factor A character string specifying the continuous covariate of
+#'   interest included in the model.
+#' @param significance_threshold A numeric value indicating the threshold for
+#'   highlighting significant cell-groups. Defaults to 0.05.
+#' @param remove_unwanted_effects A logical value indicating whether to remove
+#'   unwanted variation from the data before plotting. Defaults to `FALSE`.
+#'
+#' @return A `ggplot` object representing the scatterplot of cell proportions
+#'   against the continuous covariate, faceted by cell group.
+#'
+#' @export
+#'
+#' @examples
+#'
+#' print("cmdstanr is needed to run this example.")
+#' # Note: Before running the example, ensure that the 'cmdstanr' package is installed:
+#' # install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
+#'
+#' \donttest{
+#' if (instantiate::stan_cmdstan_exists()) {
+#'     data("seurat_obj")
+#'
+#'     estimate <- sccomp_estimate(
+#'       seurat_obj,
+#'       formula_composition = ~ continuous_covariate,
+#'       formula_variability = ~ 1,
+#'       sample = "sample",
+#'       cell_group = "cell_group",
+#'       cores = 1
+#'     ) |>
+#'     sccomp_test()
+#'
+#'     # Plot proportions against the continuous covariate
+#'     sccomp_scatterplot(
+#'         .data = estimate,
+#'         factor = "continuous_covariate",
+#'         significance_threshold = 0.05
+#'     )
+#' }
+#' }
+#' @references
+#' S. Mangiola, A.J. Roth-Schulze, M. Trussart, E. Zozaya-Valdés, M. Ma, Z. Gao, A.F. Rubin, T.P. Speed, H. Shim, & A.T. Papenfuss, sccomp: Robust differential composition and variability analysis for single-cell data, Proc. Natl. Acad. Sci. U.S.A. 120 (33) e2203828120, https://doi.org/10.1073/pnas.2203828120 (2023).
+sccomp_scatterplot = function(
+    .data,
+    factor,
+    significance_threshold = 0.05,
+    remove_unwanted_effects = FALSE
+){
+
+  .cell_group = attr(.data, ".cell_group")
+  .count = attr(.data, ".count")
+  .sample = attr(.data, ".sample")
+
+  # Check if test have been done
+  if(.data |> select(ends_with("FDR")) |> ncol() |> equals(0))
+    stop("sccomp says: to produce plots, you need to run the function sccomp_test() on your estimates.")
+
+  .data_filtered =
+    subset_results_by_factor(.data, factor, keep_intercept = FALSE)
+
+  data_proportion =
+    .data_filtered |>
+
+    # Otherwise does not work
+    select(-`factor`) |>
+
+    pivot_wider(names_from = parameter, values_from = c(contains("c_"), contains("v_"))) |>
+    left_join(
+      attr(.data, "count_data") ,
+      by = quo_name(.cell_group)
+    ) |>
+    with_groups(!!.sample, ~ mutate(.x, proportion = (!!.count)/sum(!!.count)) )
+
+  if(remove_unwanted_effects){
+    .data_adjusted =
+      .data |>
+      sccomp_remove_unwanted_effects(formula_composition_keep = as.formula("~ " |> paste(factor))) |>
+      rename(proportion = adjusted_proportion)
+
+    data_proportion =
+      data_proportion |>
+      select(-proportion) |>
+      left_join(.data_adjusted, by = join_by(!!.cell_group, !!.sample))
+  }
+  else
+    message( "sccomp says: When visualising proportions, especially for complex models, consider setting `remove_unwanted_effects=TRUE`. This will adjust the proportions, preserving only the observed effect.")
+
+  # If I don't have outliers add them
+  if(.data |> attr("outliers") |> is.null() |> not())
+    data_proportion =
+    data_proportion |>
+    left_join(.data |> attr("outliers"), by = c(quo_name(.sample), quo_name(.cell_group)))
+ else
+  data_proportion = data_proportion |> mutate(outlier = FALSE)
+
+  plot_scatterplot(
+    .data = .data,
+    data_proportion = data_proportion,
+    factor_of_interest = factor,
+    significance_threshold = significance_threshold,
+    my_theme = sccomp_theme()
+  ) +
+    ggtitle(sprintf("Grouped by %s (for multi-factor models, associations could be hardly observable with unidimensional data stratification)", factor))
+
+}
 
 #' Plot Scatterplot of Cell-group Proportion
 #'
diff --git a/inst/NEWS.rd b/inst/NEWS.rd
index 92708f94..9f4f15bd 100644
--- a/inst/NEWS.rd
+++ b/inst/NEWS.rd
@@ -1,6 +1,12 @@
 \name{NEWS}
 \title{News for Package \pkg{sccomp}}
 
+\section{News in version 2.1.34}{
+\itemize{
+    \item Exported \code{sccomp_scatterplot()} for visualising cell-group proportions against a continuous covariate, complementing \code{sccomp_boxplot()} for discrete factors. The function accepts \code{.data}, \code{factor}, \code{significance_threshold}, and \code{remove_unwanted_effects}, and is used by the \code{plot()} method for numeric covariates.
+    \item Documented \code{sccomp_scatterplot()} in the introduction vignette (Continuous Factor Analysis section).
+}}
+
 \section{News in version 2.1.33}{
 \itemize{
     \item **New feature: smooth terms (splines) in formulas.** Added support for \code{s()} and \code{t2()} smooth terms in \code{sccomp} formulas, enabling flexible non-linear modelling of continuous covariates. Introduced utility functions for parsing smooth terms and a dedicated vignette demonstrating their use.
diff --git a/man/sccomp_scatterplot.Rd b/man/sccomp_scatterplot.Rd
new file mode 100644
index 00000000..f82d1ee0
--- /dev/null
+++ b/man/sccomp_scatterplot.Rd
@@ -0,0 +1,69 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_scatterplot.R
+\name{sccomp_scatterplot}
+\alias{sccomp_scatterplot}
+\title{sccomp_scatterplot}
+\usage{
+sccomp_scatterplot(
+  .data,
+  factor,
+  significance_threshold = 0.05,
+  remove_unwanted_effects = FALSE
+)
+}
+\arguments{
+\item{.data}{A tibble containing the results from \code{sccomp_estimate} and
+\code{sccomp_test}, including the columns: cell_group name, sample name, read
+counts, factor(s), p-values, and significance indicators.}
+
+\item{factor}{A character string specifying the continuous covariate of
+interest included in the model.}
+
+\item{significance_threshold}{A numeric value indicating the threshold for
+highlighting significant cell-groups. Defaults to 0.05.}
+
+\item{remove_unwanted_effects}{A logical value indicating whether to remove
+unwanted variation from the data before plotting. Defaults to \code{FALSE}.}
+}
+\value{
+A \code{ggplot} object representing the scatterplot of cell proportions
+against the continuous covariate, faceted by cell group.
+}
+\description{
+Creates a scatterplot of cell-group proportions against a continuous covariate,
+optionally overlaying the posterior predictive trend from the fitted model.
+Use this visualisation when the factor of interest is numeric; for discrete
+factors prefer \code{\link[=sccomp_boxplot]{sccomp_boxplot()}}.
+}
+\examples{
+
+print("cmdstanr is needed to run this example.")
+# Note: Before running the example, ensure that the 'cmdstanr' package is installed:
+# install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
+
+\donttest{
+if (instantiate::stan_cmdstan_exists()) {
+    data("seurat_obj")
+
+    estimate <- sccomp_estimate(
+      seurat_obj,
+      formula_composition = ~ continuous_covariate,
+      formula_variability = ~ 1,
+      sample = "sample",
+      cell_group = "cell_group",
+      cores = 1
+    ) |>
+    sccomp_test()
+
+    # Plot proportions against the continuous covariate
+    sccomp_scatterplot(
+        .data = estimate,
+        factor = "continuous_covariate",
+        significance_threshold = 0.05
+    )
+}
+}
+}
+\references{
+S. Mangiola, A.J. Roth-Schulze, M. Trussart, E. Zozaya-Valdés, M. Ma, Z. Gao, A.F. Rubin, T.P. Speed, H. Shim, & A.T. Papenfuss, sccomp: Robust differential composition and variability analysis for single-cell data, Proc. Natl. Acad. Sci. U.S.A. 120 (33) e2203828120, https://doi.org/10.1073/pnas.2203828120 (2023).
+}
diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R
index d7235cbb..0658327c 100644
--- a/tests/testthat/test-plot.R
+++ b/tests/testthat/test-plot.R
@@ -423,6 +423,24 @@ test_that("sccomp_plot_intervals_2D includes regression line from prec parameter
   )
 })
 
+test_that("sccomp_scatterplot works for continuous covariates", {
+  skip_cmdstan()
+
+  my_estimate |>
+    sccomp_test() |>
+    sccomp_scatterplot(factor = "continuous_covariate", significance_threshold = 0.025) |>
+    expect_s3_class("ggplot")
+
+  my_estimate |>
+    sccomp_test() |>
+    sccomp_scatterplot(
+      factor = "continuous_covariate",
+      significance_threshold = 0.025,
+      remove_unwanted_effects = TRUE
+    ) |>
+    expect_s3_class("ggplot")
+})
+
 test_that("sccomp_boxplot can accept additional ggplot layers", {
    skip_cmdstan()
 
diff --git a/vignettes/introduction.Rmd b/vignettes/introduction.Rmd
index 14452791..90953cc2 100644
--- a/vignettes/introduction.Rmd
+++ b/vignettes/introduction.Rmd
@@ -261,7 +261,7 @@ sccomp_result =
 
 ## Visualization and Summary Plots
 
-A plot of group proportions, faceted by groups. The blue boxplots represent the posterior predictive check. If the model is descriptively adequate for the data, the blue boxplots should roughly overlay the black boxplots, which represent the observed data. The outliers are coloured in red. A boxplot will be returned for every (discrete) covariate present in formula_composition. The colour coding represents the significant associations for composition and/or variability.
+A plot of group proportions, faceted by groups. The blue boxplots represent the posterior predictive check. If the model is descriptively adequate for the data, the blue boxplots should roughly overlay the black boxplots, which represent the observed data. The outliers are coloured in red. A boxplot will be returned for every (discrete) covariate present in formula_composition; continuous covariates are visualised with `sccomp_scatterplot()` instead (see Continuous Factor Analysis below). The colour coding represents the significant associations for composition and/or variability.
 
 ```{r boxplot-basic, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
 sccomp_result |> 
@@ -363,11 +363,19 @@ res =
       cell_group = "cell_group",
       cores = 1, verbose=FALSE,
       max_sampling_iterations = 2000
-    )
+    ) |>
+    sccomp_test()
 
 res
 ```
 
+For continuous covariates, use `sccomp_scatterplot()` instead of `sccomp_boxplot()`. Points are observed proportions; the black line is a smooth of the data and the blue ribbon is the posterior predictive trend.
+
+```{r continuous-scatterplot, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+res |>
+  sccomp_scatterplot(factor = "continuous_covariate")
+```
+
 ## Random Effect Modeling (Mixed-Effect Modeling)
 
 `sccomp` supports multilevel modeling by allowing the inclusion of random effects in the compositional and variability formulas. This is particularly useful when your data has hierarchical or grouped structures, such as measurements nested within subjects, batches, or experimental units. By incorporating random effects, sccomp can account for variability at different levels of your data, improving model fit and inference accuracy.