From 5c37cdfa829ac5532eccd58d98be32a578cfe48f Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 12:54:21 +0200
Subject: [PATCH 01/25] ppd_dens_overlay

---
 R/helpers-ppc.R          |  8 ++++-
 R/ppd-distributions.R    | 74 +++++++++++++++++++++++++++-------------
 man/PPD-distributions.Rd |  3 ++
 3 files changed, 60 insertions(+), 25 deletions(-)

diff --git a/R/helpers-ppc.R b/R/helpers-ppc.R
index e268b315..3b9dfd07 100644
--- a/R/helpers-ppc.R
+++ b/R/helpers-ppc.R
@@ -594,4 +594,10 @@ 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(show_marginal = FALSE) {
+  if (isTRUE(show_marginal)) {
+    expression(PPD, italic(y)[pred])
+  } else {
+    expression(italic(y)[pred])
+  }
+}
diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 70c4e5a6..803d324c 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 yreps.
 #'
 #' @template details-binomial
 #' @template return-ggplot-or-data
@@ -37,25 +38,46 @@ 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 = "marginal"),
+          linewidth = 1,
+          alpha = alpha,
           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(show_marginal = TRUE),
+          values = get_color(c("d", "m")),
+          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()
+    }
+
   }
 
 
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index 45bbbe74..e6c3ec49 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,
@@ -77,6 +78,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 yreps.}
+
 \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

From 34bac15b7f44e55c986f50170f7d684eeab695ea Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 12:56:50 +0200
Subject: [PATCH 02/25] Update ppd-distributions.R

---
 R/ppd-distributions.R | 1 -
 1 file changed, 1 deletion(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 803d324c..989c2bb9 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -77,7 +77,6 @@ ppd_dens_overlay <-
         overlay_ppd_densities(
           mapping = aes(color = "marginal"),
           linewidth = 1,
-          alpha = alpha,
           trim = trim,
           bw = bw,
           adjust = adjust,

From 7b9b4134237a23e2e773ef9ee7622a716fa5216b Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 12:58:30 +0200
Subject: [PATCH 03/25] ppd_ecdf_overlay

---
 R/ppd-distributions.R | 35 +++++++++++++++++++++++++++--------
 1 file changed, 27 insertions(+), 8 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 989c2bb9..7ae7bf13 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -107,6 +107,7 @@ ppd_dens_overlay <-
 #' @export
 ppd_ecdf_overlay <-
   function(ypred,
+           show_marginal= FALSE,
            ...,
            discrete = FALSE,
            pad = TRUE,
@@ -115,7 +116,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),
@@ -129,16 +130,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(show_marginal = TRUE),
+          values = get_color(c("d", "m")),
+          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()
+    }
   }
 
 

From 3a2c6e5369b5d46acaf11cdc222da873a0a828c5 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 13:07:20 +0200
Subject: [PATCH 04/25] Update PPD-distributions.Rd

---
 man/PPD-distributions.Rd | 1 +
 1 file changed, 1 insertion(+)

diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index e6c3ec49..3a27d5d0 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -31,6 +31,7 @@ ppd_dens_overlay(
 
 ppd_ecdf_overlay(
   ypred,
+  show_marginal = FALSE,
   ...,
   discrete = FALSE,
   pad = TRUE,

From 68612d7d53217b326b56429ec0ac042b9ed44905 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 13:10:33 +0200
Subject: [PATCH 05/25] ppd_dens

---
 R/ppd-distributions.R    | 63 ++++++++++++++++++++++++++++------------
 man/PPD-distributions.Rd | 10 ++++++-
 2 files changed, 54 insertions(+), 19 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 7ae7bf13..b968a53a 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -165,6 +165,7 @@ ppd_ecdf_overlay <-
 #' @export
 ppd_dens <-
     function(ypred,
+             show_marginal = FALSE,
              ...,
              trim = FALSE,
              size = 0.5,
@@ -174,30 +175,56 @@ 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)
+
+      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
+          ) +
+          scale_color_ppd(
+            labels = ypred_label(show_marginal = TRUE),
+            values = get_color(c("d", "m")),
+            guide = guide_legend(
+              override.aes = list(size = 2 * size, alpha = 1))
+          ) +
+          scale_fill_ppd(
+            labels = ypred_label(show_marginal = TRUE),
+            values = get_color(c("d", "m")),
+            guide = guide_legend(
+              override.aes = list(size = 2 * size, alpha = 1))
+          )
+      } else {
+        p +
+          scale_color_ppd() +
+          scale_fill_ppd() +
+          legend_none()
+      }
+    }
 
 
 #' @rdname PPD-distributions
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index 3a27d5d0..634a9376 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -39,7 +39,15 @@ 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)
 

