diff --git a/DESCRIPTION b/DESCRIPTION
index 60309a335..6c5d28f34 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -113,6 +113,7 @@ Suggests:
     coxme,
     cplm,
     dbscan,
+    distributional,
     domir (>= 0.2.0),
     drc,
     DRR,
@@ -221,4 +222,4 @@ Config/testthat/edition: 3
 Config/testthat/parallel: true
 Config/Needs/website: easystats/easystatstemplate
 Config/rcmdcheck/ignore-inconsequential-notes: true
-Remotes: easystats/insight
+Remotes: easystats/insight, easystats/bayestestR
diff --git a/R/equivalence_test.R b/R/equivalence_test.R
index 86371a0dc..3d2b44ea5 100644
--- a/R/equivalence_test.R
+++ b/R/equivalence_test.R
@@ -10,15 +10,14 @@ bayestestR::equivalence_test
 #'
 #' @param x A statistical model.
 #' @param range The range of practical equivalence of an effect. May be
-#'   `"default"`, to automatically define this range based on properties of the
-#'   model's data.
+#' `"default"`, to automatically define this range based on properties of the
+#' model's data.
 #' @param ci Confidence Interval (CI) level. Default to `0.95` (`95%`).
 #' @param rule Character, indicating the rules when testing for practical
-#'   equivalence. Can be `"bayes"`, `"classic"` or `"cet"`. See
-#'   'Details'.
+#' equivalence. Can be `"bayes"`, `"classic"` or `"cet"`. See 'Details'.
 #' @param test Hypothesis test for computing contrasts or pairwise comparisons.
-#'   See [`?ggeffects::test_predictions`](https://strengejacke.github.io/ggeffects/reference/test_predictions.html)
-#'   for details.
+#' See [`?ggeffects::test_predictions`](https://strengejacke.github.io/ggeffects/reference/test_predictions.html)
+#' for details.
 #' @param verbose Toggle warnings and messages.
 #' @param ... Arguments passed to or from other methods.
 #' @inheritParams model_parameters.merMod
@@ -28,18 +27,17 @@ bayestestR::equivalence_test
 #' readings can be found in the references. See also [`p_significance()`] for
 #' a unidirectional equivalence test.
 #'
-#' @details
-#' In classical null hypothesis significance testing (NHST) within a frequentist
-#' framework, it is not possible to accept the null hypothesis, H0 - unlike
-#' in Bayesian statistics, where such probability statements are possible.
-#' "[...] one can only reject the null hypothesis if the test
+#' @details In classical null hypothesis significance testing (NHST) within a
+#' frequentist framework, it is not possible to accept the null hypothesis, H0 -
+#' unlike in Bayesian statistics, where such probability statements are
+#' possible. "[...] one can only reject the null hypothesis if the test
 #' statistics falls into the critical region(s), or fail to reject this
 #' hypothesis. In the latter case, all we can say is that no significant effect
 #' was observed, but one cannot conclude that the null hypothesis is true."
-#' (_Pernet 2017_). One way to address this issues without Bayesian methods
-#' is *Equivalence Testing*, as implemented in `equivalence_test()`.
-#' While you either can reject the null hypothesis or claim an inconclusive result
-#' in NHST, the equivalence test - according to _Pernet_ - adds a third category,
+#' (_Pernet 2017_). One way to address this issues without Bayesian methods is
+#' *Equivalence Testing*, as implemented in `equivalence_test()`. While you
+#' either can reject the null hypothesis or claim an inconclusive result in
+#' NHST, the equivalence test - according to _Pernet_ - adds a third category,
 #' *"accept"*. Roughly speaking, the idea behind equivalence testing in a
 #' frequentist framework is to check whether an estimate and its uncertainty
 #' (i.e. confidence interval) falls within a region of "practical equivalence".
@@ -95,25 +93,25 @@ bayestestR::equivalence_test
 #' ## p-Values
 #' The equivalence p-value is the area of the (cumulative) confidence
 #' distribution that is outside of the region of equivalence. It can be
-#' interpreted as p-value for *rejecting* the alternative hypothesis
-#' and *accepting* the "null hypothesis" (i.e. assuming practical
-#' equivalence). That is, a high p-value means we reject the assumption of
-#' practical equivalence and accept the alternative hypothesis.
+#' interpreted as p-value for *rejecting* the alternative hypothesis and
+#' *accepting* the "null hypothesis" (i.e. assuming practical equivalence). That
+#' is, a high p-value means we reject the assumption of practical equivalence
+#' and accept the alternative hypothesis.
 #'
 #' ## Second Generation p-Value (SGPV)
 #' Second generation p-values (SGPV) were proposed as a statistic that
 #' represents _the proportion of data-supported hypotheses that are also null
 #' hypotheses_ _(Blume et al. 2018, Lakens and Delacre 2020)_. It represents the
