diff --git a/DESCRIPTION b/DESCRIPTION
index 350cffc6b..49b658367 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,7 +23,7 @@ URL: https://insightsengineering.github.io/teal.modules.general/,
 BugReports:
     https://github.com/insightsengineering/teal.modules.general/issues
 Depends:
-    ggplot2 (>= 3.4.0),
+    ggplot2 (>= 3.5.0),
     R (>= 4.1),
     shiny (>= 1.8.1),
     teal (>= 1.0.0.9003),
@@ -32,7 +32,7 @@ Imports:
     bslib (>= 0.8.0),
     checkmate (>= 2.1.0),
     colourpicker (>= 1.3.0),
-    dplyr (>= 1.0.5),
+    dplyr (>= 1.1.0),
     DT (>= 0.13),
     forcats (>= 1.0.0),
     generics (>= 0.1.3),
diff --git a/NAMESPACE b/NAMESPACE
index 10d634628..9a8bb56e3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -11,6 +11,7 @@ S3method(create_sparklines,numeric)
 S3method(teal.reporter::to_rmd,markdown_internal)
 S3method(tools::toHTML,markdown_internal)
 export(add_facet_labels)
+export(geom_mosaic)
 export(get_scatterplotmatrix_stats)
 export(tm_a_pca)
 export(tm_a_regression)
@@ -33,4 +34,5 @@ import(shiny)
 import(teal)
 import(teal.transform)
 importFrom(dplyr,"%>%")
+importFrom(dplyr,.data)
 importFrom(lifecycle,deprecated)
diff --git a/NEWS.md b/NEWS.md
index fcde995b0..ec081aaba 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,6 +3,7 @@
 ### Enhancements
 
 - Modules now return a `teal_report` object that contains the data, code and reporter. All the reporter buttons were removed from the modules' UI.
+- Support case when both variables are categorical in association and bivariate plots.
 
 # teal.modules.general 0.5.1
 