From c07c0d07db41820789848e4dc86c26ce7988b1b0 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 13:17:28 +0200
Subject: [PATCH 06/25] ppd_hist

---
 R/ppd-distributions.R    | 43 ++++++++++++++++++++++++++++++++--------
 man/PPD-distributions.Rd | 10 +++++++++-
 2 files changed, 44 insertions(+), 9 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index b968a53a..a985e2d3 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -231,6 +231,7 @@ ppd_dens <-
 #' @export
 ppd_hist <-
   function(ypred,
+           show_marginal = FALSE,
            ...,
            binwidth = NULL,
            bins = NULL,
@@ -239,29 +240,55 @@ ppd_hist <-
     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",
+            y = after_stat(density)),
         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)
+
+    if (isTRUE(show_marginal)) {
+      data2 <- transform(data, rep_label = "PPD")
+
+      p +
+        geom_histogram(
+          aes(color = "PPD",
+              fill = "PPD",
+              y = after_stat(density)),
+          linewidth = 1,
+          binwidth = binwidth,
+          bins = bins,
+          breaks = breaks,
+          data = data2
+        ) +
+        scale_color_ppd(
+          labels = ypred_label(show_marginal = TRUE),
+          values = get_color(c("d", "m"))
+        ) +
+        scale_fill_ppd(
+          labels = ypred_label(show_marginal = TRUE),
+          values = get_color(c("d", "m"))
+        )
+
+    } else {
+      p +
+        scale_color_ppd() +
+        scale_fill_ppd() +
+        legend_none()
+    }
   }
 
 #' @rdname PPD-distributions
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index 634a9376..f02ec607 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -49,7 +49,15 @@ ppd_dens(
   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 = TRUE
+)
 
 ppd_dots(ypred, ..., binwidth = NA, quantiles = 100, freq = TRUE)
 

From 00cabc633b9499930098d043ff692c0c8c2f2b74 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 13:22:13 +0200
Subject: [PATCH 07/25] ppd_dots

---
 R/ppd-distributions.R    | 41 ++++++++++++++++++++++++++++++++--------
 man/PPD-distributions.Rd |  9 ++++++++-
 2 files changed, 41 insertions(+), 9 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index a985e2d3..1a023e4a 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -295,6 +295,7 @@ ppd_hist <-
 #' @export
 ppd_dots <-
   function(ypred,
+           show_marginal = FALSE,
            ...,
            binwidth = NA,
            quantiles = 100,
@@ -304,28 +305,52 @@ 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)
+
+    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,
+          ...
+        ) +
+        scale_color_ppd(
+          labels = ypred_label(show_marginal = TRUE),
+          values = get_color(c("d", "m"))
+        ) +
+        scale_fill_ppd(
+          labels = ypred_label(show_marginal = TRUE),
+          values = get_color(c("d", "m"))
+        )
+
+    } else {
+      p +
+        scale_color_ppd() +
+        scale_fill_ppd() +
+        legend_none()
+    }
   }
 
 
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index f02ec607..98e22994 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -59,7 +59,14 @@ ppd_hist(
   freq = TRUE
 )
 
