diff --git a/NEWS.md b/NEWS.md
index b3e8922c..4412a35c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,6 @@
 # bayesplot (development version)
 
+* New `show_marginal` argument to `ppd_*()` functions to show the PPD - the marginal predictive distribution by @mattansb (#425)
 * Documentation added for all exported `*_data()` functions (#209)
 * Improved documentation for `binwidth`, `bins`, and `breaks` arguments to clarify they are passed to `ggplot2::geom_area()` and `ggdist::stat_dots()` in addition to `ggplot2::geom_histogram()`
 * Improved documentation for `freq` argument to clarify it applies to frequency polygons in addition to histograms
diff --git a/R/helpers-gg.R b/R/helpers-gg.R
index 953a51d4..a53de233 100644
--- a/R/helpers-gg.R
+++ b/R/helpers-gg.R
@@ -125,22 +125,26 @@ scale_fill_ppc <-
 
 scale_color_ppd <-
   function(name = NULL,
-           values = get_color("mh"),
-           labels = ypred_label(),
+           values = NULL,
+           labels = NULL,
            ...) {
-    scale_color_ppc(name = name,
-                    values = values,
-                    labels = labels,
-                    ...)
+    scale_color_ppc(
+      name = name,
+      values = values %||% setNames(get_color(c("dh", "mh")), nm = c("PPD", "ypred")),
+      labels = labels %||% ypred_label(),
+      ...
+    )
   }
 
 scale_fill_ppd <-
   function(name = NULL,
-           values = get_color("m"),
-           labels = ypred_label(),
+           values = NULL,
+           labels = NULL,
            ...) {
-    scale_fill_ppc(name = name,
-                   values = values,
-                   labels = labels,
-                   ...)
+    scale_fill_ppc(
+      name = name,
+      values = values %||% setNames(get_color(c("d", "m")), nm = c("PPD", "ypred")),
+      labels = labels %||% ypred_label(),
+      ...
+    )
   }
diff --git a/R/helpers-ppc.R b/R/helpers-ppc.R
index e268b315..267f1126 100644
--- a/R/helpers-ppc.R
+++ b/R/helpers-ppc.R
@@ -594,4 +594,9 @@ u_scale <- function(x) {
 create_rep_ids <- function(ids) paste('italic(y)[rep] (', ids, ")")
 y_label <- function() expression(italic(y))
 yrep_label <- function() expression(italic(y)[rep])
-ypred_label <- function() expression(italic(y)[pred])
+ypred_label <- function() {
+  c(
+    PPD = "PPD",
+    ypred = expression(italic(y)[pred])
+  )
+}
diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 70c4e5a6..ebd5eb43 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -10,6 +10,7 @@
 #'
 #' @template args-ypred
 #' @inheritParams PPC-distributions
+#' @param show_marginal Plot the marginal PPD along with the `ypred`s.
 #'
 #' @template details-binomial
 #' @template return-ggplot-or-data
@@ -19,6 +20,7 @@
 #' color_scheme_set("brightblue")
 #' preds <- example_yrep_draws()
 #' ppd_dens_overlay(ypred = preds[1:50, ])
+#' ppd_dens_overlay(ypred = preds[1:50, ], show_marginal = TRUE)
 #' ppc_dens_overlay(y = example_y_data(), yrep = preds[1:50, ])
 #'
 NULL
@@ -37,25 +39,45 @@ ppd_data <- function(ypred, group = NULL) {
 #' @rdname PPD-distributions
 #' @export
 ppd_dens_overlay <-
-    function(ypred,
-             ...,
-             size = 0.25,
-             alpha = 0.7,
-             trim = FALSE,
-             bw = "nrd0",
-             adjust = 1,
-             kernel = "gaussian",
-             bounds = NULL,
-             n_dens = 1024) {
-      check_ignored_arguments(...)
-      bounds <- validate_density_bounds(bounds)
+  function(ypred,
+           show_marginal = FALSE,
+           ...,
+           size = 0.25,
+           alpha = 0.7,
+           trim = FALSE,
+           bw = "nrd0",
+           adjust = 1,
+           kernel = "gaussian",
+           bounds = NULL,
+           n_dens = 1024) {
+    check_ignored_arguments(...)
+    bounds <- validate_density_bounds(bounds)
 
-      data <- ppd_data(ypred)
-      ggplot(data, mapping = aes(x = .data$value)) +
-        overlay_ppd_densities(
+    data <- ppd_data(ypred)
+    p <- ggplot(data, mapping = aes(x = .data$value)) +
+      overlay_ppd_densities(
         mapping = aes(group = .data$rep_id, color = "ypred"),
         linewidth = size,
         alpha = alpha,
+        trim = trim,
+        bw = bw,
+        adjust = adjust,
+        kernel = kernel,
+        bounds = bounds,
+        n = n_dens
+      ) +
+      bayesplot_theme_get() +
+      dont_expand_axes() +
+      yaxis_title(FALSE) +
+      xaxis_title(FALSE) +
+      yaxis_text(FALSE) +
+      yaxis_ticks(FALSE)
+
+    if (isTRUE(show_marginal)) {
+      p +
+        overlay_ppd_densities(
+          mapping = aes(color = "PPD"),
+          linewidth = 1,
           trim = trim,
           bw = bw,
           adjust = adjust,
@@ -63,18 +85,22 @@ ppd_dens_overlay <-
           bounds = bounds,
           n = n_dens
         ) +
-      scale_color_ppd(
+        scale_color_ppd(
+          labels = ypred_label(),
+          values = setNames(get_color(c("d", "m")), nm = c("PPD", "ypred")),
+          guide = guide_legend(
+            override.aes = list(size = 2 * size, alpha = 1))
+        )
+    } else {
+      p + scale_color_ppd(
         values = get_color("m"),
-        guide = guide_legend( # in case user turns legend back on
+        # in case user turns legend back on
+        guide = guide_legend(
           override.aes = list(size = 2 * size, alpha = 1))
       ) +
-      bayesplot_theme_get() +
-      dont_expand_axes() +
-      yaxis_title(FALSE) +
-      xaxis_title(FALSE) +
-      yaxis_text(FALSE) +
-      yaxis_ticks(FALSE) +
-      legend_none()
+        legend_none()
+    }
+
   }
 
 
@@ -82,6 +108,7 @@ ppd_dens_overlay <-
 #' @export
 ppd_ecdf_overlay <-
   function(ypred,
+           show_marginal= FALSE,
            ...,
            discrete = FALSE,
            pad = TRUE,
@@ -90,7 +117,7 @@ ppd_ecdf_overlay <-
     check_ignored_arguments(...)
 
     data <- ppd_data(ypred)
-    ggplot(data, mapping = aes(x = .data$value)) +
+    p <- ggplot(data, mapping = aes(x = .data$value)) +
       hline_at(
         c(0, 0.5, 1),
         linewidth = c(0.2, 0.1, 0.2),
@@ -104,16 +131,34 @@ ppd_ecdf_overlay <-
         alpha = alpha,
         pad = pad
       ) +
-      scale_color_ppd(
-        values = get_color("m"),
-        guide = guide_legend( # in case user turns legend back on
-          override.aes = list(linewidth = 2 * size, alpha = 1))
-      ) +
       scale_y_continuous(breaks = c(0, 0.5, 1)) +
       bayesplot_theme_get() +
       yaxis_title(FALSE) +
-      xaxis_title(FALSE) +
-      legend_none()
+      xaxis_title(FALSE)
+
+    if (isTRUE(show_marginal)) {
+      p +
+        stat_ecdf(
+          mapping = aes(color = "PPD"),
+          geom = if (discrete) "step" else "line",
+          linewidth = 1,
+          pad = pad
+        ) +
+        scale_color_ppd(
+          labels = ypred_label(),
+          values = setNames(get_color(c("d", "m")), nm = c("PPD", "ypred")),
+          guide = guide_legend(
+            override.aes = list(size = 2 * size, alpha = 1))
+        )
+    } else {
+      p +
+        scale_color_ppd(
+          values = get_color("m"),
+          guide = guide_legend( # in case user turns legend back on
+            override.aes = list(linewidth = 2 * size, alpha = 1))
+        ) +
+        legend_none()
+    }
   }
 
 
@@ -121,6 +166,7 @@ ppd_ecdf_overlay <-
 #' @export
 ppd_dens <-
     function(ypred,
+             show_marginal = FALSE,
              ...,
              trim = FALSE,
              size = 0.5,
@@ -130,73 +176,105 @@ ppd_dens <-
       bounds <- validate_density_bounds(bounds)
 
       data <- ppd_data(ypred)
-      ggplot(data, mapping = aes(
-        x = .data$value,
-        color = "ypred",
-        fill = "ypred"
-      )) +
+      p <- ggplot(data, mapping = aes(.data$value)) +
         geom_density(
+          aes(color = "ypred",
+              fill = "ypred"),
           linewidth = size,
           alpha = alpha,
           trim = trim,
           bounds = bounds
         ) +
-      scale_color_ppd() +
-      scale_fill_ppd() +
-      bayesplot_theme_get() +
-      facet_wrap_parsed("rep_label") +
-      force_axes_in_facets() +
-      dont_expand_y_axis() +
-      legend_none() +
-      yaxis_text(FALSE) +
-      yaxis_title(FALSE) +
-      yaxis_ticks(FALSE) +
-      xaxis_title(FALSE) +
-      facet_text(FALSE)
-  }
+        bayesplot_theme_get() +
+        facet_wrap_parsed("rep_label") +
+        force_axes_in_facets() +
+        dont_expand_y_axis() +
+        yaxis_text(FALSE) +
+        yaxis_title(FALSE) +
+        yaxis_ticks(FALSE) +
+        xaxis_title(FALSE) +
+        facet_text(FALSE) +
+        scale_color_ppd() +
+        scale_fill_ppd()
+
+      if (isTRUE(show_marginal)) {
+        data2 <- transform(data, rep_label = "PPD")
+
+        p +
+          geom_density(
+            aes(color = "PPD",
+                fill = "PPD"),
+            linewidth = 1,
+            trim = trim,
+            bounds = bounds,
+            data = data2
+          )
+      } else {
+        p +
+          legend_none()
+      }
+    }
 
 
 #' @rdname PPD-distributions
 #' @export
 ppd_hist <-
   function(ypred,
+           show_marginal = FALSE,
            ...,
            binwidth = NULL,
            bins = NULL,
            breaks = NULL,
-           freq = TRUE) {
+           freq = !show_marginal) {
     check_ignored_arguments(...)
 
     data <- ppd_data(ypred)
-    ggplot(data, mapping = set_hist_aes(
-      freq,
-      color = "ypred",
-      fill = "ypred"
-    )) +
+    p <- ggplot(data, mapping = set_hist_aes(freq)) +
       geom_histogram(
+        aes(color = "ypred",
+            fill = "ypred"),
         linewidth = 0.25,
         binwidth = binwidth,
         bins = bins,
         breaks = breaks
       ) +
-      scale_color_ppd() +
-      scale_fill_ppd() +
       bayesplot_theme_get() +
       facet_wrap_parsed("rep_label") +
       force_axes_in_facets() +
       dont_expand_y_axis() +
-      legend_none() +
       yaxis_text(FALSE) +
       yaxis_title(FALSE) +
       yaxis_ticks(FALSE) +
       xaxis_title(FALSE) +
-      facet_text(FALSE)
+      facet_text(FALSE) +
+      scale_color_ppd() +
+      scale_fill_ppd()
+
+    if (isTRUE(show_marginal)) {
+      data2 <- transform(data, rep_label = "PPD")
+
+      p +
+        geom_histogram(
+          aes(color = "PPD",
+              fill = "PPD"),
+          linewidth = 1,
+          binwidth = binwidth,
+          bins = bins,
+          breaks = breaks,
+          data = data2
+        )
+
+    } else {
+      p +
+        legend_none()
+    }
   }
 
 #' @rdname PPD-distributions
 #' @export
 ppd_dots <-
   function(ypred,
+           show_marginal = FALSE,
            ...,
            binwidth = NA,
            quantiles = 100,
@@ -206,28 +284,44 @@ ppd_dots <-
     suggested_package("ggdist")
 
     data <- ppd_data(ypred)
-    ggplot(data, mapping = set_hist_aes(
-      freq,
-      color = "ypred",
-      fill = "ypred"
-    )) +
+
+    p <- ggplot(data, mapping = set_hist_aes(freq)) +
       ggdist::stat_dots(
+        aes(color = "ypred",
+            fill = "ypred"),
         binwidth = binwidth,
         quantiles = quantiles,
         ...
       ) +
-      scale_color_ppd() +
-      scale_fill_ppd() +
       bayesplot_theme_get() +
       facet_wrap_parsed("rep_label") +
       force_axes_in_facets() +
       dont_expand_y_axis() +
-      legend_none() +
       yaxis_text(FALSE) +
       yaxis_title(FALSE) +
       yaxis_ticks(FALSE) +
       xaxis_title(FALSE) +
-      facet_text(FALSE)
+      facet_text(FALSE) +
+      scale_color_ppd() +
+      scale_fill_ppd()
+
+    if (isTRUE(show_marginal)) {
+      data2 <- transform(data, rep_label = "PPD")
+
+      p +
+        ggdist::stat_dots(
+          aes(color = "PPD",
+              fill = "PPD"),
+          data = data2,
+          binwidth = binwidth,
+          quantiles = quantiles,
+          ...
+        )
+
+    } else {
+      p +
+        legend_none()
+    }
   }
 
 
@@ -235,10 +329,11 @@ ppd_dots <-
 #' @export
 ppd_freqpoly <-
   function(ypred,
+           show_marginal = FALSE,
            ...,
            binwidth = NULL,
            bins = NULL,
-           freq = TRUE,
+           freq = !show_marginal,
            size = 0.5,
            alpha = 1) {
 
@@ -249,12 +344,10 @@ ppd_freqpoly <-
     }
 
     data <- ppd_data(ypred, group = dots$group)
-    ggplot(data, mapping = set_hist_aes(
-      freq,
-      color = "ypred",
-      fill = "ypred"
-    )) +
+    p <- ggplot(data, mapping = set_hist_aes(freq)) +
       geom_area(
+        aes(color = "ypred",
+            fill = "ypred"),
         stat = "bin",
         binwidth = binwidth,
         bins = bins,
@@ -262,8 +355,6 @@ ppd_freqpoly <-
         alpha = alpha
       ) +
       facet_wrap_parsed("rep_label") +
-      scale_color_ppd() +
-      scale_fill_ppd() +
       bayesplot_theme_get() +
       force_axes_in_facets() +
       dont_expand_y_axis() +
@@ -272,7 +363,28 @@ ppd_freqpoly <-
       yaxis_ticks(FALSE) +
       xaxis_title(FALSE) +
       facet_text(FALSE) +
-      legend_none()
+      scale_color_ppd() +
+      scale_fill_ppd()
+
+
+    if (isTRUE(show_marginal)) {
+      data2 <- transform(data, rep_label = "PPD")
+
+      p +
+        geom_area(
+          aes(color = "PPD",
+              fill = "PPD"),
+          data = data2,
+          stat = "bin",
+          binwidth = binwidth,
+          bins = bins,
+          linewidth = 1,
+        )
+
+    } else {
+      p +
+        legend_none()
+    }
   }
 
 
@@ -281,10 +393,11 @@ ppd_freqpoly <-
 ppd_freqpoly_grouped <-
   function(ypred,
            group,
+           show_marginal = FALSE,
            ...,
            binwidth = NULL,
            bins = NULL,
-           freq = TRUE,
+           freq = !show_marginal,
            size = 0.5,
            alpha = 1) {
 
@@ -307,6 +420,7 @@ ppd_freqpoly_grouped <-
 #' @export
 ppd_boxplot <-
   function(ypred,
+           show_marginal = FALSE,
            ...,
            notch = TRUE,
            size = 0.5,
@@ -314,13 +428,13 @@ ppd_boxplot <-
     check_ignored_arguments(...)
 
     data <- ppd_data(ypred)
-    ggplot(data, mapping = aes(
+    p <- ggplot(data, mapping = aes(
       x = .data$rep_label,
-      y = .data$value,
-      color = "ypred",
-      fill = "ypred"
+      y = .data$value
     )) +
       geom_boxplot(
+        aes(color = "ypred",
+            fill = "ypred"),
         notch = notch,
         linewidth = size,
         alpha = alpha,
@@ -335,8 +449,24 @@ ppd_boxplot <-
       yaxis_title(FALSE) +
       xaxis_ticks(FALSE) +
       xaxis_text(FALSE) +
-      xaxis_title(FALSE) +
-      legend_none()
+      xaxis_title(FALSE)
+
+    if (isTRUE(show_marginal)) {
+      p +
+        geom_boxplot(
+          aes(x = "PPD",
+              color = "PPD",
+              fill = "PPD"),
+          notch = notch,
+          linewidth = 1,
+          outlier.color = get_color("lh"),
+          outlier.alpha = 2/3,
+          outlier.size = 1
+        )
+
+    } else {
+      p + legend_none()
+    }
   }
 
 
diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index d7146405..412d6f1e 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -13,6 +13,7 @@
 #'
 #' @template args-ypred
 #' @inheritParams PPC-test-statistics
+#' @inheritParams PPD-distributions
 #'
 #' @template details-binomial
 #' @template return-ggplot-or-data
@@ -22,11 +23,13 @@
 #' yrep <- example_yrep_draws()
 #' ppd_stat(yrep)
 #' ppd_stat(yrep, stat = "sd") + legend_none()
+#' ppd_stat(yrep, show_marginal = TRUE)
 #'
 #' # use your own function for the 'stat' argument
 #' color_scheme_set("brightblue")
 #' q25 <- function(y) quantile(y, 0.25)
 #' ppd_stat(yrep, stat = "q25") # legend includes function name
+#' ppd_stat(yrep, stat = "q25", show_marginal = TRUE)
 NULL
 
 #' @rdname PPD-test-statistics
@@ -34,12 +37,13 @@ NULL
 ppd_stat <-
   function(ypred,
            stat = "mean",
+           show_marginal = FALSE,
            ...,
            discrete = FALSE,
            binwidth = NULL,
            bins = NULL,
            breaks = NULL,
-           freq = TRUE) {
+           freq = !show_marginal) {
     stopifnot(length(stat) == 1)
     dots <- list(...)
     if (!from_grouped(dots)) {
@@ -50,34 +54,49 @@ ppd_stat <-
     data <- ppd_stat_data(
       ypred = ypred,
       group = dots$group,
-      stat = match.fun(stat)
+      stat = match.fun(stat),
+      show_marginal = show_marginal
     )
+    data$type <- ifelse(grepl("ypred", data$variable), "ypred", "PPD")
+
     graph <- ggplot(data, mapping = set_hist_aes(
       freq,
-      color = "ypred",
-      fill = "ypred"
+      color = .data$type,
+      fill = .data$type
     ))
-    graph <- graph + if (discrete) {
-      geom_bar(
-        color = get_color("lh"),
-        linewidth = 0.25,
-        na.rm = TRUE,
-        position = "identity",
-      )
+
+    if (discrete) {
+      graph <- graph +
+        geom_bar(
+          linewidth = 0.25,
+          na.rm = TRUE,
+          position = "identity",
+        )
+    } else {
+      graph <- graph +
+        geom_histogram(
+          data = data[data$type != "PPD",],
+          linewidth = 0.25,
+          na.rm = TRUE,
+          binwidth = binwidth,
+          bins = bins,
+          breaks = breaks
+        )
+      if (isTRUE(show_marginal)) {
+        graph <- graph +
+          geom_vline(
+            aes(xintercept = .data$value, color = .data$type),
+            data = data[data$type == "PPD",],
+            key_glyph = "rect",
+            linewidth = 2
+          )
+      }
     }
-    else {
-      geom_histogram(
-        linewidth = 0.25,
-        na.rm = TRUE,
-        binwidth = binwidth,
-        bins = bins,
-        breaks = breaks
-      ) }
-      graph +
-      scale_color_ppd(guide = "none") +
-      scale_fill_ppd(labels = Typred_label(), guide = guide_legend(
-        title = stat_legend_title(stat, deparse(substitute(stat)))
-      )) +
+    graph +
+      scale_fill_ppd(guide = "none") +
+      scale_color_ppd(labels = Typred_label(), guide = guide_legend(
+        title = stat_legend_title(stat, deparse(substitute(stat))),
+        override.aes = list(fill = get_color(c(if (show_marginal) "d", "m"))))) +
       bayesplot_theme_get() +
       dont_expand_y_axis() +
       xaxis_title(FALSE) +
@@ -93,6 +112,7 @@ ppd_stat_grouped <-
   function(ypred,
            group,
            stat = "mean",
+           show_marginal = FALSE,
            ...,
            discrete = FALSE,
            facet_args = list(),
@@ -114,6 +134,7 @@ ppd_stat_grouped <-
 ppd_stat_freqpoly <-
   function(ypred,
            stat = "mean",
+           show_marginal = FALSE,
            ...,
            facet_args = list(),
            binwidth = NULL,
@@ -129,15 +150,18 @@ ppd_stat_freqpoly <-
     data <- ppd_stat_data(
       ypred = ypred,
       group = dots$group,
-      stat = match.fun(stat)
+      stat = match.fun(stat),
+      show_marginal = show_marginal
     )
-    ggplot(data, mapping = set_hist_aes(freq)) +
+    data$type <- ifelse(grepl("ypred", data$variable), "ypred", "PPD")
+
+    p <- ggplot(data, mapping = set_hist_aes(freq, color = .data$type)) +
       geom_freqpoly(
-        aes(color = "ypred"),
         linewidth = 0.5,
         na.rm = TRUE,
         binwidth = binwidth,
-        bins = bins
+        bins = bins,
+        data = data[data$type != "PPD",]
       ) +
       scale_color_ppd(
         name = stat_legend_title(stat, deparse(substitute(stat))),
@@ -149,6 +173,18 @@ ppd_stat_freqpoly <-
       yaxis_text(FALSE) +
       yaxis_ticks(FALSE) +
       yaxis_title(FALSE)
+
+    if (isTRUE(show_marginal)) {
+      p <- p +
+        geom_vline(
+          aes(xintercept = .data$value, color = .data$type),
+          data = data[data$type == "PPD",],
+          key_glyph = "path",
+          linewidth = 2
+        )
+    }
+
+    p
   }
 
 
@@ -158,6 +194,7 @@ ppd_stat_freqpoly_grouped <-
   function(ypred,
            group,
            stat = "mean",
+           show_marginal = FALSE,
            ...,
            facet_args = list(),
            binwidth = NULL,
@@ -177,6 +214,7 @@ ppd_stat_freqpoly_grouped <-
 ppd_stat_2d <-
   function(ypred,
            stat = c("mean", "sd"),
+           show_marginal = FALSE,
            ...,
            size = 2.5,
            alpha = 0.7) {
@@ -196,20 +234,25 @@ ppd_stat_2d <-
     data <- ppd_stat_data(
       ypred = ypred,
       group = NULL,
-      stat = c(match.fun(stat[[1]]), match.fun(stat[[2]]))
+      stat = c(match.fun(stat[[1]]), match.fun(stat[[2]])),
+      show_marginal = show_marginal
     )
+    data$type <- ifelse(grepl("ypred", data$variable), "ypred", "PPD")
+
     ggplot(data) +
       geom_point(
         mapping = aes(
           x = .data$value,
           y = .data$value2,
-          fill = "ypred",
-          color = "ypred"
+          fill = .data$type,
+          color = .data$type,
+          shape = .data$type
         ),
-        shape = 21,
         size = size,
         alpha = alpha
       ) +
+      scale_shape_manual(lgnd_title, labels = Typred_label(),
+                         values = c(ypred = 21, PPD = 23)) +
       scale_fill_ppd(lgnd_title, labels = Typred_label()) +
       scale_color_ppd(lgnd_title, labels = Typred_label()) +
       labs(x = stat_labs[1], y = stat_labs[2]) +
@@ -219,7 +262,7 @@ ppd_stat_2d <-
 
 #' @rdname PPD-test-statistics
 #' @export
-ppd_stat_data <- function(ypred, group = NULL, stat) {
+ppd_stat_data <- function(ypred, group = NULL, stat, show_marginal = FALSE) {
   if (!(length(stat) %in% 1:2)) {
     abort("'stat' must have length 1 or 2.")
   }
@@ -239,7 +282,8 @@ ppd_stat_data <- function(ypred, group = NULL, stat) {
     predictions = ypred,
     y = NULL,
     group = group,
-    stat = stat
+    stat = stat,
+    show_marginal = show_marginal
   )
 }
 
@@ -263,7 +307,7 @@ ppd_stat_data <- function(ypred, group = NULL, stat) {
 #' ppc_stat_data(y, yrep, group, stat = "median")
 #'
 #' @importFrom dplyr group_by ungroup summarise rename
-.ppd_stat_data <- function(predictions, y = NULL, group = NULL, stat) {
+.ppd_stat_data <- function(predictions, y = NULL, group = NULL, stat, show_marginal = FALSE) {
   stopifnot(length(stat) %in% c(1,2))
   if (length(stat) == 1) {
     stopifnot(is.function(stat)) # sanity check, should already be validated
@@ -292,10 +336,10 @@ ppd_stat_data <- function(ypred, group = NULL, stat) {
   )
   colnames(d) <- gsub(".", "_", colnames(d), fixed = TRUE)
   molten_d <- reshape2::melt(d, id.vars = "group")
-  molten_d <- group_by(molten_d, .data$group, .data$variable)
 
   data <-
     molten_d %>%
+    group_by(.data$group, .data$variable) %>%
     summarise(
       value1 = stat1(.data$value),
       value2 = if (!is.null(stat2))
@@ -304,6 +348,22 @@ ppd_stat_data <- function(ypred, group = NULL, stat) {
     rename(value = "value1") %>%
     ungroup()
 
+  if (isTRUE(show_marginal)) {
+    data_ppd <- molten_d %>%
+      dplyr::filter(.data$variable != "y") %>%
+      group_by(.data$group) %>%
+      summarise(
+        variable = "PPD",
+        value1 = stat1(.data$value),
+        value2 = if (!is.null(stat2))
+          stat2(.data$value) else NA
+      ) %>%
+      rename(value = "value1") %>%
+      ungroup()
+
+    data <- rbind(data, data_ppd)
+  }
+
   if (is.null(stat2)) {
     data$value2 <- NULL
   }
@@ -328,4 +388,7 @@ stat_group_facets <- function(facet_args, scales_default = "free") {
   do.call("facet_wrap", facet_args)
 }
 
-Typred_label <- function() expression(italic(T)(italic(y)[pred]))
+Typred_label <- function() {
+  expression(PPD = italic(T)(PPD),
+             ypred = italic(T)(italic(y)[pred]))
+}
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index a53c62e4..9c82f45f 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -17,6 +17,7 @@ ppd_data(ypred, group = NULL)
 
 ppd_dens_overlay(
   ypred,
+  show_marginal = FALSE,
   ...,
   size = 0.25,
   alpha = 0.7,
@@ -30,6 +31,7 @@ ppd_dens_overlay(
 
 ppd_ecdf_overlay(
   ypred,
+  show_marginal = FALSE,
   ...,
   discrete = FALSE,
   pad = TRUE,
@@ -37,18 +39,42 @@ ppd_ecdf_overlay(
   alpha = 0.7
 )
 
-ppd_dens(ypred, ..., trim = FALSE, size = 0.5, alpha = 1, bounds = NULL)
+ppd_dens(
+  ypred,
+  show_marginal = FALSE,
+  ...,
+  trim = FALSE,
+  size = 0.5,
+  alpha = 1,
+  bounds = NULL
+)
 
-ppd_hist(ypred, ..., binwidth = NULL, bins = NULL, breaks = NULL, freq = TRUE)
+ppd_hist(
+  ypred,
+  show_marginal = FALSE,
+  ...,
+  binwidth = NULL,
+  bins = NULL,
+  breaks = NULL,
+  freq = !show_marginal
+)
 
-ppd_dots(ypred, ..., binwidth = NA, quantiles = 100, freq = TRUE)
+ppd_dots(
+  ypred,
+  show_marginal = FALSE,
+  ...,
+  binwidth = NA,
+  quantiles = 100,
+  freq = TRUE
+)
 
 ppd_freqpoly(
   ypred,
+  show_marginal = FALSE,
   ...,
   binwidth = NULL,
   bins = NULL,
-  freq = TRUE,
+  freq = !show_marginal,
   size = 0.5,
   alpha = 1
 )
@@ -56,15 +82,23 @@ ppd_freqpoly(
 ppd_freqpoly_grouped(
   ypred,
   group,
+  show_marginal = FALSE,
   ...,
   binwidth = NULL,
   bins = NULL,
-  freq = TRUE,
+  freq = !show_marginal,
   size = 0.5,
   alpha = 1
 )
 
-ppd_boxplot(ypred, ..., notch = TRUE, size = 0.5, alpha = 1)
+ppd_boxplot(
+  ypred,
+  show_marginal = FALSE,
+  ...,
+  notch = TRUE,
+  size = 0.5,
+  alpha = 1
+)
 }
 \arguments{
 \item{ypred}{An \code{S} by \code{N} matrix of draws from the posterior (or prior)
@@ -77,6 +111,8 @@ Will be coerced to \link[base:factor]{factor} if not already a factor.
 Each value in \code{group} is interpreted as the group level pertaining
 to the corresponding observation.}
 
+\item{show_marginal}{Plot the marginal PPD along with the \code{ypred}s.}
+
 \item{...}{For dot plots, optional additional arguments to pass to \code{\link[ggdist:stat_dots]{ggdist::stat_dots()}}.}
 
 \item{size, alpha}{Passed to the appropriate geom to control the appearance of
@@ -143,6 +179,7 @@ the input contains the "success" \emph{proportions} (not discrete
 color_scheme_set("brightblue")
 preds <- example_yrep_draws()
 ppd_dens_overlay(ypred = preds[1:50, ])
+ppd_dens_overlay(ypred = preds[1:50, ], show_marginal = TRUE)
 ppc_dens_overlay(y = example_y_data(), yrep = preds[1:50, ])
 
 }
diff --git a/man/PPD-test-statistics.Rd b/man/PPD-test-statistics.Rd
index 7c8b41f3..21e789e5 100644
--- a/man/PPD-test-statistics.Rd
+++ b/man/PPD-test-statistics.Rd
@@ -14,18 +14,20 @@
 ppd_stat(
   ypred,
   stat = "mean",
+  show_marginal = FALSE,
   ...,
   discrete = FALSE,
   binwidth = NULL,
   bins = NULL,
   breaks = NULL,
-  freq = TRUE
+  freq = !show_marginal
 )
 
 ppd_stat_grouped(
   ypred,
   group,
   stat = "mean",
+  show_marginal = FALSE,
   ...,
   discrete = FALSE,
   facet_args = list(),
@@ -38,6 +40,7 @@ ppd_stat_grouped(
 ppd_stat_freqpoly(
   ypred,
   stat = "mean",
+  show_marginal = FALSE,
   ...,
   facet_args = list(),
   binwidth = NULL,
@@ -49,6 +52,7 @@ ppd_stat_freqpoly_grouped(
   ypred,
   group,
   stat = "mean",
+  show_marginal = FALSE,
   ...,
   facet_args = list(),
   binwidth = NULL,
@@ -56,9 +60,16 @@ ppd_stat_freqpoly_grouped(
   freq = TRUE
 )
 
-ppd_stat_2d(ypred, stat = c("mean", "sd"), ..., size = 2.5, alpha = 0.7)
+ppd_stat_2d(
+  ypred,
+  stat = c("mean", "sd"),
+  show_marginal = FALSE,
+  ...,
+  size = 2.5,
+  alpha = 0.7
+)
 
-ppd_stat_data(ypred, group = NULL, stat)
+ppd_stat_data(ypred, group = NULL, stat, show_marginal = FALSE)
 }
 \arguments{
 \item{ypred}{An \code{S} by \code{N} matrix of draws from the posterior (or prior)
@@ -73,6 +84,8 @@ statistic. If specified as a string (or strings) then the legend will
 display the function name(s). If specified as a function (or functions)
 then generic naming is used in the legend.}
 
+\item{show_marginal}{Plot the marginal PPD along with the \code{ypred}s.}
+
 \item{...}{Currently unused.}
 
 \item{discrete}{For \code{ppc_stat()} and \code{ppc_stat_grouped()}, if \code{TRUE} then a
@@ -130,11 +143,13 @@ the input contains the "success" \emph{proportions} (not discrete
 yrep <- example_yrep_draws()
 ppd_stat(yrep)
 ppd_stat(yrep, stat = "sd") + legend_none()
+ppd_stat(yrep, show_marginal = TRUE)
 
 # use your own function for the 'stat' argument
 color_scheme_set("brightblue")
 q25 <- function(y) quantile(y, 0.25)
 ppd_stat(yrep, stat = "q25") # legend includes function name
+ppd_stat(yrep, stat = "q25", show_marginal = TRUE)
 }
 \references{
 Gabry, J. , Simpson, D. , Vehtari, A. , Betancourt, M. and