diff --git a/DESCRIPTION b/DESCRIPTION
index 7cac1e67d..7bc5d1687 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -98,10 +98,11 @@ Suggests:
     rtdists,
     RWiener,
     sandwich,
+    scales,
     see (>= 0.11.0),
     survival,
     testthat (>= 3.2.1),
-    tinyplot,
+    tinyplot (>= 0.6.0),
     tinytable,
     vdiffr,
     withr
diff --git a/R/tinyplot.R b/R/tinyplot.R
index b1984bba0..e4f1ad455 100644
--- a/R/tinyplot.R
+++ b/R/tinyplot.R
@@ -1,33 +1,55 @@
 #' @rdname visualisation_recipe.estimate_predicted
-#' @param theme A character string specifying the theme to use for the plot.
-#' Defaults to `"tufte"`. For other options please see [`tinyplot::tinytheme()`].
-#' Use `NULL` if no theme should be applied.
+#' @param type The type of `tinyplot` visualization. It is recommended that
+#' users leave as `NULL` (the default), in which case the plot type will be
+#' determined automatically by the underlying `modelbased` object.
+#' @param dodge Dodge value for grouped plots. If `NULL` (the default), then
+#' the dodging behavior is determined by the number of groups and
+#' `getOption("modelbased_tinyplot_dodge")`.
+#' @param ... Other arguments passed to \code{\link[tinyplot]{tinyplot}}.
 #'
 #' @examplesIf all(insight::check_if_installed(c("tinyplot", "marginaleffects"), quietly = TRUE))
 #' # ==============================================
 #' # tinyplot
 #' # ==============================================
 #' \donttest{
+#' library(tinyplot)
 #' data(efc, package = "modelbased")
 #' efc <- datawizard::to_factor(efc, c("e16sex", "c172code", "e42dep"))
 #' m <- lm(neg_c_7 ~ e16sex + c172code + barthtot, data = efc)
 #'
 #' em <- estimate_means(m, "c172code")
-#' tinyplot::plt(em)
+#' plt(em)
 #'
+#' # pass additional tinyplot arguments for customization, e.g.
+#' plt(em, theme = "classic")
+#' plt(em, theme = "classic", flip = TRUE)
+#' # etc.
+#'
+#' # Aside: use tinyplot::tinytheme() to set a persistent theme
+#' tinytheme("classic")
+#'
+#' # continuous variable example
 #' em <- estimate_means(m, "barthtot")
-#' tinyplot::plt(em)
+#' plt(em)
 #'
+#' # grouped example
 #' m <- lm(neg_c_7 ~ e16sex * c172code + e42dep, data = efc)
 #' em <- estimate_means(m, c("e16sex", "c172code"))