-ppd_dots(ypred, ..., binwidth = NA, quantiles = 100, freq = TRUE)
+ppd_dots(
+  ypred,
+  show_marginal = FALSE,
+  ...,
+  binwidth = NA,
+  quantiles = 100,
+  freq = TRUE
+)
 
 ppd_freqpoly(
   ypred,

From a7a34aaf01dd2a07577043bfac96487862335dc9 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 13:36:42 +0200
Subject: [PATCH 08/25] fix hist

---
 R/ppd-distributions.R    | 8 +++-----
 man/PPD-distributions.Rd | 2 +-
 2 files changed, 4 insertions(+), 6 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 1a023e4a..855a82f6 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -236,15 +236,14 @@ ppd_hist <-
            binwidth = NULL,
            bins = NULL,
            breaks = NULL,
-           freq = TRUE) {
+           freq = !show_marginal) {
     check_ignored_arguments(...)
 
     data <- ppd_data(ypred)
     p <- ggplot(data, mapping = set_hist_aes(freq)) +
       geom_histogram(
         aes(color = "ypred",
-            fill = "ypred",
-            y = after_stat(density)),
+            fill = "ypred"),
         linewidth = 0.25,
         binwidth = binwidth,
         bins = bins,
@@ -266,8 +265,7 @@ ppd_hist <-
       p +
         geom_histogram(
           aes(color = "PPD",
-              fill = "PPD",
-              y = after_stat(density)),
+              fill = "PPD"),
           linewidth = 1,
           binwidth = binwidth,
           bins = bins,
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index 98e22994..4fa8cb27 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -56,7 +56,7 @@ ppd_hist(
   binwidth = NULL,
   bins = NULL,
   breaks = NULL,
-  freq = TRUE
+  freq = !show_marginal
 )
 
 ppd_dots(

From 24b5000b1afbb1512f3369aedffaefcc1cf5db6a Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 13:37:57 +0200
Subject: [PATCH 09/25] ppd_freqpoly

---
 R/ppd-distributions.R    | 48 +++++++++++++++++++++++++++++++---------
 man/PPD-distributions.Rd |  3 ++-
 2 files changed, 40 insertions(+), 11 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 855a82f6..23660a52 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -356,10 +356,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) {
 
@@ -370,12 +371,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,
@@ -383,8 +382,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() +
@@ -392,8 +389,39 @@ ppd_freqpoly <-
       yaxis_title(FALSE) +
       yaxis_ticks(FALSE) +
       xaxis_title(FALSE) +
-      facet_text(FALSE) +
-      legend_none()
+      facet_text(FALSE)
+
+
+    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,
+        ) +
+        scale_color_ppd(
+          labels = ypred_label(show_marginal = TRUE),
+          values = get_color(c("d", "m")),
+          guide = guide_legend(override.aes = list(size = 2 * size, alpha = 1))
+        ) +
+        scale_fill_ppd(
+          labels = ypred_label(show_marginal = TRUE),
+          values = get_color(c("d", "m")),
+          guide = guide_legend(override.aes = list(size = 2 * size, alpha = 1))
+        )
+
+    } else {
+      p +
+        scale_color_ppd() +
+        scale_fill_ppd() +
+        legend_none()
+    }
   }
 
 
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index 4fa8cb27..62da866c 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -70,10 +70,11 @@ ppd_dots(
 
 ppd_freqpoly(
   ypred,
+  show_marginal = FALSE,
   ...,
   binwidth = NULL,
   bins = NULL,
-  freq = TRUE,
+  freq = !show_marginal,
   size = 0.5,
   alpha = 1
 )

From 58b3463e9bec471a86029dad9bd81ab4b6d7d241 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 13:40:28 +0200
Subject: [PATCH 10/25] fix colors

---
 R/ppd-distributions.R | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 23660a52..99fbb45b 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -208,7 +208,7 @@ ppd_dens <-
           ) +
           scale_color_ppd(
             labels = ypred_label(show_marginal = TRUE),
-            values = get_color(c("d", "m")),
+            values = get_color(c("dh", "mh")),
             guide = guide_legend(
               override.aes = list(size = 2 * size, alpha = 1))
           ) +