-#' proportion of the _full_ confidence interval range (assuming a normally
-#' distributed, equal-tailed interval) that is inside the ROPE. The SGPV ranges
-#' from zero to one. Higher values indicate that the effect is more likely to be
-#' practically equivalent ("not of interest").
+#' proportion of the _full_ confidence interval range (assuming a normally or
+#' t-distributed, equal-tailed interval, based on the model) that is inside the
+#' ROPE. The SGPV ranges from zero to one. Higher values indicate that the
+#' effect is more likely to be practically equivalent ("not of interest").
 #'
 #' Note that the assumed interval, which is used to calculate the SGPV, is an
 #' estimation of the _full interval_ based on the chosen confidence level. For
 #' example, if the 95% confidence interval of a coefficient ranges from -1 to 1,
-#' the underlying _full (normally distributed) interval_ approximately ranges
-#' from -1.9 to 1.9, see also following code:
+#' the underlying _full (normally or t-distributed) interval_ approximately
+#' ranges from -1.9 to 1.9, see also following code:
 #'
 #' ```
 #' # simulate full normal distribution
@@ -390,6 +388,7 @@ equivalence_test.ggeffects <- function(x,
   focal <- attributes(x)$original.terms
   obj_name <- attributes(x)$model.name
   ci <- attributes(x)$ci.lvl
+  dof <- attributes(x)$df
 
   x <- .get_ggeffects_model(x)
 
@@ -419,6 +418,7 @@ equivalence_test.ggeffects <- function(x,
         ci_narrow,
         range_rope = range,
         rule = rule,
+        dof = dof,
         verbose = verbose
       )
     }, conf_int, conf_int2
@@ -490,6 +490,18 @@ equivalence_test.ggeffects <- function(x,
   }
 
 
+  # ==== check degrees of freedom ====
+
+  df_column <- grep("(df|df_error)", colnames(x))
+  if (length(df_column) > 0) {
+    dof <- unique(x[[df_column]])
+    if (length(dof) > 1) {
+      dof <- Inf
+    }
+  } else {
+    dof <- Inf
+  }
+
   # ==== requested confidence intervals ====
 
   params <- conf_int <- .ci_generic(x, ci = ci)
@@ -513,6 +525,7 @@ equivalence_test.ggeffects <- function(x,
         ci_narrow,
         range_rope = range,
         rule = rule,
+        dof = dof,
         verbose = verbose
       )
     }, conf_int, conf_int2