-#' tinyplot::plt(em)
+#' plt(em)
+#'
+#' # use plt_add (alias tinyplot_add) to add layers
+#' plt_add(type = "l", lty = 2)
+#'
+#' # Reset to default theme
+#' tinytheme()
 #' }
 #' @exportS3Method tinyplot::tinyplot
 tinyplot.estimate_means <- function(
   x,
+  type = NULL,
+  dodge = NULL,
   show_data = FALSE,
   numeric_as_discrete = NULL,
-  theme = "tufte",
   ...
 ) {
   insight::check_if_installed("tinyplot")
@@ -48,9 +70,8 @@ tinyplot.estimate_means <- function(
   data <- aes$data
   aes <- aes$aes
 
-  # save additional arguments, once for theming and once for the plot
+  # save additional arguments, will pass via do.call to tinyplot
   dots <- list(...)
-  theme_dots <- dots
 
   # preparation of settings / arguments ----------------------------------
 
@@ -73,6 +94,11 @@ tinyplot.estimate_means <- function(
     }
   }
 
+  # type placeholder
+  if (!is.null(type)) {
+    aes$type <- type
+  }
+
   # handle non-standard plot types -------------------------------
 
   if (aes$type == "grouplevel") {
@@ -110,18 +136,49 @@ tinyplot.estimate_means <- function(
   # Set dodge value for grouped point or pointrange plots.
   # The value 0.07 was chosen to reduce overlap in this context; adjust via
   # option if needed.
-  dodge_value <- getOption("modelbased_tinyplot_dodge", 0.07)
-  if (!is.null(aes$color) && aes$type %in% c("pointrange", "point")) {
+
+  dodge_value <- if (!is.null(dodge)) {
+    dodge
+  } else {
+    getOption("modelbased_tinyplot_dodge", 0.07)
+  }
+  if (
+    !is.null(aes$color) &&
+      aes$type %in% c("pointrange", "point", "l", "errorbar", "ribbon")
+  ) {
     dots$dodge <- dodge_value
   }
 
-  ## TODO: legend labels?
   ## TODO: show residuals?
 
   # x/y labels --------------------------------
   dots$xlab <- aes$labs$x
   dots$ylab <- aes$labs$y
 
+  # legend labels --------------------------------
+
+  # we also need to account for custom legend options passed through dots
+  if (is.null(dots$legend)) {
+    dots$legend = list(title = aes$labs$colour)
+  } else if (inherits(dots$legend, "list")) {
+    if (!("title" %in% names(dots$legend))) {
+      dots$legend = utils::modifyList(
+        dots$legend,
+        list(title = aes$labs$colour),
+        keep.null = TRUE
+      )
+    }
+  } else if (!isFALSE(dots$legend)) {
+    dots$legend = tryCatch(
+      utils::modifyList(
+        as.list(dots$legend),
+        list(title = aes$labs$colour),
+        keep.null = TRUE
+      ),
+      error = function(e) dots$legend
+    )
+  }
+
   # add aesthetics to the plot description
   plot_args <- insight::compact_list(c(
     list(plot_description, data = data, type = aes$type),
@@ -129,12 +186,6 @@ tinyplot.estimate_means <- function(
     dots
   ))
 
-  # default theme
-  if (!is.null(theme)) {
-    theme_dots[c(elements, "facet", "xlab", "ylab", "flip")] <- NULL
-    do.call(tinyplot::tinytheme, c(list(theme = theme), theme_dots))
-  }
-
   # add data points if requested --------------------------------
 
   if (show_data) {
diff --git a/R/visualisation_recipe.R b/R/visualisation_recipe.R
index 86336be32..20ad906f7 100644
--- a/R/visualisation_recipe.R
+++ b/R/visualisation_recipe.R
@@ -45,7 +45,7 @@
 #' @param point,line,pointrange,ribbon,facet,grid Additional
 #' aesthetics and parameters for the geoms (see customization example).
 #' @param ... Arguments passed from `plot()` to `visualisation_recipe()`, or
-#' to `tinyplot()` and `tinytheme()` if you use that method.
+#' to `tinyplot()` if you use that method.
 #'
 #' @details There are two options to remove the confidence bands or errors bars
 #' from the plot. To remove error bars, simply set the `pointrange` geom to
diff --git a/man/visualisation_recipe.estimate_predicted.Rd b/man/visualisation_recipe.estimate_predicted.Rd
index 64d73aabf..ca63f7e5f 100644
--- a/man/visualisation_recipe.estimate_predicted.Rd
+++ b/man/visualisation_recipe.estimate_predicted.Rd
@@ -15,9 +15,10 @@
 
 \method{tinyplot}{estimate_means}(
   x,
+  type = NULL,
+  dodge = NULL,
   show_data = FALSE,
   numeric_as_discrete = NULL,
-  theme = "tufte",
   ...
 )
 
@@ -60,7 +61,15 @@
 \item{x}{A modelbased object.}
 
 \item{...}{Arguments passed from \code{plot()} to \code{visualisation_recipe()}, or
-to \code{tinyplot()} and \code{tinytheme()} if you use that method.}
+to \code{tinyplot()} if you use that method.}
+
+\item{type}{The type of \code{tinyplot} visualization. It is recommended that
+users leave as \code{NULL} (the default), in which case the plot type will be
+determined automatically by the underlying \code{modelbased} object.}
+
+\item{dodge}{Dodge value for grouped plots. If \code{NULL} (the default), then
+the dodging behavior is determined by the number of groups and
+\code{getOption("modelbased_tinyplot_dodge")}.}
 
 \item{show_data}{Logical, if \code{TRUE}, display the "raw" data as a background
 to the model-based estimation. This argument will be ignored for plotting
@@ -76,10 +85,6 @@ predictor. Use \code{FALSE} to always use continuous color scales for numeric
 predictors. It is possible to set a global default value using \code{options()},
 e.g. \code{options(modelbased_numeric_as_discrete = 10)}.}
 
-\item{theme}{A character string specifying the theme to use for the plot.
-Defaults to \code{"tufte"}. For other options please see \code{\link[tinyplot:tinytheme]{tinyplot::tinytheme()}}.
-Use \code{NULL} if no theme should be applied.}
-
 \item{show_residuals}{Logical, if \code{TRUE}, display residuals of the model
 as a background to the model-based estimation. Residuals will be computed
 for the predictors in the data grid, using \code{\link[=residualize_over_grid]{residualize_over_grid()}}.}
@@ -148,19 +153,36 @@ when using \code{tinyplot::plt()}. Should be a number between \code{0} and \code
 # tinyplot
 # ==============================================
 \donttest{
+library(tinyplot)
 data(efc, package = "modelbased")
 efc <- datawizard::to_factor(efc, c("e16sex", "c172code", "e42dep"))
 m <- lm(neg_c_7 ~ e16sex + c172code + barthtot, data = efc)
 
 em <- estimate_means(m, "c172code")
-tinyplot::plt(em)
+plt(em)
+
+# pass additional tinyplot arguments for customization, e.g.
+plt(em, theme = "classic")
+plt(em, theme = "classic", flip = TRUE)
+# etc.
 
+# Aside: use tinyplot::tinytheme() to set a persistent theme
+tinytheme("classic")
+
+# continuous variable example
 em <- estimate_means(m, "barthtot")
-tinyplot::plt(em)
+plt(em)
 
+# grouped example
 m <- lm(neg_c_7 ~ e16sex * c172code + e42dep, data = efc)
 em <- estimate_means(m, c("e16sex", "c172code"))
-tinyplot::plt(em)
+plt(em)
+
+# use plt_add (alias tinyplot_add) to add layers
+plt_add(type = "l", lty = 2)
+
+# Reset to default theme
+tinytheme()
 }
 \dontshow{\}) # examplesIf}
 \dontshow{if (all(insight::check_if_installed(c("marginaleffects", "see", "ggplot2"), quietly = TRUE)) && getRversion() >= "4.1.0") withAutoprint(\{ # examplesIf}
diff --git a/vignettes/plotting_tinyplot.Rmd b/vignettes/plotting_tinyplot.Rmd
index 6981b72b7..138818c6a 100644
--- a/vignettes/plotting_tinyplot.Rmd
+++ b/vignettes/plotting_tinyplot.Rmd
@@ -25,7 +25,7 @@ knitr::opts_chunk$set(
 options(modelbased_join_dots = FALSE)
 options(modelbased_select = "minimal")
 
-pkgs <- c("marginaleffects", "tinyplot")
+pkgs <- c("marginaleffects", "tinyplot", "scales")
 if (!all(insight::check_if_installed(pkgs, quietly = TRUE))) {
   knitr::opts_chunk$set(eval = FALSE)
 }
@@ -43,6 +43,7 @@ The simplest case is possibly plotting one categorical predictor. Predicted valu
 ```{r}
 library(modelbased)
 library(tinyplot)
+tinytheme("classic", palette.qualitative = "Tableau 10")
 
 data(efc, package = "modelbased")
 efc <- datawizard::to_factor(efc, c("e16sex", "c172code", "e42dep"))
@@ -51,6 +52,21 @@ m <- lm(neg_c_7 ~ e16sex + c172code + barthtot, data = efc)
 estimate_means(m, "c172code") |> plt()
 ```
 
+In general, plots can be further modified using functions or arguments from the **tinyplot** package.
+
+```{r}
+estimate_means(m, "c172code") |>
+  plt(type = "errorbar", flip = TRUE)
+plt_add(type = "l", lty = 2)
+```
+
+**Pro-tip:** You can pass a labeling function to wrap long axis labels. Here we use one from the `scales` package.
+
+```{r}
+estimate_means(m, "c172code") |>
+  plt(xaxl = scales::label_wrap(20), flip = TRUE)
+```
+
 ## One predictor - numeric
 
 For numeric predictors, the range of predictions at different values of the focal predictor are plotted, the uncertainty is displayed as confidence band.
@@ -68,6 +84,13 @@ m <- lm(neg_c_7 ~ e16sex * c172code + e42dep, data = efc)
 estimate_means(m, c("e16sex", "c172code")) |> plt()
 ```
 
+Again, you can layer on top of this plot using standard **tinyplot** functions and arguments.
+
+```{r}
+estimate_means(m, c("e16sex", "c172code")) |> plt()
+plt_add(type = "l", lty = 2)
+```
+
 ## Two predictors - numeric * categorical
 
 For two predictors, where the first is numeric and the second categorical, range of predictions including confidence bands are shown, with the different levels of the second (categorical) predictor mapped to colors again.
@@ -77,14 +100,14 @@ m <- lm(neg_c_7 ~ barthtot * c172code + e42dep, data = efc)
 estimate_means(m, c("barthtot", "c172code")) |> plt()
 ```
 
-In general, plots can be further modified using functions or arguments from the **tinyplot** package. Thereby, other themes, color scales, faceting and so on, can be applied.
+One potentially useful customization for these numeric * categorical cases is mapping predictive groups by facets (in addition to, or instead of colors). Below we make an additional tweak by wrapping the facet names, since these are rather long.
 
 ```{r}
 estimate_means(m, c("barthtot", "c172code")) |>
-  plt(facet = ~c172code)
-
-estimate_means(m, c("barthtot", "c172code")) |>
-  plt(palette = "okabe")
+  within({
+    c172code = gsub(" of education$", "\nof education", c172code)
+  }) |> 
+  plt(facet = "by", legend = FALSE)
 ```
 
 ## Two predictors - categorical * numeric
@@ -122,10 +145,17 @@ estimate_means(
     "barthtot = c(30, 50, 80)"
   )
 ) |> plt()
-
 estimate_means(m, c("c172code", "barthtot = [fivenum]")) |> plt()
 ```
 
+One aesthetic issue in the preceding plots is the fact that the middle x-axis category is missing (hidden) due to space limitations. Again, this is a common problem when we have several discrete categories and long label strings. One option is to increase the horizontal spacing by moving the legend (e.g., `legend = bottom!`). But a more sure-fire way to ensure that all tick labels are printed is by using a labelling function that wraps the
+text and/or flipping the plot.
+
+```{r}
+estimate_means(m, c("c172code", "barthtot = [fivenum]")) |>
+  plt(xaxl = scales::label_wrap(20), flip = TRUE)
+```
+
 ## Three numeric predictors
 
 The default plot-setting for three numeric predictors can be rather confusing.
@@ -135,12 +165,19 @@ m <- lm(neg_c_7 ~ c12hour * barthtot * c160age, data = efc)
 estimate_means(m, c("c12hour", "barthtot", "c160age")) |> plt()
 ```
 
-Instead, it is recommended to use `length`, create a "reference grid", or again specify meaningful values directly in the `by` argument.
+Instead, it is recommended to use `length`, create a "reference grid", or again specify meaningful values directly in the `by` argument. Note that this will have the ancillary effect of generating facets by the third variable (here: "cs160age", i.e. carer age).
 
 ```{r}
-estimate_means(m, c("c12hour", "barthtot", "c160age"), length = 2) |> plt()
-
-estimate_means(m, c("c12hour", "barthtot", "c160age"), range = "grid") |> plt()
+estimate_means(m, c("c12hour", "barthtot", "c160age"), length = 2) |>
+  plt(
+    main = "Effect of care on elder outcomes",
+    sub = "Facets denote representative carer ages"
+  )
+estimate_means(m, c("c12hour", "barthtot", "c160age"), range = "grid") |>
+  plt(
+    main = "Effect of care on elder outcomes",
+    sub = "Facets denote representative carer ages (mean +/- 1 SD)"
+  )
 ```
 
 ## Three categorical predictors
@@ -152,11 +189,19 @@ m <- lm(neg_c_7 ~ e16sex * c172code * e42dep, data = efc)
 estimate_means(m, c("e16sex", "c172code", "e42dep")) |> plt()
 ```
 
+Again, though we can improve the final aesthetic with a few tweaks.
+
+```{r}
+estimate_means(m, c("e16sex", "c172code", "e42dep")) |>
+  plt(dodge = 0.02, theme = "clean2")
+```
+
 ## Smooth plots
 
 Remember that by default a range of ten values is chosen for numeric focal predictors. While this mostly works well for plotting linear relationships, plots may look less smooth for certain models that involve quadratic or cubic terms, or splines, or for instance if you have GAMs.
 
 ```{r}
+tinytheme("classic", palette.qualitative = "Tableau 10")
 m <- lm(neg_c_7 ~ e16sex * c12hour + e16sex * I(c12hour^2), data = efc)
 estimate_means(m, c("c12hour", "e16sex")) |> plt()
 ```
@@ -167,6 +212,11 @@ In this case, simply increase the number of representative values by setting `le
 estimate_means(m, c("c12hour", "e16sex"), length = 200) |> plt()
 ```
 
+```{r}
+# reset theme
+tinytheme()
+```
+
 ```{r echo=FALSE}
 # reset options
 options(