diff --git a/DESCRIPTION b/DESCRIPTION
index d0d8d57..8f3e8c1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: dynamite
 Title: Bayesian Modeling and Causal Inference for Multivariate
     Longitudinal Data
-Version: 1.5.0
+Version: 1.5.1
 Authors@R: c(
     person("Santtu", "Tikka", email = "santtuth@gmail.com", 
            role = c("aut", "cre"), comment = c(ORCID = "0000-0003-4039-4342")),
diff --git a/NEWS.md b/NEWS.md
index 76f3d6a..28de9f7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,11 @@
+# dynamite 1.5.1
+
+  * The `type` argument of `coef()` and `plot()` has been replaced by `types` accepting multiple types simultaneously, similar to `as.data.table()` and `as.data.frame()`.
+  * The functions `plot_betas()`, `plot_deltas()`, `plot_nus()`, `plot_lambdas()` and `plot_psis()` have been deprecated and are now provided via the default plot method by selecting the appropriate `types`.
+  * A new argument `plot_type` has been added to control what type of plot will be drawn by the `plot()` method. The default value `"default"` draws the posterior means and posterior intervals of all parameters. The old functionality of drawing posterior densities and traceplots is provided by the option `"trace"`.
+  * The `plot()` method has gained the argument `n_params` to limit the amount of parameters drawn at once (per parameter type).
+  * Both time-varying and time-invariant parameters can now be plotted simultaneously.
+
 # dynamite 1.5.0
 
   * Estimation of dynamic multivariate panel models with multiple imputation is now available via the function `dynamice()` which uses the `mice` package.
@@ -8,7 +16,6 @@
   * Added a vignette on `get_code()` and `get_data()` functions and how they can be used to modify the generated Stan code and perform variational Bayes inference.
   * Contemporaneous dependencies are now allowed between different components of multivariate distributions, e.g., `obs(c(y, x) ~ x | 1, family = "mvgaussian")`.
   * Ordered probit and logit regressions are now available via `obs(., family = "cumulative", link = "probit")` and `obs(., family = "cumulative", link = "logit")`, respectively.
-  * The functions `plot_betas()`, `plot_deltas()`, `plot_nus()`, `plot_lambdas()` and `plot_psis()` have been deprecated and are now provided via the default plot method via the argument `plot_type` instead.
 
 # dynamite 1.4.11
 
diff --git a/R/as_data_frame.R b/R/as_data_frame.R
index e86b317..8fcf615 100644
--- a/R/as_data_frame.R
+++ b/R/as_data_frame.R
@@ -44,19 +44,19 @@
 #' @param x \[`dynamitefit`]\cr The model fit object.
 #' @param row.names Ignored.
 #' @param optional Ignored.
+#' @param types \[`character()`]\cr Type(s) of the parameters for which the
+#'   samples should be extracted. See details of possible values. Default is
+#'   all values listed in details except spline coefficients `omega`.
+#'   This argument is mutually exclusive with `parameters`.
 #' @param parameters \[`character()`]\cr Parameter(s) for which the samples
 #'   should be extracted. Possible options can be found with function
-#'   `get_parameter_names()`. Default is all parameters of specific type for all
-#'   responses.
+#'   `get_parameter_names()`. Default is all parameters of specific type for
+#'   all responses. This argument is mutually exclusive with `types`.
 #' @param responses \[`character()`]\cr Response(s) for which the samples
 #'   should be extracted. Possible options are elements of
 #'   `unique(x$priors$response)`, and the default is this entire vector.
 #'    Ignored if the argument `parameters` is supplied.
-#' @param types \[`character()`]\cr Type(s) of the parameters for which the
-#'   samples should be extracted. See details of possible values. Default is
-#'   all values listed in details except spline coefficients `omega`,
 #'   `omega_alpha`, and `omega_psi`. See also [dynamite::get_parameter_types()].
-#'    Ignored if the argument `parameters` is supplied.
 #' @param times \[`double()`]\cr Time point(s) to keep. If `NULL`
 #'   (the default), all time points are kept.
 #' @param groups \[`character()`] Group name(s) to keep. If `NULL`
@@ -107,23 +107,24 @@
 #' )
 #'
 as.data.frame.dynamitefit <- function(x, row.names = NULL, optional = FALSE,
-                                      parameters = NULL, responses = NULL,
-                                      types = NULL, times = NULL, groups = NULL,
+                                      types = NULL, parameters = NULL,
+                                      responses = NULL,
+                                      times = NULL, groups = NULL,
                                       summary = FALSE, probs = c(0.05, 0.95),
                                       include_fixed = TRUE, ...) {
   out <- as.data.table.dynamitefit(
-    x,
+    x = x,
     keep.rownames = FALSE,
-    row.names,
-    optional,
-    parameters,
-    responses,
-    types,
-    times,
-    groups,
-    summary,
-    probs,
-    include_fixed,
+    row.names = row.names,
+    optional = optional,
+    parameters = parameters,
+    responses = responses,
+    types = types,
+    times = times,
+    groups = groups,
+    summary = summary,
+    probs = probs,
+    include_fixed = include_fixed,
     ...
   )
   tibble::tibble(data.table::setDF(out))
diff --git a/R/as_data_table.R b/R/as_data_table.R
index 65bcce9..b31176d 100644
--- a/R/as_data_table.R
+++ b/R/as_data_table.R
@@ -23,8 +23,8 @@
 #'
 as.data.table.dynamitefit <- function(x, keep.rownames = FALSE,
                                       row.names = NULL, optional = FALSE,
-                                      parameters = NULL,
-                                      responses = NULL, types = NULL,
+                                      types = NULL, parameters = NULL,
+                                      responses = NULL,
                                       times = NULL, groups = NULL,
                                       summary = FALSE, probs = c(0.05, 0.95),
                                       include_fixed = TRUE, ...) {
@@ -115,12 +115,6 @@ as.data.table.dynamitefit <- function(x, keep.rownames = FALSE,
       "Model does not contain response variable{?s} {.var {responses[!z]}}."
     )
   }
-  all_types <- c(
-    "alpha", "beta", "delta", "tau", "tau_alpha", "xi",
-    "sigma_nu", "sigma", "phi", "nu", "lambda", "sigma_lambda",
-    "psi", "tau_psi", "corr", "corr_psi", "corr_nu",
-    "omega", "omega_alpha", "omega_psi", "cutpoints"
-  )
   if (is.null(types)) {
     types <- ifelse_(
       is.null(parameters),
@@ -534,7 +528,11 @@ as_data_table_omega <- function(x, draws, n_draws, response, category, ...) {
     names(get_channel(x, response)$J_varying)
   )
   k <- length(var_names)
-  data.table::data.table(
+  params_ord <- rep(
+    paste0(rep(var_names, each = D), "_d", rep(seq_len(D), k)),
+    each = n_draws
+  )
+  tmp <- data.table::data.table(
     parameter = rep(
       paste0(var_names, "_d", rep(seq_len(D), each = k)),
       each = n_draws
@@ -542,6 +540,7 @@ as_data_table_omega <- function(x, draws, n_draws, response, category, ...) {
     value = c(draws),
     category = category
   )
+  tmp[order(match(parameter, params_ord))]
 }
 
 #' @describeIn as_data_table_default Data Table for a "omega_alpha" Parameter
@@ -722,3 +721,67 @@ as_data_table_cutpoints <- function(x, draws, response,
   }
 }
 
+# Parameter types ---------------------------------------------------------
+
+all_types <- c(
+  "alpha",
+  "beta",
+  "corr",
+  "corr_nu",
+  "corr_psi",
+  "cutpoints",
+  "delta",
+  "lambda",
+  "nu",
+  "omega",
+  "omega_alpha",
+  "omega_psi",
+  "phi",
+  "psi",
+  "sigma_nu",
+  "sigma",
+  "sigma_lambda",
+  "tau",
+  "tau_alpha",
+  "tau_psi",
+  "xi"
+)
+
+fixed_types <- c(
+  "alpha",
+  "beta",
+  "corr",
+  "corr_nu",
+  "corr_psi",
+  "cutpoints",
+  "lambda",
+  "nu",
+  "omega",
+  "omega_alpha",
+  "omega_psi",
+  "phi",
+  "sigma",
+  "sigma_lambda",
+  "sigma_nu",
+  "tau",
+  "tau_alpha",
+  "tau_psi",
+  "xi"
+)
+
+varying_types <- c(
+  "alpha",
+  "cutpoints",
+  "delta",
+  "psi"
+)
+
+default_types <- c(
+  "alpha",
+  "beta",
+  "cutpoints",
+  "delta",
+  "lambda",
+  "nu",
+  "psi"
+)
diff --git a/R/coef.R b/R/coef.R
index fd8ea8b..05bf3b3 100644
--- a/R/coef.R
+++ b/R/coef.R
@@ -5,17 +5,6 @@
 #' @export
 #' @family output
 #' @param object \[`dynamitefit`]\cr The model fit object.
-#' @param type  \[`character(1)`]\cr Either `beta` (the default) for
-#'   time-invariant coefficients, `delta` for time-varying coefficients,
-#'   `nu` for random effects, `lambda` for factor loadings, or `psi` for
-#'   latent factor. Ignored if the argument `parameters` is supplied.
-#' @param include_alpha \[`logical(1)`]\cr If `TRUE` (default), extracts also
-#'   time-invariant intercept term alpha if time-invariant parameters beta are
-#'   extracted, and time-varying alpha if time-varying delta are extracted.
-#'   Ignored if the argument `parameters` is supplied.
-#' @param include_cutpoints \[`logical(1)`]\cr If `TRUE` (default), plots also
-#'   the cutpoints if such parameters exists in the model (either time-varying
-#'   or time-invariant).
 #' @param summary \[`logical(1)`]\cr If `TRUE` (default), returns posterior
 #'   mean, standard deviation, and posterior quantiles (as defined by the
 #'   `probs` argument) for all parameters. If `FALSE`, returns the
@@ -31,12 +20,9 @@
 #' deltas <- coef(gaussian_example_fit, type = "delta")
 #'
 coef.dynamitefit <- function(object,
-                             parameters = NULL,
-                             type = c("beta", "delta", "nu", "lambda", "psi"),
+                             types = NULL, parameters = NULL,
                              responses = NULL, times = NULL, groups = NULL,
-                             summary = TRUE, probs = c(0.05, 0.95),
-                             include_alpha = TRUE,
-                             include_cutpoints = TRUE, ...) {
+                             summary = TRUE, probs = c(0.05, 0.95), ...) {
   stopifnot_(
     !missing(object),
     "Argument {.arg object} is missing."
@@ -46,58 +32,24 @@ coef.dynamitefit <- function(object,
     "Argument {.arg object} must be a {.cls dynamitefit} object."
   )
   if (is.null(parameters)) {
-    type <- onlyif(is.character(type), tolower(type))
-    type <- try(
-      match.arg(type, c("beta", "delta", "nu", "lambda", "psi")),
-      silent = TRUE
-    )
-    stopifnot_(
-      !inherits(type, "try-error"),
-      "Argument {.arg type} must be either {.val beta}, {.val delta},
-      {.val nu}, {.val lambda}, or {.val psi}."
-    )
-    stopifnot_(
-      checkmate::test_flag(x = include_alpha),
-      "Argument {.arg include_alpha} must be single {.cls logical} value."
-    )
+    types <- onlyif(is.character(types), tolower(types))
+    invalid_types <- types[!types %in% all_types]
     stopifnot_(
-      checkmate::test_flag(x = include_cutpoints),
-      "Argument {.arg include_cutpoints} must be a single {.cls logical} value."
+      length(invalid_types) == 0L,
+      c(
+        "Argument {.arg types} must contain only valid types.",
+        `x` = "Invalid types found: {invalid_types}"
+      )
     )
   }
-  if (is.null(parameters) && type %in% c("beta", "delta")) {
-    types <- c(
-      type,
-      onlyif(include_alpha, "alpha"),
-      onlyif(include_cutpoints, "cutpoints")
-    )
-    out <- as.data.frame.dynamitefit(
-      object,
-      parameters = parameters,
-      responses = responses,
-      types = types,
-      times = times,
-      groups = groups,
-      summary = summary,
-      probs = probs
-    )
-    # remove extra alphas
-    if (identical(type, "delta")) {
-      out <- out[!is.na(out$time), ]
-    } else {
-      out <- out[is.na(out$time), ]
-    }
-  } else {
-    out <- as.data.frame.dynamitefit(
-      object,
-      parameters = parameters,
-      responses = responses,
-      types = type,
-      times = times,
-      groups = groups,
-      summary = summary,
-      probs = probs
-    )
-  }
-  out
+  as.data.frame.dynamitefit(
+    x = object,
+    types = types,
+    parameters = parameters,
+    responses = responses,
+    times = times,
+    groups = groups,
+    summary = summary,
+    probs = probs
+  )
 }
diff --git a/R/deprecated.R b/R/deprecated.R
index 75cdaf2..a3d2c9d 100644
--- a/R/deprecated.R
+++ b/R/deprecated.R
@@ -19,29 +19,29 @@
 #' @param ... Not used.
 #' @section Details:
 #'
-#'   * `plot_betas` is now called via `plot(., plot_type = "beta")`
-#'   * `plot_deltas` is now called via `plot(., plot_type = "delta")`
-#'   * `plot_nus` is now called via `plot(., plot_type = "nu")`
-#'   * `plot_lambdas` is now called via `plot(., plot_type = "lambda")`
-#'   * `plot_psis` is now called via `plot(., plot_type = "psi")`
+#'   * `plot_betas` is now called via `plot(., types = "beta")`
+#'   * `plot_deltas` is now called via `plot(., types = "delta")`
+#'   * `plot_nus` is now called via `plot(., types = "nu")`
+#'   * `plot_lambdas` is now called via `plot(., types = "lambda")`
+#'   * `plot_psis` is now called via `plot(., types = "psi")`
 #'
 plot_betas <- function(x, ...) {
   .Deprecated("plot")
-  plot(x, plot_type = "beta", ...)
+  plot(x, types = "beta", ...)
 }
 plot_deltas <- function(x, ...) {
   .Deprecated("plot")
-  plot(x, plot_type = "delta", ...)
+  plot(x, types = "delta", ...)
 }
 plot_nus <- function(x, ...) {
   .Deprecated("plot")
-  plot(x, plot_type = "nu", ...)
+  plot(x, types = "nu", ...)
 }
 plot_lambdas <- function(x, ...) {
   .Deprecated("plot")
-  plot(x, plot_type = "lambda", ...)
+  plot(x, types = "lambda", ...)
 }
 plot_psis <- function(x, ...) {
   .Deprecated("plot")
-  plot(x, plot_type = "psi", ...)
+  plot(x, types = "psi", ...)
 }
diff --git a/R/plot.R b/R/plot.R
index 0e6cdbb..114a2f4 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -238,24 +238,19 @@ plot_dynamiteformula_ggplot <- function(g, vertex_size, label_size) {
 #' @family plotting
 #' @param x \[`dynamitefit`]\cr The model fit object.
 #' @param plot_type \[`character(1)`]\cr What type of plot to draw? The default
-#'   is `"default"` which draws posterior densities and traceplots of the
-#'   selected parameters. Other options are `"beta"`, `"delta"`, `"nu"`,
-#'   `"lambda"` and `"psi"`. Option `"beta"` draws the posterior intervals
-#'   and means of the `beta` type parameters. Option `"delta"` draws the
-#'   posterior means and intervals as function of time for the `delta` type
-#'   parameters. Option `"nu"` draws the posterior means and intervals for the
-#'   group-specific random intercepts. Option `"lambda"` draws the posterior
-#'   means and intervals for the factor loading, and option `"psi"` draws
-#'   the means and posterior intervals for the latent factors as a function of
-#'   time.
-#' @param type \[`character(1)`]\cr Type of the parameter for which the plots
+#'   is `"default"` which draws posterior means and intervals of the parameters
+#'   selected by `types` or `parameters`. If both `"types"` and
+#'   `parameters` are `NULL`, all parameters are drawn up to the maximum
+#'   specified by `n_params`. Option `"trace"` instead draws posterior
+#'   densities and traceplots of the parameters.
+#' @param types \[`character(1)`]\cr Types of the parameter for which the plots
 #'   should be drawn. Possible options can be found with the function
 #'   [dynamite::get_parameter_types()]. Ignored if the argument `parameters`
-#'   is supplied or if `plot_type` is not `"default"`.
+#'   is supplied.
 #' @param parameters \[`charecter()`]\cr Parameter name(s) for which the plots
 #'   should be drawn. Possible options can be found with the function
-#'   [dynamite::get_parameter_names()]. The default is all parameters of a
-#'   specific type, limited by `n_params`.
+#'   [dynamite::get_parameter_names()]. The default is all parameters,
+#'   limited by `n_params`.
 #' @param responses \[`character()`]\cr Response(s) for which the plots should
 #'   be drawn. Possible options are `unique(x$priors$response)`. Default is
 #'   all responses. Ignored if the argument `parameters` is supplied.
@@ -264,26 +259,21 @@ plot_dynamiteformula_ggplot <- function(g, vertex_size, label_size) {
 #'   included, up to the maximum number of parameters specified by `n_params`
 #'   starting from the first non-fixed time point.
 #' @param groups \[`character(1)`]\cr Group name(s) for which the plots
-#'   should be drawn when `plot_type` is `"nu"` or `"lambda"`.
-#'   The default is all groups.
+#'   should be drawn for group-specific parameters.
 #' @param level \[`numeric(1)`]\cr Level for posterior intervals.
 #'   Default is 0.05, leading to 90% intervals.
 #' @param alpha \[`numeric(1)`]\cr Opacity level for `geom_ribbon`.
 #'   Default is 0.5.
 #' @param scales \[`character(1)`]\cr Should y-axis of the panels be `"fixed"`
 #'   (the default) or `"free"`? See [ggplot2::facet_wrap()].
-#' @param include_alpha \[`logical(1)`]\cr If `TRUE` (default), plots also
-#'   the alphas if such parameters exists in the model (either time-varying or
-#'   time-invariant)
-#' @param include_cutpoints \[`logical(1)`]\cr If `TRUE` (default), plots also
-#'   the cutpoints if such parameters exists in the model (either time-varying
-#'   or time-invariant).
-#' @param n_params \[`integer()`]\cr Maximum number of parameters to plot.
-#'   The default value is set by `plot_type`: 5 for `"default"`, 50 for
-#'   `"beta"`, `"nu"` and `"lambda"`, and 3 for `"delta"` and `"psi"`. This
-#'   argument is intended to prevent accidental plots that may be very large
-#'   and time consuming to plot. Please use the `parameters`, `times`, and
-#'   `groups` arguments to fine-tune which parameters to plot.
+#' @param n_params \[`integer()`]\cr A single value or a vector of length 2
+#'   specifying the maximum number of parameters to plot. If a single value
+#'   is provided, the same limit is used for all parameters. If a vector is
+#'   supplied, the first element defines the maximum number of time-invariant
+#'   parameters to plot and the second the maximum number of time-varying
+#'   parameters to plot. The defaults values are 50 for time-invariant
+#'   parameters and 3 for time-varying parameters. The default value is 5
+#'   for `plot_type == "trace"`.
 #' @param ... Not used..
 #' @return A `ggplot` object.
 #' @srrstats {BS6.1, RE6.0, RE6.1, BS6.2, BS6.3, BS6.5} Implements the `plot()`
@@ -293,13 +283,11 @@ plot_dynamiteformula_ggplot <- function(g, vertex_size, label_size) {
 #' data.table::setDTthreads(1) # For CRAN
 #' plot(gaussian_example_fit, type = "beta")
 #'
-plot.dynamitefit <- function(x, plot_type = c(
-                               "default", "beta", "delta", "nu", "lambda", "psi"
-                             ),
-                             type = NULL, parameters = NULL, responses = NULL,
-                             groups = NULL, times = NULL, level = 0.05,
-                             alpha = 0.5, scales = c("fixed", "free"),
-                             include_alpha = TRUE, include_cutpoints = TRUE,
+plot.dynamitefit <- function(x, plot_type = c("default", "trace"),
+                             types = NULL, parameters = NULL,
+                             responses = NULL, groups = NULL, times = NULL,
+                             level = 0.05, alpha = 0.5,
+                             scales = c("fixed", "free"),
                              n_params = NULL, ...) {
   stopifnot_(
     !missing(x),
@@ -309,23 +297,14 @@ plot.dynamitefit <- function(x, plot_type = c(
     is.dynamitefit(x),
     "Argument {.arg x} must be a {.cls dynamitefit} object."
   )
+  plot_type <- onlyif(is.character(plot_type), tolower(plot_type))
   plot_type <- try(
-    match.arg(
-      plot_type,
-      c("default", "beta", "delta", "nu", "lambda", "psi")
-    ),
+    match.arg(plot_type, c("default", "trace")),
     silent = TRUE
   )
   stopifnot_(
     !inherits(plot_type, "try-error"),
-    "Argument {.arg type} must one one of
-    {.val default}, {.val beta}, {.val delta}, {.val nu}, {.val lambda},
-    or {.val psi}."
-  )
-  stopifnot_(
-    !identical(plot_type, "default") || !is.null(parameters) || !is.null(type),
-    "Either {.arg parameters} or {.arg type} must be provided when
-     {.arg plot_type} is {.val default}."
+    "Argument {.arg type} must be either {.val default} or {.val trace}."
   )
   stopifnot_(
     checkmate::test_character(
@@ -337,12 +316,13 @@ plot.dynamitefit <- function(x, plot_type = c(
     "Argument {.arg parameters} must be a {.cls character} vector."
   )
   stopifnot_(
-    checkmate::test_string(
-      x = type,
-      na.ok = FALSE,
+    checkmate::test_character(
+      x = types,
+      any.missing = FALSE,
+      min.len = 1L,
       null.ok = TRUE
     ),
-    "Argument {.arg type} must be a single {.cls character} string."
+    "Argument {.arg types} must be a {.cls character} vector."
   )
   stopifnot_(
     checkmate::test_number(
@@ -370,48 +350,94 @@ plot.dynamitefit <- function(x, plot_type = c(
     !inherits(scales, "try-error"),
     "Argument {.arg scales} must be either {.val fixed} or {.val free}."
   )
-  switch(
-    plot_type,
-    default = plot_default(x, parameters, type, responses, times, n_params),
-    beta = plot_beta(
+  stopifnot_(
+    checkmate::test_integerish(
+      x = n_params,
+      lower = 1.0,
+      any.missing = FALSE,
+      null.ok = TRUE,
+      min.len = 1L,
+      max.len = 2L
+    ),
+    "Argument {.arg n_params} must a
+     positive integer vector of at most length 2."
+  )
+  n_params <- ifelse_(
+    !is.null(n_params) && length(n_params) == 1L,
+    c(n_params, n_params),
+    n_params
+  )
+  if (identical(plot_type, "trace")) {
+    plot_trace(
       x,
+      types,
       parameters,
       responses,
-      level,
-      include_alpha,
-      include_cutpoints,
+      times,
+      groups,
       n_params
-    ),
-    delta = plot_delta(
+    )
+  } else {
+    p_fixed <- plot_fixed(
       x,
+      types,
       parameters,
       responses,
+      times,
+      groups,
       level,
       alpha,
       scales,
-      include_alpha,
-      include_cutpoints,
+      n_params[1L]
+    )
+    p_varying <- plot_varying(
+      x,
+      types,
+      parameters,
+      responses,
       times,
-      n_params
-    ),
-    nu = plot_nu(x, parameters, responses, level, groups, n_params),
-    lambda = plot_lambda(x, responses, level, groups, n_params),
-    psi = plot_psi(x, responses, level, alpha, scales, n_params)
-  )
+      groups,
+      level,
+      alpha,
+      scales,
+      n_params[2L]
+    )
+    stopifnot_(
+      !is.null(p_fixed) || !is.null(p_varying),
+      "The model does not contain any of the specified parameters."
+    )
+    if (is.null(p_fixed)) {
+      p_varying
+    } else if (is.null(p_varying)) {
+      p_fixed
+    } else {
+      patchwork::wrap_plots(
+        p_fixed,
+        patchwork::free(p_varying),
+        heights = c(0.33, 0.66),
+        ncol = 1L,
+        nrow = 2L,
+        byrow = TRUE,
+        guides = "keep"
+      )
+    }
+  }
 }
 
 #' Traceplots and Density Plots for a `dynamitefit` Object
 #'
 #' @inheritParams plot.dynamitefit
 #' @noRd
-plot_default <- function(x, parameters, type, responses, times, n_params) {
+plot_trace <- function(x, types, parameters, responses,
+                       times, groups, n_params) {
   out <- suppressWarnings(
     as_draws_df.dynamitefit(
       x,
       parameters = parameters,
       responses = responses,
-      types = type,
-      times = times
+      types = types,
+      times = times,
+      groups = groups
     )
   )
   # avoid NSE notes from R CMD check
@@ -448,47 +474,72 @@ plot_default <- function(x, parameters, type, responses, times, n_params) {
   )
 }
 
-#' Plot Time-invariant Regression Coefficients of a Dynamite Model
+#' Plot Time-invariant Parameters of a Dynamite Model
 #'
 #' @inheritParams plot.dynamitefit
 #' @noRd
-plot_beta <- function(x, parameters, responses, level,
-                      include_alpha, include_cutpoints, n_params) {
-
+plot_fixed <- function(x, types, parameters, responses,
+                       times, groups, level, alpha, scales, n_params) {
   if (!is.null(parameters)) {
-    beta_names <- get_parameter_names(x, types = "beta")
-    found_pars <- parameters %in% beta_names
+    fixed_names <- get_parameter_names(x, types = fixed_types)
+    found_pars <- parameters %in% fixed_names
     stopifnot_(
       all(found_pars),
       c(
        "Parameter{?s} {.var {parameters[!found_pars]}} not found or
         {?it is/they are} of wrong type:",
-        `i` = 'Use {.fun get_parameter_names} with {.arg types = "beta"} to
+        `i` = 'Use {.fun get_parameter_names} to
         check suitable parameter names.'
       )
     )
   }
+  types <- ifelse_(
+    is.null(types),
+    default_types[default_types %in% fixed_types],
+    types[types %in% fixed_types]
+  )
+  if (length(types) == 0L) {
+    return(NULL)
+  }
   coefs <- coef.dynamitefit(
     x,
+    types = types,
     parameters = parameters,
     responses = responses,
-    type = "beta",
-    probs = c(level, 1 - level),
-    include_alpha = include_alpha,
-    include_cutpoints = include_cutpoints
-  )
-  stopifnot_(
-    nrow(coefs) > 0L,
-    "The model does not contain fixed coefficients beta."
+    times = times,
+    groups = groups,
+    probs = c(level, 1 - level)
   )
+  coefs <- coefs[is.na(coefs$time), ]
+  if (nrow(coefs) == 0L) {
+    return(NULL)
+  }
   coefs <- filter_params(coefs, n_params, 50)
-  title <- paste0(
-    "Posterior mean and ",
-    100 * (1 - 2 * level),
-    "% intervals of the time-invariant coefficients"
+  n_coefs <- nrow(coefs)
+  coefs$parameter <- glue::glue(
+    "{coefs$parameter}_{coefs$category}_{coefs$group}"
+  )
+  # remove NAs from parameters which are not category specific
+  coefs$parameter <- gsub("_NA", "", coefs$parameter)
+  coefs$parameter <- factor(coefs$parameter, levels = coefs$parameter)
+  title_spec <- "time-invariant parameters"
+  if (n_unique(types) == 1L) {
+    title_spec <- switch(
+      types[1L],
+      alpha = "time-invariant intercepts",
+      beta = "time-invariant regression coefficients",
+      cutpoints = "time-invariant cutpoints",
+      nu = "random intercepts",
+      lambda = "latent factor loadings",
+      "time-invariant parameters"
+    )
+  }
+  title <- glue::glue(
+    "Posterior means and {100 * (1 - 2 * level)} ",
+    "% intervals of the {title_spec}"
   )
   # avoid NSE notes from R CMD check
-  time <- mean <- category <- parameter <- NULL
+  time <- mean <- category <- parameter <- type <- NULL
   if (any(!is.na(coefs$category))) {
     p <- ggplot2::ggplot(
       coefs,
@@ -503,6 +554,7 @@ plot_beta <- function(x, parameters, responses, level,
     p <- ggplot2::ggplot(coefs, ggplot2::aes(mean, parameter))
   }
   p +
+    ggplot2::facet_wrap(~type, scales = "free") +
     ggplot2::geom_pointrange(
       ggplot2::aes(
         xmin = !!rlang::sym(paste0("q", 100 * level)),
@@ -510,49 +562,67 @@ plot_beta <- function(x, parameters, responses, level,
       ),
       position = ggplot2::position_dodge(0.5)
     ) +
-    ggplot2::labs(title = title, x = "Value", y = "Parameter")
+    ggplot2::labs(x = "Value", y = "Parameter") +
+    ggplot2::ggtitle(title)
 }
 
 
-#' Plot Time-varying Regression Coefficients of a Dynamite Model
+#' Plot Time-varying Parameters of a Dynamite Model
 #'
 #' @inheritParams plot.dynamitefit
 #' @noRd
-plot_delta <- function(x, parameters, responses, level, alpha,
-                       scales, include_alpha, include_cutpoints,
-                       times, n_params) {
+plot_varying <- function(x, types, parameters, responses,
+                         times, groups, level, alpha, scales, n_params) {
   if (!is.null(parameters)) {
-    delta_names <- get_parameter_names(x, types = "delta")
-    found_pars <- parameters %in% delta_names
+    varying_names <- get_parameter_names(x, types = varying_types)
+    found_pars <- parameters %in% varying_names
     stopifnot_(
       all(found_pars),
       c(
         "Parameter{?s} {.var {parameters[!found_pars]}} not found or
         {?it is/they are} of wrong type:",
-        `i` = 'Use {.fun get_parameter_names} with {.arg types = "delta"} to
+        `i` = 'Use {.fun get_parameter_names} to
         check suitable parameter names.'
       )
     )
   }
+  types <- ifelse_(
+    is.null(types),
+    default_types[default_types %in% varying_types],
+    types[types %in% varying_types]
+  )
+  if (length(types) == 0L) {
+    return(NULL)
+  }
   coefs <- coef.dynamitefit(
     x,
+    types = types,
     parameters = parameters,
     responses = responses,
-    type = "delta",
     times = times,
-    probs = c(level, 1 - level),
-    include_alpha = include_alpha,
-    include_cutpoints = include_cutpoints
-  )
-  stopifnot_(
-    nrow(coefs) > 0L,
-    "The model does not contain varying coefficients delta."
+    groups = groups,
+    probs = c(level, 1 - level)
   )
-  coefs <- filter_params(coefs, n_params, 3L)
-  title <- paste0(
-    "Posterior mean and ",
-    100 * (1 - 2 * level),
-    "% intervals of the time-varying coefficients"
+  coefs <- coefs[!is.na(coefs$time), ]
+  if (nrow(coefs) == 0L) {
+    return(NULL)
+  }
+  coefs <- filter_params(coefs, n_params, 3)
+  n_coefs <- nrow(coefs)
+  title_spec <- "time-varying parameters"
+  if (n_unique(types) == 1L) {
+    title_spec <- switch(
+      types[1L],
+      alpha = "time-varying intercepts",
+      cutpoints = "time-invariant cutpoints",
+      delta = "time-invariant regression coefficients",
+      psi = "latent factors",
+      "time-varying parameters"
+    )
+  }
+  title <- glue::glue(
+    "Posterior means and {100 * (1 - 2 * level)} %",
+    "intervals of the {title_spec}"
   )
   # avoid NSE notes from R CMD check
   time <- mean <- category <- parameter <- NULL
@@ -581,149 +651,8 @@ plot_delta <- function(x, parameters, responses, level, alpha,
     ggplot2::geom_line(na.rm = TRUE) +
     ggplot2::scale_x_continuous(limits = range(coefs$time)) +
     ggplot2::facet_wrap("parameter", scales = scales) +
-    ggplot2::labs(title = title, x = "Time", y = "Value")
-}
-
-#' Plot Random effects of a Dynamite Model
-#'
-#' @inheritParams plot.dynamitefit
-#' @noRd
-plot_nu <- function(x, parameters, responses, level, groups, n_params) {
-  if (!is.null(parameters)) {
-    nu_names <- get_parameter_names(x, types = "nu")
-    found_pars <- parameters %in% nu_names
-    stopifnot_(
-      all(found_pars),
-      c(
-        "Parameter{?s} {.var {parameters[!found_pars]}} not found or
-        {?it is/they are} of wrong type:",
-        `i` = 'Use {.fun get_parameter_names} with {.arg types = "nu"} to
-        check suitable parameter names.'
-      )
-    )
-  }
-  coefs <- try(
-    coef.dynamitefit(
-      x,
-      parameters = parameters,
-      type = "nu",
-      groups = groups,
-      responses = responses,
-      probs = c(level, 1 - level)
-    ),
-    silent = TRUE
-  )
-  stopifnot_(
-    !inherits(coefs, "try-error"),
-    "The model does not contain random effects nu."
-  )
-  coefs <- filter_params(coefs, n_params, 50L)
-  # avoid NSE notes from R CMD check
-  mean <- parameter <- NULL
-  coefs$parameter <- glue::glue(
-    "{coefs$parameter}_{coefs$category}_{coefs$group}"
-  )
-  # remove NAs from parameters which are not category specific
-  coefs$parameter <- gsub("_NA", "", coefs$parameter)
-  coefs$parameter <- factor(coefs$parameter, levels = coefs$parameter)
-  title <- paste0(
-    "Posterior mean and ",
-    100 * (1 - 2 * level),
-    "% intervals of the random intercepts"
-  )
-  ggplot2::ggplot(coefs, ggplot2::aes(mean, parameter)) +
-    ggplot2::geom_pointrange(
-      ggplot2::aes(
-        xmin = !!rlang::sym(paste0("q", 100 * level)),
-        xmax = !!rlang::sym(paste0("q", 100 * (1 - level)))
-      )
-    ) +
-    ggplot2::labs(title = title, x = "Value", y = "Parameter")
-}
-
-#' Plot Factor Loadings of a Dynamite Model
-#'
-#' @inheritParams plot.dynamitefit
-#' @noRd
-plot_lambda <- function(x, responses, level, groups, n_params) {
-  coefs <- try(
-    coef.dynamitefit(
-      x,
-      type = "lambda",
-      responses = responses,
-      probs = c(level, 1 - level)
-    ),
-    silent = TRUE
-  )
-  stopifnot_(
-    !inherits(coefs, "try-error"),
-    "The model does not contain latent factor psi."
-  )
-  coefs <- filter_params(coefs, n_params, 50L)
-  # avoid NSE notes from R CMD check
-  time <- mean <- parameter <- NULL
-  coefs$parameter <- glue::glue("{coefs$parameter}_{coefs$group}")
-  coefs$parameter <- factor(coefs$parameter, levels = coefs$parameter)
-  title <- paste0(
-    "Posterior mean and ",
-    100 * (1 - 2 * level),
-    "% intervals of the latent factor loadings"
-  )
-  ggplot2::ggplot(coefs, ggplot2::aes(mean, parameter)) +
-    ggplot2::geom_pointrange(ggplot2::aes(
-      xmin = !!rlang::sym(paste0("q", 100 * level)),
-      xmax = !!rlang::sym(paste0("q", 100 * (1 - level)))
-    )) +
-    ggplot2::labs(title = title, x = "Value", y = "Parameter")
-}
-
-#' Plot Latent Factors of a Dynamite Model
-#'
-#' @inheritParams plot.dynamitefit
-#' @noRd
-plot_psi <- function(x, responses, level, alpha, scales, times, n_params) {
-  coefs <- coef.dynamitefit(
-    x,
-    type = "psi",
-    responses = responses,
-    probs = c(level, 1 - level)
-  )
-  stopifnot_(
-    nrow(coefs) > 0L,
-    "The model does not contain latent factor psi."
-  )
-  coefs <- filter_params(coefs, n_params, 3L)
-  title <- paste0(
-    "Posterior mean and ",
-    100 * (1 - 2 * level),
-    "% intervals of the latent factors"
-  )
-  # avoid NSE notes from R CMD check
-  time <- mean <- category <- NULL
-  if (any(!is.na(coefs$category))) {
-    p <- ggplot2::ggplot(
-      coefs,
-      ggplot2::aes(
-        time,
-        mean,
-        colour = category,
-        fill = category
-      )
-    )
-  } else {
-    p <- ggplot2::ggplot(coefs, ggplot2::aes(time, mean))
-  }
-  p +
-    ggplot2::geom_ribbon(
-      ggplot2::aes(
-        ymin = !!rlang::sym(paste0("q", 100 * level)),
-        ymax = !!rlang::sym(paste0("q", 100 * (1 - level)))
-      ),
-      alpha = alpha
-    ) +
-    ggplot2::geom_line() +
-    ggplot2::facet_wrap("parameter", scales = scales) +
-    ggplot2::labs(title = title, x = "Time", y = "Value")
+    ggplot2::labs(x = "Time", y = "Value") +
+    ggplot2::ggtitle(title)
 }
 
 #' Select only specific number of parameters for plotting
@@ -747,23 +676,39 @@ filter_params <- function(x, n_params, n_params_default) {
     n_vars <- length(vars)
     keep_params <- logical(n_vars)
     keep_params[seq_len(min(n_vars, n_params))] <- TRUE
+    onlyif(
+      !n_params_set && !all(keep_params),
+      warning_(c(
+        "Number of parameters to be plotted ({n_u_params}) exceeds the
+       maximum number of parameters ({n_params}). The remaining parameters will
+       not be plotted.",
+       `i` = "Please increase {.arg n_params} to plot more parameters."
+      ))
+    )
   } else {
     params <- glue::glue("{x$parameter}_{x$group}")
     params <- gsub("_NA", "", params)
-    u_params <- unique(params)
-    n_u_params <- length(u_params)
-    n_params <- min(n_params, n_u_params)
-    keep_params <- params %in% u_params[seq_len(n_params)]
+    keep_params <- logical(nrow(x))
+    for (type in unique(x$type)) {
+      type_params <- params[x$type == type]
+      u_params <- unique(type_params)
+      n_u_params <- length(u_params)
+      max_params <- min(n_params, n_u_params)
+      keep_type_params <- type_params %in% u_params[seq_len(n_params)]
+      keep_params <- keep_params |
+        params %in% u_params[seq_len(n_params)]
+      onlyif(
+        !n_params_set && !all(keep_type_params),
+        warning_(c(
+          "Number of parameters to be plotted ({n_u_params}) exceeds the
+         maximum number of parameters ({max_params}) for parameters
+         of type {.var type}. The remaining parameters of this type will
+         not be plotted.",
+         `i` = "Please increase {.arg n_params} to plot more parameters."
+        ))
+      )
+    }
   }
-  onlyif(
-    !n_params_set && !all(keep_params),
-    warning_(c(
-      "Number of parameters to be plotted ({n_u_params}) exceeds the
-       maximum number of parameters ({n_params}). The remaining parameters will
-       not be plotted.",
-      `i` = "Please increase {.arg n_params} to plot more parameters."
-    ))
-  )
   ifelse_(
     is_draws,
     vars[keep_params],
diff --git a/README.Rmd b/README.Rmd
index c28b55a..913b066 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -34,7 +34,7 @@ The `dynamite` [R](https://www.r-project.org/) package provides an easy-to-use i
 
 * Support for regular time-invariant effects, group-level random effects, and time-varying effects modeled via Bayesian P-splines.
 * Joint modeling of multiple measurements per individual (multiple channels) based directly on the assumed data generating process. Individual channels can be univariate or multivariate.
-* Support for various distributions: Currently Gaussian, Multivariate Gaussian, Student t, Categorical, Multinomial, Poisson, Bernoulli, Binomial, Negative Binomial, Gamma, Exponential, and Beta distributions are available, and these can be combined arbitrarily in multichannel models.
+* Support for various distributions: Currently Gaussian, Multivariate Gaussian, Student t, Categorical, Ordered, Multinomial, Poisson, Bernoulli, Binomial, Negative Binomial, Gamma, Exponential, and Beta distributions are available, and these can be combined arbitrarily in multichannel models.
 * Allows evaluating realistic long-term counterfactual predictions which take into account the dynamic structure of the model by posterior predictive distribution simulation.
 * Transparent quantification of parameter and predictive uncertainty due to a fully Bayesian approach.
 * User-friendly and efficient R interface with state-of-the-art estimation via Stan. Both `rstan` and `cmdstanr` backends are supported.
@@ -86,17 +86,17 @@ print(gaussian_example_fit)
 
 Posterior estimates of time-varying effects:
 ```{r, fig.width = 9, fig.height = 4}
-plot(gaussian_example_fit, plot_type = "delta", scales = "free")
+plot(gaussian_example_fit, types = c("alpha", "delta"), scales = "free")
 ```
 
-And group-specific intercepts:
+And group-specific intercepts (for first 10 groups):
 ```{r, fig.width = 9, fig.height = 4}
-plot(gaussian_example_fit, plot_type = "nu", groups = 1:10)
+plot(gaussian_example_fit, types = "nu", groups = 1:10)
 ```
 
 Traceplots and density plots for time-invariant parameters:
 ```{r, fig.width = 9, fig.height = 4}
-plot(gaussian_example_fit, type = "beta")
+plot(gaussian_example_fit, plot_type = "trace", types = "beta")
 ```
 
 Posterior predictive samples for the first 4 groups (samples based on the posterior distribution of model parameters and observed data on first time point):
diff --git a/README.md b/README.md
index 365431b..176b43a 100644
--- a/README.md
+++ b/README.md
@@ -32,8 +32,8 @@ methodology from many other approaches are:
   channels) based directly on the assumed data generating process.
   Individual channels can be univariate or multivariate.
 - Support for various distributions: Currently Gaussian, Multivariate
-  Gaussian, Student t, Categorical, Multinomial, Poisson, Bernoulli,
-  Binomial, Negative Binomial, Gamma, Exponential, and Beta
+  Gaussian, Student t, Categorical, Ordered, Multinomial, Poisson,
+  Bernoulli, Binomial, Negative Binomial, Gamma, Exponential, and Beta
   distributions are available, and these can be combined arbitrarily in
   multichannel models.
 - Allows evaluating realistic long-term counterfactual predictions which
@@ -105,8 +105,8 @@ print(gaussian_example_fit)
 #> 
 #> Elapsed time (seconds):
 #>         warmup sample
-#> chain:1  5.511  3.354
-#> chain:2  5.504  3.462
+#> chain:1  5.801  3.542
+#> chain:2  5.658  3.544
 #> 
 #> Summary statistics of the time- and group-invariant parameters:
 #> # A tibble: 6 × 10
@@ -123,15 +123,15 @@ print(gaussian_example_fit)
 Posterior estimates of time-varying effects:
 
 ``` r
-plot(gaussian_example_fit, plot_type = "delta", scales = "free")
+plot(gaussian_example_fit, types = c("alpha", "delta"), scales = "free")
 ```
 
 <img src="man/figures/README-unnamed-chunk-7-1.png" style="display: block; margin: auto;" />
 
-And group-specific intercepts:
+And group-specific intercepts (for first 10 groups):
 
 ``` r
-plot(gaussian_example_fit, plot_type = "nu", groups = 1:10)
+plot(gaussian_example_fit, types = "nu", groups = 1:10)
 ```
 
 <img src="man/figures/README-unnamed-chunk-8-1.png" style="display: block; margin: auto;" />
@@ -139,7 +139,7 @@ plot(gaussian_example_fit, plot_type = "nu", groups = 1:10)
 Traceplots and density plots for time-invariant parameters:
 
 ``` r
-plot(gaussian_example_fit, type = "beta")
+plot(gaussian_example_fit, plot_type = "trace", types = "beta")
 ```
 
 <img src="man/figures/README-unnamed-chunk-9-1.png" style="display: block; margin: auto;" />
diff --git a/man/as.data.frame.dynamitefit.Rd b/man/as.data.frame.dynamitefit.Rd
index 30ce9bf..a054306 100644
--- a/man/as.data.frame.dynamitefit.Rd
+++ b/man/as.data.frame.dynamitefit.Rd
@@ -8,9 +8,9 @@
   x,
   row.names = NULL,
   optional = FALSE,
+  types = NULL,
   parameters = NULL,
   responses = NULL,
-  types = NULL,
   times = NULL,
   groups = NULL,
   summary = FALSE,
@@ -26,21 +26,21 @@
 
 \item{optional}{Ignored.}
 
+\item{types}{[\code{character()}]\cr Type(s) of the parameters for which the
+samples should be extracted. See details of possible values. Default is
+all values listed in details except spline coefficients \code{omega}.
+This argument is mutually exclusive with \code{parameters}.}
+
 \item{parameters}{[\code{character()}]\cr Parameter(s) for which the samples
 should be extracted. Possible options can be found with function
-\code{get_parameter_names()}. Default is all parameters of specific type for all
-responses.}
+\code{get_parameter_names()}. Default is all parameters of specific type for
+all responses. This argument is mutually exclusive with \code{types}.}
 
 \item{responses}{[\code{character()}]\cr Response(s) for which the samples
 should be extracted. Possible options are elements of
 \code{unique(x$priors$response)}, and the default is this entire vector.
-Ignored if the argument \code{parameters} is supplied.}
-
-\item{types}{[\code{character()}]\cr Type(s) of the parameters for which the
-samples should be extracted. See details of possible values. Default is
-all values listed in details except spline coefficients \code{omega},
-\code{omega_alpha}, and \code{omega_psi}. See also \code{\link[=get_parameter_types]{get_parameter_types()}}.
-Ignored if the argument \code{parameters} is supplied.}
+Ignored if the argument \code{parameters} is supplied.
+\code{omega_alpha}, and \code{omega_psi}. See also \code{\link[=get_parameter_types]{get_parameter_types()}}.}
 
 \item{times}{[\code{double()}]\cr Time point(s) to keep. If \code{NULL}
 (the default), all time points are kept.}
diff --git a/man/as.data.table.dynamitefit.Rd b/man/as.data.table.dynamitefit.Rd
index 61e8e2c..af9f52f 100644
--- a/man/as.data.table.dynamitefit.Rd
+++ b/man/as.data.table.dynamitefit.Rd
@@ -10,9 +10,9 @@
   keep.rownames = FALSE,
   row.names = NULL,
   optional = FALSE,
+  types = NULL,
   parameters = NULL,
   responses = NULL,
-  types = NULL,
   times = NULL,
   groups = NULL,
   summary = FALSE,
@@ -30,21 +30,21 @@
 
 \item{optional}{Ignored.}
 
+\item{types}{[\code{character()}]\cr Type(s) of the parameters for which the
+samples should be extracted. See details of possible values. Default is
+all values listed in details except spline coefficients \code{omega}.
+This argument is mutually exclusive with \code{parameters}.}
+
 \item{parameters}{[\code{character()}]\cr Parameter(s) for which the samples
 should be extracted. Possible options can be found with function
-\code{get_parameter_names()}. Default is all parameters of specific type for all
-responses.}
+\code{get_parameter_names()}. Default is all parameters of specific type for
+all responses. This argument is mutually exclusive with \code{types}.}
 
 \item{responses}{[\code{character()}]\cr Response(s) for which the samples
 should be extracted. Possible options are elements of
 \code{unique(x$priors$response)}, and the default is this entire vector.
-Ignored if the argument \code{parameters} is supplied.}
-
-\item{types}{[\code{character()}]\cr Type(s) of the parameters for which the
-samples should be extracted. See details of possible values. Default is
-all values listed in details except spline coefficients \code{omega},
-\code{omega_alpha}, and \code{omega_psi}. See also \code{\link[=get_parameter_types]{get_parameter_types()}}.
-Ignored if the argument \code{parameters} is supplied.}
+Ignored if the argument \code{parameters} is supplied.
+\code{omega_alpha}, and \code{omega_psi}. See also \code{\link[=get_parameter_types]{get_parameter_types()}}.}
 
 \item{times}{[\code{double()}]\cr Time point(s) to keep. If \code{NULL}
 (the default), all time points are kept.}
diff --git a/man/as_draws-dynamitefit.Rd b/man/as_draws-dynamitefit.Rd
index 3aa64fb..c2a9c22 100644
--- a/man/as_draws-dynamitefit.Rd
+++ b/man/as_draws-dynamitefit.Rd
@@ -24,19 +24,19 @@
 
 \item{parameters}{[\code{character()}]\cr Parameter(s) for which the samples
 should be extracted. Possible options can be found with function
-\code{get_parameter_names()}. Default is all parameters of specific type for all
-responses.}
+\code{get_parameter_names()}. Default is all parameters of specific type for
+all responses. This argument is mutually exclusive with \code{types}.}
 
 \item{responses}{[\code{character()}]\cr Response(s) for which the samples
 should be extracted. Possible options are elements of
 \code{unique(x$priors$response)}, and the default is this entire vector.
-Ignored if the argument \code{parameters} is supplied.}
+Ignored if the argument \code{parameters} is supplied.
+\code{omega_alpha}, and \code{omega_psi}. See also \code{\link[=get_parameter_types]{get_parameter_types()}}.}
 
 \item{types}{[\code{character()}]\cr Type(s) of the parameters for which the
 samples should be extracted. See details of possible values. Default is
-all values listed in details except spline coefficients \code{omega},
-\code{omega_alpha}, and \code{omega_psi}. See also \code{\link[=get_parameter_types]{get_parameter_types()}}.
-Ignored if the argument \code{parameters} is supplied.}
+all values listed in details except spline coefficients \code{omega}.
+This argument is mutually exclusive with \code{parameters}.}
 
 \item{times}{[\code{double()}]\cr Time point(s) to keep. If \code{NULL}
 (the default), all time points are kept.}
diff --git a/man/coef.dynamitefit.Rd b/man/coef.dynamitefit.Rd
index b79b2ac..dc8ceff 100644
--- a/man/coef.dynamitefit.Rd
+++ b/man/coef.dynamitefit.Rd
@@ -6,35 +6,34 @@
 \usage{
 \method{coef}{dynamitefit}(
   object,
+  types = NULL,
   parameters = NULL,
-  type = c("beta", "delta", "nu", "lambda", "psi"),
   responses = NULL,
   times = NULL,
   groups = NULL,
   summary = TRUE,
   probs = c(0.05, 0.95),
-  include_alpha = TRUE,
-  include_cutpoints = TRUE,
   ...
 )
 }
 \arguments{
 \item{object}{[\code{dynamitefit}]\cr The model fit object.}
 
+\item{types}{[\code{character()}]\cr Type(s) of the parameters for which the
+samples should be extracted. See details of possible values. Default is
+all values listed in details except spline coefficients \code{omega}.
+This argument is mutually exclusive with \code{parameters}.}
+
 \item{parameters}{[\code{character()}]\cr Parameter(s) for which the samples
 should be extracted. Possible options can be found with function
-\code{get_parameter_names()}. Default is all parameters of specific type for all
-responses.}
-
-\item{type}{[\code{character(1)}]\cr Either \code{beta} (the default) for
-time-invariant coefficients, \code{delta} for time-varying coefficients,
-\code{nu} for random effects, \code{lambda} for factor loadings, or \code{psi} for
-latent factor. Ignored if the argument \code{parameters} is supplied.}
+\code{get_parameter_names()}. Default is all parameters of specific type for
+all responses. This argument is mutually exclusive with \code{types}.}
 
 \item{responses}{[\code{character()}]\cr Response(s) for which the samples
 should be extracted. Possible options are elements of
 \code{unique(x$priors$response)}, and the default is this entire vector.
-Ignored if the argument \code{parameters} is supplied.}
+Ignored if the argument \code{parameters} is supplied.
+\code{omega_alpha}, and \code{omega_psi}. See also \code{\link[=get_parameter_types]{get_parameter_types()}}.}
 
 \item{times}{[\code{double()}]\cr Time point(s) to keep. If \code{NULL}
 (the default), all time points are kept.}
@@ -50,15 +49,6 @@ posterior samples instead.}
 \item{probs}{[\code{numeric()}]\cr Quantiles of interest. Default is
 \code{c(0.05, 0.95)}.}
 
-\item{include_alpha}{[\code{logical(1)}]\cr If \code{TRUE} (default), extracts also
-time-invariant intercept term alpha if time-invariant parameters beta are
-extracted, and time-varying alpha if time-varying delta are extracted.
-Ignored if the argument \code{parameters} is supplied.}
-
-\item{include_cutpoints}{[\code{logical(1)}]\cr If \code{TRUE} (default), plots also
-the cutpoints if such parameters exists in the model (either time-varying
-or time-invariant).}
-
 \item{...}{Ignored.}
 }
 \value{
diff --git a/man/dynamite-deprecated.Rd b/man/dynamite-deprecated.Rd
index f4e97fe..0e63fa2 100644
--- a/man/dynamite-deprecated.Rd
+++ b/man/dynamite-deprecated.Rd
@@ -30,11 +30,11 @@ package. They will eventually be completely removed.
 \section{Details}{
 
 \itemize{
-\item \code{plot_betas} is now called via \code{plot(., plot_type = "beta")}
-\item \code{plot_deltas} is now called via \code{plot(., plot_type = "delta")}
-\item \code{plot_nus} is now called via \code{plot(., plot_type = "nu")}
-\item \code{plot_lambdas} is now called via \code{plot(., plot_type = "lambda")}
-\item \code{plot_psis} is now called via \code{plot(., plot_type = "psi")}
+\item \code{plot_betas} is now called via \code{plot(., types = "beta")}
+\item \code{plot_deltas} is now called via \code{plot(., types = "delta")}
+\item \code{plot_nus} is now called via \code{plot(., types = "nu")}
+\item \code{plot_lambdas} is now called via \code{plot(., types = "lambda")}
+\item \code{plot_psis} is now called via \code{plot(., types = "psi")}
 }
 }
 
diff --git a/man/figures/README-unnamed-chunk-7-1.png b/man/figures/README-unnamed-chunk-7-1.png
index 70b0491..d540319 100644
Binary files a/man/figures/README-unnamed-chunk-7-1.png and b/man/figures/README-unnamed-chunk-7-1.png differ
diff --git a/man/figures/README-unnamed-chunk-8-1.png b/man/figures/README-unnamed-chunk-8-1.png
index 16bde4e..40bcdad 100644
Binary files a/man/figures/README-unnamed-chunk-8-1.png and b/man/figures/README-unnamed-chunk-8-1.png differ
diff --git a/man/plot.dynamitefit.Rd b/man/plot.dynamitefit.Rd
index c44e661..a7c9180 100644
--- a/man/plot.dynamitefit.Rd
+++ b/man/plot.dynamitefit.Rd
@@ -6,8 +6,8 @@
 \usage{
 \method{plot}{dynamitefit}(
   x,
-  plot_type = c("default", "beta", "delta", "nu", "lambda", "psi"),
-  type = NULL,
+  plot_type = c("default", "trace"),
+  types = NULL,
   parameters = NULL,
   responses = NULL,
   groups = NULL,
@@ -15,8 +15,6 @@
   level = 0.05,
   alpha = 0.5,
   scales = c("fixed", "free"),
-  include_alpha = TRUE,
-  include_cutpoints = TRUE,
   n_params = NULL,
   ...
 )
@@ -25,34 +23,28 @@
 \item{x}{[\code{dynamitefit}]\cr The model fit object.}
 
 \item{plot_type}{[\code{character(1)}]\cr What type of plot to draw? The default
-is \code{"default"} which draws posterior densities and traceplots of the
-selected parameters. Other options are \code{"beta"}, \code{"delta"}, \code{"nu"},
-\code{"lambda"} and \code{"psi"}. Option \code{"beta"} draws the posterior intervals
-and means of the \code{beta} type parameters. Option \code{"delta"} draws the
-posterior means and intervals as function of time for the \code{delta} type
-parameters. Option \code{"nu"} draws the posterior means and intervals for the
-group-specific random intercepts. Option \code{"lambda"} draws the posterior
-means and intervals for the factor loading, and option \code{"psi"} draws
-the means and posterior intervals for the latent factors as a function of
-time.}
+is \code{"default"} which draws posterior means and intervals of the parameters
+selected by \code{types} or \code{parameters}. If both \code{"types"} and
+\code{parameters} are \code{NULL}, all parameters are drawn up to the maximum
+specified by \code{n_params}. Option \code{"trace"} instead draws posterior
+densities and traceplots of the parameters.}
 
-\item{type}{[\code{character(1)}]\cr Type of the parameter for which the plots
+\item{types}{[\code{character(1)}]\cr Types of the parameter for which the plots
 should be drawn. Possible options can be found with the function
 \code{\link[=get_parameter_types]{get_parameter_types()}}. Ignored if the argument \code{parameters}
-is supplied or if \code{plot_type} is not \code{"default"}.}
+is supplied.}
 
 \item{parameters}{[\code{charecter()}]\cr Parameter name(s) for which the plots
 should be drawn. Possible options can be found with the function
-\code{\link[=get_parameter_names]{get_parameter_names()}}. The default is all parameters of a
-specific type, limited by \code{n_params}.}
+\code{\link[=get_parameter_names]{get_parameter_names()}}. The default is all parameters,
+limited by \code{n_params}.}
 
 \item{responses}{[\code{character()}]\cr Response(s) for which the plots should
 be drawn. Possible options are \code{unique(x$priors$response)}. Default is
 all responses. Ignored if the argument \code{parameters} is supplied.}
 
 \item{groups}{[\code{character(1)}]\cr Group name(s) for which the plots
-should be drawn when \code{plot_type} is \code{"nu"} or \code{"lambda"}.
-The default is all groups.}
+should be drawn for group-specific parameters.}
 
 \item{times}{[\code{double()}]\cr Time point(s) for which the plots should be
 drawn for time-varying parameters. By default, all time points are
@@ -68,20 +60,14 @@ Default is 0.5.}
 \item{scales}{[\code{character(1)}]\cr Should y-axis of the panels be \code{"fixed"}
 (the default) or \code{"free"}? See \code{\link[ggplot2:facet_wrap]{ggplot2::facet_wrap()}}.}
 
-\item{include_alpha}{[\code{logical(1)}]\cr If \code{TRUE} (default), plots also
-the alphas if such parameters exists in the model (either time-varying or
-time-invariant)}
-
-\item{include_cutpoints}{[\code{logical(1)}]\cr If \code{TRUE} (default), plots also
-the cutpoints if such parameters exists in the model (either time-varying
-or time-invariant).}
-
-\item{n_params}{[\code{integer()}]\cr Maximum number of parameters to plot.
-The default value is set by \code{plot_type}: 5 for \code{"default"}, 50 for
-\code{"beta"}, \code{"nu"} and \code{"lambda"}, and 3 for \code{"delta"} and \code{"psi"}. This
-argument is intended to prevent accidental plots that may be very large
-and time consuming to plot. Please use the \code{parameters}, \code{times}, and
-\code{groups} arguments to fine-tune which parameters to plot.}
+\item{n_params}{[\code{integer()}]\cr A single value or a vector of length 2
+specifying the maximum number of parameters to plot. If a single value
+is provided, the same limit is used for all parameters. If a vector is
+supplied, the first element defines the maximum number of time-invariant
+parameters to plot and the second the maximum number of time-varying
+parameters to plot. The defaults values are 50 for time-invariant
+parameters and 3 for time-varying parameters. The default value is 5
+for \code{plot_type == "trace"}.}
 
 \item{...}{Not used..}
 }
diff --git a/tests/testthat/test-errors.R b/tests/testthat/test-errors.R
index 6d0d088..4a3fed0 100644
--- a/tests/testthat/test-errors.R
+++ b/tests/testthat/test-errors.R
@@ -1404,31 +1404,14 @@ test_that("constrained prior for unconstrained parameter fails", {
 
 test_that("plot errors when the input is not a dynamitefit object", {
   expect_error(
-    plot.dynamitefit(1, type = "beta"),
+    plot.dynamitefit(1, types = "beta"),
     "Argument `x` must be a <dynamitefit> object."
   )
 })
 
-test_that("plot errors when no type is defined", {
-  expect_error(
-    plot(categorical_example_fit),
-    paste0(
-      "Either `parameters` or `type` must be provided when `plot_type` ",
-      "is \"default\"\\."
-    )
-  )
-})
-
-test_that("plot errors when type is vector", {
-  expect_error(
-    plot(gaussian_example_fit, type = c("beta", "delta")),
-    "Argument `type` must be a single <character> string."
-  )
-})
-
 test_that("plot errors when no variable is found ", {
   expect_error(
-    plot(categorical_example_fit, type = "delta"),
+    plot(categorical_example_fit, types = "delta"),
     paste0(
       "No parameters of type `delta` found for any of the response ",
       "channels `x` and `y`."
@@ -1436,69 +1419,12 @@ test_that("plot errors when no variable is found ", {
   )
 })
 
-test_that("plotting deltas errors when the model does not contain deltas", {
-  expect_error(
-    plot(categorical_example_fit, plot_type = "delta"),
-    "The model does not contain varying coefficients delta."
-  )
-})
-
 test_that("plotting nus errors when the model does not contain nus", {
   expect_error(
-    plot(categorical_example_fit, plot_type = "nu"),
-    "The model does not contain random effects nu."
-  )
-})
-
-test_that("plotting betas errors when incorrect parameters are supplied", {
-  expect_error(
-    plot(
-      gaussian_example_fit,
-      plot_type = "beta",
-      parameters = "delta_y_x"
-    ),
+    plot(categorical_example_fit, types = "nu"),
     paste0(
-      "Parameter `delta_y_x` not found or it is of wrong type:\n",
-      'i Use `get_parameter_names\\(\\)` with `types = "beta"` to check ',
-      'suitable parameter names\\.'
-    )
-  )
-})
-
-test_that("plotting deltas errors when incorrect parameters are supplied", {
-  expect_error(
-    plot(
-      gaussian_example_fit,
-      plot_type = "delta",
-      parameters = c("a", "delta_y_x")
-    ),
-    paste0(
-      "Parameter `a` not found or it is of wrong type:\n",
-      'i Use `get_parameter_names\\(\\)` with `types = "delta"` to check ',
-      'suitable parameter names\\.'
-    )
-  )
-})
-
-test_that("plotting nus errors when incorrect parameters are supplied", {
-  expect_error(
-    plot(
-      gaussian_example_fit,
-      plot_type = "nu",
-      parameters = 5
-    ),
-    "Argument `parameters` must be a <character> vector\\."
-  )
-  expect_error(
-    plot(
-      gaussian_example_fit,
-      plot_type = "nu",
-      parameters = "test"
-    ),
-    paste0(
-      "Parameter `test` not found or it is of wrong type:\n",
-      "i Use `get_parameter_names\\(\\)` with `types = \"nu\"` to check ",
-      "suitable parameter names\\."
+      "No parameters of type `nu` found for any of the response ",
+      "channels `x` and `y`."
     )
   )
 })
diff --git a/tests/testthat/test-output.R b/tests/testthat/test-output.R
index e8fda77..793846a 100644
--- a/tests/testthat/test-output.R
+++ b/tests/testthat/test-output.R
@@ -171,25 +171,25 @@ test_that("tikz formula plot works", {
 
 test_that("betas can be plotted", {
   expect_error(
-    plot(gaussian_example_fit, plot_type = "beta"),
+    plot(gaussian_example_fit, types = "beta"),
     NA
   )
   expect_error(
-    plot(categorical_example_fit, plot_type = "beta"),
+    plot(categorical_example_fit, types = "beta"),
     NA
   )
 })
 
 test_that("deltas can be plotted", {
   expect_error(
-    plot(gaussian_example_fit, plot_type = "delta"),
+    plot(gaussian_example_fit, types = "delta"),
     NA
   )
 })
 
 test_that("nus can be plotted", {
   expect_error(
-    plot(gaussian_example_fit, plot_type = "nu"),
+    plot(gaussian_example_fit, types = "nu"),
     NA
   )
 })
diff --git a/vignettes/dynamite.Rmd b/vignettes/dynamite.Rmd
index eb068c6..681a12d 100644
--- a/vignettes/dynamite.Rmd
+++ b/vignettes/dynamite.Rmd
@@ -497,13 +497,13 @@ For `"dynamitefit"` objects, the `summary()` method is a shortcut for `as.data.f
 The default `plot()` method for `"dynamitefit"` objects produces the marginal posterior densities and traceplots of the MCMC chains for the desired parameters by using the `ggplot2` package to produce the plot. This method has the same `parameters` and `responses` arguments as the `as.data.frame()` method, and the argument `type`, which is analogous to the `types` argument of `as.data.frame()`, but accepts only a single character string for plotting parameters of a specific type. Figure \@ref(fig:parameterposteriorplot) shows the posterior density and traceplot of the `beta_y_z` parameter, which is the time-invariant regression coefficient of `y` in the `gaussian_example_fit` model.
 
 ```{r parameterposteriorplot, echo=TRUE, eval=TRUE, warning=FALSE, fig.width=7, fig.height=3.5, fig.cap="Posterior density and traceplots for the time-invariant regression coefficient of `z` in the model for the response variable `y` (parameter `beta_y_z`) in the `gaussian_example_fit` model."}
-plot(gaussian_example_fit, responses = "y", type = "beta")
+plot(gaussian_example_fit, responses = "y", types = "beta")
 ```
 
 There are also customized plotting functions available for specific parameters types such as time-invariant regression coefficients, time-varying regression coefficients, group-specific random effects, latent factor loadings, and latent factors. For instance, we can obtain the posteriors of the time-varying coefficients for the `gaussian_example_fit` object to see how the parameter values change over time, as shown in Figure \@ref(fig:gaussiandeltas).
 
 ```{r gaussiandeltas, echo=TRUE, eval=TRUE, warning=FALSE, fig.width=7, fig.height=3.5, fig.cap="Posterior mean and 90\\% intervals of the time-varying coefficients for the `gaussian_example_fit` model. The panels from left to right show the time-varying intercept for `y`, the time-varying effect of `x` on `y`, and the time-varying effect of `lag(y)` (the previous time-point) on `y`."}
-plot(gaussian_example_fit, plot_type = "delta") + labs(title = "")
+plot(gaussian_example_fit, types = "delta") + labs(title = "")
 ```
 
 The generated Stan code of the model can be extracted with the method `get_code()` as a `"character"` string. This feature of geared towards advanced users who may for example need to make slight modifications to the generated code in order to adapt the model to a specific scenario that cannot be accomplished with the `dynamite` model syntax. The generated code also contains helpful annotations describing the model blocks, parameters, and complicated code sections. Using the argument `blocks`, we can extract only specific blocks of the full model code. To illustrate, we extract the parameters block of the `gaussian_example_fit` model code as the full model code is too large to display.