diff --git a/R/geom_mosaic.R b/R/geom_mosaic.R
new file mode 100644
index 000000000..6ecfddb4b
--- /dev/null
+++ b/R/geom_mosaic.R
@@ -0,0 +1,226 @@
+# minimal implementation of ggplot2 mosaic after ggmosaic was archived in CRAN
+#
+# This was heavily inspired by github.com/haleyjeppson/ggmosaic package but
+# simplified to only support 2 categorical variables
+
+#' Mosaic Rectangles Layer for ggplot2
+#'
+#' Adds a mosaic-style rectangles layer to a ggplot, visualizing the
+#' joint distribution of categorical variables.
+#' Each rectangle's size reflects the proportion of observations for
+#' combinations of `x` and `fill`.
+#'
+#' @param mapping Set of aesthetic mappings created by `aes()`. Must specify `x` and `fill`.
+#' @param data The data to be displayed in this layer.
+#' @param stat The statistical transformation to use on the data. Defaults to `"rects"`.
+#' @param position Position adjustment. Defaults to `"identity"`.
+#' @param ... Other arguments passed to `layer()`.
+#' @param na.rm Logical. Should missing values be removed?
+#' @param show.legend Logical. Should this layer be included in the legends?
+#' @param inherit.aes Logical. If `FALSE`, overrides default aesthetics.
+#'
+#' @return A ggplot2 layer that adds mosaic rectangles to the plot.
+#'
+#' @examples
+#' df <- data.frame(RACE = c("Black", "White", "Black", "Asian"), SEX = c("M", "M", "F", "F"))
+#' library(ggplot2)
+#' ggplot(df) +
+#'   geom_mosaic(aes(x = RACE, fill = SEX))
+#' @export
+geom_mosaic <- function(mapping = NULL, data = NULL,
+                        stat = "mosaic", position = "identity",
+                        ...,
+                        na.rm = FALSE, # nolint: object_name_linter.
+                        show.legend = TRUE, # nolint: object_name_linter.
+                        inherit.aes = TRUE) { # nolint: object_name_linter.
+
+  aes_x <- mapping$x
+  if (!is.null(aes_x)) {
+    aes_x <- list(rlang::quo_get_expr(mapping$x))
+    var_x <- paste0("x__", as.character(aes_x))
+    mapping[[var_x]] <- mapping$x
+  }
+
+  aes_fill <- mapping$fill
+  if (!is.null(aes_fill)) {
+    aes_fill <- rlang::quo_text(mapping$fill)
+  }
+
+  mapping$x <- structure(1L)
+
+  layer <- ggplot2::layer(
+    geom = GeomMosaic,
+    stat = "mosaic",
+    data = data,
+    mapping = mapping,
+    position = position,
+    show.legend = show.legend,
+    inherit.aes = inherit.aes,
+    check.aes = FALSE,
+    params = list(na.rm = na.rm, ...)
+  )
+  list(layer, .scale_x_mosaic())
+}
+
+#' @keywords internal
+GeomMosaic <- ggplot2::ggproto( # nolint: object_name_linter.
+  "GeomMosaic", ggplot2::GeomRect,
+  default_aes = ggplot2::aes(
+    colour = NA, linewidth = 0.5, linetype = 1, alpha = 1, fill = "grey30"
+  ),
+  draw_panel = function(data, panel_params, coord) {
+    if (all(is.na(data$colour))) data$colour <- scales::alpha(data$fill, data$alpha)
+    ggplot2::GeomRect$draw_panel(data, panel_params, coord)
+  },
+  required_aes = c("xmin", "xmax", "ymin", "ymax")
+)
+
+#' @keywords internal
+StatMosaic <- ggplot2::ggproto( # nolint: object_name_linter.
+  "StatMosaic", ggplot2::Stat,
+  required_aes = c("x", "fill"),
+  compute_group = function(data, scales) data,
+  compute_panel = function(data, scales) {
+    data$x <- data[, grepl("x__", colnames(data))]
+    result <- .calculate_coordinates(data)
+
+    results_non_zero <- result[result$.n != 0, ]
+    breaks <- unique(with(results_non_zero, (xmin + xmax) / 2))
+    labels <- unique(results_non_zero$x)
+    result$x <- list(list2env(list(breaks = breaks[breaks != 0], labels = labels[breaks != 0])))
+
+    result$group <- 1
+    result$PANEL <- unique(data$PANEL)
+    result
+  }
+)
+
+#' Determining scales for mosaics
+#'
+#' @param breaks,labels,minor_breaks One of:
+#' - `NULL` for no breaks / labels.
+#' - [ggplot2::waiver()] for the default breaks / labels computed by the scale.
+#' - A numeric / character vector giving the positions of the breaks / labels.
+#' - A function.
+#' See [ggplot2::scale_x_continuous()] for more details.
+#' @param na.value The value to be used for `NA` values.
+#' @param position For position scales, The position of the axis.
+#' left or right for y axes, top or bottom for x axes.
+#' @param ... other arguments passed to `continuous_scale()`.
+#' @keywords internal
+.scale_x_mosaic <- function(breaks = unique,
+                            minor_breaks = NULL,
+                            labels = unique,
+                            na.value = NA_real_, # nolint: object_name_linter.
+                            position = "bottom",
+                            ...) {
+  ggplot2::continuous_scale(
+    aesthetics = c(
+      "x", "xmin", "xmax", "xend", "xintercept", "xmin_final", "xmax_final",
+      "xlower", "xmiddle", "xupper"
+    ),
+    palette = identity,
+    breaks = breaks,
+    minor_breaks = minor_breaks,
+    labels = labels,
+    na.value = na.value,
+    position = position,
+    super = ScaleContinuousMosaic, ,
+    guide = ggplot2::waiver(),
+    ...
+  )
+}
+
+#' @keywords internal
+ScaleContinuousMosaic <- ggplot2::ggproto( # nolint: object_name_linter.
+  "ScaleContinuousMosaic", ggplot2::ScaleContinuousPosition,
+  train = function(self, x) {
+    if (length(x) == 0) {
+      return()
+    }
+    if (is.list(x)) {
+      scale_x <- x[[1]]
+      # re-assign the scale values now that we have the information - but only if necessary
+      if (is.function(self$breaks)) self$breaks <- scale_x$breaks
+      if (is.function(self$labels)) self$labels <- as.vector(scale_x$labels)
+      return(NULL)
+    }
+    if (is_discrete(x)) {
+      self$range$train(x = c(0, 1))
+      return(NULL)
+    }
+    self$range$train(x, call = self$call)
+  },
+  map = function(self, x, limits = self$get_limits()) {
+    if (is_discrete(x)) {
+      return(x)
+    }
+    if (is.list(x)) {
+      return(0)
+    } # need a number
+    scaled <- as.numeric(self$oob(x, limits))
+    ifelse(!is.na(scaled), scaled, self$na.value)
+  },
+  dimension = function(self, expand = c(0, 0)) c(-0.05, 1.05)
+)
+
+#' @noRd
+is_discrete <- function(x) is.factor(x) || is.character(x) || is.logical(x)
+
+#' @describeIn geom_mosaic
+#' Computes the coordinates for rectangles in a mosaic plot based
+#' on combinations of `x` and `fill` variables.
+#' For each unique `x` and `fill`, calculates the proportional
+#' widths and heights, stacking rectangles within each `x` group.
+#'
+#' ### Value
+#'
+#' A data frame with columns: `x`, `fill`, `xmin`, `xmax`, `ymin`, `ymax`,
+#' representing the position and size of each rectangle.
+#'
+#' @keywords internal
+.calculate_coordinates <- function(data) {
+  # Example: compute rectangles from x and y
+  result <- data |>
+    # Count combinations of X and Y
+    dplyr::count(.data$x, .data$fill, .drop = FALSE) |>
+    # Compute total for each X group
+    dplyr::mutate(
+      .by = .data$x,
+      x_total = sum(.data$n),
+      prop = .data$n / .data$x_total,
+      prop = dplyr::if_else(is.nan(.data$prop), 0, .data$prop)
+    ) |>
+    dplyr::arrange(dplyr::desc(.data$x_total), .data$x, .data$fill) |>
+    # Compute total sample size to turn counts into widths
+    dplyr::mutate(
+      N_total = dplyr::n(),
+      x_width = .data$x_total / .data$N_total
+    ) |>
+    # Convert counts to x widths
+    dplyr::mutate(
+      .by = .data$x,
+      x_width_last = dplyr::if_else(dplyr::row_number() == dplyr::n(), .data$x_width, 0)
+    ) |>
+    # Compute x-min/x-max for each group
+    dplyr::mutate(
+      xmin = cumsum(dplyr::lag(.data$x_width_last, default = 0)),
+      xmax = .data$xmin + .data$x_width
+    ) |>
+    # Compute y-min/y-max for stacked proportions
+    dplyr::mutate(
+      .by = .data$x,
+      ymin = c(0, utils::head(cumsum(.data$prop), -1)),
+      ymax = cumsum(.data$prop)
+    ) |>
+    dplyr::mutate(
+      xmin = .data$xmin / max(.data$xmax),
+      xmax = .data$xmax / max(.data$xmax),
+      xmin = dplyr::if_else(.data$n == 0, 0, .data$xmin + 0.005),
+      xmax = dplyr::if_else(.data$n == 0, 0, .data$xmax - 0.005),
+      ymin = dplyr::if_else(.data$n == 0, 0, .data$ymin + 0.005),
+      ymax = dplyr::if_else(.data$n == 0, 0, .data$ymax - 0.005)
+    ) |>
+    dplyr::select(.data$x, .data$fill, .data$xmin, .data$xmax, .data$ymin, .data$ymax, .n = .data$n)
+  result
+}
diff --git a/R/teal.modules.general.R b/R/teal.modules.general.R
index 07e7cf03c..d733f4fdc 100644
--- a/R/teal.modules.general.R
+++ b/R/teal.modules.general.R
@@ -7,8 +7,7 @@
 #' @import shiny
 #' @import teal
 #' @import teal.transform