@@ -274,7 +274,7 @@ ppd_hist <-
         ) +
         scale_color_ppd(
           labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m"))
+          values = get_color(c("dh", "mh"))
         ) +
         scale_fill_ppd(
           labels = ypred_label(show_marginal = TRUE),
@@ -336,7 +336,7 @@ ppd_dots <-
         ) +
         scale_color_ppd(
           labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m"))
+          values = get_color(c("dh", "mh"))
         ) +
         scale_fill_ppd(
           labels = ypred_label(show_marginal = TRUE),
@@ -407,7 +407,7 @@ ppd_freqpoly <-
         ) +
         scale_color_ppd(
           labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m")),
+          values = get_color(c("dh", "mh")),
           guide = guide_legend(override.aes = list(size = 2 * size, alpha = 1))
         ) +
         scale_fill_ppd(

From e6fcb6ae2c9fa022a16b832a3b938a94993ce970 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 15:17:27 +0200
Subject: [PATCH 11/25] simplify code

---
 R/helpers-gg.R        |  4 +--
 R/helpers-ppc.R       | 11 +++----
 R/ppd-distributions.R | 72 ++++++++++---------------------------------
 3 files changed, 24 insertions(+), 63 deletions(-)

diff --git a/R/helpers-gg.R b/R/helpers-gg.R
index 953a51d4..1ac16af9 100644
--- a/R/helpers-gg.R
+++ b/R/helpers-gg.R
@@ -125,7 +125,7 @@ scale_fill_ppc <-
 
 scale_color_ppd <-
   function(name = NULL,
-           values = get_color("mh"),
+           values = setNames(get_color(c("dh", "mh")), nm = c("PPD", "ypred")),
            labels = ypred_label(),
            ...) {
     scale_color_ppc(name = name,
@@ -136,7 +136,7 @@ scale_color_ppd <-
 
 scale_fill_ppd <-
   function(name = NULL,
-           values = get_color("m"),
+           values = setNames(get_color(c("d", "m")), nm = c("PPD", "ypred")),
            labels = ypred_label(),
            ...) {
     scale_fill_ppc(name = name,
diff --git a/R/helpers-ppc.R b/R/helpers-ppc.R
index 3b9dfd07..267f1126 100644
--- a/R/helpers-ppc.R
+++ b/R/helpers-ppc.R
@@ -594,10 +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(show_marginal = FALSE) {
-  if (isTRUE(show_marginal)) {
-    expression(PPD, italic(y)[pred])
-  } else {
-    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 99fbb45b..94f94ff6 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -75,7 +75,7 @@ ppd_dens_overlay <-
     if (isTRUE(show_marginal)) {
       p +
         overlay_ppd_densities(
-          mapping = aes(color = "marginal"),
+          mapping = aes(color = "PPD"),
           linewidth = 1,
           trim = trim,
           bw = bw,
@@ -85,8 +85,8 @@ ppd_dens_overlay <-
           n = n_dens
         ) +
         scale_color_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m")),
+          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))
         )
@@ -144,8 +144,8 @@ ppd_ecdf_overlay <-
           pad = pad
         ) +
         scale_color_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m")),
+          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))
         )
@@ -192,7 +192,9 @@ ppd_dens <-
         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")
@@ -205,23 +207,9 @@ ppd_dens <-
             trim = trim,
             bounds = bounds,
             data = data2
-          ) +
-          scale_color_ppd(
-            labels = ypred_label(show_marginal = TRUE),
-            values = get_color(c("dh", "mh")),
-            guide = guide_legend(
-              override.aes = list(size = 2 * size, alpha = 1))
-          ) +
-          scale_fill_ppd(
-            labels = ypred_label(show_marginal = TRUE),
-            values = get_color(c("d", "m")),
-            guide = guide_legend(
-              override.aes = list(size = 2 * size, alpha = 1))
           )
       } else {
         p +
-          scale_color_ppd() +
-          scale_fill_ppd() +
           legend_none()
       }
     }
