diff --git a/DESCRIPTION b/DESCRIPTION
index f34136fd9..1bfba562d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.21.5.6
+Version: 0.21.5.7
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -182,7 +182,7 @@ Suggests:
     PMCMRplus,
     poorman,
     posterior,
-    PROreg,
+    PROreg (>= 1.3.0),
     pscl,
     psych,
     pvclust,
@@ -199,7 +199,7 @@ Suggests:
     survival,
     testthat (>= 3.2.1),
     tidyselect,
-    tinytable,
+    tinytable (>= 0.1.0),
     TMB,
     truncreg,
     VGAM,
diff --git a/NEWS.md b/NEWS.md
index 970877c72..3e9bfe8f9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -9,6 +9,13 @@
 * `include_reference` can now directly be set to `TRUE` in `model_parameters()`
   and doesn't require a call to `print()` anymore.
 
+* `compare_parameters()` gains a `include_reference` argument, to add the
+  reference category of categorical predictors to the parameters table.
+
+* `print_md()` for `compare_parameters()` now by default uses the *tinytable*
+  package to create markdown tables. This allows better control for column
+  heading spanning over multiple columns.
+
 ## Bug fixes
 
 * Fixed issue with parameter names for `model_parameters()` and objects from
diff --git a/R/compare_parameters.R b/R/compare_parameters.R
index cf6285cb7..2c6fa7943 100644
--- a/R/compare_parameters.R
+++ b/R/compare_parameters.R
@@ -83,6 +83,8 @@ compare_parameters <- function(...,
                                coefficient_names = NULL,
                                keep = NULL,
                                drop = NULL,
+                               include_reference = FALSE,
+                               groups = NULL,
                                verbose = TRUE) {
   models <- list(...)
 
@@ -186,6 +188,7 @@ compare_parameters <- function(...,
         keep = keep,
         drop = drop,
         wb_component = FALSE,
+        include_reference = include_reference,
         verbose = verbose
       )
     }