-#' @importFrom dplyr %>%
-#'
+#' @importFrom dplyr %>% .data
 #'
 #' @name teal.modules.general
 #' @keywords internal
diff --git a/R/tm_g_association.R b/R/tm_g_association.R
index aadaee5e2..cbce628df 100644
--- a/R/tm_g_association.R
+++ b/R/tm_g_association.R
@@ -506,12 +506,12 @@ srv_tm_g_association <- function(id,
           substitute(
             expr = {
               plots <- plot_calls
-              plot <- gridExtra::arrangeGrob(plots[[1]], plots[[2]], ncol = 1)
+              plot <- gridExtra::arrangeGrob(grobs = plots, ncol = 1)
             },
             env = list(
               plot_calls = do.call(
                 "call",
-                c(list("list", ref_call), var_calls),
+                c(list("list", ref_call), unname(var_calls)),
                 quote = TRUE
               )
             )
diff --git a/R/tm_g_bivariate.R b/R/tm_g_bivariate.R
index e363858df..55e1b508c 100644
--- a/R/tm_g_bivariate.R
+++ b/R/tm_g_bivariate.R
@@ -724,7 +724,9 @@ srv_g_bivariate <- function(id,
       })
     )
 
-    plot_r <- reactive(req(decorated_output_q_facets())[["plot"]])
+    plot_r <- reactive({
+      req(decorated_output_q_facets())[["plot"]]
+    })
 
     pws <- teal.widgets::plot_with_settings_srv(
       id = "myplot",
@@ -768,7 +770,7 @@ bivariate_plot_call <- function(data_name,
     y <- if (is.call(y)) y else as.name(y)
   }
 
-  cl <- bivariate_ggplot_call(
+  bivariate_ggplot_call(
     x_class = x_class,
     y_class = y_class,
     freq = freq,
@@ -927,7 +929,13 @@ bivariate_ggplot_call <- function(x_class,
     )
     # Factor and character plots
   } else if (x_class == "factor" && y_class == "factor") {
-    stop("Categorical variables 'x' and 'y' are currently not supported.")
+    plot_call <- reduce_plot_call(
+      plot_call,
+      substitute(
+        teal.modules.general::geom_mosaic(ggplot2::aes(x = xval, fill = yval)),
+        env = list(xval = x, yval = y)
+      )
+    )
   } else {
     stop("x y type combination not allowed")
   }
diff --git a/man/dot-scale_x_mosaic.Rd b/man/dot-scale_x_mosaic.Rd
new file mode 100644
index 000000000..610143b36
--- /dev/null
+++ b/man/dot-scale_x_mosaic.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/geom_mosaic.R
+\name{.scale_x_mosaic}
+\alias{.scale_x_mosaic}
+\title{Determining scales for mosaics}
+\usage{
+.scale_x_mosaic(
+  breaks = unique,
+  minor_breaks = NULL,
+  labels = unique,
+  na.value = NA_real_,
+  position = "bottom",
+  ...
+)
+}
+\arguments{
+\item{breaks, labels, minor_breaks}{One of:
+\itemize{
+\item \code{NULL} for no breaks / labels.
+\item \code{\link[ggplot2:waiver]{ggplot2::waiver()}} for the default breaks / labels computed by the scale.
+\item A numeric / character vector giving the positions of the breaks / labels.
+\item A function.
+See \code{\link[ggplot2:scale_continuous]{ggplot2::scale_x_continuous()}} for more details.
+}}
+
+\item{na.value}{The value to be used for \code{NA} values.}
+
+\item{position}{For position scales, The position of the axis.
+left or right for y axes, top or bottom for x axes.}
+
+\item{...}{other arguments passed to \code{continuous_scale()}.}
+}
+\description{
+Determining scales for mosaics
+}
+\keyword{internal}
diff --git a/man/geom_mosaic.Rd b/man/geom_mosaic.Rd
new file mode 100644
index 000000000..dadbb4902
--- /dev/null
+++ b/man/geom_mosaic.Rd
@@ -0,0 +1,66 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/geom_mosaic.R
+\name{geom_mosaic}
+\alias{geom_mosaic}
+\alias{.calculate_coordinates}
+\title{Mosaic Rectangles Layer for ggplot2}
+\usage{
+geom_mosaic(
+  mapping = NULL,
+  data = NULL,
+  stat = "mosaic",
+  position = "identity",
+  ...,
+  na.rm = FALSE,
+  show.legend = TRUE,
+  inherit.aes = TRUE
+)
+
+.calculate_coordinates(data)
+}
+\arguments{
+\item{mapping}{Set of aesthetic mappings created by \code{aes()}. Must specify \code{x} and \code{fill}.}
+
+\item{data}{The data to be displayed in this layer.}
+
+\item{stat}{The statistical transformation to use on the data. Defaults to \code{"rects"}.}
+
+\item{position}{Position adjustment. Defaults to \code{"identity"}.}
+
+\item{...}{Other arguments passed to \code{layer()}.}
+
+\item{na.rm}{Logical. Should missing values be removed?}
+
+\item{show.legend}{Logical. Should this layer be included in the legends?}
+
+\item{inherit.aes}{Logical. If \code{FALSE}, overrides default aesthetics.}
+}
+\value{
+A ggplot2 layer that adds mosaic rectangles to the plot.
+}
+\description{
+Adds a mosaic-style rectangles layer to a ggplot, visualizing the
+joint distribution of categorical variables.
+Each rectangle's size reflects the proportion of observations for
+combinations of \code{x} and \code{fill}.
+}
+\section{Functions}{
+\itemize{
+\item \code{.calculate_coordinates()}: Computes the coordinates for rectangles in a mosaic plot based
+on combinations of \code{x} and \code{fill} variables.
+For each unique \code{x} and \code{fill}, calculates the proportional
+widths and heights, stacking rectangles within each \code{x} group.
+\subsection{Value}{
+
+A data frame with columns: \code{x}, \code{fill}, \code{xmin}, \code{xmax}, \code{ymin}, \code{ymax},
+representing the position and size of each rectangle.
+}
+
+}}
+\examples{
+df <- data.frame(RACE = c("Black", "White", "Black", "Asian"), SEX = c("M", "M", "F", "F"))
+library(ggplot2)
+ggplot(df) +
+  geom_mosaic(aes(x = RACE, fill = SEX))
+}
+\keyword{internal}
diff --git a/tests/testthat/test_bivariate_ggplot_call.R b/tests/testthat/test_bivariate_ggplot_call.R
index c85d4fd1b..618ada9e2 100644
--- a/tests/testthat/test_bivariate_ggplot_call.R
+++ b/tests/testthat/test_bivariate_ggplot_call.R
@@ -18,35 +18,24 @@ testthat::test_that("bivariate_ggplot_call with numerics", {
   )
 })
 
-testthat::test_that("bivariate_ggplot_call with factor, char, logical", {
-  error_message <- "Categorical variables 'x' and 'y' are currently not supported."
-  testthat::expect_error(
-    bivariate_ggplot_call("factor", "factor") %>% deparse(width.cutoff = 300),
-    error_message
-  )
-  testthat::expect_error(
-    bivariate_ggplot_call("logical", "factor") %>% deparse(width.cutoff = 300),
-    error_message
-  )
-  testthat::expect_error(
-    bivariate_ggplot_call("character", "factor") %>% deparse(width.cutoff = 300),
-    error_message
-  )
-  testthat::expect_error(
-    bivariate_ggplot_call("logical", "character") %>% deparse(width.cutoff = 300),
-    error_message
-  )
-  testthat::expect_error(
-    bivariate_ggplot_call("character", "logical") %>% deparse(width.cutoff = 300),
-    error_message
-  )
-  testthat::expect_error(
-    bivariate_ggplot_call("logical", "logical") %>% deparse(width.cutoff = 300),
-    error_message
-  )
-  testthat::expect_error(
-    bivariate_ggplot_call("character", "character") %>% deparse(width.cutoff = 300),
-    error_message
+testthat::describe("bivariate_ggplot_call with arguments:", {
+  possible_classes <- c("factor", "logical", "character")
+  comb <- expand.grid(a = possible_classes, b = possible_classes, stringsAsFactors = FALSE)
+  apply(
+    comb,
+    1,
+    function(x) {
+      it(sprintf("%s and %s", x[[1]], x[[2]]), {
+        testthat::expect_match(
+          deparse(
+            bivariate_ggplot_call(x[[1]], x[[2]], data_name = "ANL", x = "x", y = "y"),
+            width.cutoff = 300
+          ),
+          "teal.modules.general::geom_mosaic",
+          all = FALSE
+        )
+      })
+    }
   )
 })