@@ -257,7 +245,9 @@ ppd_hist <-
       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")
@@ -271,20 +261,10 @@ ppd_hist <-
           bins = bins,
           breaks = breaks,
           data = data2
-        ) +
-        scale_color_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("dh", "mh"))
-        ) +
-        scale_fill_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m"))
         )
 
     } else {
       p +
-        scale_color_ppd() +
-        scale_fill_ppd() +
         legend_none()
     }
   }
@@ -320,7 +300,9 @@ ppd_dots <-
       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")
@@ -333,20 +315,10 @@ ppd_dots <-
           binwidth = binwidth,
           quantiles = quantiles,
           ...
-        ) +
-        scale_color_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("dh", "mh"))
-        ) +
-        scale_fill_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m"))
         )
 
     } else {
       p +
-        scale_color_ppd() +
-        scale_fill_ppd() +
         legend_none()
     }
   }
@@ -389,7 +361,9 @@ ppd_freqpoly <-
       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)) {
@@ -404,22 +378,10 @@ ppd_freqpoly <-
           binwidth = binwidth,
           bins = bins,
           linewidth = 1,
-        ) +
-        scale_color_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("dh", "mh")),
-          guide = guide_legend(override.aes = list(size = 2 * size, alpha = 1))
-        ) +
-        scale_fill_ppd(
-          labels = ypred_label(show_marginal = TRUE),
-          values = get_color(c("d", "m")),
-          guide = guide_legend(override.aes = list(size = 2 * size, alpha = 1))
         )
 
     } else {
       p +
-        scale_color_ppd() +
-        scale_fill_ppd() +
         legend_none()
     }
   }

From a17a0e5b5712f141175dda67b3433a5c17ac7825 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 15:19:28 +0200
Subject: [PATCH 12/25] ppd_freqpoly

---
 R/ppd-distributions.R    | 3 ++-
 man/PPD-distributions.Rd | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 94f94ff6..47ca1d38 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -392,10 +392,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) {
 
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index 62da866c..6f0687de 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -82,10 +82,11 @@ ppd_freqpoly(
 ppd_freqpoly_grouped(
   ypred,
   group,
+  show_marginal = FALSE,
   ...,
   binwidth = NULL,
   bins = NULL,
-  freq = TRUE,
+  freq = !show_marginal,
   size = 0.5,
   alpha = 1
 )

From 5d9d505e16f17398791b4a655a08b768d5a891d1 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 15:23:56 +0200
Subject: [PATCH 13/25] ppd_boxplot

---
 R/ppd-distributions.R    | 29 +++++++++++++++++++++++------
 man/PPD-distributions.Rd |  9 ++++++++-
 2 files changed, 31 insertions(+), 7 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 47ca1d38..047c8d42 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -419,6 +419,7 @@ ppd_freqpoly_grouped <-
 #' @export
 ppd_boxplot <-
   function(ypred,
+           show_marginal = FALSE,
            ...,
            notch = TRUE,
            size = 0.5,
@@ -426,13 +427,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,
@@ -447,8 +448,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("mh"),
+          outlier.alpha = 2/3,
+          outlier.size = 1
+        )
+
+    } else {
+      p + legend_none()
+    }
   }
 
 
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index 6f0687de..b9224426 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -91,7 +91,14 @@ ppd_freqpoly_grouped(
   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)

From 314a864f6e0cf62869fa576bdaa62f298cfdb5a5 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 20:17:49 +0200
Subject: [PATCH 14/25] ppd_stat

---
 R/ppd-test-statistics.R    | 94 ++++++++++++++++++++++++++------------
 man/PPD-test-statistics.Rd |  3 +-
 2 files changed, 67 insertions(+), 30 deletions(-)

diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index d7146405..35b51b09 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -34,12 +34,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 +51,48 @@ 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 = subset(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 = subset(data, type == "PPD"),
+            key_glyph = "rect"
+          )
+      }
     }