@@ -253,6 +256,7 @@ compare_parameters <- function(...,
   attr(all_models, "model_names") <- gsub("\"", "", unlist(lapply(model_names, insight::safe_deparse)), fixed = TRUE)
   attr(all_models, "output_style") <- select
   attr(all_models, "all_attributes") <- object_attributes
+  attr(all_models, "parameter_groups") <- groups
   class(all_models) <- c("compare_parameters", "see_compare_parameters", unique(class(all_models)))
 
   all_models
diff --git a/R/display.R b/R/display.R
index deb95ca05..ba45e6a66 100644
--- a/R/display.R
+++ b/R/display.R
@@ -53,7 +53,8 @@
 #' experimental, thus, only a fixed layout-style is available at the moment
 #' (columns for estimates, confidence intervals and p-values). However, it
 #' is possible to include other model components, like zero-inflation, or random
-#' effects in the table. See 'Examples'.
+#' effects in the table. See 'Examples'. An alternative is to set `engine = "tt"`
+#' in `print_html()` to use the _tinytable_ package for creating HTML tables.
 #'
 #' @seealso [print.parameters_model()]
 #'
diff --git a/R/format.R b/R/format.R
index c12e20a86..790489abf 100644
--- a/R/format.R
+++ b/R/format.R
@@ -253,6 +253,7 @@ format.compare_parameters <- function(x,
                                       zap_small = FALSE,
                                       format = NULL,
                                       groups = NULL,
+                                      engine = NULL,
                                       ...) {
   m_class <- attributes(x)$model_class
   x$Method <- NULL
@@ -393,7 +394,7 @@ format.compare_parameters <- function(x,
   # group parameters - this function find those parameters that should be
   # grouped, reorders parameters into groups and indents lines that belong
   # to one group, adding a header for each group
-  if (!is.null(groups)) {
+  if (!is.null(groups) && !identical(engine, "tt")) {
     out <- .parameter_groups(out, groups)
   }
   indent_groups <- attributes(x)$indent_groups
diff --git a/R/print.compare_parameters.R b/R/print.compare_parameters.R
index 2bcdab704..8845d0701 100644
--- a/R/print.compare_parameters.R
+++ b/R/print.compare_parameters.R
@@ -31,6 +31,9 @@ print.compare_parameters <- function(x,
   if (missing(select)) {
     select <- attributes(x)$output_style
   }
+  if (missing(groups)) {
+    groups <- attributes(x)$parameter_groups
+  }
 
   formatted_table <- format(
     x,
diff --git a/R/print.parameters_model.R b/R/print.parameters_model.R
index b707ed422..422791258 100644
--- a/R/print.parameters_model.R
+++ b/R/print.parameters_model.R
@@ -137,6 +137,10 @@
 #' default for the `select` argument. See argument's documentation for available
 #' options.
 #'
+#' - `easystats_html_engine`: `options(easystats_html_engine = "gt")` will set
+#' the default HTML engine for tables to `gt`, i.e. the _gt_ package is used to
+#' create HTML tables. If set to `tt`, the _tinytable_ package is used.
+#'
 #' @details `summary()` is a convenient shortcut for
 #'   `print(object, select = "minimal", show_sigma = TRUE, show_formula = TRUE)`.
 #'
@@ -308,10 +312,10 @@ print.parameters_model <- function(x,
     footer <- ""
   }
   if (!identical(footer, "")) {
-    if (!is.null(footer)) {
-      footer <- paste0("\n", footer, "\n", footer_stats)
-    } else {
+    if (is.null(footer)) {
       footer <- footer_stats
+    } else {
+      footer <- paste0("\n", footer, "\n", footer_stats)
     }
   }
 
@@ -412,7 +416,7 @@ print.parameters_random <- function(x, digits = 2, ...) {
                           show_formula = FALSE,
                           format = "text") {
   # get attributes
-  sigma <- attributes(x)$sigma
+  model_sigma <- attributes(x)$sigma
   show_summary <- isTRUE(attributes(x)$show_summary)
   verbose <- .additional_arguments(x, "verbose", TRUE)
 
@@ -423,7 +427,7 @@ print.parameters_random <- function(x, digits = 2, ...) {
   show_r2 <- .additional_arguments(x, "show_summary", FALSE)
 
   # set defaults, if necessary
-  if (is.null(sigma)) {
+  if (is.null(model_sigma)) {
     show_sigma <- FALSE
   }
 
@@ -522,10 +526,10 @@ print.parameters_random <- function(x, digits = 2, ...) {
   random_params$Term[is.na(random_params$Term)] <- ""
   random_params$SD[is.na(random_params$SD)] <- ""
 
-  non_empty <- random_params$Term != "" & random_params$Type != ""
+  non_empty <- random_params$Term != "" & random_params$Type != "" # nolint
   random_params$Line[non_empty] <- sprintf("%s (%s)", random_params$Type[non_empty], random_params$Term[non_empty])
 
-  non_empty <- random_params$Term != "" & random_params$Type == ""
+  non_empty <- random_params$Term != "" & random_params$Type == "" # nolint
   random_params$Line[non_empty] <- sprintf("%s", random_params$Term[non_empty])
 
   # final fix, indentions
@@ -567,11 +571,11 @@ print.parameters_random <- function(x, digits = 2, ...) {
   col_width <- rep(NA, length(shared_cols))
   for (i in seq_along(shared_cols)) {
     col_width[i] <- max(unlist(lapply(formatted_table, function(j) {
-      col <- j[[shared_cols[i]]]
-      if (!is.null(col)) {
-        max(nchar(col))
-      } else {
+      column <- j[[shared_cols[i]]]
+      if (is.null(column)) {
         NA
+      } else {
+        max(nchar(column))
       }
     })))
   }
diff --git a/R/print_html.R b/R/print_html.R
index 9b21d4180..d71aa194b 100644
--- a/R/print_html.R
+++ b/R/print_html.R
@@ -46,6 +46,9 @@ print_html.parameters_model <- function(x,
   if (missing(select) || is.null(select)) {
     select <- attributes(x)$output_style
   }
+  if (missing(groups)) {
+    groups <- attributes(x)$parameter_groups
+  }
 
   # we need glue-like syntax right now...
   if (!is.null(select)) {
@@ -150,6 +153,7 @@ print_html.compare_parameters <- function(x,
                                           font_size = "100%",
                                           line_padding = 4,
                                           column_labels = NULL,
+                                          engine = "gt",
                                           ...) {
   # check if user supplied digits attributes
   if (missing(digits)) {
@@ -169,6 +173,29 @@ print_html.compare_parameters <- function(x,
     select <- attributes(x)$output_style
   }
 
+  # markdown engine?
+  engine <- match.arg(getOption("easystats_html_engine", engine), c("gt", "default", "tt"))
+
+  # for tiny table, we can just call print_md()
+  if (engine == "tt") {
+    return(print_md(
+      x,
+      digits = digits,
+      ci_digits = ci_digits,
+      p_digits = p_digits,
+      caption = caption,
+      subtitle = subtitle,
+      footer = footer,
+      select = select,
+      split_components = TRUE,
+      ci_brackets = ci_brackets,
+      zap_small = zap_small,
+      groups = groups,
+      engine = "tt",
+      outformat = "html"
+    ))
+  }
+
   # we need glue-like syntax right now...
   select <- .convert_to_glue_syntax(style = select, "<br>")
 
diff --git a/R/print_md.R b/R/print_md.R
index 7d2237d3b..98a7adbec 100644
--- a/R/print_md.R
+++ b/R/print_md.R
@@ -91,10 +91,10 @@ print_md.parameters_model <- function(x,
     footer <- ""
   }
   if (!identical(footer, "")) {
-    if (!is.null(footer)) {
-      footer <- paste0("\n", footer, "\n", footer_stats)
-    } else {
+    if (is.null(footer)) {
       footer <- footer_stats
+    } else {
+      footer <- paste0("\n", footer, "\n", footer_stats)
     }
   }
 
@@ -136,6 +136,7 @@ print_md.compare_parameters <- function(x,
                                         ci_brackets = c("(", ")"),
                                         zap_small = FALSE,
                                         groups = NULL,
+                                        engine = "tt",
                                         ...) {
   # check if user supplied digits attributes
   if (missing(digits)) {
@@ -152,6 +153,12 @@ print_md.compare_parameters <- function(x,
   if (missing(select) || is.null(select)) {
     select <- attributes(x)$output_style
   }
+  if (missing(groups)) {
+    groups <- attributes(x)$parameter_groups
+  }
+
+  # markdown engine?
+  engine <- match.arg(engine, c("tt", "default"))
 
   formatted_table <- format(
     x,
@@ -164,16 +171,35 @@ print_md.compare_parameters <- function(x,
     ci_brackets = ci_brackets,
     format = "markdown",
     zap_small = zap_small,
-    groups = groups
+    groups = groups,
+    engine = engine
   )
 
-  insight::export_table(
-    formatted_table,
-    format = "markdown",
-    caption = caption,
-    subtitle = subtitle,
-    footer = footer
-  )
+  if (identical(engine, "tt")) {
+    # retrieve output format - print_md() may be called from print_html()
+    dots <- list(...)
+    if (identical(dots$outformat, "html")) {
+      outformat <- "html"
+    } else {
+      outformat <- "markdown"
+    }
+    .export_table_tt(
+      x,
+      formatted_table,
+      groups,
+      caption = caption,
+      footer = footer,
+      outformat = outformat
+    )
+  } else {
+    insight::export_table(
+      formatted_table,
+      format = "markdown",
+      caption = caption,
+      subtitle = subtitle,
+      footer = footer
+    )
+  }
 }
 
 
@@ -262,14 +288,12 @@ print_md.equivalence_test_lm <- function(x, digits = 2, ci_brackets = c("(", ")"
   rule <- attributes(x)$rule
   rope <- attributes(x)$rope
 
-  if (!is.null(rule)) {
-    if (rule == "cet") {
-      table_caption <- "Conditional Equivalence Testing"
-    } else if (rule == "classic") {
-      table_caption <- "TOST-test for Practical Equivalence"
-    } else {
-      table_caption <- "Test for Practical Equivalence"
-    }
+  if (is.null(rule)) {
+    table_caption <- "Test for Practical Equivalence"
+  } else if (rule == "cet") {
+    table_caption <- "Conditional Equivalence Testing"
+  } else if (rule == "classic") {
+    table_caption <- "TOST-test for Practical Equivalence"
   } else {
     table_caption <- "Test for Practical Equivalence"
   }
@@ -327,3 +351,100 @@ print_md.parameters_distribution <- function(x, digits = 2, ci_brackets = c("(",
 
   insight::export_table(formatted_table, format = "markdown", align = "firstleft", ...)
 }
+
+
+# helper -----------------------
+
+.export_table_tt <- function(x, formatted_table, groups, caption = NULL, footer = NULL, outformat = "markdown") {
+  insight::check_if_installed("tinytable", minimum_version = "0.1.0")
+  row_groups <- NULL
+  # check if we have a list of tables
+  if (!is.data.frame(formatted_table) && is.list(formatted_table) && length(formatted_table) > 1) {
+    # sanity check - cannot combine multiple tables when we have groups
+    if (!is.null(groups)) {
+      insight::format_error("Cannot combine multiple tables when groups are present.")
+    }
+    # add table caption as group variable, and bind tables
+    # we then extract row headers based on values in the group indices
+    formatted_table <- lapply(formatted_table, function(i) {
+      i$group <- attr(i, "table_caption")
+      i
+    })
+    # bind tables
+    formatted_table <- do.call(rbind, formatted_table)
+    # find positions for sub headers
+    row_groups <- as.list(which(!duplicated(formatted_table$group)))
+    names(row_groups) <- formatted_table$group[unlist(row_groups)]
+    # remove no longer needed group variable
+    formatted_table$group <- NULL
+  }
+  # we need to find out which columns refer to which model, in order to
+  # add a column heading for each model
+  models <- attributes(x)$model_names
+  col_names <- gsub("(.*) \\((.*)\\)$", "\\2", colnames(formatted_table))
+  col_groups <- sapply(models, function(i) which(i == col_names), simplify = FALSE)
+  # clean column names. These still contain the model name
+  colnames(formatted_table) <- gsub("(.*) \\((.*)\\)$", "\\1", colnames(formatted_table))
+  # check if we have column spans at all?
+  if (all(lengths(col_groups) == 1)) {
+    col_groups <- NULL
+  }
+  # group rows?
+  if (!is.null(groups)) {
+    # make sure we have numeric indices for groups
+    groups <- lapply(groups, function(g) {
+      if (is.character(g)) {
+        # if groups were provided as parameter names, we find the row position
+        # by matching the parameter name
+        match(g, formatted_table$Parameter)
+      } else {
+        # else, we assume that the group is a row position
+        g
+      }
+    })
+    # sanity check - do all rows match a parameter?
+    group_indices <- unlist(groups, use.names = FALSE)
+    if (anyNA(group_indices) || any(group_indices < 1) || any(group_indices > nrow(formatted_table))) {
+      insight::format_error("Some group indices do not match any parameter.")
+    }
+    # if row indices are not sorted, we need to resort the parameters data frame
+    if (is.unsorted(unlist(groups))) {
+      new_rows <- c(unlist(groups), setdiff(seq_len(nrow(formatted_table)), unlist(groups)))
+      formatted_table <- formatted_table[new_rows, ]
+      # we need to update indices in groups as well. Therefore, we need to convert
+      # list of row indices into a vector with row indices, then subtract the
+      # differences of old and new row positions, and then split that vector into
+      # a list again
+      groups <- stats::setNames(unlist(groups), rep(names(groups), lengths(groups)))
+      groups <- groups - (unlist(groups) - sort(unlist(groups)))
+      groups <- split(unname(groups), factor(names(groups), levels = unique(names(groups))))
+    }
+    # find matching rows for groups
+    row_groups <- lapply(seq_along(groups), function(i) {
+      g <- groups[[i]]
+      if (is.character(g)) {
+        # if groups were provided as parameter names, we find the row position
+        # by matching the parameter name
+        g <- match(g, formatted_table$Parameter)[1]
+      } else {
+        # else, we assume that the group is a row position
+        g <- g[1]
+      }
+      g
+    })
+    # set element names
+    names(row_groups) <- names(groups)
+    if (identical(outformat, "markdown")) {
+      # for markdown, format italic
+      names(row_groups) <- paste0("*", names(row_groups), "*")
+    }
+  }
+  # replace NA in formatted_table by ""
+  formatted_table[is.na(formatted_table)] <- ""
+  # create base table
+  out <- tinytable::tt(formatted_table, notes = footer, caption = caption)
+  # insert sub header rows and column spans
+  out <- tinytable::group_tt(out, i = row_groups, j = col_groups)
+  out@output <- outformat
+  out
+}
diff --git a/R/print_table.R b/R/print_table.R
index 029f40c3a..d86399c09 100644
--- a/R/print_table.R
+++ b/R/print_table.R
@@ -154,20 +154,18 @@ print_table <- function(x, digits = 2, p_digits = 3, theme = "default", ...) {
   }
 
   # base table
-  out <- tinytable::tt(x, caption = NULL, notes = NULL, ...)
+  out <- tinytable::tt(as.data.frame(x), caption = NULL, notes = NULL, ...)
   # add subheaders, if any
-  if (!is.null(row_header_labels)) {
+  if (is.null(row_header_labels)) {
+    out <- tinytable::group_tt(out, j = col_groups)
+  } else {
     out <- tinytable::group_tt(out, i = row_header_labels, j = col_groups)
     out <- tinytable::style_tt(out, i = row_header_pos, italic = TRUE)
-  } else {
-    out <- tinytable::group_tt(out, j = col_groups)
   }
   # style table
   out <- insight::apply_table_theme(out, x, theme = theme, sub_header_positions = row_header_pos)
-  # workaround, to make sure HTML is default output
-  m <- attr(out, "tinytable_meta")
-  m$output <- "html"
-  attr(out, "tinytable_meta") <- m
+  # make sure HTML is default output
+  out@output <- "html"
 
   out
 }
diff --git a/R/utils_format.R b/R/utils_format.R
index 974e4da14..aba6a8381 100644
--- a/R/utils_format.R
+++ b/R/utils_format.R
@@ -53,7 +53,9 @@
     # here we either have "<br>" or " " as line breaks, followed by empty "()"
     i <- gsub("<br>()", "", i, fixed = TRUE)
     i <- gsub(" ()", "", i, fixed = TRUE)
+    i <- gsub(" (, )", "", i, fixed = TRUE)
     i[i == "()"] <- ""
+    i[i == "(, )"] <- ""
     # remove other non-matched patterns
     i <- gsub("{stars}", "", i, fixed = TRUE)
     i <- gsub("{rhat}", "", i, fixed = TRUE)
diff --git a/inst/WORDLIST b/inst/WORDLIST
index d31a439cc..7637348ab 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -313,6 +313,7 @@ subscales
 systemfit
 th
 tidymodels
+tinytable
 tseries
 unicode
 varEST
diff --git a/man/compare_parameters.Rd b/man/compare_parameters.Rd
index db1dbf328..0218c9bb0 100644
--- a/man/compare_parameters.Rd
+++ b/man/compare_parameters.Rd
@@ -20,6 +20,8 @@ compare_parameters(
   coefficient_names = NULL,
   keep = NULL,
   drop = NULL,
+  include_reference = FALSE,
+  groups = NULL,
   verbose = TRUE
 )
 
@@ -38,6 +40,8 @@ compare_models(
   coefficient_names = NULL,
   keep = NULL,
   drop = NULL,
+  include_reference = FALSE,
+  groups = NULL,
   verbose = TRUE
 )
 }
@@ -179,6 +183,22 @@ names.}
 
 \item{drop}{See \code{keep}.}
 
+\item{include_reference}{Logical, if \code{TRUE}, the reference level of factors will
+be added to the parameters table. This is only relevant for models with
+categorical predictors. The coefficient for the reference level is always
+\code{0} (except when \code{exponentiate = TRUE}, then the coefficient will be \code{1}),
+so this is just for completeness.}
+
+\item{groups}{Named list, can be used to group parameters in the printed output.
+List elements may either be character vectors that match the name of those
+parameters that belong to one group, or list elements can be row numbers
+of those parameter rows that should belong to one group. The names of the
+list elements will be used as group names, which will be inserted as "header
+row". A possible use case might be to emphasize focal predictors and control
+variables, see 'Examples'. Parameters will be re-ordered according to the
+order used in \code{groups}, while all non-matching parameters will be added
+to the end.}
+
 \item{verbose}{Toggle warnings and messages.}
 }
 \value{
diff --git a/man/display.parameters_model.Rd b/man/display.parameters_model.Rd
index 1a36422f7..f351e0c63 100644
--- a/man/display.parameters_model.Rd
+++ b/man/display.parameters_model.Rd
@@ -305,7 +305,8 @@ which prints the output as a formatted HTML table. It is still somewhat
 experimental, thus, only a fixed layout-style is available at the moment
 (columns for estimates, confidence intervals and p-values). However, it
 is possible to include other model components, like zero-inflation, or random
-effects in the table. See 'Examples'.
+effects in the table. See 'Examples'. An alternative is to set \code{engine = "tt"}
+in \code{print_html()} to use the \emph{tinytable} package for creating HTML tables.
 }
 \examples{
 \dontshow{if (require("gt", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
diff --git a/man/model_parameters.Rd b/man/model_parameters.Rd
index 70e012556..d124d2cc9 100644
--- a/man/model_parameters.Rd
+++ b/man/model_parameters.Rd
@@ -374,6 +374,9 @@ will replace the interaction mark (by default, \code{*}) with the related charac
 \item \code{parameters_select}: \verb{options(parameters_select = <value>)} will set the
 default for the \code{select} argument. See argument's documentation for available
 options.
+\item \code{easystats_html_engine}: \code{options(easystats_html_engine = "gt")} will set
+the default HTML engine for tables to \code{gt}, i.e. the \emph{gt} package is used to
+create HTML tables. If set to \code{tt}, the \emph{tinytable} package is used.
 }
 }
 
diff --git a/man/print.parameters_model.Rd b/man/print.parameters_model.Rd
index 5b602b1d9..18bec59c7 100644
--- a/man/print.parameters_model.Rd
+++ b/man/print.parameters_model.Rd
@@ -185,6 +185,9 @@ will replace the interaction mark (by default, \code{*}) with the related charac
 \item \code{parameters_select}: \verb{options(parameters_select = <value>)} will set the
 default for the \code{select} argument. See argument's documentation for available
 options.
+\item \code{easystats_html_engine}: \code{options(easystats_html_engine = "gt")} will set
+the default HTML engine for tables to \code{gt}, i.e. the \emph{gt} package is used to
+create HTML tables. If set to \code{tt}, the \emph{tinytable} package is used.
 }
 }
 
diff --git a/tests/testthat/_snaps/compare_parameters.md b/tests/testthat/_snaps/compare_parameters.md
index dac4e1410..2e807d8ab 100644
--- a/tests/testthat/_snaps/compare_parameters.md
+++ b/tests/testthat/_snaps/compare_parameters.md
@@ -44,3 +44,127 @@
       --------------------------------------------------------
       SD (Intercept: persons) |    |    |  1.08 ( 0.49,  2.37)
 
+# compare_parameters, print_md
+
+    Code
+      print(out)
+    Output
+      +----------------+-----------------------+--------+----------------------+--------+
+      |                | lm1                            | lm2                           |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Parameter      | Estimate (ci)         | p      | Estimate (ci)        | p      |
+      +================+=======================+========+======================+========+
+      | *Groups*                                                                        |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | grp (2)        | -4.31 (-15.95, 7.32)  | 0.465  | 0.32 (-22.56, 23.20) | 0.978  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | grp (3)        | -1.31 (-13.47, 10.84) | 0.831  | 3.77 (-19.72, 27.26) | 0.752  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | *Interactions*                                                                  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Days * grp (2) |                       |        | -1.01 (-5.35, 3.32)  | 0.645  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Days * grp (3) |                       |        | -1.11 (-5.53, 3.31)  | 0.621  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | *Controls*                                                                      |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Days           | 10.44 (8.84, 12.03)   | <0.001 | 11.23 (7.87, 14.60)  | <0.001 |
+      +----------------+-----------------------+--------+----------------------+--------+
+      |                |                       |        |                      |        |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Observations   | 180                   |        | 180                  |        |
+      +----------------+-----------------------+--------+----------------------+--------+ 
+
+---
+
+    Code
+      print_md(cp)
+    Output
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      |                         | lm1                            | lm2                           |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | Parameter               | Estimate (ci)         | p      | Estimate (ci)        | p      |
+      +=========================+=======================+========+======================+========+
+      | Fixed Effects                                                                            |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | Days                    | 10.44 (8.84, 12.03)   | <0.001 | 11.23 (7.87, 14.60)  | <0.001 |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | grp (2)                 | -4.31 (-15.95, 7.32)  | 0.465  | 0.32 (-22.56, 23.20) | 0.978  |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | grp (3)                 | -1.31 (-13.47, 10.84) | 0.831  | 3.77 (-19.72, 27.26) | 0.752  |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | Days * grp (2)          |                       |        | -1.01 (-5.35, 3.32)  | 0.645  |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | Days * grp (3)          |                       |        | -1.11 (-5.53, 3.31)  | 0.621  |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | Random Effects                                                                           |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | SD (Intercept: Subject) | 37.06 (25.85, 53.13)  |        | 37.08 (25.85, 53.19) |        |
+      +-------------------------+-----------------------+--------+----------------------+--------+
+      | SD (Residual)           | 31.13 (27.89, 34.75)  |        | 31.30 (28.02, 34.96) |        |
+      +-------------------------+-----------------------+--------+----------------------+--------+ 
+
+---
+
+    Code
+      print(out)
+    Output
+      +----------------+-----------------------+--------+----------------------+--------+
+      |                | lm1                            | lm2                           |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Parameter      | Estimate (ci)         | p      | Estimate (ci)        | p      |
+      +================+=======================+========+======================+========+
+      | *Groups*                                                                        |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | grp (1)        | 0.00                  |        | 0.00                 |        |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | grp (2)        | -4.31 (-15.95, 7.32)  | 0.465  | 0.32 (-22.56, 23.20) | 0.978  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | grp (3)        | -1.31 (-13.47, 10.84) | 0.831  | 3.77 (-19.72, 27.26) | 0.752  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | *Interactions*                                                                  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Days * grp (2) |                       |        | -1.01 (-5.35, 3.32)  | 0.645  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Days * grp (3) |                       |        | -1.11 (-5.53, 3.31)  | 0.621  |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | *Controls*                                                                      |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Days           | 10.44 (8.84, 12.03)   | <0.001 | 11.23 (7.87, 14.60)  | <0.001 |
+      +----------------+-----------------------+--------+----------------------+--------+
+      |                |                       |        |                      |        |
+      +----------------+-----------------------+--------+----------------------+--------+
+      | Observations   | 180                   |        | 180                  |        |
+      +----------------+-----------------------+--------+----------------------+--------+ 
+
+---
+
+    Code
+      print(out)
+    Output
+      +----------------+-----------------------+----------------------+
+      | Parameter      | lm1                   | lm2                  |
+      +================+=======================+======================+
+      | *Groups*                                                      |
+      +----------------+-----------------------+----------------------+
+      | grp (1)        | 0.00                  | 0.00                 |
+      +----------------+-----------------------+----------------------+
+      | grp (2)        | -4.31 (-15.95, 7.32)  | 0.32 (-22.56, 23.20) |
+      +----------------+-----------------------+----------------------+
+      | grp (3)        | -1.31 (-13.47, 10.84) | 3.77 (-19.72, 27.26) |
+      +----------------+-----------------------+----------------------+
+      | *Interactions*                                                |
+      +----------------+-----------------------+----------------------+
+      | Days * grp (2) |                       | -1.01 (-5.35, 3.32)  |
+      +----------------+-----------------------+----------------------+
+      | Days * grp (3) |                       | -1.11 (-5.53, 3.31)  |
+      +----------------+-----------------------+----------------------+
+      | *Controls*                                                    |
+      +----------------+-----------------------+----------------------+
+      | Days           | 10.44 (8.84, 12.03)   | 11.23 (7.87, 14.60)  |
+      +----------------+-----------------------+----------------------+
+      |                |                       |                      |
+      +----------------+-----------------------+----------------------+
+      | Observations   | 180                   | 180                  |
+      +----------------+-----------------------+----------------------+ 
+
diff --git a/tests/testthat/test-compare_parameters.R b/tests/testthat/test-compare_parameters.R
index 4c37549ff..fe261d8d9 100644
--- a/tests/testthat/test-compare_parameters.R
+++ b/tests/testthat/test-compare_parameters.R
@@ -1,117 +1,211 @@
+skip_if_not_installed("withr")
+
 # make sure we have the correct interaction mark for tests
-options(parameters_interaction = "*")
-
-data(iris)
-m1 <- lm(Sepal.Length ~ Species, data = iris)
-m2 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)
-counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
-outcome <- gl(3, 1, 9)
-treatment <- gl(3, 3)
-m3 <- glm(counts ~ outcome + treatment, family = poisson())
-
-x <- compare_parameters(m1, m2, m3)
-test_that("compare_parameters, default", {
-  expect_identical(
-    colnames(x),
-    c(
-      "Parameter", "Component", "Effects", "Coefficient.m1", "SE.m1", "CI.m1",
-      "CI_low.m1", "CI_high.m1", "t.m1", "df_error.m1", "p.m1", "Coefficient.m2",
-      "SE.m2", "CI.m2", "CI_low.m2", "CI_high.m2", "t.m2", "df_error.m2",
-      "p.m2", "Log-Mean.m3", "SE.m3", "CI.m3", "CI_low.m3", "CI_high.m3",
-      "z.m3", "df_error.m3", "p.m3"
-    )
-  )
-  out <- capture.output(x)
-  expect_length(out, 14)
-  out <- format(x, select = "ci")
-  expect_identical(colnames(out), c("Parameter", "m1", "m2", "m3"))
-  expect_identical(
-    out$Parameter,
-    c(
-      "(Intercept)", "Species (versicolor)", "Species (virginica)",
-      "Petal Length", "Species (versicolor) * Petal Length",
-      "Species (virginica) * Petal Length", "outcome (2)", "outcome (3)",
-      "treatment (2)", "treatment (3)", NA, "Observations"
-    )
-  )
-})
-
-
-x <- compare_parameters(m1, m2, m3, select = "se_p2")
-test_that("compare_parameters, se_p2", {
-  expect_identical(
-    colnames(x),
-    c(
-      "Parameter", "Component", "Effects", "Coefficient.m1", "SE.m1", "CI.m1",
-      "CI_low.m1", "CI_high.m1", "t.m1", "df_error.m1", "p.m1", "Coefficient.m2",
-      "SE.m2", "CI.m2", "CI_low.m2", "CI_high.m2", "t.m2", "df_error.m2",
-      "p.m2", "Log-Mean.m3", "SE.m3", "CI.m3", "CI_low.m3", "CI_high.m3",
-      "z.m3", "df_error.m3", "p.m3"
-    )
-  )
-  out <- capture.output(x)
-  expect_length(out, 14)
-  out <- format(x, select = "se_p2")
-  expect_identical(
-    colnames(out),
-    c(
-      "Parameter", "Estimate (SE) (m1)", "p (m1)", "Estimate (SE) (m2)",
-      "p (m2)", "Estimate (SE) (m3)", "p (m3)"
-    )
-  )
-  expect_identical(
-    out$Parameter,
-    c(
-      "(Intercept)", "Species (versicolor)", "Species (virginica)",
-      "Petal Length", "Species (versicolor) * Petal Length",
-      "Species (virginica) * Petal Length", "outcome (2)", "outcome (3)",
-      "treatment (2)", "treatment (3)", NA, "Observations"
-    )
-  )
-})
-
-
-data(mtcars)
-m1 <- lm(mpg ~ wt, data = mtcars)
-m2 <- glm(vs ~ wt + cyl, data = mtcars, family = "binomial")
-
-test_that("compare_parameters, column name with escaping regex characters", {
-  out <- utils::capture.output(compare_parameters(m1, m2, column_names = c("linear model (m1)", "logistic reg. (m2)")))
-  expect_identical(out[1], "Parameter    |    linear model (m1) |   logistic reg. (m2)")
-})
-
-
-data(mtcars)
-m1 <- lm(mpg ~ hp, mtcars)
-m2 <- lm(mpg ~ hp, mtcars)
-test_that("compare_parameters, proper printing for CI=NULL #820", {
-  expect_snapshot(compare_parameters(m1, m2, ci = NULL))
-})
-
-
-skip_on_cran()
-
-
-test_that("compare_parameters, correct random effects", {
-  skip_if_not_installed("glmmTMB")
-  skip_if_not(getRversion() >= "4.0.0")
-
-  data("fish")
-  m0 <- glm(count ~ child + camper, data = fish, family = poisson())
-
-  m1 <- glmmTMB::glmmTMB(
-    count ~ child + camper + (1 | persons) + (1 | ID),
-    data = fish,
-    family = poisson()
-  )
-
-  m2 <- glmmTMB::glmmTMB(
-    count ~ child + camper + zg + (1 | ID),
-    ziformula = ~ child + (1 | persons),
-    data = fish,
-    family = glmmTMB::truncated_poisson()
-  )
-
-  cp <- compare_parameters(m0, m1, m2, effects = "all", component = "all")
-  expect_snapshot(cp)
-})
+withr::with_options(
+  list(parameters_interaction = "*"),
+  {
+    data(iris)
+    m1 <- lm(Sepal.Length ~ Species, data = iris)
+    m2 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)
+    counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
+    outcome <- gl(3, 1, 9)
+    treatment <- gl(3, 3)
+    m3 <- glm(counts ~ outcome + treatment, family = poisson())
+
+    x <- compare_parameters(m1, m2, m3)
+    test_that("compare_parameters, default", {
+      expect_identical(
+        colnames(x),
+        c(
+          "Parameter", "Component", "Effects", "Coefficient.m1", "SE.m1", "CI.m1",
+          "CI_low.m1", "CI_high.m1", "t.m1", "df_error.m1", "p.m1", "Coefficient.m2",
+          "SE.m2", "CI.m2", "CI_low.m2", "CI_high.m2", "t.m2", "df_error.m2",
+          "p.m2", "Log-Mean.m3", "SE.m3", "CI.m3", "CI_low.m3", "CI_high.m3",
+          "z.m3", "df_error.m3", "p.m3"
+        )
+      )
+      out <- capture.output(x)
+      expect_length(out, 14)
+      out <- format(x, select = "ci")
+      expect_identical(colnames(out), c("Parameter", "m1", "m2", "m3"))
+      expect_identical(
+        out$Parameter,
+        c(
+          "(Intercept)", "Species (versicolor)", "Species (virginica)",
+          "Petal Length", "Species (versicolor) * Petal Length",
+          "Species (virginica) * Petal Length", "outcome (2)", "outcome (3)",
+          "treatment (2)", "treatment (3)", NA, "Observations"
+        )
+      )
+    })
+
+
+    x <- compare_parameters(m1, m2, m3, select = "se_p2")
+    test_that("compare_parameters, se_p2", {
+      expect_identical(
+        colnames(x),
+        c(
+          "Parameter", "Component", "Effects", "Coefficient.m1", "SE.m1", "CI.m1",
+          "CI_low.m1", "CI_high.m1", "t.m1", "df_error.m1", "p.m1", "Coefficient.m2",
+          "SE.m2", "CI.m2", "CI_low.m2", "CI_high.m2", "t.m2", "df_error.m2",
+          "p.m2", "Log-Mean.m3", "SE.m3", "CI.m3", "CI_low.m3", "CI_high.m3",
+          "z.m3", "df_error.m3", "p.m3"
+        )
+      )
+      out <- capture.output(x)
+      expect_length(out, 14)
+      out <- format(x, select = "se_p2")
+      expect_identical(
+        colnames(out),
+        c(
+          "Parameter", "Estimate (SE) (m1)", "p (m1)", "Estimate (SE) (m2)",
+          "p (m2)", "Estimate (SE) (m3)", "p (m3)"
+        )
+      )
+      expect_identical(
+        out$Parameter,
+        c(
+          "(Intercept)", "Species (versicolor)", "Species (virginica)",
+          "Petal Length", "Species (versicolor) * Petal Length",
+          "Species (virginica) * Petal Length", "outcome (2)", "outcome (3)",
+          "treatment (2)", "treatment (3)", NA, "Observations"
+        )
+      )
+    })
+
+
+    data(mtcars)
+    m1 <- lm(mpg ~ wt, data = mtcars)
+    m2 <- glm(vs ~ wt + cyl, data = mtcars, family = "binomial")
+
+    test_that("compare_parameters, column name with escaping regex characters", {
+      out <- utils::capture.output(compare_parameters(m1, m2, column_names = c("linear model (m1)", "logistic reg. (m2)")))
+      expect_identical(out[1], "Parameter    |    linear model (m1) |   logistic reg. (m2)")
+    })
+
+
+    data(mtcars)
+    m1 <- lm(mpg ~ hp, mtcars)
+    m2 <- lm(mpg ~ hp, mtcars)
+    test_that("compare_parameters, proper printing for CI=NULL #820", {
+      expect_snapshot(compare_parameters(m1, m2, ci = NULL))
+    })
+
+
+    skip_on_cran()
+
+
+    test_that("compare_parameters, correct random effects", {
+      skip_if_not_installed("glmmTMB")
+      skip_if_not(getRversion() >= "4.0.0")
+
+      data("fish")
+      m0 <- glm(count ~ child + camper, data = fish, family = poisson())
+
+      m1 <- glmmTMB::glmmTMB(
+        count ~ child + camper + (1 | persons) + (1 | ID),
+        data = fish,
+        family = poisson()
+      )
+
+      m2 <- glmmTMB::glmmTMB(
+        count ~ child + camper + zg + (1 | ID),
+        ziformula = ~ child + (1 | persons),
+        data = fish,
+        family = glmmTMB::truncated_poisson()
+      )
+
+      cp <- compare_parameters(m0, m1, m2, effects = "all", component = "all")
+      expect_snapshot(cp)
+    })
+
+
+    test_that("compare_parameters, print_md", {
+      skip_if_not_installed("lme4")
+      data(sleepstudy, package = "lme4")
+      set.seed(1234)
+      sleepstudy$grp <- as.factor(sample.int(3, nrow(sleepstudy), replace = TRUE))
+      lm1 <- lme4::lmer(Reaction ~ Days + grp + (1 | Subject), data = sleepstudy)
+      lm2 <- lme4::lmer(Reaction ~ Days * grp + (1 | Subject), data = sleepstudy)
+
+      cp <- compare_parameters(lm1, lm2, select = "{estimate} ({ci})|{p}", drop = "^\\(Intercept")
+      out <- print_md(cp, groups = list(
+        Groups = c("grp (2)", "grp (3)"),
+        Interactions = c("Days * grp (2)", "Days * grp (3)"),
+        Controls = "Days"
+      ))
+      expect_snapshot(print(out))
+
+      cp <- compare_parameters(lm1, lm2, select = "{estimate} ({ci})|{p}", drop = "^\\(Intercept", effects = "all")
+      expect_snapshot(print_md(cp))
+
+      # error
+      expect_error(
+        print_md(cp, groups = list(
+          Groups = c("grp (2)", "grp (3)"),
+          Interactions = c("Days * grp (2)", "Days * grp (3)"),
+          Controls = "Days"
+        )),
+        regex = "Cannot combine"
+      )
+
+      # with reference level
+      cp <- compare_parameters(lm1, lm2, select = "{estimate} ({ci})|{p}", drop = "^\\(Intercept", include_reference = TRUE)
+      out <- print_md(cp, groups = list(
+        Groups = 2:4,
+        Interactions = 5:6,
+        Controls = 1
+      ))
+      expect_snapshot(print(out))
+
+      # with reference level
+      cp <- compare_parameters(lm1, lm2, drop = "^\\(Intercept", include_reference = TRUE)
+      out <- print_md(cp, groups = list(
+        Groups = 2:4,
+        Interactions = 5:6,
+        Controls = 1
+      ))
+      expect_snapshot(print(out))
+
+      # error
+      cp <- compare_parameters(lm1, lm2, select = "{estimate} ({ci})|{p}", drop = "^\\(Intercept")
+      expect_error(
+        print_md(cp, groups = list(
+          Groups = c("grp (2)", "grp (3)"),
+          Interactions = c("Days * grp (2)", "Days * grp (3)"),
+          Controls = "XDays"
+        )),
+        regex = "Some group indices"
+      )
+      expect_error(
+        print_md(cp, groups = list(
+          Groups = 1:2,
+          Interactions = 4:5,
+          Controls = 10
+        )),
+        regex = "Some group indices"
+      )
+
+      # output identical for both calls
+      cp1 <- compare_parameters(lm1, lm2, select = "{estimate} ({ci})|{p}", drop = "^\\(Intercept")
+      out1 <- capture.output(print_md(cp1, groups = list(
+        Groups = c("grp (2)", "grp (3)"),
+        Interactions = c("Days * grp (2)", "Days * grp (3)"),
+        Controls = "Days"
+      )))
+      cp2 <- compare_parameters(
+        lm1,
+        lm2,
+        select = "{estimate} ({ci})|{p}",
+        drop = "^\\(Intercept",
+        groups = list(
+          Groups = c("grp (2)", "grp (3)"),
+          Interactions = c("Days * grp (2)", "Days * grp (3)"),
+          Controls = "Days"
+        )
+      )
+      out2 <- capture.output(print_md(cp2))
+      expect_identical(out1, out2)
+    })
+  }
+)
diff --git a/tests/testthat/test-model_parameters_df.R b/tests/testthat/test-model_parameters_df.R
index 5ad60ce95..da8b44508 100644
--- a/tests/testthat/test-model_parameters_df.R
+++ b/tests/testthat/test-model_parameters_df.R
@@ -28,7 +28,7 @@ test_that("model_parameters.glm", {
 
 
 test_that("model_parameters.BBmm", {
-  skip_if_not_installed("PROreg")
+  skip_if_not_installed("PROreg", minimum_version = "1.3.0")
   set.seed(1234)
 
   # defining the parameters
@@ -63,17 +63,17 @@ test_that("model_parameters.BBmm", {
   }))
   params <- suppressWarnings(model_parameters(model))
   expect_equal(params$df_error, c(96, 96), tolerance = 1e-3)
-  expect_equal(params$CI_low, c(0.26363, -1.46645), tolerance = 1e-3)
-  expect_equal(params$p, c(0.00811, 0), tolerance = 1e-3)
+  expect_equal(params$CI_low, c(0.26366, -1.46628), tolerance = 1e-3)
+  expect_equal(params$p, c(0.00814, 0), tolerance = 1e-3)
 
   params <- suppressWarnings(model_parameters(model, ci_method = "normal"))
   expect_equal(params$df_error, c(Inf, Inf), tolerance = 1e-3)
-  expect_equal(params$CI_low, c(0.27359, -1.46136), tolerance = 1e-3)
-  expect_equal(params$p, c(0.00811, 0), tolerance = 1e-3)
+  expect_equal(params$CI_low, c(0.27313, -1.46119), tolerance = 1e-3)
+  expect_equal(params$p, c(0.00814, 0), tolerance = 1e-3)
 })
 
 test_that("model_parameters.BBreg", {
-  skip_if_not_installed("PROreg")
+  skip_if_not_installed("PROreg", minimum_version = "1.3.0")
   set.seed(18)
   # we simulate a covariate, fix the paramters of the beta-binomial
   # distribution and simulate a response variable.