diff --git a/DESCRIPTION b/DESCRIPTION
index b5f2b0b5b..b592a6a4f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: bayestestR
 Title: Understand and Describe Bayesian Models and Posterior Distributions
-Version: 0.14.0.5
+Version: 0.14.0.6
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
@@ -67,7 +67,7 @@ Depends:
     R (>= 3.6)
 Imports:
     insight (>= 0.20.4.2),
-    datawizard (>= 0.10.0),
+    datawizard (>= 0.12.3.1),
     graphics,
     methods,
     stats,
@@ -127,4 +127,4 @@ Config/testthat/parallel: true
 Config/rcmdcheck/ignore-inconsequential-notes: true
 Config/Needs/website: easystats/easystatstemplate
 Config/Needs/check: stan-dev/cmdstanr
-Remotes: easystats/insight
+Remotes: easystats/insight, easystats/datawizard
diff --git a/NEWS.md b/NEWS.md
index 3eebdf519..db28c58f5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,10 @@
 
 ## Changes
 
+* Support for `posterior::rvar`-type column in data frames. 
+  For example, a data frame `df` with an `rvar` column `".pred"` can now be 
+  called directly via `p_direction(df, rvar_col = ".pred")`.
+
 * Added support for `{marginaleffects}`
 
 * Results from objects generated by `{emmeans}` (`emmGrid`/`emm_list`) now
@@ -17,9 +21,6 @@
   - `p_direction()` gets a `remove_na` argument, which defaults to `TRUE`, to
     remove `NA` values from the input before calculating the pd-values.
 
-  - The `data.frame` method for `p_direction()` gets an `rvar_col` argument, to
-    specify the column that contains the `rvar` objects.
-
   - Besides the existing `as.numeric()` method, `p_direction()` now also has an
     `as.vector()` method.
 
diff --git a/R/bayesfactor_parameters.R b/R/bayesfactor_parameters.R
index f11f9c665..84127a40c 100644
--- a/R/bayesfactor_parameters.R
+++ b/R/bayesfactor_parameters.R
@@ -85,26 +85,22 @@
 #' (Note that by default, `brms::brm()` uses flat priors for fixed-effects;
 #' See example below.)
 #' \cr\cr
-#' It is important to provide the correct `prior` for meaningful results.
+#' It is important to provide the correct `prior` for meaningful results,
+#' to match the `posterior`-type input:
 #'
-#'   - When `posterior` is a numerical vector, `prior` should also be a numerical vector.
-#'   - When `posterior` is a `data.frame`, `prior` should also be a `data.frame`, with matching column order.
-#'   - When `posterior` is a `stanreg`, `brmsfit` or other supported Bayesian model:
-#'       - `prior` can be set to `NULL`, in which case prior samples are drawn internally.
-#'       - `prior` can also be a model equivalent to `posterior` but with samples from
-#'         the priors *only*. See [unupdate()].
-#'       - **Note:** When `posterior` is a `brmsfit_multiple` model, `prior` **must** be provided.
-#'   - When `posterior` is an output from a `{marginaleffects}` function, `prior` should also be an an output
-#'   from a `{marginaleffects}` function equivalent to `posterior` but created
-#'   with a model of priors samples *only*.
-#'   - When `posterior` is an `emmGrid` / `emm_list` object:
-#'       - `prior` should also be an `emmGrid` / `emm_list` object equivalent to `posterior` but
-#'        created with a model of priors samples *only*. See [unupdate()].
-#'       - `prior` can also be the original (posterior) *model*. If so, the function will try to
-#'        update the `emmGrid` / `emm_list` to use the [unupdate()]d prior-model.
-#'        (*This cannot be done for `brmsfit` models.*)
-#'       - **Note**: When the `emmGrid` has undergone any transformations (`"log"`, `"response"`, etc.),
-#'         or `regrid`ing, then `prior` must be an `emmGrid` object, as stated above.
+#' - **A numeric vector** - `prior` should also be a _numeric vector_, representing the prior-estimate.
+#' - **A data frame** - `prior` should also be a _data frame_, representing the prior-estimates, in matching column order.
+#'   - If `rvar_col` is specified, `prior` should be _the name of an `rvar` column_ that represents the prior-estimates.
+#' - **Supported Bayesian model (`stanreg`, `brmsfit`, etc.)**
+#'   - `prior` should be _a model an equivalent model with MCMC samples from the priors **only**_. See [unupdate()].
+#'   - If `prior` is set to `NULL`, [unupdate()] is called internally (not supported for `brmsfit_multiple` model).
+#' - **Output from a `{marginaleffects}` function** - `prior` should also be _an equivalent output_ from a `{marginaleffects}` function based on a prior-model
+#'  (See [unupdate()]).
+#' - **Output from an `{emmeans}` function**
+#'   - `prior` should also be _an equivalent output_ from an `{emmeans}` function based on a prior-model (See [unupdate()]).
+#'   - `prior` can also be _the original (posterior) model_, in which case the function
+#'      will try to "unupdate" the estimates (not supported if the estimates have undergone
+#'      any transformations -- `"log"`, `"response"`, etc. -- or any `regrid`ing).
 #'
 #' @section Interpreting Bayes Factors:
 #' A Bayes factor greater than 1 can be interpreted as evidence against the
@@ -193,8 +189,8 @@ bayesfactor_parameters <- function(posterior,
                                    prior = NULL,
                                    direction = "two-sided",
                                    null = 0,
-                                   verbose = TRUE,
-                                   ...) {
+                                   ...,
+                                   verbose = TRUE) {
   UseMethod("bayesfactor_parameters")
 }
 