-    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 +108,7 @@ ppd_stat_grouped <-
   function(ypred,
            group,
            stat = "mean",
+           show_marginal = FALSE,
            ...,
            discrete = FALSE,
            facet_args = list(),
@@ -219,7 +235,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 +255,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 +280,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) {
   stopifnot(length(stat) %in% c(1,2))
   if (length(stat) == 1) {
     stopifnot(is.function(stat)) # sanity check, should already be validated
@@ -292,10 +309,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 +321,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 +361,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-test-statistics.Rd b/man/PPD-test-statistics.Rd
index 168f2bb3..9ddb0335 100644
--- a/man/PPD-test-statistics.Rd
+++ b/man/PPD-test-statistics.Rd
@@ -14,6 +14,7 @@
 ppd_stat(
   ypred,
   stat = "mean",
+  show_marginal = FALSE,
   ...,
   discrete = FALSE,
   binwidth = NULL,
@@ -58,7 +59,7 @@ ppd_stat_freqpoly_grouped(
 
 ppd_stat_2d(ypred, stat = c("mean", "sd"), ..., size = 2.5, alpha = 0.7)
 
-ppd_stat_data(ypred, group = NULL, stat)
+ppd_stat_data(ypred, group = NULL, stat, show_marginal)
 }
 \arguments{
 \item{ypred}{An \code{S} by \code{N} matrix of draws from the posterior (or prior)

From ff2de90868e7c0ef9aa7d4fcc80d9f9781f91581 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 20:19:00 +0200
Subject: [PATCH 15/25] Update PPD-test-statistics.Rd

---
 man/PPD-test-statistics.Rd | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/man/PPD-test-statistics.Rd b/man/PPD-test-statistics.Rd
index 9ddb0335..cfaee44f 100644
--- a/man/PPD-test-statistics.Rd
+++ b/man/PPD-test-statistics.Rd
@@ -20,13 +20,14 @@ ppd_stat(
   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(),
@@ -59,7 +60,7 @@ ppd_stat_freqpoly_grouped(
 
 ppd_stat_2d(ypred, stat = c("mean", "sd"), ..., size = 2.5, alpha = 0.7)
 
-ppd_stat_data(ypred, group = NULL, stat, show_marginal)
+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)

From 239f5e0d22ecbb186871f7a2b935f72833b18294 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 20:26:29 +0200
Subject: [PATCH 16/25] ppd_stat_freqpoly

---
 R/ppd-test-statistics.R    | 24 ++++++++++++++++++++----
 man/PPD-test-statistics.Rd |  2 ++
 2 files changed, 22 insertions(+), 4 deletions(-)

diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index 35b51b09..70e9408b 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -130,6 +130,7 @@ ppd_stat_grouped <-
 ppd_stat_freqpoly <-
   function(ypred,
            stat = "mean",
+           show_marginal = FALSE,
            ...,
            facet_args = list(),
            binwidth = NULL,
@@ -145,15 +146,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 = subset(data, type != "PPD")
       ) +
       scale_color_ppd(
         name = stat_legend_title(stat, deparse(substitute(stat))),
@@ -165,6 +169,17 @@ 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 = subset(data, type == "PPD"),
+          key_glyph = "path"
+        )
+    }
+
+    p
   }
 
 
@@ -174,6 +189,7 @@ ppd_stat_freqpoly_grouped <-
   function(ypred,
            group,
            stat = "mean",
+           show_marginal = FALSE,
            ...,
            facet_args = list(),
            binwidth = NULL,
diff --git a/man/PPD-test-statistics.Rd b/man/PPD-test-statistics.Rd
index cfaee44f..c5d2fb2f 100644
--- a/man/PPD-test-statistics.Rd
+++ b/man/PPD-test-statistics.Rd
@@ -40,6 +40,7 @@ ppd_stat_grouped(
 ppd_stat_freqpoly(
   ypred,
   stat = "mean",
+  show_marginal = FALSE,
   ...,
   facet_args = list(),
   binwidth = NULL,
@@ -51,6 +52,7 @@ ppd_stat_freqpoly_grouped(
   ypred,
   group,
   stat = "mean",
+  show_marginal = FALSE,
   ...,
   facet_args = list(),
   binwidth = NULL,

From 6e901b85d2c4418b41ca963932ea4bb10dab18aa Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 20:31:48 +0200
Subject: [PATCH 17/25] ppd_stat_2d

---
 R/ppd-test-statistics.R    | 14 ++++++++++----
 man/PPD-test-statistics.Rd |  9 ++++++++-
 2 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index 70e9408b..330d9a80 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -209,6 +209,7 @@ ppd_stat_freqpoly_grouped <-
 ppd_stat_2d <-
   function(ypred,
            stat = c("mean", "sd"),
+           show_marginal = FALSE,
            ...,
            size = 2.5,
            alpha = 0.7) {
@@ -228,20 +229,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]) +
diff --git a/man/PPD-test-statistics.Rd b/man/PPD-test-statistics.Rd
index c5d2fb2f..0492e342 100644
--- a/man/PPD-test-statistics.Rd
+++ b/man/PPD-test-statistics.Rd
@@ -60,7 +60,14 @@ 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, show_marginal = FALSE)
 }

From 926537332ddd8879afa07277de727148ab6e7fb3 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 20:31:50 +0200
Subject: [PATCH 18/25] Update NEWS.md

---
 NEWS.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/NEWS.md b/NEWS.md
index 25c323e3..aaec1cd9 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)
 * New functions `mcmc_dots` and `mcmc_dots_by_chain` for dot plots of MCMC draws by @behramulukir (#402)
 * Default to `quantiles=100` for all dot plots by @behramulukir (#402)
 

From 70f6d7215ffd9da1859a29350fe9fe746ac10694 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 20:41:19 +0200
Subject: [PATCH 19/25] make line thicker

---
 R/ppd-test-statistics.R | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index 330d9a80..a6e32f50 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -84,7 +84,8 @@ ppd_stat <-
           geom_vline(
             aes(xintercept = .data$value, color = .data$type),
             data = subset(data, type == "PPD"),
-            key_glyph = "rect"
+            key_glyph = "rect",
+            linewidth = 2
           )
       }
     }
@@ -175,7 +176,8 @@ ppd_stat_freqpoly <-
         geom_vline(
           aes(xintercept = .data$value, color = .data$type),
           data = subset(data, type == "PPD"),
-          key_glyph = "path"
+          key_glyph = "path",
+          linewidth = 2
         )
     }
 

From 360447f25ddfcf8fb709ece74d8103bc489aa7a1 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 20:41:33 +0200
Subject: [PATCH 20/25] lighter outliers

---
 R/ppd-distributions.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index 047c8d42..ff1726d2 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -458,7 +458,7 @@ ppd_boxplot <-
               fill = "PPD"),
           notch = notch,
           linewidth = 1,
-          outlier.color = get_color("mh"),
+          outlier.color = get_color("lh"),
           outlier.alpha = 2/3,
           outlier.size = 1
         )

From 0cfb6286ddef66b95af9b92806de397a20308540 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Sun, 8 Mar 2026 23:11:03 +0200
Subject: [PATCH 21/25] docs and examples

---
 R/ppd-distributions.R      | 3 ++-
 R/ppd-test-statistics.R    | 2 ++
 man/PPD-distributions.Rd   | 3 ++-
 man/PPD-test-statistics.Rd | 2 ++
 4 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/R/ppd-distributions.R b/R/ppd-distributions.R
index ff1726d2..ebd5eb43 100644
--- a/R/ppd-distributions.R
+++ b/R/ppd-distributions.R
@@ -10,7 +10,7 @@
 #'
 #' @template args-ypred
 #' @inheritParams PPC-distributions
-#' @param show_marginal Plot the marginal PPD along with the yreps.
+#' @param show_marginal Plot the marginal PPD along with the `ypred`s.
 #'
 #' @template details-binomial
 #' @template return-ggplot-or-data
@@ -20,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
diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index a6e32f50..5c0e3bf9 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -22,11 +22,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
diff --git a/man/PPD-distributions.Rd b/man/PPD-distributions.Rd
index b9224426..e8a2411b 100644
--- a/man/PPD-distributions.Rd
+++ b/man/PPD-distributions.Rd
@@ -111,7 +111,7 @@ 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 yreps.}
+\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()}}.}
 
