diff --git a/DESCRIPTION b/DESCRIPTION
index 564af64cb..230fe9905 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: bayestestR
 Title: Understand and Describe Bayesian Models and Posterior Distributions
-Version: 0.15.0.3
+Version: 0.15.0.4
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/NEWS.md b/NEWS.md
index 1af4ab5a5..bd534d59e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,10 @@
+# bayestestR 0.15.0.4
+
+## Bug fixes
+
+* Fix to `emmeans` / `marginaleffects` / `data.frame(<rvar>)` methods when using multiple credible levels (#688).
+
+
 # bayestestR 0.15.0
 
 ## Changes
diff --git a/R/utils.R b/R/utils.R
index 9bf719549..206a9c951 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -230,7 +230,7 @@
   datagrid <- insight::get_datagrid(object)
   grid_names <- colnames(datagrid)
 
-  if (long) {
+  if (long || nrow(datagrid) < nrow(results)) {
     datagrid$Parameter <- unique(results$Parameter)
     results <- datawizard::data_merge(datagrid, results, by = "Parameter")
     results$Parameter <- NULL
@@ -270,7 +270,7 @@
   grid_names <- colnames(object)[!is_rvar]
   datagrid <- data.frame(object[, grid_names, drop = FALSE])
 
-  if (long) {
+  if (long || nrow(datagrid) < nrow(results)) {
     datagrid$Parameter <- unique(results$Parameter)
     results <- datawizard::data_merge(datagrid, results, by = "Parameter")
     results$Parameter <- NULL
diff --git a/tests/testthat/test-data.frame-with-rvar.R b/tests/testthat/test-data.frame-with-rvar.R
index 4399c75b0..5d59300f7 100644
--- a/tests/testthat/test-data.frame-with-rvar.R
+++ b/tests/testthat/test-data.frame-with-rvar.R
@@ -45,6 +45,9 @@ test_that("data.frame w/ rvar_col descrive_posterior etc", {
   res <- eti(dfx, rvar_col = "my_rvar", ci = c(0.8, 0.95))
   res.ref <- eti(dfx$my_rvar, ci = c(0.8, 0.95))
   expect_true(all(c("mu", "sigma") %in% colnames(res)))
+  expect_identical(as.data.frame(res[c("mu", "sigma")]),
+                   data.frame(mu = c(0, 0, 0.5, 0.5, 1, 1),
+                              sigma = c(1, 1, 0.5, 0.5, 0.25, 0.25)))
   expect_identical(nrow(format(res)), 3L)
   expect_identical(ncol(format(res)), 4L)
   expect_equal(res[setdiff(colnames(res), c("mu", "sigma"))],
diff --git a/tests/testthat/test-emmGrid.R b/tests/testthat/test-emmGrid.R
index ed05ad940..0b300e460 100644
--- a/tests/testthat/test-emmGrid.R
+++ b/tests/testthat/test-emmGrid.R
@@ -31,6 +31,11 @@ test_that("emmGrid hdi", {
 
   xhdi2 <- hdi(emc_, ci = 0.95)
   expect_equal(xhdi$CI_low, xhdi2$CI_low)
+
+  xhdi3 <- hdi(all_, ci = c(0.9, 0.95))
+  expect_equal(as.data.frame(xhdi3[1:2]),
+               data.frame(group = c("1", "1", "2", "2", ".", "."),
+                          contrast = c(".", ".", ".", ".", "group1 - group2", "group1 - group2")))
 })
 
 test_that("emmGrid point_estimate", {
diff --git a/tests/testthat/test-marginaleffects.R b/tests/testthat/test-marginaleffects.R
index 8fb654499..1f5fe7642 100644
--- a/tests/testthat/test-marginaleffects.R
+++ b/tests/testthat/test-marginaleffects.R
@@ -31,6 +31,19 @@ withr::with_environment(
       ignore_attr = TRUE
     )
 
+    # multi ci levels
+    res <- hdi(mfx, ci = c(0.8, 0.9))
+    expect_equal(as.data.frame(res[1:3]),
+                 data.frame(term = c("am", "am", "am", "am", "cyl", "cyl",
+                                     "cyl", "cyl", "cyl", "cyl", "cyl", "cyl",
+                                     "hp", "hp", "hp", "hp"),
+                            contrast = c("1 - 0", "1 - 0", "1 - 0", "1 - 0",
+                                         "6 - 4", "6 - 4", "8 - 4", "8 - 4",
+                                         "6 - 4", "6 - 4", "8 - 4", "8 - 4",
+                                         "dY/dX", "dY/dX", "dY/dX", "dY/dX"),
+                            am = c(0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0,
+                                   1, 1)))
+
     # estimate_density
     mfx <- marginaleffects::comparisons(mod,
       variables = "cyl",