@@ -623,6 +636,7 @@ equivalence_test.ggeffects <- function(x,
                                       ci_narrow,
                                       range_rope,
                                       rule,
+                                      dof = Inf,
                                       verbose) {
   final_ci <- NULL
 
@@ -674,7 +688,7 @@ equivalence_test.ggeffects <- function(x,
   data.frame(
     CI_low = final_ci[1],
     CI_high = final_ci[2],
-    SGPV = .rope_coverage(ci = ci, range_rope, ci_range = final_ci),
+    SGPV = .rope_coverage(ci = ci, range_rope, ci_range = final_ci, dof = dof),
     ROPE_low = range_rope[1],
     ROPE_high = range_rope[2],
     ROPE_Equivalence = decision,
@@ -726,8 +740,8 @@ equivalence_test.ggeffects <- function(x,
 # same range / limits as the confidence interval, thus indeed representing a
 # normally distributed confidence interval. We then calculate the probability
 # mass of this interval that is inside the ROPE.
-.rope_coverage <- function(ci = 0.95, range_rope, ci_range) {
-  out <- .generate_posterior_from_ci(ci, ci_range)
+.rope_coverage <- function(ci = 0.95, range_rope, ci_range, dof = Inf) {
+  out <- .generate_posterior_from_ci(ci, ci_range, dof = dof)
   # compare: ci_range and range(out)
   # The SGPV refers to the proportion of the confidence interval inside the
   # full ROPE - thus, we set ci = 1 here
@@ -736,23 +750,39 @@ equivalence_test.ggeffects <- function(x,
 }
 
 
-.generate_posterior_from_ci <- function(ci = 0.95, ci_range, precision = 10000) {
+.generate_posterior_from_ci <- function(ci = 0.95, ci_range, dof = Inf, precision = 10000) {
   # this function creates an approximate normal distribution that covers the
   # CI-range, i.e. we "simulate" a posterior distribution from a frequentist CI
 
-  # first we need the z-values of thq quantiles from the CI range
-  z_value <- stats::qnorm((1 + ci) / 2)
-  # then we need the range of the CI (in units), also to calculate the mean of
+  # sanity check - dof argument
+  if (is.null(dof)) {
+    dof <- Inf
+  }
+  # first we need the range of the CI (in units), also to calculate the mean of
   # the normal distribution
   diff_ci <- abs(diff(ci_range))
   mean_dist <- ci_range[2] - (diff_ci / 2)
+  # then we need the critical values of the quantiles from the CI range
+  z_value <- stats::qt((1 + ci) / 2, df = dof)
   # the range of Z-scores (from lower to upper quantile) gives us the range of
   # the provided interval in terms of standard deviations. now we divide the
   # known range of the provided CI (in units) by the z-score-range, which will
   # give us the standard deviation of the distribution.
   sd_dist <- diff_ci / diff(c(-1 * z_value, z_value))
-  # we now know all parameters (mean and sd) to simulate a normal distribution
-  bayestestR::distribution_normal(n = precision, mean = mean_dist, sd = sd_dist)
+  # generate normal-distribution if we don't have t-distribution, or if
+  # we don't have necessary packages installed
+  if (is.infinite(dof) || !insight::check_if_installed("distributional", quietly = TRUE)) {
+    # tell user to install "distributional"
+    if (!is.infinite(dof)) {
+      insight::format_alert("For models with only few degrees of freedom, install the {distributional} package to increase accuracy of `p_direction()`, `p_significance()` and `equivalence_test()`.") # nolint
+    }
+    # we now know all parameters (mean and sd) to simulate a normal distribution
+    bayestestR::distribution_normal(n = precision, mean = mean_dist, sd = sd_dist)
+  } else {
+    insight::check_if_installed("distributional")
+    out <- distributional::dist_student_t(df = dof, mu = mean_dist, sigma = sd_dist)
+    sort(unlist(distributional::generate(out, times = precision), use.names = FALSE))
+  }
 }
 
 
diff --git a/R/p_significance.R b/R/p_significance.R
index 17124d0fe..2629d793f 100644
--- a/R/p_significance.R
+++ b/R/p_significance.R
@@ -35,21 +35,22 @@ bayestestR::p_significance
 #' for functions related to checking effect existence and significance.
 #'
 #' @details `p_significance()` returns the proportion of the _full_ confidence
-#' interval range (assuming a normally distributed, equal-tailed interval) that
-#' is outside a certain range (the negligible effect, or ROPE, see argument
-#' `threshold`). If there are values of the distribution both below and above
-#' the ROPE, `p_significance()` returns the higher probability of a value being
-#' outside the ROPE. Typically, this value should be larger than 0.5 to indicate
-#' practical significance. However, if the range of the negligible effect is
-#' rather large compared to the range of the confidence interval,
-#' `p_significance()` will be less than 0.5, which indicates no clear practical
-#' significance.
+#' interval range (assuming a normally or t-distributed, equal-tailed interval,
+#' based on the model) that is outside a certain range (the negligible effect,
+#' or ROPE, see argument `threshold`). If there are values of the distribution
+#' both below and above the ROPE, `p_significance()` returns the higher
+#' probability of a value being outside the ROPE. Typically, this value should
+#' be larger than 0.5 to indicate practical significance. However, if the range
+#' of the negligible effect is rather large compared to the range of the
+#' confidence interval, `p_significance()` will be less than 0.5, which
+#' indicates no clear practical significance.
 #'
 #' Note that the assumed interval, which is used to calculate the practical
 #' significance, is an estimation of the _full interval_ based on the chosen
 #' confidence level. For example, if the 95% confidence interval of a
-#' coefficient ranges from -1 to 1, the underlying _full (normally distributed)
-#' interval_ approximately ranges from -1.9 to 1.9, see also following code:
+#' coefficient ranges from -1 to 1, the underlying _full (normally or
+#' t-distributed) interval_ approximately ranges from -1.9 to 1.9, see also
+#' following code:
 #'
 #' ```
 #' # simulate full normal distribution
@@ -178,11 +179,12 @@ p_significance.lm <- function(x, threshold = "default", ci = 0.95, verbose = TRU
 .posterior_ci <- function(x, ci, ...) {
   # first, we need CIs
   out <- ci(x, ci = ci, ...)
+  dof <- .safe(insight::get_df(x, type = "wald"), Inf)
   # we now iterate all confidence intervals and create an approximate normal
   # distribution that covers the CI-range.
   posterior <- as.data.frame(lapply(seq_len(nrow(out)), function(i) {
     ci_range <- as.numeric(out[i, c("CI_low", "CI_high")])
-    .generate_posterior_from_ci(ci, ci_range)
+    .generate_posterior_from_ci(ci, ci_range, dof = dof)
   }))
   colnames(posterior) <- out$Parameter
 
diff --git a/man/equivalence_test.lm.Rd b/man/equivalence_test.lm.Rd
index 91f8b2b93..8796960d9 100644
--- a/man/equivalence_test.lm.Rd
+++ b/man/equivalence_test.lm.Rd
@@ -44,8 +44,7 @@ model's data.}
 \item{ci}{Confidence Interval (CI) level. Default to \code{0.95} (\verb{95\%}).}
 
 \item{rule}{Character, indicating the rules when testing for practical
-equivalence. Can be \code{"bayes"}, \code{"classic"} or \code{"cet"}. See
-'Details'.}
+equivalence. Can be \code{"bayes"}, \code{"classic"} or \code{"cet"}. See 'Details'.}
 
 \item{verbose}{Toggle warnings and messages.}
 
@@ -67,17 +66,17 @@ A data frame.
 Compute the (conditional) equivalence test for frequentist models.
 }
 \details{
-In classical null hypothesis significance testing (NHST) within a frequentist
-framework, it is not possible to accept the null hypothesis, H0 - unlike
-in Bayesian statistics, where such probability statements are possible.
-"\link{...} one can only reject the null hypothesis if the test
+In classical null hypothesis significance testing (NHST) within a
+frequentist framework, it is not possible to accept the null hypothesis, H0 -
+unlike in Bayesian statistics, where such probability statements are
+possible. "\link{...} one can only reject the null hypothesis if the test
 statistics falls into the critical region(s), or fail to reject this
 hypothesis. In the latter case, all we can say is that no significant effect
 was observed, but one cannot conclude that the null hypothesis is true."
-(\emph{Pernet 2017}). One way to address this issues without Bayesian methods
-is \emph{Equivalence Testing}, as implemented in \code{equivalence_test()}.
-While you either can reject the null hypothesis or claim an inconclusive result
-in NHST, the equivalence test - according to \emph{Pernet} - adds a third category,
+(\emph{Pernet 2017}). One way to address this issues without Bayesian methods is
+\emph{Equivalence Testing}, as implemented in \code{equivalence_test()}. While you
+either can reject the null hypothesis or claim an inconclusive result in
+NHST, the equivalence test - according to \emph{Pernet} - adds a third category,
 \emph{"accept"}. Roughly speaking, the idea behind equivalence testing in a
 frequentist framework is to check whether an estimate and its uncertainty
 (i.e. confidence interval) falls within a region of "practical equivalence".
@@ -138,10 +137,10 @@ only uses the regular intervals.
 
 The equivalence p-value is the area of the (cumulative) confidence
 distribution that is outside of the region of equivalence. It can be
-interpreted as p-value for \emph{rejecting} the alternative hypothesis
-and \emph{accepting} the "null hypothesis" (i.e. assuming practical
-equivalence). That is, a high p-value means we reject the assumption of
-practical equivalence and accept the alternative hypothesis.
+interpreted as p-value for \emph{rejecting} the alternative hypothesis and
+\emph{accepting} the "null hypothesis" (i.e. assuming practical equivalence). That
+is, a high p-value means we reject the assumption of practical equivalence
+and accept the alternative hypothesis.
 }
 
 \subsection{Second Generation p-Value (SGPV)}{
@@ -149,16 +148,16 @@ practical equivalence and accept the alternative hypothesis.
 Second generation p-values (SGPV) were proposed as a statistic that
 represents \emph{the proportion of data-supported hypotheses that are also null
 hypotheses} \emph{(Blume et al. 2018, Lakens and Delacre 2020)}. It represents the
-proportion of the \emph{full} confidence interval range (assuming a normally
-distributed, equal-tailed interval) that is inside the ROPE. The SGPV ranges
-from zero to one. Higher values indicate that the effect is more likely to be
-practically equivalent ("not of interest").
+proportion of the \emph{full} confidence interval range (assuming a normally or
+t-distributed, equal-tailed interval, based on the model) that is inside the
+ROPE. The SGPV ranges from zero to one. Higher values indicate that the
+effect is more likely to be practically equivalent ("not of interest").
 
 Note that the assumed interval, which is used to calculate the SGPV, is an
 estimation of the \emph{full interval} based on the chosen confidence level. For
 example, if the 95\% confidence interval of a coefficient ranges from -1 to 1,
-the underlying \emph{full (normally distributed) interval} approximately ranges
-from -1.9 to 1.9, see also following code:
+the underlying \emph{full (normally or t-distributed) interval} approximately
+ranges from -1.9 to 1.9, see also following code:
 
 \if{html}{\out{<div class="sourceCode">}}\preformatted{# simulate full normal distribution
 out <- bayestestR::distribution_normal(10000, 0, 0.5)
diff --git a/man/p_direction.lm.Rd b/man/p_direction.lm.Rd
index cdf13d24c..c057861be 100644
--- a/man/p_direction.lm.Rd
+++ b/man/p_direction.lm.Rd
@@ -58,7 +58,7 @@ with Bayesian statistics (Makowski et al., 2019).
 
 In most cases, it seems that the \emph{pd} has a direct correspondence with the
 frequentist one-sided \emph{p}-value through the formula (for two-sided \emph{p}):
-\deqn{p = 2 \times (1 - p_d)}{p = 2 * (1 - pd)}
+\ifelse{html}{\out{p = 2 * (1 - p<sub>d</sub>)}}{\eqn{p = 2 \times (1 - p_d)}}
 Thus, a two-sided p-value of respectively \code{.1}, \code{.05}, \code{.01} and \code{.001} would
 correspond approximately to a \emph{pd} of \verb{95\%}, \verb{97.5\%}, \verb{99.5\%} and \verb{99.95\%}.
 See \code{\link[bayestestR:pd_to_p]{pd_to_p()}} for details.
diff --git a/man/p_significance.lm.Rd b/man/p_significance.lm.Rd
index c7f49e7a3..4c598e487 100644
--- a/man/p_significance.lm.Rd
+++ b/man/p_significance.lm.Rd
@@ -53,21 +53,22 @@ point-null by default, \code{p_significance()} compares to a range-null.
 }
 \details{
 \code{p_significance()} returns the proportion of the \emph{full} confidence
-interval range (assuming a normally distributed, equal-tailed interval) that
-is outside a certain range (the negligible effect, or ROPE, see argument
-\code{threshold}). If there are values of the distribution both below and above
-the ROPE, \code{p_significance()} returns the higher probability of a value being
-outside the ROPE. Typically, this value should be larger than 0.5 to indicate
-practical significance. However, if the range of the negligible effect is
-rather large compared to the range of the confidence interval,
-\code{p_significance()} will be less than 0.5, which indicates no clear practical
-significance.
+interval range (assuming a normally or t-distributed, equal-tailed interval,
+based on the model) that is outside a certain range (the negligible effect,
+or ROPE, see argument \code{threshold}). If there are values of the distribution
+both below and above the ROPE, \code{p_significance()} returns the higher
+probability of a value being outside the ROPE. Typically, this value should
+be larger than 0.5 to indicate practical significance. However, if the range
+of the negligible effect is rather large compared to the range of the
+confidence interval, \code{p_significance()} will be less than 0.5, which
+indicates no clear practical significance.
 
 Note that the assumed interval, which is used to calculate the practical
 significance, is an estimation of the \emph{full interval} based on the chosen
 confidence level. For example, if the 95\% confidence interval of a
-coefficient ranges from -1 to 1, the underlying \emph{full (normally distributed)
-interval} approximately ranges from -1.9 to 1.9, see also following code:
+coefficient ranges from -1 to 1, the underlying \emph{full (normally or
+t-distributed) interval} approximately ranges from -1.9 to 1.9, see also
+following code:
 
 \if{html}{\out{<div class="sourceCode">}}\preformatted{# simulate full normal distribution
 out <- bayestestR::distribution_normal(10000, 0, 0.5)
diff --git a/tests/testthat/_snaps/p_direction.md b/tests/testthat/_snaps/p_direction.md
index 14e9c1fa9..8d61719ae 100644
--- a/tests/testthat/_snaps/p_direction.md
+++ b/tests/testthat/_snaps/p_direction.md
@@ -8,8 +8,8 @@
       Parameter   |         95% CI |     pd
       -------------------------------------
       (Intercept) | [24.44, 48.94] |   100%
-      gear        | [-1.69,  2.41] | 63.55%
-      wt          | [-4.77, -1.28] | 99.97%
-      cyl         | [-2.17,  0.55] | 87.91%
-      hp          | [-0.05,  0.01] | 90.61%
+      gear        | [-1.69,  2.41] | 63.59%
+      wt          | [-4.77, -1.28] | 99.92%
+      cyl         | [-2.17,  0.55] | 88.20%
+      hp          | [-0.05,  0.01] | 91.17%
 
diff --git a/tests/testthat/_snaps/p_significance.md b/tests/testthat/_snaps/p_significance.md
index 0678b1b85..b80688ac1 100644
--- a/tests/testthat/_snaps/p_significance.md
+++ b/tests/testthat/_snaps/p_significance.md
@@ -8,8 +8,8 @@
       Parameter   |         95% CI |     ps
       -------------------------------------
       (Intercept) | [24.44, 48.94] |   100%
-      gear        | [-1.69,  2.41] | 40.93%
-      wt          | [-4.77, -1.28] | 99.67%
+      gear        | [-1.69,  2.41] | 39.83%
+      wt          | [-4.77, -1.28] | 99.59%
       cyl         | [-2.17,  0.55] | 61.88%
       hp          | [-0.05,  0.01] |  0.00%
 
diff --git a/tests/testthat/test-model_parameters_robust.R b/tests/testthat/test-model_parameters_robust.R
index 2d4564cf5..c5a18e5bb 100644
--- a/tests/testthat/test-model_parameters_robust.R
+++ b/tests/testthat/test-model_parameters_robust.R
@@ -8,41 +8,31 @@ model <- lm(mpg ~ wt * am + cyl + gear, data = mtcars)
 test_that("model_parameters, robust CL", {
   params1 <- model_parameters(
     model,
-    robust = TRUE,
-    vcov_estimation = "CL",
-    vcov_type = "HC1",
-    verbose = FALSE
-  )
-  params2 <- model_parameters(
-    model,
-    robust = TRUE,
-    vcov_estimation = "CL",
+    vcov = "CL",
     vcov_args = list(type = "HC1"),
     verbose = FALSE
   )
   robust_se <- unname(sqrt(diag(sandwich::vcovCL(model))))
   expect_equal(params1$SE, robust_se, tolerance = 1e-3)
   expect_equal(params1$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
-  expect_equal(params2$SE, robust_se, tolerance = 1e-3)
-  expect_equal(params2$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
 })
 
 test_that("model_parameters, robust", {
-  params <- model_parameters(model, robust = TRUE, verbose = FALSE)
+  params <- model_parameters(model, vcov = "HC", verbose = FALSE)
   robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
   expect_equal(params$SE, robust_se, tolerance = 1e-3)
   expect_equal(params$p, c(0, 0.0259, 0.01478, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
 })
 
 test_that("ci, robust", {
-  params <- ci(model, robust = TRUE, verbose = FALSE)
+  params <- ci(model, vcov = "HC", verbose = FALSE)
   robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
   upper_ci <- as.vector(coef(model) + qt(0.975, df.residual(model)) * robust_se)
   expect_equal(params$CI_high, upper_ci, tolerance = 1e-3, ignore_attr = TRUE)
 })
 
 test_that("model_parameters, robust CL", {
-  params <- model_parameters(model, robust = TRUE, vcov_estimation = "CL", vcov_type = "HC1", verbose = FALSE)
+  params <- model_parameters(model, vcov = "vcovCL", verbose = FALSE)
   robust_se <- unname(sqrt(diag(sandwich::vcovCL(model))))
   expect_equal(params$SE, robust_se, tolerance = 1e-3)
   expect_equal(params$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
@@ -51,7 +41,7 @@ test_that("model_parameters, robust CL", {
 model2 <- lm(mpg ~ wt * am + cyl + gear, data = datawizard::standardize(mtcars))
 
 test_that("model_parameters, robust", {
-  params <- model_parameters(model, standardize = "refit", robust = TRUE, verbose = FALSE)
+  params <- model_parameters(model, standardize = "refit", vcov = "HC", verbose = FALSE)
   robust_se <- unname(sqrt(diag(sandwich::vcovHC(model2))))
   expect_equal(params$SE, robust_se, tolerance = 1e-3)
   expect_equal(params$p, c(0.28624, 0.0259, 0.43611, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
@@ -67,8 +57,7 @@ test_that("model_parameters, robust CR", {
   params <- model_parameters(
     model,
     robust = TRUE,
-    vcov_estimation = "CR",
-    vcov_type = "CR1",
+    vcov = "CR1",
     vcov_args = list(cluster = iris$cluster),
     verbose = FALSE
   )
@@ -86,57 +75,55 @@ data(mtcars)
 mtcars$am <- as.factor(mtcars$am)
 model <- lm(mpg ~ wt * am + cyl + gear, data = mtcars)
 
-if (packageVersion("parameters") >= "0.17.0") {
-  test_that("model_parameters, robust", {
-    expect_warning(expect_warning(expect_warning(model_parameters(model, robust = TRUE))))
-    params <- model_parameters(model, vcov = "HC3")
-    robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
-    expect_equal(params$SE, robust_se, tolerance = 1e-3)
-    expect_equal(params$p, c(0, 0.0259, 0.01478, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
-  })
-
-  test_that("ci, robust", {
-    params <- ci(model, vcov = "HC3")
-    robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
-    upper_ci <- as.vector(coef(model) + qt(0.975, df.residual(model)) * robust_se)
-    expect_equal(params$CI_high, upper_ci, tolerance = 1e-3, ignore_attr = TRUE)
-  })
-
-  test_that("model_parameters, robust CL", {
-    params <- model_parameters(model, vcov = "vcovCL", vcov_args = list(type = "HC1"))
-    robust_se <- unname(sqrt(diag(sandwich::vcovCL(model))))
-    expect_equal(params$SE, robust_se, tolerance = 1e-3)
-    expect_equal(params$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
-  })
-
-  d <- datawizard::standardize(mtcars)
-  model2 <- lm(mpg ~ wt * am + cyl + gear, data = d)
-
-  test_that("model_parameters, robust", {
-    params <- model_parameters(model, standardize = "refit", vcov = "HC3")
-    robust_se <- unname(sqrt(diag(sandwich::vcovHC(model2))))
-    expect_equal(params$SE, robust_se, tolerance = 1e-3)
-    expect_equal(params$p, c(0.28624, 0.0259, 0.43611, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
-  })
-
-
-  # cluster-robust standard errors, using clubSandwich
-  data(iris)
-  model <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
-  iris$cluster <- factor(rep(LETTERS[1:8], length.out = nrow(iris)))
-
-  test_that("model_parameters, robust CR", {
-    params <- model_parameters(model, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$cluster))
-    robust_se <- unname(sqrt(diag(clubSandwich::vcovCR(model, type = "CR1", cluster = iris$cluster))))
-    expect_equal(params$SE, robust_se, tolerance = 1e-3)
-    expect_equal(params$p, c(0.01246, 0.04172, 0.18895, 0.57496, 0, 0), tolerance = 1e-3)
-  })
-
-  test_that("model_parameters, normal", {
-    params <- model_parameters(model)
-    expect_equal(params$p, c(0.13267, 0.21557, 0.36757, 0.77012, 3e-05, 0), tolerance = 1e-3)
-  })
-}
+test_that("model_parameters, robust", {
+  expect_warning(expect_warning(expect_warning(model_parameters(model, robust = TRUE))))
+  params <- model_parameters(model, vcov = "HC3")
+  robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
+  expect_equal(params$SE, robust_se, tolerance = 1e-3)
+  expect_equal(params$p, c(0, 0.0259, 0.01478, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
+})
+
+test_that("ci, robust", {
+  params <- ci(model, vcov = "HC3")
+  robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
+  upper_ci <- as.vector(coef(model) + qt(0.975, df.residual(model)) * robust_se)
+  expect_equal(params$CI_high, upper_ci, tolerance = 1e-3, ignore_attr = TRUE)
+})
+
+test_that("model_parameters, robust CL", {
+  params <- model_parameters(model, vcov = "vcovCL", vcov_args = list(type = "HC1"))
+  robust_se <- unname(sqrt(diag(sandwich::vcovCL(model))))
+  expect_equal(params$SE, robust_se, tolerance = 1e-3)
+  expect_equal(params$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
+})
+
+d <- datawizard::standardize(mtcars)
+model2 <- lm(mpg ~ wt * am + cyl + gear, data = d)
+
+test_that("model_parameters, robust", {
+  params <- model_parameters(model, standardize = "refit", vcov = "HC3")
+  robust_se <- unname(sqrt(diag(sandwich::vcovHC(model2))))
+  expect_equal(params$SE, robust_se, tolerance = 1e-3)
+  expect_equal(params$p, c(0.28624, 0.0259, 0.43611, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
+})
+
+
+# cluster-robust standard errors, using clubSandwich
+data(iris)
+model <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
+iris$cluster <- factor(rep(LETTERS[1:8], length.out = nrow(iris)))
+
+test_that("model_parameters, robust CR", {
+  params <- model_parameters(model, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$cluster))
+  robust_se <- unname(sqrt(diag(clubSandwich::vcovCR(model, type = "CR1", cluster = iris$cluster))))
+  expect_equal(params$SE, robust_se, tolerance = 1e-3)
+  expect_equal(params$p, c(0.01246, 0.04172, 0.18895, 0.57496, 0, 0), tolerance = 1e-3)
+})
+
+test_that("model_parameters, normal", {
+  params <- model_parameters(model)
+  expect_equal(params$p, c(0.13267, 0.21557, 0.36757, 0.77012, 3e-05, 0), tolerance = 1e-3)
+})
 
 test_that("ci_ml1, robust", {
   skip("TODO: this one actually is not correct.")
diff --git a/tests/testthat/test-model_parameters_std_mixed.R b/tests/testthat/test-model_parameters_std_mixed.R
index a4640037e..7aa2518a9 100644
--- a/tests/testthat/test-model_parameters_std_mixed.R
+++ b/tests/testthat/test-model_parameters_std_mixed.R
@@ -266,9 +266,8 @@ test_that("model_parameters, standardize-refit robust", {
     standardize = "refit",
     effects = "fixed",
     robust = TRUE,
-    vcov_estimation = "CR",
-    vcov_type = "CR1",
-    vcov_args = list(cluster = iris$grp),
+    vcov = "CR",
+    vcov_args = list(type = "CR1", cluster = iris$grp),
     verbose = FALSE
   )
   expect_equal(c(nrow(params), ncol(params)), c(7, 10))
diff --git a/tests/testthat/test-p_direction.R b/tests/testthat/test-p_direction.R
index dff61e296..0eb3d4902 100644
--- a/tests/testthat/test-p_direction.R
+++ b/tests/testthat/test-p_direction.R
@@ -1,5 +1,6 @@
 skip_on_cran()
 skip_if_not_installed("bayestestR")
+skip_if_not_installed("distributional")
 
 test_that("p_direction", {
   data(mtcars)
@@ -12,11 +13,11 @@ test_that("p_direction", {
 
   set.seed(123)
   x <- p_direction(m, ci = 0.8)
-  expect_equal(x$pd, c(1, 0.6382, 0.9997, 0.884, 0.9107), tolerance = 1e-4)
+  expect_equal(x$pd, c(1, 0.6359, 0.9992, 0.882, 0.9117), tolerance = 1e-3)
 
   set.seed(123)
   x <- p_direction(m, null = 0.2)
-  expect_equal(x$pd, c(1, 0.5617, 0.9999, 0.9276, 1), tolerance = 1e-4)
+  expect_equal(x$pd, c(1, 0.5567, 0.9997, 0.9309, 1), tolerance = 1e-3)
 })
 
 test_that("p_direction, glmmTMB", {
diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R
index c95756e01..7ec769851 100644
--- a/tests/testthat/test-p_significance.R
+++ b/tests/testthat/test-p_significance.R
@@ -12,11 +12,11 @@ test_that("p_significance", {
 
   set.seed(123)
   x <- p_significance(m, ci = 0.8)
-  expect_equal(x$ps, c(1, 0.4074, 0.9972, 0.6212, 0), tolerance = 1e-4)
+  expect_equal(x$ps, c(1, 0.3983, 0.9959, 0.6188, 0), tolerance = 1e-3)
 
   set.seed(123)
   x <- p_significance(m, threshold = 0.5)
-  expect_equal(x$ps, c(1, 0.4478, 0.9977, 0.6737, 0), tolerance = 1e-4)
+  expect_equal(x$ps, c(1, 0.4393, 0.9969, 0.6803, 0), tolerance = 1e-4)
 })
 
 test_that("p_significance, glmmTMB", {
diff --git a/vignettes/model_parameters_robust.Rmd b/vignettes/model_parameters_robust.Rmd
index 86fbc1638..8718d01f5 100644
--- a/vignettes/model_parameters_robust.Rmd
+++ b/vignettes/model_parameters_robust.Rmd
@@ -1,6 +1,6 @@
 ---
 title: "Robust Estimation of Standard Errors, Confidence Intervals, and p-values"
-output: 
+output:
   rmarkdown::html_vignette:
     toc: true
     fig_width: 10.08
@@ -10,7 +10,7 @@ vignette: >
   %\VignetteIndexEntry{Robust Estimation of Standard Errors, Confidence Intervals, and p-values}
   \usepackage[utf8]{inputenc}
   %\VignetteEngine{knitr::rmarkdown}
-editor_options: 
+editor_options:
   chunk_output_type: console
 bibliography: bibliography.bib
 ---
@@ -72,12 +72,13 @@ There are two arguments (see
 [`?standard_error`](https://easystats.github.io/parameters/reference/standard_error.html)
 for further details) that allow for choosing different methods and options of
 robust estimation:
+
 - `vcov`
 - `vcov_args`
 
 Let us start with a simple example, which uses a heteroskedasticity-consistent
 covariance matrix estimation with estimation-type "HC3" (i.e.
-`sandwich::vcovHC(type = "HC3")`). 
+`sandwich::vcovHC(type = "HC3")`).
 
 First let's create a simple linear regression model, which we know violates
 homoscedasticity assumption, and thus robust estimation methods are to be
@@ -92,8 +93,8 @@ check_heteroscedasticity(model)
 ```
 
 We would extract model parameters both with and without robust estimation to
-highlight the difference it makes to standard errors, confidence intervals, 
-*t*-statistic, and *p*-values. Also, note that the coefficient estimate remains 
+highlight the difference it makes to standard errors, confidence intervals,
+*t*-statistic, and *p*-values. Also, note that the coefficient estimate remains
 unchanged.
 
 ```{r}
@@ -115,7 +116,7 @@ If a different type of covariance matrix estimation is required, use the
 `vcov`-argument. This argument accepts the name of a function from the
 *sandwich* or *clubSandwich* packages as a string, such as `"vcovCL"` (or just
 its suffix `"CL"`). *parameters* will then call the corresponding function with
-the content of `vcov_args` as arguments. 
+the content of `vcov_args` as arguments.
 
 The specific estimation type can be controlled by passing a `type` argument via `vcov_args`. See `?sandwich::vcovCL` for information about the different types of covariance matrices that this function can produce (`HC0` to `HC3`). In the next
 example, we use a clustered covariance matrix estimation with `HC1`-estimation
@@ -231,7 +232,7 @@ model_parameters(
 ```
 
 Notice that robust estimation returns different standard errors, confidence
-intervals, test statistic and *p*-values compared to the standard estimation. 
+intervals, test statistic and *p*-values compared to the standard estimation.
 Also, note that the coefficient estimate remains unchanged.
 
 ### Robust Covariance Matrix Estimation on Standardized Mixed Model Parameters