@@ -204,8 +200,8 @@ bayesfactor_pointnull <- function(posterior,
                                   prior = NULL,
                                   direction = "two-sided",
                                   null = 0,
-                                  verbose = TRUE,
-                                  ...) {
+                                  ...,
+                                  verbose = TRUE) {
   if (length(null) > 1L && verbose) {
     insight::format_alert("`null` is a range - computing a ROPE based Bayes factor.")
   }
@@ -226,8 +222,8 @@ bayesfactor_rope <- function(posterior,
                              prior = NULL,
                              direction = "two-sided",
                              null = rope_range(posterior, verbose = FALSE),
-                             verbose = TRUE,
-                             ...) {
+                             ...,
+                             verbose = TRUE) {
   if (length(null) < 2 && verbose) {
     insight::format_alert("'null' is a point - computing a Savage-Dickey (point null) Bayes factor.")
   }
@@ -260,8 +256,8 @@ bayesfactor_parameters.numeric <- function(posterior,
                                            prior = NULL,
                                            direction = "two-sided",
                                            null = 0,
-                                           verbose = TRUE,
-                                           ...) {
+                                           ...,
+                                           verbose = TRUE) {
   # nm <- insight::safe_deparse(substitute(posterior)
 
   if (is.null(prior)) {
@@ -293,11 +289,11 @@ bayesfactor_parameters.stanreg <- function(posterior,
                                            prior = NULL,
                                            direction = "two-sided",
                                            null = 0,
-                                           verbose = TRUE,
                                            effects = c("fixed", "random", "all"),
                                            component = c("conditional", "location", "smooth_terms", "sigma", "zi", "zero_inflated", "all"),
                                            parameters = NULL,
-                                           ...) {
+                                           ...,
+                                           verbose = TRUE) {
   cleaned_parameters <- insight::clean_parameters(posterior)
   effects <- match.arg(effects)
   component <- match.arg(component)
@@ -339,8 +335,8 @@ bayesfactor_parameters.blavaan <- function(posterior,
                                            prior = NULL,
                                            direction = "two-sided",
                                            null = 0,
-                                           verbose = TRUE,
-                                           ...) {
+                                           ...,
+                                           verbose = TRUE) {
   cleaned_parameters <- insight::clean_parameters(posterior)
 
   samps <- .clean_priors_and_posteriors(posterior, prior,
@@ -372,8 +368,8 @@ bayesfactor_parameters.emmGrid <- function(posterior,
                                            prior = NULL,
                                            direction = "two-sided",
                                            null = 0,
-                                           verbose = TRUE,
-                                           ...) {
+                                           ...,
+                                           verbose = TRUE) {
   samps <- .clean_priors_and_posteriors(
     posterior,
     prior,
@@ -406,13 +402,33 @@ bayesfactor_parameters.comparisons <- bayesfactor_parameters.emmGrid
 
 
 #' @rdname bayesfactor_parameters
+#' @inheritParams p_direction
 #' @export
 bayesfactor_parameters.data.frame <- function(posterior,
                                               prior = NULL,
                                               direction = "two-sided",
                                               null = 0,
-                                              verbose = TRUE,
-                                              ...) {
+                                              rvar_col = NULL,
+                                              ...,
+                                              verbose = TRUE) {
+  x_rvar <- .possibly_extract_rvar_col(posterior, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::bayesfactor_parameters
+    cl$posterior <- x_rvar
+    cl$rvar_col <- NULL
+    prior_rvar <- .possibly_extract_rvar_col(posterior, prior)
+    if (length(prior_rvar) > 0L) {
+      cl$prior <- prior_rvar
+    }
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(posterior))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, posterior))
+  }
+
   # find direction
   direction <- .get_direction(direction)
 
@@ -469,11 +485,11 @@ bayesfactor_parameters.draws <- function(posterior,
                                          prior = NULL,
                                          direction = "two-sided",
                                          null = 0,
-                                         verbose = TRUE,
-                                         ...) {
+                                         ...,
+                                         verbose = TRUE) {
   bayesfactor_parameters(
     .posterior_draws_to_df(posterior),
-    prior = prior,
+    prior = if (!is.null(prior)) .posterior_draws_to_df(prior),
     direction = direction,
     null = null,
     verbose = verbose,
diff --git a/R/bayesfactor_restricted.R b/R/bayesfactor_restricted.R
index c5ddf4f46..1be2a7298 100644
--- a/R/bayesfactor_restricted.R
+++ b/R/bayesfactor_restricted.R
@@ -117,7 +117,7 @@
 #' Retrieved from https://richarddmorey.org/category/order-restrictions/.
 #'
 #' @export
-bayesfactor_restricted <- function(posterior, hypothesis, prior = NULL, verbose = TRUE, ...) {
+bayesfactor_restricted <- function(posterior, ...) {
   UseMethod("bayesfactor_restricted")
 }
 
@@ -195,7 +195,23 @@ bayesfactor_restricted.predictions <- bayesfactor_restricted.emmGrid
 bayesfactor_restricted.comparisons <- bayesfactor_restricted.emmGrid
 
 #' @export
-bayesfactor_restricted.data.frame <- function(posterior, hypothesis, prior = NULL, ...) {
+#' @rdname bayesfactor_restricted
+#' @inheritParams p_direction
+bayesfactor_restricted.data.frame <- function(posterior, hypothesis, prior = NULL, rvar_col = NULL, ...) {
+  x_rvar <- .possibly_extract_rvar_col(posterior, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::bayesfactor_restricted
+    cl$posterior <- x_rvar
+    cl$rvar_col <- NULL
+    prior_rvar <- .possibly_extract_rvar_col(posterior, prior)
+    if (length(prior_rvar) > 0L) {
+      cl$prior <- prior_rvar
+    }
+    return(eval.parent(cl))
+  }
+
+
   p_hypothesis <- parse(text = hypothesis)
 
   if (is.null(prior)) {
@@ -251,7 +267,11 @@ bayesfactor_restricted.data.frame <- function(posterior, hypothesis, prior = NUL
 
 #' @export
 bayesfactor_restricted.draws <- function(posterior, hypothesis, prior = NULL, ...) {
-  bayesfactor_restricted(.posterior_draws_to_df(posterior), hypothesis = hypothesis, prior = prior, ...)
+  bayesfactor_restricted(.posterior_draws_to_df(posterior),
+    hypothesis = hypothesis,
+    prior = if (!is.null(prior)) .posterior_draws_to_df(prior),
+    ...
+  )
 }
 
 #' @export
diff --git a/R/bci.R b/R/bci.R
index 190c229c6..b5e61ccf4 100644
--- a/R/bci.R
+++ b/R/bci.R
@@ -42,10 +42,26 @@ bci.numeric <- function(x, ci = 0.95, verbose = TRUE, ...) {
 
 
 #' @rdname bci
+#' @inheritParams p_direction
 #' @export
-bci.data.frame <- function(x, ci = 0.95, verbose = TRUE, ...) {
+bci.data.frame <- function(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...) {
+  obj_name <- insight::safe_deparse_symbol(substitute(x))
+
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::bci
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x, long = length(ci) > 1L))
+  }
+
   dat <- .compute_interval_dataframe(x = x, ci = ci, verbose = verbose, fun = "bci")
-  attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
+  attr(dat, "object_name") <- obj_name
   dat
 }
 
@@ -168,7 +184,7 @@ bci.sim <- function(x, ci = 0.95, parameters = NULL, verbose = TRUE, ...) {
 bci.emmGrid <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xdf <- insight::get_parameters(x)
   dat <- bci(xdf, ci = ci, verbose = verbose, ...)
-  dat <- .append_datagrid(dat, x)
+  dat <- .append_datagrid(dat, x, long = length(ci) > 1L)
   attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   dat
 }
@@ -181,7 +197,7 @@ bci.emm_list <- bci.emmGrid
 bci.slopes <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xrvar <- .get_marginaleffects_draws(x)
   dat <- bci(xrvar, ci = ci, verbose = verbose, ...)
-  dat <- .append_datagrid(dat, x)
+  dat <- .append_datagrid(dat, x, long = length(ci) > 1L)
   attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   dat
 }
diff --git a/R/ci.R b/R/ci.R
index be5e91f05..0011fcf62 100644
--- a/R/ci.R
+++ b/R/ci.R
@@ -157,8 +157,25 @@ ci.numeric <- function(x, ci = 0.95, method = "ETI", verbose = TRUE, BF = 1, ...
 
 
 #' @rdname ci
+#' @inheritParams p_direction
 #' @export
-ci.data.frame <- ci.numeric
+ci.data.frame <- function(x, ci = 0.95, method = "ETI", BF = 1, rvar_col = NULL, verbose = TRUE, ...) {
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::ci
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(x))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x, long = length(ci) > 1L))
+  }
+
+  .ci_bayesian(x, ci = ci, method = method, verbose = verbose, BF = BF, ...)
+}
 
 
 #' @export
@@ -181,7 +198,7 @@ ci.emmGrid <- function(x, ci = NULL, ...) {
   if (is.null(ci)) ci <- 0.95
   xdf <- insight::get_parameters(x)
   out <- ci(xdf, ci = ci, ...)
-  out <- .append_datagrid(out, x)
+  out <- .append_datagrid(out, x, long = length(ci) > 1L)
   out
 }
 
@@ -200,7 +217,7 @@ ci.slopes <- function(x, ci = NULL, ...) {
   if (is.null(ci)) ci <- 0.95
   xrvar <- .get_marginaleffects_draws(x)
   out <- ci(xrvar, ci = ci, ...)
-  out <- .append_datagrid(out, x)
+  out <- .append_datagrid(out, x, long = length(ci) > 1L)
   out
 }
 
diff --git a/R/describe_posterior.R b/R/describe_posterior.R
index 56b3c8432..d238c23d9 100644
--- a/R/describe_posterior.R
+++ b/R/describe_posterior.R
@@ -383,10 +383,16 @@ describe_posterior.default <- function(posterior, ...) {
       dot_args <- list(...)
       dot_args$verbose <- !"rope" %in% test
       test_equi <- .prepare_output(
-        equivalence_test(x_df,
-          range = rope_range,
-          ci = rope_ci,
-          dot_args
+        do.call(
+          equivalence_test,
+          c(
+            dot_args,
+            list(
+              x = x_df,
+              range = rope_range,
+              ci = rope_ci
+            )
+          )
         ),
         cleaned_parameters,
         is_stanmvreg
@@ -459,7 +465,7 @@ describe_posterior.default <- function(posterior, ...) {
     test_psig$.rowid <- seq_len(nrow(test_psig))
   } else if (!all(is.na(test_rope$Parameter))) {
     test_rope$.rowid <- seq_len(nrow(test_rope))
-  } else if (!all(is.na(test_bf$Parameter))) {
+  } else if (!all(is.na(test_bf$Parameter))) { # nolint
     test_bf$.rowid <- seq_len(nrow(test_bf))
   } else {
     estimates$.rowid <- seq_len(nrow(estimates))
@@ -568,7 +574,60 @@ describe_posterior.double <- describe_posterior.numeric
 
 
 #' @export
-describe_posterior.data.frame <- describe_posterior.numeric
+#' @rdname describe_posterior
+#' @inheritParams p_direction
+describe_posterior.data.frame <- function(posterior,
+                                          centrality = "median",
+                                          dispersion = FALSE,
+                                          ci = 0.95,
+                                          ci_method = "eti",
+                                          test = c("p_direction", "rope"),
+                                          rope_range = "default",
+                                          rope_ci = 0.95,
+                                          keep_iterations = FALSE,
+                                          bf_prior = NULL,
+                                          BF = 1,
+                                          rvar_col = NULL,
+                                          verbose = TRUE,
+                                          ...) {
+  x_rvar <- .possibly_extract_rvar_col(posterior, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::describe_posterior
+    cl$posterior <- x_rvar
+    cl$rvar_col <- NULL
+    prior_rvar <- .possibly_extract_rvar_col(posterior, bf_prior)
+    if (length(prior_rvar) > 0L) {
+      cl$bf_prior <- prior_rvar
+    }
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(posterior))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, posterior))
+  }
+
+
+  out <- .describe_posterior(
+    posterior,
+    centrality = centrality,
+    dispersion = dispersion,
+    ci = ci,
+    ci_method = ci_method,
+    test = test,
+    rope_range = rope_range,
+    rope_ci = rope_ci,
+    keep_iterations = keep_iterations,
+    bf_prior = bf_prior,
+    BF = BF,
+    verbose = verbose,
+    ...
+  )
+
+  class(out) <- unique(c("describe_posterior", "see_describe_posterior", class(out)))
+  out
+}
 
 
 #' @export
diff --git a/R/equivalence_test.R b/R/equivalence_test.R
index 163db680d..57d4b2058 100644
--- a/R/equivalence_test.R
+++ b/R/equivalence_test.R
@@ -145,8 +145,24 @@ equivalence_test.numeric <- function(x, range = "default", ci = 0.95, verbose =
 
 
 #' @rdname equivalence_test
+#' @inheritParams p_direction
 #' @export
-equivalence_test.data.frame <- function(x, range = "default", ci = 0.95, verbose = TRUE, ...) {
+equivalence_test.data.frame <- function(x, range = "default", ci = 0.95, rvar_col = NULL, verbose = TRUE, ...) {
+  obj_name <- insight::safe_deparse_symbol(substitute(x))
+
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::equivalence_test
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x))
+  }
+
   l <- insight::compact_list(lapply(
     x,
     equivalence_test,
@@ -163,7 +179,7 @@ equivalence_test.data.frame <- function(x, range = "default", ci = 0.95, verbose
   )
   row.names(out) <- NULL
 
-  attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
+  attr(out, "object_name") <- obj_name
   class(out) <- unique(c("equivalence_test", "see_equivalence_test_df", class(out)))
 
   out
diff --git a/R/estimate_density.R b/R/estimate_density.R
index a8307c37c..9ee4c424f 100644
--- a/R/estimate_density.R
+++ b/R/estimate_density.R
@@ -221,6 +221,7 @@ estimate_density.numeric <- function(x,
 
 
 #' @rdname estimate_density
+#' @inheritParams p_direction
 #' @export
 estimate_density.data.frame <- function(x,
                                         method = "kernel",
@@ -232,7 +233,25 @@ estimate_density.data.frame <- function(x,
                                         select = NULL,
                                         by = NULL,
                                         at = NULL,
+                                        rvar_col = NULL,
                                         ...) {
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::estimate_density
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(x))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    out <- .append_datagrid(out, x, long = TRUE)
+    class(out) <- .set_density_class(out)
+    return(out)
+  }
+
+
   # Sanity
   if (!is.null(at)) {
     insight::format_warning(paste0(
@@ -375,11 +394,7 @@ estimate_density.emmGrid <- function(x,
     bw = bw, ...
   )
 
-  # This doesn't use .append_datagrid because we get a non-grid output
-  dgrid <- insight::get_datagrid(x)
-  dgrid$Parameter <- unique(out$Parameter)
-  out <- datawizard::data_join(dgrid, out, by = "Parameter")
-  out$Parameter <- NULL
+  out <- .append_datagrid(out, x, long = TRUE)
   class(out) <- .set_density_class(out)
   out
 }
@@ -403,11 +418,7 @@ estimate_density.slopes <- function(x,
     bw = bw, ...
   )
 
-  # This doesn't use .append_datagrid because we get a non-grid output
-  dgrid <- insight::get_datagrid(x)
-  dgrid$Parameter <- unique(out$Parameter)
-  out <- datawizard::data_join(dgrid, out, by = "Parameter")
-  out$Parameter <- NULL
+  out <- .append_datagrid(out, x, long = TRUE)
   class(out) <- .set_density_class(out)
   out
 }
diff --git a/R/eti.R b/R/eti.R
index 9e2b64935..4bcf75af1 100644
--- a/R/eti.R
+++ b/R/eti.R
@@ -66,9 +66,27 @@ eti.numeric <- function(x, ci = 0.95, verbose = TRUE, ...) {
 
 
 #' @export
-eti.data.frame <- function(x, ci = 0.95, verbose = TRUE, ...) {
+#' @rdname eti
+#' @inheritParams p_direction
+eti.data.frame <- function(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...) {
+  obj_name <- insight::safe_deparse_symbol(substitute(x))
+
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::eti
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x, long = length(ci) > 1L))
+  }
+
+
   dat <- .compute_interval_dataframe(x = x, ci = ci, verbose = verbose, fun = "eti")
-  attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
+  attr(dat, "object_name") <- obj_name
   dat
 }
 
@@ -175,7 +193,7 @@ eti.sim <- function(x, ci = 0.95, parameters = NULL, verbose = TRUE, ...) {
 eti.emmGrid <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xdf <- insight::get_parameters(x)
   dat <- eti(xdf, ci = ci, verbose = verbose, ...)
-  dat <- .append_datagrid(dat, x)
+  dat <- .append_datagrid(dat, x, long = length(ci) > 1L)
   attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   dat
 }
@@ -187,7 +205,7 @@ eti.emm_list <- eti.emmGrid
 eti.slopes <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xrvar <- .get_marginaleffects_draws(x)
   dat <- eti(xrvar, ci = ci, verbose = verbose, ...)
-  dat <- .append_datagrid(dat, x)
+  dat <- .append_datagrid(dat, x, long = length(ci) > 1L)
   attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   dat
 }
diff --git a/R/hdi.R b/R/hdi.R
index bf1aa4244..0662a2c52 100644
--- a/R/hdi.R
+++ b/R/hdi.R
@@ -143,10 +143,26 @@ hdi.numeric <- function(x, ci = 0.95, verbose = TRUE, ...) {
 
 
 #' @rdname hdi
+#' @inheritParams p_direction
 #' @export
-hdi.data.frame <- function(x, ci = 0.95, verbose = TRUE, ...) {
+hdi.data.frame <- function(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...) {
+  obj_name <- insight::safe_deparse_symbol(substitute(x))
+
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::hdi
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x, long = length(ci) > 1L))
+  }
+
   dat <- .compute_interval_dataframe(x = x, ci = ci, verbose = verbose, fun = "hdi")
-  attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
+  attr(dat, "object_name") <- obj_name
   dat
 }
 
@@ -264,7 +280,7 @@ hdi.sim <- function(x, ci = 0.95, parameters = NULL, verbose = TRUE, ...) {
 hdi.emmGrid <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xdf <- insight::get_parameters(x)
   out <- hdi(xdf, ci = ci, verbose = verbose, ...)
-  out <- .append_datagrid(out, x)
+  out <- .append_datagrid(out, x, long = length(ci) > 1L)
   attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   out
 }
@@ -276,7 +292,7 @@ hdi.emm_list <- hdi.emmGrid
 hdi.slopes <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xrvar <- .get_marginaleffects_draws(x)
   out <- hdi(xrvar, ci = ci, verbose = verbose, ...)
-  out <- .append_datagrid(out, x)
+  out <- .append_datagrid(out, x, long = length(ci) > 1L)
   attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   out
 }
diff --git a/R/map_estimate.R b/R/map_estimate.R
index 68bfafae0..aea57b96c 100644
--- a/R/map_estimate.R
+++ b/R/map_estimate.R
@@ -148,8 +148,23 @@ map_estimate.brmsfit <- function(x, precision = 2^10, method = "kernel", effects
 
 
 #' @rdname map_estimate
+#' @inheritParams p_direction
 #' @export
-map_estimate.data.frame <- function(x, precision = 2^10, method = "kernel", ...) {
+map_estimate.data.frame <- function(x, precision = 2^10, method = "kernel", rvar_col = NULL, ...) {
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::map_estimate
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(x))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x))
+  }
+
   .map_estimate_models(x, precision = precision, method = method)
 }
 
diff --git a/R/p_direction.R b/R/p_direction.R
index 12e8e05b5..c5b1f84e0 100644
--- a/R/p_direction.R
+++ b/R/p_direction.R
@@ -122,6 +122,7 @@
 #' df <- data.frame(replicate(4, rnorm(100)))
 #' p_direction(df)
 #' p_direction(df, method = "kernel")
+#'
 #' \donttest{
 #' # rstanarm models
 #' # -----------------------------------------------
@@ -148,6 +149,14 @@
 #' p_direction(bf)
 #' p_direction(bf, method = "kernel")
 #' }
+#'
+#' @examplesIf requireNamespace("posterior", quietly = TRUE)
+#' # Using "rvar_col"
+#' x <- data.frame(mu = c(0, 0.5, 1), sigma = c(1, 0.5, 0.25))
+#' x$my_rvar <- posterior::rvar_rng(rnorm, 3, mean = x$mu, sd = x$sigma)
+#' x
+#' p_direction(x, rvar_col = "my_rvar")
+#'
 #' @export
 p_direction <- function(x, ...) {
   UseMethod("p_direction")
@@ -187,6 +196,8 @@ p_direction.numeric <- function(x,
 
 
 #' @rdname p_direction
+#' @param rvar_col A single character - the name of an `rvar` column in the data
+#'   frame to be processed. See example in [p_direction()].
 #' @export
 p_direction.data.frame <- function(x,
                                    method = "direct",
@@ -196,48 +207,18 @@ p_direction.data.frame <- function(x,
                                    rvar_col = NULL,
                                    ...) {
   obj_name <- insight::safe_deparse_symbol(substitute(x))
-
-  if (is.null(rvar_col)) {
-    return(.p_direction_df(
-      x,
-      method = method,
-      null = null,
-      as_p = as_p,
-      remove_na = remove_na,
-      obj_name = obj_name,
-      ...
-    ))
-  }
-
-  if (length(rvar_col) != 1L && !rvar_col %in% colnames(x)) {
-    insight::format_error("The `rvar_col` argument must be a single, valid column name.")
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::p_direction
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+    return(.append_datagrid(out, x))
   }
 
-  out <- p_direction(
-    x[[rvar_col]],
-    method = method,
-    null = null,
-    as_p = as_p,
-    remove_na = remove_na,
-    ...
-  )
-
-  x[["pd"]] <- out[["pd"]]
-  attr(x, "object_name") <- obj_name
-  attr(x, "as_p") <- as_p
 
-  x
-}
-
-
-#' @keywords internal
-.p_direction_df <- function(x,
-                            method = "direct",
-                            null = 0,
-                            as_p = FALSE,
-                            remove_na = TRUE,
-                            obj_name = NULL,
-                            ...) {
   x <- .select_nums(x)
 
   if (ncol(x) == 1) {
@@ -282,6 +263,9 @@ p_direction.data.frame <- function(x,
 }
 
 
+
+
+
 #' @export
 p_direction.draws <- function(x,
                               method = "direct",
diff --git a/R/p_map.R b/R/p_map.R
index 60efe1580..a0267e867 100644
--- a/R/p_map.R
+++ b/R/p_map.R
@@ -124,7 +124,23 @@ p_map.get_predicted <- function(x,
 
 
 #' @export
-p_map.data.frame <- function(x, null = 0, precision = 2^10, method = "kernel", ...) {
+#' @rdname p_map
+#' @inheritParams p_direction
+p_map.data.frame <- function(x, null = 0, precision = 2^10, method = "kernel", rvar_col = NULL, ...) {
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::p_map
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(x))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x))
+  }
+
   x <- .select_nums(x)
 
   if (ncol(x) == 1) {
diff --git a/R/p_rope.R b/R/p_rope.R
index 9b2522228..3a8ae940f 100644
--- a/R/p_rope.R
+++ b/R/p_rope.R
@@ -38,7 +38,27 @@ p_rope.numeric <- function(x, range = "default", verbose = TRUE, ...) {
 
 
 #' @export
-p_rope.data.frame <- p_rope.numeric
+#' @rdname p_rope
+#' @inheritParams p_direction
+p_rope.data.frame <- function(x, range = "default", rvar_col = NULL, verbose = TRUE, ...) {
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::p_rope
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(x))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x))
+  }
+
+  out <- .p_rope(rope(x, range = range, ci = 1, verbose = verbose, ...))
+  attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
+  out
+}
 
 
 #' @export
@@ -99,7 +119,7 @@ p_rope.stanreg <- function(x,
                              "auxiliary"
                            ),
                            parameters = NULL,
-                           verbose = verbose,
+                           verbose = TRUE,
                            ...) {
   out <- .p_rope(rope(
     x,
diff --git a/R/p_significance.R b/R/p_significance.R
index 68a179fa2..46c8bd1c5 100644
--- a/R/p_significance.R
+++ b/R/p_significance.R
@@ -109,8 +109,25 @@ p_significance.get_predicted <- function(x,
 
 
 #' @export
-p_significance.data.frame <- function(x, threshold = "default", ...) {
+#' @rdname p_significance
+#' @inheritParams p_direction
+p_significance.data.frame <- function(x, threshold = "default", rvar_col = NULL, ...) {
   obj_name <- insight::safe_deparse_symbol(substitute(x))
+
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::p_significance
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x))
+  }
+
+
   threshold <- .select_threshold_ps(threshold = threshold)
   x <- .select_nums(x)
 
@@ -253,18 +270,18 @@ p_significance.stanreg <- function(x,
   component <- match.arg(component)
   threshold <- .select_threshold_ps(model = x, threshold = threshold, verbose = verbose)
 
-  data <- p_significance(
+  result <- p_significance(
     insight::get_parameters(x, effects = effects, component = component, parameters = parameters),
     threshold = threshold
   )
 
   cleaned_parameters <- insight::clean_parameters(x)
-  out <- .prepare_output(data, cleaned_parameters, inherits(x, "stanmvreg"))
+  out <- .prepare_output(result, cleaned_parameters, inherits(x, "stanmvreg"))
 
   attr(out, "clean_parameters") <- cleaned_parameters
   attr(out, "threshold") <- threshold
   attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
-  class(out) <- class(data)
+  class(out) <- class(result)
 
   out
 }
@@ -289,18 +306,18 @@ p_significance.brmsfit <- function(x,
   component <- match.arg(component)
   threshold <- .select_threshold_ps(model = x, threshold = threshold, verbose = verbose)
 
-  data <- p_significance(
+  result <- p_significance(
     insight::get_parameters(x, effects = effects, component = component, parameters = parameters),
     threshold = threshold
   )
 
   cleaned_parameters <- insight::clean_parameters(x)
-  out <- .prepare_output(data, cleaned_parameters)
+  out <- .prepare_output(result, cleaned_parameters)
 
   attr(out, "clean_parameters") <- cleaned_parameters
   attr(out, "threshold") <- threshold
   attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
-  class(out) <- class(data)
+  class(out) <- class(result)
 
   out
 }
diff --git a/R/point_estimate.R b/R/point_estimate.R
index 15b305766..c11a6040a 100644
--- a/R/point_estimate.R
+++ b/R/point_estimate.R
@@ -83,7 +83,7 @@ point_estimate.numeric <- function(x, centrality = "all", dispersion = FALSE, th
     estimate_list <- centrality
   }
 
-  out <- data.frame(".temp" = 0)
+  out <- data.frame(.temp = 0)
 
   # Median
   if ("median" %in% estimate_list) {
@@ -130,7 +130,29 @@ point_estimate.numeric <- function(x, centrality = "all", dispersion = FALSE, th
 
 
 #' @export
-point_estimate.data.frame <- function(x, centrality = "all", dispersion = FALSE, threshold = 0.1, ...) {
+#' @rdname point_estimate
+#' @inheritParams p_direction
+point_estimate.data.frame <- function(x,
+                                      centrality = "all",
+                                      dispersion = FALSE,
+                                      threshold = 0.1,
+                                      rvar_col = NULL,
+                                      ...) {
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::point_estimate
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(x))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x))
+  }
+
+
   x <- .select_nums(x)
 
   if (ncol(x) == 1) {
@@ -140,7 +162,7 @@ point_estimate.data.frame <- function(x, centrality = "all", dispersion = FALSE,
     estimates <- do.call(rbind, estimates)
   }
 
-  out <- cbind(data.frame("Parameter" = names(x), stringsAsFactors = FALSE), estimates)
+  out <- cbind(data.frame(Parameter = names(x), stringsAsFactors = FALSE), estimates)
   rownames(out) <- NULL
   attr(out, "data") <- x
   attr(out, "centrality") <- centrality
@@ -193,7 +215,11 @@ point_estimate.BGGM <- point_estimate.bcplm
 
 
 #' @export
-point_estimate.bamlss <- function(x, centrality = "all", dispersion = FALSE, component = c("conditional", "location", "all"), ...) {
+point_estimate.bamlss <- function(x,
+                                  centrality = "all",
+                                  dispersion = FALSE,
+                                  component = c("conditional", "location", "all"),
+                                  ...) {
   component <- match.arg(component)
   out <- point_estimate(
     insight::get_parameters(x, component = component),
diff --git a/R/print.equivalence_test.R b/R/print.equivalence_test.R
index 978253fc7..044d4554b 100644
--- a/R/print.equivalence_test.R
+++ b/R/print.equivalence_test.R
@@ -35,7 +35,7 @@ print.equivalence_test <- function(x, digits = 2, ...) {
 
   ci <- unique(x$CI)
   keep.columns <- c(
-    attr(x, "grid_cols"), "Parameter", "Effects", "Component",
+    attr(x, "idvars"), "Parameter", "Effects", "Component",
     "ROPE_Equivalence", "ROPE_Percentage", "CI", "HDI"
   )
 
@@ -84,7 +84,7 @@ print.equivalence_test <- function(x, digits = 2, ...) {
 
 
 .dynGet <- function(x,
-                    ifnotfound = stop(gettextf("%s not found", sQuote(x)), domain = NA),
+                    ifnotfound = stop(gettextf("%s not found", sQuote(x)), domain = NA, call. = FALSE),
                     minframe = 1L,
                     inherits = FALSE) {
   x <- insight::safe_deparse(x)
diff --git a/R/print.rope.R b/R/print.rope.R
index 43cc9a2c2..d56d5369f 100644
--- a/R/print.rope.R
+++ b/R/print.rope.R
@@ -28,7 +28,7 @@ print.rope <- function(x, digits = 2, ...) {
 
   # These are the base columns we want to print
   cols <- c(
-    attr(x, "grid_cols"), "Parameter", "ROPE_Percentage", "Effects", "Component",
+    attr(x, "idvars"), "Parameter", "ROPE_Percentage", "Effects", "Component",
     if (is_multivariate) c("ROPE_low", "ROPE_high")
   )
 
diff --git a/R/rope.R b/R/rope.R
index c9ba83f41..cdb692584 100644
--- a/R/rope.R
+++ b/R/rope.R
@@ -214,7 +214,24 @@ rope.get_predicted <- function(x,
 
 
 #' @export
-rope.data.frame <- function(x, range = "default", ci = 0.95, ci_method = "ETI", verbose = TRUE, ...) {
+#' @rdname rope
+#' @inheritParams p_direction
+rope.data.frame <- function(x, range = "default", ci = 0.95, ci_method = "ETI", rvar_col = NULL, verbose = TRUE, ...) {
+  obj_name <- insight::safe_deparse_symbol(substitute(x))
+
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::rope
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x))
+  }
+
   out <- .prepare_rope_df(x, range, ci, ci_method, verbose)
   HDI_area_attributes <- insight::compact_list(out$HDI_area)
   dat <- data.frame(
@@ -225,7 +242,7 @@ rope.data.frame <- function(x, range = "default", ci = 0.95, ci_method = "ETI",
   row.names(dat) <- NULL
 
   attr(dat, "HDI_area") <- HDI_area_attributes
-  attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
+  attr(dat, "object_name") <- obj_name
 
   class(dat) <- c("rope", "see_rope", "data.frame")
   dat
diff --git a/R/si.R b/R/si.R
index 497a60f74..69bb9bde3 100644
--- a/R/si.R
+++ b/R/si.R
@@ -168,7 +168,7 @@ si.emmGrid <- function(posterior, prior = NULL,
     BF = BF, verbose = verbose, ...
   )
 
-  out <- .append_datagrid(out, posterior)
+  out <- .append_datagrid(out, posterior, long = length(BF) > 1L)
   attr(out, "ci_method") <- "SI"
   attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(posterior))
   out
@@ -221,8 +221,27 @@ si.get_predicted <- function(posterior, prior = NULL, BF = 1, use_iterations = F
 
 
 #' @rdname si
+#' @inheritParams p_direction
 #' @export
-si.data.frame <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, ...) {
+si.data.frame <- function(posterior, prior = NULL, BF = 1, rvar_col = NULL, verbose = TRUE, ...) {
+  x_rvar <- .possibly_extract_rvar_col(posterior, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::si
+    cl$posterior <- x_rvar
+    cl$rvar_col <- NULL
+    prior_rvar <- .possibly_extract_rvar_col(posterior, prior)
+    if (length(prior_rvar) > 0L) {
+      cl$prior <- prior_rvar
+    }
+    out <- eval.parent(cl)
+
+    obj_name <- insight::safe_deparse_symbol(substitute(posterior))
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, posterior, long = length(BF) > 1L))
+  }
+
   if (is.null(prior)) {
     prior <- posterior
     insight::format_warning(
@@ -255,7 +274,10 @@ si.data.frame <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, ...)
 
 #' @export
 si.draws <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, ...) {
-  si(.posterior_draws_to_df(posterior), prior = prior, BF = BF, verbose = verbose, ...)
+  si(.posterior_draws_to_df(posterior),
+    prior = if (!is.null(prior)) .posterior_draws_to_df(prior),
+    BF = BF, verbose = verbose, ...
+  )
 }
 
 #' @export
@@ -276,7 +298,7 @@ si.rvar <- si.draws
     )
   }
 
-  out <- data.frame(
+  data.frame(
     Parameter = colnames(posterior),
     CI = BF,
     CI_low = sis[, 1],
@@ -310,12 +332,12 @@ si.rvar <- si.draws
 
   f_prior <- .logspline(prior, ...)
   f_posterior <- .logspline(posterior, ...)
-  d_prior <- logspline::dlogspline(x_axis, f_prior)
-  d_posterior <- logspline::dlogspline(x_axis, f_posterior)
+  d_prior <- logspline::dlogspline(x_axis, f_prior, log = TRUE)
+  d_posterior <- logspline::dlogspline(x_axis, f_posterior, log = TRUE)
 
-  relative_d <- d_posterior / d_prior
+  relative_d <- d_posterior - d_prior
 
-  crit <- relative_d >= BF
+  crit <- relative_d >= log(BF)
 
   cp <- rle(stats::na.omit(crit))
   if (length(cp$lengths) > 3 && verbose) {
diff --git a/R/spi.R b/R/spi.R
index 1e38c9104..53429f90f 100644
--- a/R/spi.R
+++ b/R/spi.R
@@ -66,9 +66,26 @@ spi.numeric <- function(x, ci = 0.95, verbose = TRUE, ...) {
 
 
 #' @export
-spi.data.frame <- function(x, ci = 0.95, verbose = TRUE, ...) {
+#' @rdname spi
+#' @inheritParams p_direction
+spi.data.frame <- function(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...) {
+  obj_name <- insight::safe_deparse_symbol(substitute(x))
+
+  x_rvar <- .possibly_extract_rvar_col(x, rvar_col)
+  if (length(x_rvar) > 0L) {
+    cl <- match.call()
+    cl[[1]] <- bayestestR::spi
+    cl$x <- x_rvar
+    cl$rvar_col <- NULL
+    out <- eval.parent(cl)
+
+    attr(out, "object_name") <- sprintf('%s[["%s"]]', obj_name, rvar_col)
+
+    return(.append_datagrid(out, x, long = length(ci) > 1L))
+  }
+
   dat <- .compute_interval_dataframe(x = x, ci = ci, verbose = verbose, fun = "spi")
-  attr(dat, "object_name") <- insight::safe_deparse_symbol(substitute(x))
+  attr(dat, "object_name") <- obj_name
   dat
 }
 
@@ -140,7 +157,7 @@ spi.sim <- function(x, ci = 0.95, parameters = NULL, verbose = TRUE, ...) {
 spi.emmGrid <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xdf <- insight::get_parameters(x)
   out <- spi(xdf, ci = ci, verbose = verbose, ...)
-  out <- .append_datagrid(out, x)
+  out <- .append_datagrid(out, x, long = length(ci) > 1L)
   attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   out
 }
@@ -152,7 +169,7 @@ spi.emm_list <- spi.emmGrid
 spi.slopes <- function(x, ci = 0.95, verbose = TRUE, ...) {
   xrvar <- .get_marginaleffects_draws(x)
   out <- spi(xrvar, ci = ci, verbose = verbose, ...)
-  out <- .append_datagrid(out, x)
+  out <- .append_datagrid(out, x, long = length(ci) > 1L)
   attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x))
   out
 }
@@ -335,11 +352,7 @@ spi.get_predicted <- function(x, ci = 0.95, use_iterations = FALSE, verbose = TR
   }
 
   # output
-  data.frame(
-    "CI" = ci,
-    "CI_low" = x.l,
-    "CI_high" = x.u
-  )
+  data.frame(CI = ci, CI_low = x.l, CI_high = x.u)
 }
 
 .spi_lower <- function(bw, n.sims, k, l, dens, x) {
@@ -416,7 +429,7 @@ spi.get_predicted <- function(x, ci = 0.95, use_iterations = FALSE, verbose = TR
   w.l <- quadprog::solve.QP(D.l, d.l, A.l, c(1, rep(0, range_ll_lu + 2)), range_ll_lu)
   x.l <- w.l$solution %*% x[l.l:l.u]
 
-  return(x.l)
+  x.l
 }
 
 .spi_upper <- function(bw, n.sims, ui, u, dens, x) {
diff --git a/R/utils.R b/R/utils.R
index 374a04fdc..876e7debd 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -174,30 +174,100 @@
 }
 
 #' @keywords internal
-.append_datagrid <- function(results, object) {
+.append_datagrid <- function(results, object, long = FALSE) {
+  UseMethod(".append_datagrid", object = object)
+}
+
+#' @keywords internal
+.append_datagrid.emmGrid <- function(results, object, long = FALSE) {
   # results is assumed to be a data frame with "Parameter" column
   # object is an emmeans / marginalefeects that results is based on
 
   all_attrs <- attributes(results) # save attributes for later
+  all_class <- class(results)
+
+  datagrid <- insight::get_datagrid(object)
+  grid_names <- colnames(datagrid)
 
-  grid <- insight::get_datagrid(object)
-  grid_names <- colnames(grid)
+  if (long) {
+    datagrid$Parameter <- unique(results$Parameter)
+    results <- datawizard::data_merge(datagrid, results, by = "Parameter")
+    results$Parameter <- NULL
+    class(results) <- all_class
+  } else {
+    results[colnames(datagrid)] <- datagrid
+    results$Parameter <- NULL
+    results <- results[, c(grid_names, setdiff(colnames(results), grid_names)), drop = FALSE]
+
+    # add back attributes
+    most_attrs <- all_attrs[setdiff(names(all_attrs), names(attributes(datagrid)))]
+    attributes(results)[names(most_attrs)] <- most_attrs
+  }
 
-  results[colnames(grid)] <- grid
-  results$Parameter <- NULL
-  results <- results[, c(grid_names, setdiff(colnames(results), grid_names)), drop = FALSE]
 
-  # add back attributes
-  most_attrs <- all_attrs[setdiff(names(all_attrs), names(attributes(grid)))]
-  attributes(results)[names(most_attrs)] <- most_attrs
 
-  attr(results, "grid_cols") <- grid_names
+  attr(results, "idvars") <- grid_names
   results
 }
 
+.append_datagrid.emm_list <- .append_datagrid.emmGrid
+
+.append_datagrid.slopes <- .append_datagrid.emmGrid
+
+.append_datagrid.predictions <- .append_datagrid.emmGrid
+
+.append_datagrid.comparisons <- .append_datagrid.emmGrid
+
+.append_datagrid.data.frame <- function(results, object, long = FALSE) {
+  # results is assumed to be a data frame with "Parameter" column
+  # object is a data frame with an rvar column that results is based on
+
+  all_attrs <- attributes(results) # save attributes for later
+  all_class <- class(results)
+
+  is_rvar <- vapply(object, inherits, FUN.VALUE = logical(1), "rvar")
+  grid_names <- colnames(object)[!is_rvar]
+  datagrid <- data.frame(object[, grid_names, drop = FALSE])
+
+  if (long) {
+    datagrid$Parameter <- unique(results$Parameter)
+    results <- datawizard::data_merge(datagrid, results, by = "Parameter")
+    results$Parameter <- NULL
+    class(results) <- all_class
+  } else {
+    results[grid_names] <- object[grid_names]
+    results$Parameter <- NULL
+    results <- results[, c(grid_names, setdiff(colnames(results), grid_names)), drop = FALSE]
+
+    # add back attributes
+    most_attrs <- all_attrs[setdiff(names(all_attrs), names(attributes(object)))]
+    attributes(results)[names(most_attrs)] <- most_attrs
+  }
+
+  attr(results, "idvars") <- grid_names
+  results
+}
+
+
 #' @keywords internal
 .get_marginaleffects_draws <- function(object) {
   # errors and checks are handled by marginaleffects
   insight::check_if_installed("marginaleffects")
   data.frame(marginaleffects::posterior_draws(object, shape = "DxP"))
 }
+
+#' @keywords internal
+.possibly_extract_rvar_col <- function(df, rvar_col) {
+  if (missing(rvar_col) || is.null(rvar_col)) {
+    return(NULL)
+  }
+
+  if (is.character(rvar_col) &&
+    length(rvar_col) == 1L &&
+    rvar_col %in% colnames(df) &&
+    inherits(df[[rvar_col]], "rvar")) {
+    return(df[[rvar_col]])
+  }
+
+  insight::format_error("The `rvar_col` argument must be a single, valid column name.")
+}
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 09e32ff64..91e5ec922 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -200,6 +200,7 @@ treedepth
 tweedie
 un
 underbrace
+unupdate
 versicolor
 versicolors
 virginica
diff --git a/man/bayesfactor_parameters.Rd b/man/bayesfactor_parameters.Rd
index a29ee8d37..fbc97e667 100644
--- a/man/bayesfactor_parameters.Rd
+++ b/man/bayesfactor_parameters.Rd
@@ -19,8 +19,8 @@ bayesfactor_parameters(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 bayesfactor_pointnull(
@@ -28,8 +28,8 @@ bayesfactor_pointnull(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 bayesfactor_rope(
@@ -37,8 +37,8 @@ bayesfactor_rope(
   prior = NULL,
   direction = "two-sided",
   null = rope_range(posterior, verbose = FALSE),
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 bf_parameters(
@@ -46,8 +46,8 @@ bf_parameters(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 bf_pointnull(
@@ -55,8 +55,8 @@ bf_pointnull(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 bf_rope(
@@ -64,8 +64,8 @@ bf_rope(
   prior = NULL,
   direction = "two-sided",
   null = rope_range(posterior, verbose = FALSE),
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 \method{bayesfactor_parameters}{numeric}(
@@ -73,8 +73,8 @@ bf_rope(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 \method{bayesfactor_parameters}{stanreg}(
@@ -82,12 +82,12 @@ bf_rope(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
   effects = c("fixed", "random", "all"),
   component = c("conditional", "location", "smooth_terms", "sigma", "zi",
     "zero_inflated", "all"),
   parameters = NULL,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 \method{bayesfactor_parameters}{brmsfit}(
@@ -95,12 +95,12 @@ bf_rope(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
   effects = c("fixed", "random", "all"),
   component = c("conditional", "location", "smooth_terms", "sigma", "zi",
     "zero_inflated", "all"),
   parameters = NULL,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 \method{bayesfactor_parameters}{blavaan}(
@@ -108,8 +108,8 @@ bf_rope(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
-  ...
+  ...,
+  verbose = TRUE
 )
 
 \method{bayesfactor_parameters}{data.frame}(
@@ -117,8 +117,9 @@ bf_rope(
   prior = NULL,
   direction = "two-sided",
   null = 0,
-  verbose = TRUE,
-  ...
+  rvar_col = NULL,
+  ...,
+  verbose = TRUE
 )
 }
 \arguments{
@@ -135,11 +136,11 @@ tailed) or \code{1}, \code{"right"} (right tailed).}
 \item{null}{Value of the null, either a scalar (for point-null) or a range
 (for a interval-null).}
 
-\item{verbose}{Toggle off warnings.}
-
 \item{...}{Arguments passed to and from other methods. (Can be used to pass
 arguments to internal \code{\link[logspline:logspline]{logspline::logspline()}}.)}
 
+\item{verbose}{Toggle off warnings.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
@@ -152,6 +153,9 @@ that should be returned. Meta-parameters (like \code{lp__} or \code{prior_}) are
 filtered by default, so only parameters that typically appear in the
 \code{summary()} are returned. Use \code{parameters} to select specific parameters
 for the output.}
+
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
 }
 \value{
 A data frame containing the (log) Bayes factor representing evidence
@@ -228,29 +232,27 @@ prior.
 (Note that by default, \code{brms::brm()} uses flat priors for fixed-effects;
 See example below.)
 \cr\cr
-It is important to provide the correct \code{prior} for meaningful results.
+It is important to provide the correct \code{prior} for meaningful results,
+to match the \code{posterior}-type input:
+\itemize{
+\item \strong{A numeric vector} - \code{prior} should also be a \emph{numeric vector}, representing the prior-estimate.
+\item \strong{A data frame} - \code{prior} should also be a \emph{data frame}, representing the prior-estimates, in matching column order.
 \itemize{
-\item When \code{posterior} is a numerical vector, \code{prior} should also be a numerical vector.
-\item When \code{posterior} is a \code{data.frame}, \code{prior} should also be a \code{data.frame}, with matching column order.
-\item When \code{posterior} is a \code{stanreg}, \code{brmsfit} or other supported Bayesian model:
+\item If \code{rvar_col} is specified, \code{prior} should be \emph{the name of an \code{rvar} column} that represents the prior-estimates.
+}
+\item \strong{Supported Bayesian model (\code{stanreg}, \code{brmsfit}, etc.)}
 \itemize{
-\item \code{prior} can be set to \code{NULL}, in which case prior samples are drawn internally.
-\item \code{prior} can also be a model equivalent to \code{posterior} but with samples from
-the priors \emph{only}. See \code{\link[=unupdate]{unupdate()}}.
-\item \strong{Note:} When \code{posterior} is a \code{brmsfit_multiple} model, \code{prior} \strong{must} be provided.
+\item \code{prior} should be \emph{a model an equivalent model with MCMC samples from the priors \strong{only}}. See \code{\link[=unupdate]{unupdate()}}.
+\item If \code{prior} is set to \code{NULL}, \code{\link[=unupdate]{unupdate()}} is called internally (not supported for \code{brmsfit_multiple} model).
 }
-\item When \code{posterior} is an output from a \code{{marginaleffects}} function, \code{prior} should also be an an output
-from a \code{{marginaleffects}} function equivalent to \code{posterior} but created
-with a model of priors samples \emph{only}.
-\item When \code{posterior} is an \code{emmGrid} / \code{emm_list} object:
+\item \strong{Output from a \code{{marginaleffects}} function} - \code{prior} should also be \emph{an equivalent output} from a \code{{marginaleffects}} function based on a prior-model
+(See \code{\link[=unupdate]{unupdate()}}).
+\item \strong{Output from an \code{{emmeans}} function}
 \itemize{
-\item \code{prior} should also be an \code{emmGrid} / \code{emm_list} object equivalent to \code{posterior} but
-created with a model of priors samples \emph{only}. See \code{\link[=unupdate]{unupdate()}}.
-\item \code{prior} can also be the original (posterior) \emph{model}. If so, the function will try to
-update the \code{emmGrid} / \code{emm_list} to use the \code{\link[=unupdate]{unupdate()}}d prior-model.
-(\emph{This cannot be done for \code{brmsfit} models.})
-\item \strong{Note}: When the \code{emmGrid} has undergone any transformations (\code{"log"}, \code{"response"}, etc.),
-or \code{regrid}ing, then \code{prior} must be an \code{emmGrid} object, as stated above.
+\item \code{prior} should also be \emph{an equivalent output} from an \code{{emmeans}} function based on a prior-model (See \code{\link[=unupdate]{unupdate()}}).
+\item \code{prior} can also be \emph{the original (posterior) model}, in which case the function
+will try to "unupdate" the estimates (not supported if the estimates have undergone
+any transformations -- \code{"log"}, \code{"response"}, etc. -- or any \code{regrid}ing).
 }
 }
 }
diff --git a/man/bayesfactor_restricted.Rd b/man/bayesfactor_restricted.Rd
index 6ed081598..b9ac0cfc1 100644
--- a/man/bayesfactor_restricted.Rd
+++ b/man/bayesfactor_restricted.Rd
@@ -7,18 +7,13 @@
 \alias{bayesfactor_restricted.brmsfit}
 \alias{bayesfactor_restricted.blavaan}
 \alias{bayesfactor_restricted.emmGrid}
+\alias{bayesfactor_restricted.data.frame}
 \alias{as.logical.bayesfactor_restricted}
 \title{Bayes Factors (BF) for Order Restricted Models}
 \usage{
-bayesfactor_restricted(
-  posterior,
-  hypothesis,
-  prior = NULL,
-  verbose = TRUE,
-  ...
-)
+bayesfactor_restricted(posterior, ...)
 
-bf_restricted(posterior, hypothesis, prior = NULL, verbose = TRUE, ...)
+bf_restricted(posterior, ...)
 
 \method{bayesfactor_restricted}{stanreg}(
   posterior,
@@ -56,20 +51,28 @@ bf_restricted(posterior, hypothesis, prior = NULL, verbose = TRUE, ...)
   ...
 )
 
+\method{bayesfactor_restricted}{data.frame}(
+  posterior,
+  hypothesis,
+  prior = NULL,
+  rvar_col = NULL,
+  ...
+)
+
 \method{as.logical}{bayesfactor_restricted}(x, which = c("posterior", "prior"), ...)
 }
 \arguments{
 \item{posterior}{A \code{stanreg} / \code{brmsfit} object, \code{emmGrid} or a data frame - representing
 a posterior distribution(s) from (see Details).}
 
+\item{...}{Currently not used.}
+
 \item{hypothesis}{A character vector specifying the restrictions as logical conditions (see examples below).}
 
 \item{prior}{An object representing a prior distribution (see Details).}
 
 \item{verbose}{Toggle off warnings.}
 
-\item{...}{Currently not used.}
-
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
@@ -77,6 +80,9 @@ returned? Only applies to mixed models. May be abbreviated.}
 conditional model or the zero-inflated part of the model be returned? May
 be abbreviated. Only applies to \pkg{brms}-models.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{x}{An object of class \code{bayesfactor_restricted}}
 
 \item{which}{Should the logical matrix be of the posterior or prior distribution(s)?}
@@ -119,29 +125,27 @@ prior.
 (Note that by default, \code{brms::brm()} uses flat priors for fixed-effects;
 See example below.)
 \cr\cr
-It is important to provide the correct \code{prior} for meaningful results.
+It is important to provide the correct \code{prior} for meaningful results,
+to match the \code{posterior}-type input:
 \itemize{
-\item When \code{posterior} is a numerical vector, \code{prior} should also be a numerical vector.
-\item When \code{posterior} is a \code{data.frame}, \code{prior} should also be a \code{data.frame}, with matching column order.
-\item When \code{posterior} is a \code{stanreg}, \code{brmsfit} or other supported Bayesian model:
+\item \strong{A numeric vector} - \code{prior} should also be a \emph{numeric vector}, representing the prior-estimate.
+\item \strong{A data frame} - \code{prior} should also be a \emph{data frame}, representing the prior-estimates, in matching column order.
 \itemize{
-\item \code{prior} can be set to \code{NULL}, in which case prior samples are drawn internally.
-\item \code{prior} can also be a model equivalent to \code{posterior} but with samples from
-the priors \emph{only}. See \code{\link[=unupdate]{unupdate()}}.
-\item \strong{Note:} When \code{posterior} is a \code{brmsfit_multiple} model, \code{prior} \strong{must} be provided.
-}
-\item When \code{posterior} is an output from a \code{{marginaleffects}} function, \code{prior} should also be an an output
-from a \code{{marginaleffects}} function equivalent to \code{posterior} but created
-with a model of priors samples \emph{only}.
-\item When \code{posterior} is an \code{emmGrid} / \code{emm_list} object:
+\item If \code{rvar_col} is specified, \code{prior} should be \emph{the name of an \code{rvar} column} that represents the prior-estimates.
+}
+\item \strong{Supported Bayesian model (\code{stanreg}, \code{brmsfit}, etc.)}
+\itemize{
+\item \code{prior} should be \emph{a model an equivalent model with MCMC samples from the priors \strong{only}}. See \code{\link[=unupdate]{unupdate()}}.
+\item If \code{prior} is set to \code{NULL}, \code{\link[=unupdate]{unupdate()}} is called internally (not supported for \code{brmsfit_multiple} model).
+}
+\item \strong{Output from a \code{{marginaleffects}} function} - \code{prior} should also be \emph{an equivalent output} from a \code{{marginaleffects}} function based on a prior-model
+(See \code{\link[=unupdate]{unupdate()}}).
+\item \strong{Output from an \code{{emmeans}} function}
 \itemize{
-\item \code{prior} should also be an \code{emmGrid} / \code{emm_list} object equivalent to \code{posterior} but
-created with a model of priors samples \emph{only}. See \code{\link[=unupdate]{unupdate()}}.
-\item \code{prior} can also be the original (posterior) \emph{model}. If so, the function will try to
-update the \code{emmGrid} / \code{emm_list} to use the \code{\link[=unupdate]{unupdate()}}d prior-model.
-(\emph{This cannot be done for \code{brmsfit} models.})
-\item \strong{Note}: When the \code{emmGrid} has undergone any transformations (\code{"log"}, \code{"response"}, etc.),
-or \code{regrid}ing, then \code{prior} must be an \code{emmGrid} object, as stated above.
+\item \code{prior} should also be \emph{an equivalent output} from an \code{{emmeans}} function based on a prior-model (See \code{\link[=unupdate]{unupdate()}}).
+\item \code{prior} can also be \emph{the original (posterior) model}, in which case the function
+will try to "unupdate" the estimates (not supported if the estimates have undergone
+any transformations -- \code{"log"}, \code{"response"}, etc. -- or any \code{regrid}ing).
 }
 }
 }
diff --git a/man/bci.Rd b/man/bci.Rd
index db4de886b..172d51655 100644
--- a/man/bci.Rd
+++ b/man/bci.Rd
@@ -22,7 +22,7 @@ bcai(x, ...)
 
 \method{bci}{numeric}(x, ci = 0.95, verbose = TRUE, ...)
 
-\method{bci}{data.frame}(x, ci = 0.95, verbose = TRUE, ...)
+\method{bci}{data.frame}(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...)
 
 \method{bci}{MCMCglmm}(x, ci = 0.95, verbose = TRUE, ...)
 
@@ -80,6 +80,9 @@ resemble the arguments of the \code{.numeric} or \code{.data.frame}methods.}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/ci.Rd b/man/ci.Rd
index a80f688d8..1ef0d7948 100644
--- a/man/ci.Rd
+++ b/man/ci.Rd
@@ -16,7 +16,7 @@ ci(x, ...)
 
 \method{ci}{numeric}(x, ci = 0.95, method = "ETI", verbose = TRUE, BF = 1, ...)
 
-\method{ci}{data.frame}(x, ci = 0.95, method = "ETI", verbose = TRUE, BF = 1, ...)
+\method{ci}{data.frame}(x, ci = 0.95, method = "ETI", BF = 1, rvar_col = NULL, verbose = TRUE, ...)
 
 \method{ci}{sim.merMod}(
   x,
@@ -75,6 +75,9 @@ to be estimated. Default to \code{0.95} (\verb{95\%}).}
 
 \item{BF}{The amount of support required to be included in the support interval.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/describe_posterior.Rd b/man/describe_posterior.Rd
index 1cfa19d72..bc3e0fcfe 100644
--- a/man/describe_posterior.Rd
+++ b/man/describe_posterior.Rd
@@ -3,6 +3,7 @@
 \name{describe_posterior}
 \alias{describe_posterior}
 \alias{describe_posterior.numeric}
+\alias{describe_posterior.data.frame}
 \alias{describe_posterior.stanreg}
 \alias{describe_posterior.brmsfit}
 \title{Describe Posterior Distributions}
@@ -25,6 +26,23 @@ describe_posterior(posterior, ...)
   ...
 )
 
+\method{describe_posterior}{data.frame}(
+  posterior,
+  centrality = "median",
+  dispersion = FALSE,
+  ci = 0.95,
+  ci_method = "eti",
+  test = c("p_direction", "rope"),
+  rope_range = "default",
+  rope_ci = 0.95,
+  keep_iterations = FALSE,
+  bf_prior = NULL,
+  BF = 1,
+  rvar_col = NULL,
+  verbose = TRUE,
+  ...
+)
+
 \method{describe_posterior}{stanreg}(
   posterior,
   centrality = "median",
@@ -123,6 +141,9 @@ case of models) ignored.}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{diagnostic}{Diagnostic metrics to compute.  Character (vector) or list
 with one or more of these options: \code{"ESS"}, \code{"Rhat"}, \code{"MCSE"} or \code{"all"}.}
 
diff --git a/man/equivalence_test.Rd b/man/equivalence_test.Rd
index 8f7b44326..accbbcaae 100644
--- a/man/equivalence_test.Rd
+++ b/man/equivalence_test.Rd
@@ -12,7 +12,14 @@ equivalence_test(x, ...)
 
 \method{equivalence_test}{default}(x, ...)
 
-\method{equivalence_test}{data.frame}(x, range = "default", ci = 0.95, verbose = TRUE, ...)
+\method{equivalence_test}{data.frame}(
+  x,
+  range = "default",
+  ci = 0.95,
+  rvar_col = NULL,
+  verbose = TRUE,
+  ...
+)
 
 \method{equivalence_test}{stanreg}(
   x,
@@ -55,6 +62,9 @@ model, \code{\link[=rope_range]{rope_range()}} is used.}
 \item{ci}{The Credible Interval (CI) probability, corresponding to the
 proportion of HDI, to use for the percentage in ROPE.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{verbose}{Toggle off warnings.}
 
 \item{effects}{Should results for fixed effects, random effects or both be
diff --git a/man/estimate_density.Rd b/man/estimate_density.Rd
index 6132abd14..5d5ac5c0c 100644
--- a/man/estimate_density.Rd
+++ b/man/estimate_density.Rd
@@ -18,6 +18,7 @@ estimate_density(x, ...)
   select = NULL,
   by = NULL,
   at = NULL,
+  rvar_col = NULL,
   ...
 )
 }
@@ -55,6 +56,9 @@ density estimation is performed for each group (subsets) indicated by \code{by}.
 See examples.}
 
 \item{at}{Deprecated in favour of \code{by}.}
+
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
 }
 \description{
 This function is a wrapper over different methods of density estimation. By
diff --git a/man/eti.Rd b/man/eti.Rd
index 18c7c348d..4b116beb7 100644
--- a/man/eti.Rd
+++ b/man/eti.Rd
@@ -3,6 +3,7 @@
 \name{eti}
 \alias{eti}
 \alias{eti.numeric}
+\alias{eti.data.frame}
 \alias{eti.stanreg}
 \alias{eti.brmsfit}
 \alias{eti.get_predicted}
@@ -12,6 +13,8 @@ eti(x, ...)
 
 \method{eti}{numeric}(x, ci = 0.95, verbose = TRUE, ...)
 
+\method{eti}{data.frame}(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...)
+
 \method{eti}{stanreg}(
   x,
   ci = 0.95,
@@ -49,6 +52,9 @@ resemble the arguments of the \code{.numeric} or \code{.data.frame}methods.}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/hdi.Rd b/man/hdi.Rd
index bb0003cba..927d4f3da 100644
--- a/man/hdi.Rd
+++ b/man/hdi.Rd
@@ -13,7 +13,7 @@ hdi(x, ...)
 
 \method{hdi}{numeric}(x, ci = 0.95, verbose = TRUE, ...)
 
-\method{hdi}{data.frame}(x, ci = 0.95, verbose = TRUE, ...)
+\method{hdi}{data.frame}(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...)
 
 \method{hdi}{stanreg}(
   x,
@@ -52,6 +52,9 @@ resemble the arguments of the \code{.numeric} or \code{.data.frame}methods.}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/map_estimate.Rd b/man/map_estimate.Rd
index 3ffa10290..fa1dcd338 100644
--- a/man/map_estimate.Rd
+++ b/man/map_estimate.Rd
@@ -34,7 +34,7 @@ map_estimate(x, ...)
   ...
 )
 
-\method{map_estimate}{data.frame}(x, precision = 2^10, method = "kernel", ...)
+\method{map_estimate}{data.frame}(x, precision = 2^10, method = "kernel", rvar_col = NULL, ...)
 
 \method{map_estimate}{get_predicted}(
   x,
@@ -72,6 +72,9 @@ filtered by default, so only parameters that typically appear in the
 \code{summary()} are returned. Use \code{parameters} to select specific parameters
 for the output.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{use_iterations}{Logical, if \code{TRUE} and \code{x} is a \code{get_predicted} object,
 (returned by \code{\link[insight:get_predicted]{insight::get_predicted()}}), the function is applied to the
 iterations instead of the predictions. This only applies to models that return
diff --git a/man/p_direction.Rd b/man/p_direction.Rd
index f9a2530e7..bd616918a 100644
--- a/man/p_direction.Rd
+++ b/man/p_direction.Rd
@@ -127,8 +127,8 @@ frequentist p-value using \code{\link[=pd_to_p]{pd_to_p()}}.}
 \item{remove_na}{Should missing values be removed before computation? Note
 that \code{Inf} (infinity) are \emph{not} removed.}
 
-\item{rvar_col}{Name of an \code{rvar}-type column. If \code{NULL}, each column in the
-data frame is assumed to represent draws from a posterior distribution.}
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
 
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
@@ -263,6 +263,7 @@ p_direction(posterior, method = "kernel")
 df <- data.frame(replicate(4, rnorm(100)))
 p_direction(df)
 p_direction(df, method = "kernel")
+
 \donttest{
 # rstanarm models
 # -----------------------------------------------
@@ -290,6 +291,13 @@ p_direction(bf)
 p_direction(bf, method = "kernel")
 }
 \dontshow{\}) # examplesIf}
+\dontshow{if (requireNamespace("posterior", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+# Using "rvar_col"
+x <- data.frame(mu = c(0, 0.5, 1), sigma = c(1, 0.5, 0.25))
+x$my_rvar <- posterior::rvar_rng(rnorm, 3, mean = x$mu, sd = x$sigma)
+x
+p_direction(x, rvar_col = "my_rvar")
+\dontshow{\}) # examplesIf}
 }
 \references{
 \itemize{
diff --git a/man/p_map.Rd b/man/p_map.Rd
index 61f880de2..f353c4c5b 100644
--- a/man/p_map.Rd
+++ b/man/p_map.Rd
@@ -5,6 +5,7 @@
 \alias{p_pointnull}
 \alias{p_map.numeric}
 \alias{p_map.get_predicted}
+\alias{p_map.data.frame}
 \alias{p_map.stanreg}
 \alias{p_map.brmsfit}
 \title{Bayesian p-value based on the density at the Maximum A Posteriori (MAP)}
@@ -25,6 +26,8 @@ p_pointnull(x, ...)
   ...
 )
 
+\method{p_map}{data.frame}(x, null = 0, precision = 2^10, method = "kernel", rvar_col = NULL, ...)
+
 \method{p_map}{stanreg}(
   x,
   null = 0,
@@ -72,6 +75,9 @@ iterations for predicted values (e.g., \code{brmsfit} models).}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/p_rope.Rd b/man/p_rope.Rd
index 52845be6e..8791c0224 100644
--- a/man/p_rope.Rd
+++ b/man/p_rope.Rd
@@ -3,6 +3,7 @@
 \name{p_rope}
 \alias{p_rope}
 \alias{p_rope.numeric}
+\alias{p_rope.data.frame}
 \alias{p_rope.stanreg}
 \alias{p_rope.brmsfit}
 \title{Probability of being in the ROPE}
@@ -11,6 +12,8 @@ p_rope(x, ...)
 
 \method{p_rope}{numeric}(x, range = "default", verbose = TRUE, ...)
 
+\method{p_rope}{data.frame}(x, range = "default", rvar_col = NULL, verbose = TRUE, ...)
+
 \method{p_rope}{stanreg}(
   x,
   range = "default",
@@ -18,7 +21,7 @@ p_rope(x, ...)
   component = c("location", "all", "conditional", "smooth_terms", "sigma",
     "distributional", "auxiliary"),
   parameters = NULL,
-  verbose = verbose,
+  verbose = TRUE,
   ...
 )
 
@@ -49,6 +52,9 @@ model, \code{\link[=rope_range]{rope_range()}} is used.}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/p_significance.Rd b/man/p_significance.Rd
index 4d4c944a5..30b2170e5 100644
--- a/man/p_significance.Rd
+++ b/man/p_significance.Rd
@@ -4,6 +4,7 @@
 \alias{p_significance}
 \alias{p_significance.numeric}
 \alias{p_significance.get_predicted}
+\alias{p_significance.data.frame}
 \alias{p_significance.stanreg}
 \alias{p_significance.brmsfit}
 \title{Practical Significance (ps)}
@@ -20,6 +21,8 @@ p_significance(x, ...)
   ...
 )
 
+\method{p_significance}{data.frame}(x, threshold = "default", rvar_col = NULL, ...)
+
 \method{p_significance}{stanreg}(
   x,
   threshold = "default",
@@ -66,6 +69,9 @@ iterations for predicted values (e.g., \code{brmsfit} models).}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/point_estimate.Rd b/man/point_estimate.Rd
index 210100fb8..7368ef0dc 100644
--- a/man/point_estimate.Rd
+++ b/man/point_estimate.Rd
@@ -3,6 +3,7 @@
 \name{point_estimate}
 \alias{point_estimate}
 \alias{point_estimate.numeric}
+\alias{point_estimate.data.frame}
 \alias{point_estimate.stanreg}
 \alias{point_estimate.brmsfit}
 \alias{point_estimate.BFBayesFactor}
@@ -13,6 +14,15 @@ point_estimate(x, ...)
 
 \method{point_estimate}{numeric}(x, centrality = "all", dispersion = FALSE, threshold = 0.1, ...)
 
+\method{point_estimate}{data.frame}(
+  x,
+  centrality = "all",
+  dispersion = FALSE,
+  threshold = 0.1,
+  rvar_col = NULL,
+  ...
+)
+
 \method{point_estimate}{stanreg}(
   x,
   centrality = "all",
@@ -67,6 +77,9 @@ Dispersion is not available for \code{"MAP"} or \code{"mode"} centrality indices
 the fraction (0 to 0.5) of observations to be trimmed from each end of the
 vector before the mean is computed.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/rope.Rd b/man/rope.Rd
index b4a9e1255..87036d460 100644
--- a/man/rope.Rd
+++ b/man/rope.Rd
@@ -3,6 +3,7 @@
 \name{rope}
 \alias{rope}
 \alias{rope.numeric}
+\alias{rope.data.frame}
 \alias{rope.stanreg}
 \alias{rope.brmsfit}
 \title{Region of Practical Equivalence (ROPE)}
@@ -11,6 +12,16 @@ rope(x, ...)
 
 \method{rope}{numeric}(x, range = "default", ci = 0.95, ci_method = "ETI", verbose = TRUE, ...)
 
+\method{rope}{data.frame}(
+  x,
+  range = "default",
+  ci = 0.95,
+  ci_method = "ETI",
+  rvar_col = NULL,
+  verbose = TRUE,
+  ...
+)
+
 \method{rope}{stanreg}(
   x,
   range = "default",
@@ -59,6 +70,9 @@ ROPE. Can be 'HDI' (default) or 'ETI'. See \code{\link[=ci]{ci()}}.}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/man/si.Rd b/man/si.Rd
index 45e87e7c3..0ca1492af 100644
--- a/man/si.Rd
+++ b/man/si.Rd
@@ -62,7 +62,7 @@ si(posterior, ...)
   ...
 )
 
-\method{si}{data.frame}(posterior, prior = NULL, BF = 1, verbose = TRUE, ...)
+\method{si}{data.frame}(posterior, prior = NULL, BF = 1, rvar_col = NULL, verbose = TRUE, ...)
 }
 \arguments{
 \item{posterior}{A numerical vector, \code{stanreg} / \code{brmsfit} object,
@@ -95,6 +95,9 @@ for the output.}
 (returned by \code{\link[insight:get_predicted]{insight::get_predicted()}}), the function is applied to the
 iterations instead of the predictions. This only applies to models that return
 iterations for predicted values (e.g., \code{brmsfit} models).}
+
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
 }
 \value{
 A data frame containing the lower and upper bounds of the SI.
@@ -150,29 +153,27 @@ prior.
 (Note that by default, \code{brms::brm()} uses flat priors for fixed-effects;
 See example below.)
 \cr\cr
-It is important to provide the correct \code{prior} for meaningful results.
+It is important to provide the correct \code{prior} for meaningful results,
+to match the \code{posterior}-type input:
+\itemize{
+\item \strong{A numeric vector} - \code{prior} should also be a \emph{numeric vector}, representing the prior-estimate.
+\item \strong{A data frame} - \code{prior} should also be a \emph{data frame}, representing the prior-estimates, in matching column order.
 \itemize{
-\item When \code{posterior} is a numerical vector, \code{prior} should also be a numerical vector.
-\item When \code{posterior} is a \code{data.frame}, \code{prior} should also be a \code{data.frame}, with matching column order.
-\item When \code{posterior} is a \code{stanreg}, \code{brmsfit} or other supported Bayesian model:
+\item If \code{rvar_col} is specified, \code{prior} should be \emph{the name of an \code{rvar} column} that represents the prior-estimates.
+}
+\item \strong{Supported Bayesian model (\code{stanreg}, \code{brmsfit}, etc.)}
 \itemize{
-\item \code{prior} can be set to \code{NULL}, in which case prior samples are drawn internally.
-\item \code{prior} can also be a model equivalent to \code{posterior} but with samples from
-the priors \emph{only}. See \code{\link[=unupdate]{unupdate()}}.
-\item \strong{Note:} When \code{posterior} is a \code{brmsfit_multiple} model, \code{prior} \strong{must} be provided.
-}
-\item When \code{posterior} is an output from a \code{{marginaleffects}} function, \code{prior} should also be an an output
-from a \code{{marginaleffects}} function equivalent to \code{posterior} but created
-with a model of priors samples \emph{only}.
-\item When \code{posterior} is an \code{emmGrid} / \code{emm_list} object:
+\item \code{prior} should be \emph{a model an equivalent model with MCMC samples from the priors \strong{only}}. See \code{\link[=unupdate]{unupdate()}}.
+\item If \code{prior} is set to \code{NULL}, \code{\link[=unupdate]{unupdate()}} is called internally (not supported for \code{brmsfit_multiple} model).
+}
+\item \strong{Output from a \code{{marginaleffects}} function} - \code{prior} should also be \emph{an equivalent output} from a \code{{marginaleffects}} function based on a prior-model
+(See \code{\link[=unupdate]{unupdate()}}).
+\item \strong{Output from an \code{{emmeans}} function}
 \itemize{
-\item \code{prior} should also be an \code{emmGrid} / \code{emm_list} object equivalent to \code{posterior} but
-created with a model of priors samples \emph{only}. See \code{\link[=unupdate]{unupdate()}}.
-\item \code{prior} can also be the original (posterior) \emph{model}. If so, the function will try to
-update the \code{emmGrid} / \code{emm_list} to use the \code{\link[=unupdate]{unupdate()}}d prior-model.
-(\emph{This cannot be done for \code{brmsfit} models.})
-\item \strong{Note}: When the \code{emmGrid} has undergone any transformations (\code{"log"}, \code{"response"}, etc.),
-or \code{regrid}ing, then \code{prior} must be an \code{emmGrid} object, as stated above.
+\item \code{prior} should also be \emph{an equivalent output} from an \code{{emmeans}} function based on a prior-model (See \code{\link[=unupdate]{unupdate()}}).
+\item \code{prior} can also be \emph{the original (posterior) model}, in which case the function
+will try to "unupdate" the estimates (not supported if the estimates have undergone
+any transformations -- \code{"log"}, \code{"response"}, etc. -- or any \code{regrid}ing).
 }
 }
 }
diff --git a/man/spi.Rd b/man/spi.Rd
index 995de8616..6212d9343 100644
--- a/man/spi.Rd
+++ b/man/spi.Rd
@@ -3,6 +3,7 @@
 \name{spi}
 \alias{spi}
 \alias{spi.numeric}
+\alias{spi.data.frame}
 \alias{spi.stanreg}
 \alias{spi.brmsfit}
 \alias{spi.get_predicted}
@@ -12,6 +13,8 @@ spi(x, ...)
 
 \method{spi}{numeric}(x, ci = 0.95, verbose = TRUE, ...)
 
+\method{spi}{data.frame}(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...)
+
 \method{spi}{stanreg}(
   x,
   ci = 0.95,
@@ -49,6 +52,9 @@ resemble the arguments of the \code{.numeric} or \code{.data.frame}methods.}
 
 \item{verbose}{Toggle off warnings.}
 
+\item{rvar_col}{A single character - the name of an \code{rvar} column in the data
+frame to be processed. See example in \code{\link[=p_direction]{p_direction()}}.}
+
 \item{effects}{Should results for fixed effects, random effects or both be
 returned? Only applies to mixed models. May be abbreviated.}
 
diff --git a/tests/testthat/test-data.frame-with-rvar.R b/tests/testthat/test-data.frame-with-rvar.R
new file mode 100644
index 000000000..220c534bd
--- /dev/null
+++ b/tests/testthat/test-data.frame-with-rvar.R
@@ -0,0 +1,119 @@
+test_that("data.frame w/ rvar_col descrive_posterior etc", {
+  skip_on_ci()
+  skip_on_cran()
+  skip_if_not_installed("posterior")
+
+  dfx <- data.frame(mu = c(0, 0.5, 1), sigma = c(1, 0.5, 0.25))
+  dfx$my_rvar <- posterior::rvar_rng(rnorm, 3, mean = dfx$mu, sd = dfx$sigma)
+  dfx$other_rvar <- posterior::rvar_rng(rnorm, 3, mean = dfx$mu + 0.5, sd = dfx$sigma - 0.1)
+  dfx
+
+  ## Errors
+  expect_error(p_direction(dfx, rvar_col = "mu"))
+  expect_error(p_direction(dfx, rvar_col = "my_rvarrrrrr"))
+
+
+  ## describe_posterior
+  res <- describe_posterior(dfx,
+    rvar_col = "my_rvar",
+    centrality = "MAP", ci_method = "hdi", ci = 0.8,
+    test = c("pd", "p_map", "rope", "equivalence_test"),
+    rope_ci = 1, rope_range = c(-1, 0.5)
+  )
+  res.ref <- describe_posterior(dfx$my_rvar,
+    centrality = "MAP", ci_method = "hdi", ci = 0.8,
+    test = c("pd", "p_map", "rope", "equivalence_test"),
+    rope_ci = 1, rope_range = c(-1, 0.5)
+  )
+  expect_true(all(c("mu", "sigma") %in% colnames(res)))
+  expect_equal(res[setdiff(colnames(res), c("mu", "sigma"))],
+    res.ref[setdiff(colnames(res.ref), "Parameter")],
+    ignore_attr = TRUE
+  )
+
+  ## CIs
+  res <- eti(dfx, rvar_col = "my_rvar")
+  res.ref <- eti(dfx$my_rvar)
+  expect_true(all(c("mu", "sigma") %in% colnames(res)))
+  expect_identical(nrow(format(res)), 3L)
+  expect_identical(ncol(format(res)), 3L)
+  expect_equal(res[setdiff(colnames(res), c("mu", "sigma"))],
+    res.ref[setdiff(colnames(res.ref), "Parameter")],
+    ignore_attr = TRUE
+  )
+
+  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(nrow(format(res)), 3L)
+  expect_identical(ncol(format(res)), 4L)
+  expect_equal(res[setdiff(colnames(res), c("mu", "sigma"))],
+    res.ref[setdiff(colnames(res.ref), "Parameter")],
+    ignore_attr = TRUE
+  )
+
+  ## estimate_density
+  res <- estimate_density(dfx, rvar_col = "my_rvar")
+  res.ref <- estimate_density(dfx$my_rvar)
+  expect_true(all(c("mu", "sigma") %in% colnames(res)))
+  expect_equal(res[setdiff(colnames(res), c("mu", "sigma"))],
+    res.ref[setdiff(colnames(res.ref), "Parameter")],
+    ignore_attr = TRUE
+  )
+})
+
+test_that("data.frame w/ rvar_col bayesfactors", {
+  skip_on_ci()
+  skip_on_cran()
+  skip_if_not_installed("posterior")
+  skip_if_not_installed("logspline")
+
+  dfx <- data.frame(mu = c(0, 0.5, 1), sigma = c(1, 0.5, 0.25))
+  dfx$my_rvar <- posterior::rvar_rng(rnorm, 3, mean = dfx$mu, sd = dfx$sigma)
+  dfx$other_rvar <- posterior::rvar_rng(rnorm, 3, mean = dfx$mu + 0.5, sd = dfx$sigma - 0.1)
+  dfx
+
+
+
+  ## SIs
+  res <- si(dfx, rvar_col = "my_rvar", prior = "other_rvar", verbose = FALSE)
+  res.ref <- si(dfx$my_rvar, prior = dfx$other_rvar, verbose = FALSE)
+  expect_true(all(c("mu", "sigma") %in% colnames(res)))
+  expect_identical(nrow(format(res)), 3L)
+  expect_identical(ncol(format(res)), 3L)
+  expect_equal(res[setdiff(colnames(res), c("mu", "sigma"))],
+    res.ref[setdiff(colnames(res.ref), "Parameter")],
+    ignore_attr = TRUE
+  )
+
+  res <- si(dfx,
+    rvar_col = "my_rvar", prior = "other_rvar",
+    BF = c(1, 3), verbose = FALSE
+  )
+  res.ref <- si(dfx$my_rvar,
+    prior = dfx$other_rvar,
+    BF = c(1, 3), verbose = FALSE
+  )
+  expect_true(all(c("mu", "sigma") %in% colnames(res)))
+  expect_identical(nrow(format(res)), 3L)
+  expect_identical(ncol(format(res)), 4L)
+  expect_equal(format(res[setdiff(colnames(res), c("mu", "sigma"))]),
+    format(res.ref[setdiff(colnames(res.ref), "Parameter")]),
+    ignore_attr = TRUE
+  )
+
+
+  ## bayesfactor_parameters
+  res <- bayesfactor_parameters(dfx,
+    rvar_col = "my_rvar", prior = "other_rvar",
+    verbose = FALSE
+  )
+  res.ref <- bayesfactor_parameters(dfx$my_rvar,
+    prior = dfx$other_rvar,
+    verbose = FALSE
+  )
+  expect_equal(res[setdiff(colnames(res), c("mu", "sigma"))],
+    res.ref[setdiff(colnames(res.ref), "Parameter")],
+    ignore_attr = TRUE
+  )
+})
diff --git a/tests/testthat/test-marginaleffects.R b/tests/testthat/test-marginaleffects.R
index 11be99559..a1a3dcb3f 100644
--- a/tests/testthat/test-marginaleffects.R
+++ b/tests/testthat/test-marginaleffects.R
@@ -1,4 +1,4 @@
-test_that("emmGrid descrive_posterior", {
+test_that("marginaleffects descrive_posterior", {
   skip_on_ci()
   skip_on_cran()
 
@@ -30,7 +30,7 @@ test_that("emmGrid descrive_posterior", {
   # estimate_density
   mfx <- marginaleffects::comparisons(mod,
     variables = "cyl",
-    newdata = data.frame(hp = 100, am = 0)
+    newdata = marginaleffects::datagrid(hp = 100, am = 0)
   )
   samps <- insight::get_parameters(mod)[c("cyl6", "cyl8")]
 
@@ -42,7 +42,7 @@ test_that("emmGrid descrive_posterior", {
   )
 })
 
-test_that("emmGrid bayesfactors", {
+test_that("marginaleffects bayesfactors", {
   skip_on_ci()
   skip_on_cran()
 
@@ -52,7 +52,7 @@ test_that("emmGrid bayesfactors", {
   data("mtcars")
   mtcars$cyl <- factor(mtcars$cyl)
   mod <- rstanarm::stan_glm(mpg ~ cyl + hp * am, data = mtcars, refresh = 0)
-  modp <- unupdate(mod)
+  modp <- unupdate(mod, verbose = FALSE)
 
   mfx <- marginaleffects::avg_slopes(mod, by = "am")
   mfxp <- marginaleffects::avg_slopes(modp, by = "am")