@@ -179,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 0492e342..82554f95 100644
--- a/man/PPD-test-statistics.Rd
+++ b/man/PPD-test-statistics.Rd
@@ -141,11 +141,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

From 6ad34ef0a3085d5609dd93b19ba25a8cc7cc6501 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Tue, 10 Mar 2026 13:19:36 +0200
Subject: [PATCH 22/25] fix vignette error

---
 R/ppd-test-statistics.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index 5c0e3bf9..da282366 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -306,7 +306,7 @@ ppd_stat_data <- function(ypred, group = NULL, stat, show_marginal = FALSE) {
 #' 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, show_marginal) {
+.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

From 0fda5cd471ca98c767fac30fd9c4f55d9c302252 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Tue, 10 Mar 2026 22:01:30 +0200
Subject: [PATCH 23/25] add docs to PPD-test-statistics

---
 R/ppd-test-statistics.R    | 1 +
 man/PPD-test-statistics.Rd | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index da282366..d5476ff0 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
diff --git a/man/PPD-test-statistics.Rd b/man/PPD-test-statistics.Rd
index cbd6a9b8..21e789e5 100644
--- a/man/PPD-test-statistics.Rd
+++ b/man/PPD-test-statistics.Rd
@@ -84,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

From 79479a3af97937c18b6aaec157c858bc2e7a2fd8 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Tue, 10 Mar 2026 22:04:03 +0200
Subject: [PATCH 24/25] fix binding warning

---
 R/ppd-test-statistics.R | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/ppd-test-statistics.R b/R/ppd-test-statistics.R
index d5476ff0..412d6f1e 100644
--- a/R/ppd-test-statistics.R
+++ b/R/ppd-test-statistics.R
@@ -75,7 +75,7 @@ ppd_stat <-
     } else {
       graph <- graph +
         geom_histogram(
-          data = subset(data, type != "PPD"),
+          data = data[data$type != "PPD",],
           linewidth = 0.25,
           na.rm = TRUE,
           binwidth = binwidth,
@@ -86,7 +86,7 @@ ppd_stat <-
         graph <- graph +
           geom_vline(
             aes(xintercept = .data$value, color = .data$type),
-            data = subset(data, type == "PPD"),
+            data = data[data$type == "PPD",],
             key_glyph = "rect",
             linewidth = 2
           )
@@ -161,7 +161,7 @@ ppd_stat_freqpoly <-
         na.rm = TRUE,
         binwidth = binwidth,
         bins = bins,
-        data = subset(data, type != "PPD")
+        data = data[data$type != "PPD",]
       ) +
       scale_color_ppd(
         name = stat_legend_title(stat, deparse(substitute(stat))),
@@ -178,7 +178,7 @@ ppd_stat_freqpoly <-
       p <- p +
         geom_vline(
           aes(xintercept = .data$value, color = .data$type),
-          data = subset(data, type == "PPD"),
+          data = data[data$type == "PPD",],
           key_glyph = "path",
           linewidth = 2
         )

From 4ec64c6a81086eb42a41b9649224aee8e3e407c3 Mon Sep 17 00:00:00 2001
From: "Mattan S. Ben-Shachar" <mattansb@msbstats.info>
Date: Tue, 10 Mar 2026 22:08:24 +0200
Subject: [PATCH 25/25] cleanup some code

---
 R/helpers-gg.R | 28 ++++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/R/helpers-gg.R b/R/helpers-gg.R
index 1ac16af9..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 = setNames(get_color(c("dh", "mh")), nm = c("PPD", "ypred")),
-           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 = setNames(get_color(c("d", "m")), nm = c("PPD", "ypred")),
-           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(),
+      ...
+    )
   }