diff --git a/NAMESPACE b/NAMESPACE
index cb8df574..281890ba 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -72,6 +72,7 @@ S3method(stancode,mvgam_prefit)
 S3method(standata,mvgam_prefit)
 S3method(summary,mvgam)
 S3method(summary,mvgam_prefit)
+S3method(tidy,mvgam)
 S3method(update,jsdgam)
 S3method(update,mvgam)
 S3method(variables,mvgam)
@@ -166,6 +167,7 @@ export(student_t)
 export(t2)
 export(te)
 export(ti)
+export(tidy)
 export(tweedie)
 export(variables)
 importFrom(Rcpp,evalCpp)
@@ -203,6 +205,7 @@ importFrom(brms,stancode)
 importFrom(brms,standata)
 importFrom(brms,student)
 importFrom(generics,augment)
+importFrom(generics,tidy)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_bar)
diff --git a/R/sysdata.rda b/R/sysdata.rda
index 91ec5c73..ca52d07d 100644
Binary files a/R/sysdata.rda and b/R/sysdata.rda differ
diff --git a/R/tidier_methods.R b/R/tidier_methods.R
index 23657a17..4db2c741 100644
--- a/R/tidier_methods.R
+++ b/R/tidier_methods.R
@@ -1,7 +1,332 @@
+#' @importFrom generics tidy
+#' @export
+generics::tidy
+
 #' @importFrom generics augment
 #' @export
 generics::augment
 
+
+#' Tidy an mvgam object's parameter posteriors
+#'
+#' Get parameters' posterior statistics, implementing the generic `tidy` from
+#' the package \pkg{broom}.
+#'
+#' The parameters are categorized by the column "type". For instance, the
+#' intercept of the observation model (i.e. the "formula" arg to `mvgam()`) has
+#' the "type" "observation_beta". The possible "type"s are:
+#'   * observation_family_extra_param: any extra parameters for your observation
+#'     model, e.g. sigma for a gaussian observation model. These parameters are
+#'     not directly derived from the latent trend components (continuing the
+#'     gaussian example, contrast to mu).
+#'   * observation_beta: betas from your observation model, excluding any smooths.
+#'     If your formula was `y ~ x1 + s(x2, bs='cr')`, then your intercept and
+#'     `x1`'s beta would be categorized as this.
+#'   * random_effect_group_level: Group-level random effects parameters, i.e.
+#'     the mean and sd of the distribution from which the specific random
+#'     intercepts/slopes are considered to be drawn from.
+#'   * random_effect_beta: betas for the individual random intercepts/slopes.
+#'   * trend_model_param: parameters from your `trend_model`.
+#'   * trend_beta: analog of "observation_beta", but for any `trend_formula`.
+#'   * trend_random_effect_group_level: analog of "random_effect_group_level",
+#'     but for any `trend_formula`.
+#'   * trend_random_effect_beta: analog of "random_effect_beta",
+#'     but for any `trend_formula`.
+#'
+#' Additionally, GP terms can be incorporated in several ways, leading to
+#' different "type"s (or absence!):
+#'   * `s(bs = "gp")`: No parameters returned.
+#'   * `gp()` in `formula`: "type" of "observation_param".
+#'   * `gp()` in `trend_formula`: "type" of "trend_formula_param".
+#'   * `GP()` in `trend_model`: "type" of "trend_model_param".
+#'
+#'
+#' @param x An object of class `mvgam`.
+#' @param probs The desired probability levels of the parameters' posteriors.
+#'   Defaults to `c(0.025, 0.5, 0.975)`, i.e. 2.5%, 50%, and 97.5%.
+#' @param ... Unused, included for generic consistency only.
+#' @returns A `tibble` containing:
+#'   * "parameter": The parameter in question.
+#'   * "type": The component of the model that the parameter belongs to (see details).
+#'   * "mean": The posterior mean.
+#'   * "sd": The posterior standard deviation.
+#'   * percentile(s): Any percentiles of interest from these posteriors.
+#'
+#' @family tidiers
+#'
+#' @examples
+#' \dontrun{
+#' set.seed(0)
+#' simdat <- sim_mvgam(T = 100,
+#'                     n_series = 3,
+#'                     trend_model = AR(),
+#'                     prop_trend = 0.75,
+#'                     family = gaussian())
+#' simdat$data_train$x = rnorm(nrow(simdat$data_train))
+#' simdat$data_train$year_fac = factor(simdat$data_train$year)
+#'
+#' mod <- mvgam(y ~ - 1 + s(time, by = series, bs = 'cr', k = 20) + x,
+#'              trend_formula = ~ s(year_fac, bs = 're') - 1,
+#'              trend_model = AR(cor = TRUE),
+#'              family = gaussian(),
+#'              data = simdat$data_train,
+#'              silent = 2)
+#'
+#' tidy(mod, probs = c(0.2, 0.5, 0.8))
+#' }
+#'
+#' @export
+tidy.mvgam <- function(x, probs = c(0.025, 0.5, 0.975), ...) {
+  object <- x
+  obj_vars <- variables(object)
+  digits <- 2  # TODO: Let user change?
+  partialized_mcmc_summary <- purrr::partial(mcmc_summary,
+                                             object$model_output,
+                                             ... =,
+                                             ISB = FALSE,  # Matches `x[i]`'s rather than `x`.
+                                             probs = probs,
+                                             digits = digits,
+                                             Rhat = FALSE,
+                                             n.eff = FALSE)
+  out <- tibble::tibble()
+
+  # Observation family extra parameters --------
+  xp_names_all <- obj_vars$observation_pars$orig_name
+  # no matches -> length(xp_names) == 0, even if xp_names_all is NULL
+  xp_names <- grep("vec", xp_names_all, value = TRUE, invert = TRUE)
+  if (length(xp_names) > 0) {
+    extra_params_out <- partialized_mcmc_summary(params = xp_names)
+    extra_params_out <- tibble::add_column(extra_params_out,
+                                          type = "observation_family_extra_param",
+                                          .before = 1)
+    out <- dplyr::bind_rows(out, extra_params_out)
+  }
+  # END Observation family extra parameters
+
+  # obs non-smoother betas --------
+  if (object$mgcv_model$nsdf > 0) {
+    obs_beta_name_map <- dplyr::slice_head(obj_vars$observation_betas, n = object$mgcv_model$nsdf)  # df("orig_name", "alias")
+    obs_betas_out <- partialized_mcmc_summary(params = obs_beta_name_map$orig_name)
+    row.names(obs_betas_out) <- obs_beta_name_map$alias
+    obs_betas_out <- tibble::add_column(obs_betas_out,
+                                       type = "observation_beta",
+                                       .before = 1)
+    out <- dplyr::bind_rows(out, obs_betas_out)
+  }
+  # END obs non-smoother betas
+
+  # random effects --------
+  # TODO: names for random slopes
+  re_param_name_map <- obj_vars$observation_re_params
+  if (!is.null(re_param_name_map)) {
+    re_params_out <- partialized_mcmc_summary(params = re_param_name_map$orig_name)
+    row.names(re_params_out) <- re_param_name_map$alias
+    re_params_out <- tibble::add_column(re_params_out,
+                                          type = "random_effect_group_level",
+                                          .before = 1)
+    out <- dplyr::bind_rows(out, re_params_out)
+
+    # specific betas
+    for (sp in object$mgcv_model$smooth) {
+      if (inherits(sp, "random.effect")) {
+        re_label <- sp$label
+        betas_all <- obj_vars$observation_betas
+        re_beta_idxs <- grep(re_label, betas_all$alias, fixed = TRUE)
+        re_beta_name_map <- dplyr::slice(betas_all, re_beta_idxs)
+        re_betas_out <- partialized_mcmc_summary(params = re_beta_name_map$orig_name)
+        row.names(re_betas_out) <- re_beta_name_map$alias
+        re_betas_out <- tibble::add_column(re_betas_out,
+                                           type = "random_effect_beta",
+                                           .before = 1)
+        out <- dplyr::bind_rows(out, re_betas_out)
+      }
+    }
+  }
+  # END random effects
+
+  # GPs --------
+  if (!is.null(obj_vars$trend_pars)) {
+    tm_param_names_all <- obj_vars$trend_pars$orig_name
+    gp_param_names <- grep("^alpha_gp|^rho_gp", tm_param_names_all, value = TRUE)
+    if (length(gp_param_names) > 0) {
+      gp_params_out <- partialized_mcmc_summary(params = gp_param_names)
+      # where is GP? can be in formula, trend_formula, or trend_model
+      if (grepl("^(alpha|rho)_gp_trend", gp_param_names[[1]])) {
+        param_type <- "trend_formula_param"
+      } else if (grepl("^(alpha|rho)_gp_", gp_param_names[[1]])) {  # hmph.
+        param_type <- "observation_param"
+      } else {
+        param_type <- "trend_model_param"
+      }
+      gp_params_out <- tibble::add_column(gp_params_out,
+                                          type = param_type,
+                                          .before = 1)
+      out <- dplyr::bind_rows(out, gp_params_out)
+    }
+  }
+  # END GPs
+
+  # RW, AR, CAR, VAR, ZMVN --------
+  # TODO: split out Sigma for heircor?
+  trend_model_name <- ifelse(inherits(object$trend_model, "mvgam_trend"),
+                       object$trend_model$trend_model,
+                       object$trend_model)  # str vs called obj as arg to mvgam
+  if (grepl("^VAR|^CAR|^AR|^RW|^ZMVN", trend_model_name)) {
+    # theta = MA terms
+    # alpha_cor = heirarchical corr term
+    # A = VAR auto-regressive matrix
+    # Sigma = correlated errors matrix
+    # sigma = errors
+
+    # setting up the params to extract
+    if (trend_model_name == "VAR") {
+      trend_model_params <- c("^A\\[", "^alpha_cor", "^theta", "^Sigma")
+    } else if (grepl("^CAR|^AR|^RW", trend_model_name)) {
+      cor <- inherits(object$trend_model, "mvgam_trend") && object$trend_model$cor
+      sigma_name <- ifelse(cor, "^Sigma", "^sigma")
+      trend_model_params <- c("^ar", "^alpha_cor", "^theta", sigma_name)
+    } else if (grepl("^ZMVN", trend_model_name)) {
+      trend_model_params <- c("^alpha_cor", "^Sigma")
+    }
+
+    # extracting the params
+    trend_model_params <- paste(trend_model_params, collapse = "|")
+    tm_param_names_all <- obj_vars$trend_pars$orig_name
+    tm_param_names <- grep(trend_model_params, tm_param_names_all, value = TRUE)
+    tm_params_out <- partialized_mcmc_summary(params = tm_param_names)
+    tm_params_out <- tibble::add_column(tm_params_out,
+                                          type = "trend_model_param",
+                                          .before = 1)
+    out <- dplyr::bind_rows(out, tm_params_out)
+  }
+  # END RW, AR, CAR, VAR
+
+  # 'None' trend_model with a trend_formula --------
+  if (trend_model_name == "None" && !is.null(object$trend_call)) {
+    trend_pars_names_all <- obj_vars$trend_pars$orig_name
+    trend_pars_names <- grep("sigma", trend_pars_names_all, value = TRUE)
+    if (length(trend_pars_names) > 0) {
+      trend_params_out <- partialized_mcmc_summary(params = trend_pars_names)
+      trend_params_out <- tibble::add_column(trend_params_out,
+                                             type = "trend_model_param",
+                                             .before = 1)
+      out <- dplyr::bind_rows(out, trend_params_out)
+    }
+  }
+  # END 'None' trend_model with a trend_formula
+
+  # Piecewise --------
+  # TODO: potentially lump into AR section, above; how to handle change points?
+  # to lump in, just add an
+  # `else if (grepl("^PW", trend_model_name)`, then
+  # `trend_model_params <- c("^k_trend", "^m_trend", "^delta_trend")`
+  # and change initial grep(ar car var) call
+  if (grepl("^PW", trend_model_name)) {
+    trend_model_params <- "^k_trend|^m_trend|^delta_trend"
+    tm_param_names_all <- obj_vars$trend_pars$orig_name
+    tm_param_names <- grep(trend_model_params, tm_param_names_all, value = TRUE)
+    tm_params_out <- partialized_mcmc_summary(params = tm_param_names)
+    tm_params_out <- tibble::add_column(tm_params_out,
+                                          type = "trend_model_param",
+                                          .before = 1)
+    out <- dplyr::bind_rows(out, tm_params_out)
+  }
+  # END Piecewise
+
+  # Trend formula betas --------
+  if (!is.null(object$trend_call) && object$trend_mgcv_model$nsdf > 0) {
+    trend_beta_name_map <- dplyr::slice_head(obj_vars$trend_betas,
+                                           n = object$trend_mgcv_model$nsdf)  # df("orig_name", "alias")
+    trend_betas_out <- partialized_mcmc_summary(params = trend_beta_name_map$orig_name)
+    row.names(trend_betas_out) <- trend_beta_name_map$alias
+    trend_betas_out <- tibble::add_column(trend_betas_out,
+                                        type = "trend_beta",
+                                        .before = 1)
+    out <- dplyr::bind_rows(out, trend_betas_out)
+  }
+  # END Trend formula betas
+
+  # trend random effects --------
+  trend_re_param_name_map <- obj_vars$trend_re_params
+  if (!is.null(trend_re_param_name_map)) {
+    trend_re_params_out <- partialized_mcmc_summary(params = trend_re_param_name_map$orig_name)
+    row.names(trend_re_params_out) <- trend_re_param_name_map$alias
+    trend_re_params_out <- tibble::add_column(trend_re_params_out,
+                                        type = "trend_random_effect_group_level",
+                                        .before = 1)
+    out <- dplyr::bind_rows(out, trend_re_params_out)
+
+    # specific betas
+    for (sp in object$trend_mgcv_model$smooth) {
+      if (inherits(sp, "random.effect")) {
+        trend_re_label <- sp$label
+        trend_betas_all <- obj_vars$trend_betas
+        trend_re_beta_idxs <- grep(trend_re_label, trend_betas_all$alias, fixed = TRUE)
+        trend_re_beta_name_map <- dplyr::slice(trend_betas_all, trend_re_beta_idxs)
+        trend_re_betas_out <- partialized_mcmc_summary(params = trend_re_beta_name_map$orig_name)
+        row.names(trend_re_betas_out) <- trend_re_beta_name_map$alias
+        trend_re_betas_out <- tibble::add_column(trend_re_betas_out,
+                                           type = "trend_random_effect_beta",
+                                           .before = 1)
+        out <- dplyr::bind_rows(out, trend_re_betas_out)
+      }
+    }
+  }
+  # END trend random effects
+
+  # Cleanup output --------
+  # TODO: might need to put this prior to every bind_rows to avoid rowname dups.
+  out <- tibble::rownames_to_column(out, "parameter")
+
+  # Split Sigma in case of hierarchical residual correlations
+  alpha_cor_matches <- grep("alpha_cor", out$parameter, fixed = TRUE)
+  if (length(alpha_cor_matches) > 0) {
+    out <- split_hier_Sigma(object, out)
+  }
+  # END Cleanup output
+
+  out
+}
+
+
+#' Helper function to split apart Sigma into its constituent sub-matrixes in
+#' the case of a hierarchical latent process.
+#'
+#' The default MCMC output has dummy parameters filling out Sigma to make it
+#' an nxn matrix. This removes those, and renames the remaining sub-matrixes
+#' to align with the `gr` and `subgr` sizes from `mvgam()`'s `trend_model` argument.
+#'
+#' @param object An object of class `mvgam`.
+#' @param params `tibble` The parameters that are going to be returned by
+#'   `tidy.mvgam()`. Assumed that the columns match what `tidy.mvgam()` will return.
+#'   Specifically, that there is a "parameter" column.
+#' @returns `tibble` The `params`, but with the Sigma parameters split up by `gr`.
+#' @noRd
+split_hier_Sigma <- function(object, params) {
+  params_nonSigma <- dplyr::filter(params, !grepl("^Sigma", parameter))
+  params_Sigma <- dplyr::filter(params, grepl("^Sigma", parameter))
+
+  gr <- object$trend_model$gr
+  subgr <- object$trend_model$subgr
+  gr_levels <- levels(object$obs_data[[gr]])
+  subgr_levels <- levels(object$obs_data[[subgr]])
+  n_gr <- length(gr_levels)
+  n_subgr <- length(subgr_levels)
+
+  # anything besides the dummy params should have non-zero sd
+  params_Sigma <- dplyr::filter(params_Sigma, mean != 0, sd != 0)
+  index_strs <- sub("Sigma", "", params_Sigma$parameter)[1:(n_subgr ** 2)]
+
+  # new names
+  new_names <- paste0("Sigma_",
+                     rep(seq_len(n_gr), each = n_subgr ** 2),
+                     index_strs)
+  params_Sigma["parameter"] <- new_names
+
+  dplyr::bind_rows(params_nonSigma, params_Sigma)
+}
+
+
 #' Augment an mvgam object's data
 #'
 #' Add fits and residuals to the data, implementing the generic `augment` from
@@ -27,6 +352,8 @@ generics::augment
 #'   * The residuals, along with their variability and credible bounds.
 #'
 #' @seealso \code{\link{residuals.mvgam}}, \code{\link{fitted.mvgam}}
+#' @family tidiers
+#'
 #' @examples
 #' \dontrun{
 #' set.seed(0)
diff --git a/man/augment.mvgam.Rd b/man/augment.mvgam.Rd
index 52a36440..d9b8fd43 100644
--- a/man/augment.mvgam.Rd
+++ b/man/augment.mvgam.Rd
@@ -62,4 +62,8 @@ augment(mod1, robust = TRUE, probs = c(0.25, 0.75))
 }
 \seealso{
 \code{\link{residuals.mvgam}}, \code{\link{fitted.mvgam}}
+
+Other tidiers: 
+\code{\link{tidy.mvgam}()}
 }
+\concept{tidiers}
diff --git a/man/reexports.Rd b/man/reexports.Rd
index da19e2d5..6885fb7a 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -40,6 +40,7 @@
 \alias{posterior_linpred}
 \alias{stancode}
 \alias{standata}
+\alias{tidy}
 \alias{augment}
 \title{Objects exported from other packages}
 \keyword{internal}
@@ -50,7 +51,7 @@ below to see their documentation.
 \describe{
   \item{brms}{\code{\link[brms:conditional_effects.brmsfit]{conditional_effects}}, \code{\link[brms]{gp}}, \code{\link[brms:mcmc_plot.brmsfit]{mcmc_plot}}, \code{\link[brms:set_prior]{prior}}, \code{\link[brms:set_prior]{prior_}}, \code{\link[brms:set_prior]{prior_string}}, \code{\link[brms]{set_prior}}, \code{\link[brms]{stancode}}, \code{\link[brms]{standata}}}
 
-  \item{generics}{\code{\link[generics]{augment}}}
+  \item{generics}{\code{\link[generics]{augment}}, \code{\link[generics]{tidy}}}
 
   \item{insight}{\code{\link[insight]{get_data}}}
 
diff --git a/man/tidy.mvgam.Rd b/man/tidy.mvgam.Rd
new file mode 100644
index 00000000..e0c5a969
--- /dev/null
+++ b/man/tidy.mvgam.Rd
@@ -0,0 +1,90 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tidier_methods.R
+\name{tidy.mvgam}
+\alias{tidy.mvgam}
+\title{Tidy an mvgam object's parameter posteriors}
+\usage{
+\method{tidy}{mvgam}(x, probs = c(0.025, 0.5, 0.975), ...)
+}
+\arguments{
+\item{x}{An object of class \code{mvgam}.}
+
+\item{probs}{The desired probability levels of the parameters' posteriors.
+Defaults to \code{c(0.025, 0.5, 0.975)}, i.e. 2.5\%, 50\%, and 97.5\%.}
+
+\item{...}{Unused, included for generic consistency only.}
+}
+\value{
+A \code{tibble} containing:
+\itemize{
+\item "parameter": The parameter in question.
+\item "type": The component of the model that the parameter belongs to (see details).
+\item "mean": The posterior mean.
+\item "sd": The posterior standard deviation.
+\item percentile(s): Any percentiles of interest from these posteriors.
+}
+}
+\description{
+Get parameters' posterior statistics, implementing the generic \code{tidy} from
+the package \pkg{broom}.
+}
+\details{
+The parameters are categorized by the column "type". For instance, the
+intercept of the observation model (i.e. the "formula" arg to \code{mvgam()}) has
+the "type" "observation_beta". The possible "type"s are:
+\itemize{
+\item observation_family_extra_param: any extra parameters for your observation
+model, e.g. sigma for a gaussian observation model. These parameters are
+not directly derived from the latent trend components (continuing the
+gaussian example, contrast to mu).
+\item observation_beta: betas from your observation model, excluding any smooths.
+If your formula was \code{y ~ x1 + s(x2, bs='cr')}, then your intercept and
+\code{x1}'s beta would be categorized as this.
+\item random_effect_group_level: Group-level random effects parameters, i.e.
+the mean and sd of the distribution from which the specific random
+intercepts/slopes are considered to be drawn from.
+\item random_effect_beta: betas for the individual random intercepts/slopes.
+\item trend_model_param: parameters from your \code{trend_model}.
+\item trend_beta: analog of "observation_beta", but for any \code{trend_formula}.
+\item trend_random_effect_group_level: analog of "random_effect_group_level",
+but for any \code{trend_formula}.
+\item trend_random_effect_beta: analog of "random_effect_beta",
+but for any \code{trend_formula}.
+}
+
+Additionally, GP terms can be incorporated in several ways, leading to
+different "type"s (or absence!):
+\itemize{
+\item \code{s(bs = "gp")}: No parameters returned.
+\item \code{gp()} in \code{formula}: "type" of "observation_param".
+\item \code{gp()} in \code{trend_formula}: "type" of "trend_formula_param".
+\item \code{GP()} in \code{trend_model}: "type" of "trend_model_param".
+}
+}
+\examples{
+\dontrun{
+set.seed(0)
+simdat <- sim_mvgam(T = 100,
+                    n_series = 3,
+                    trend_model = AR(),
+                    prop_trend = 0.75,
+                    family = gaussian())
+simdat$data_train$x = rnorm(nrow(simdat$data_train))
+simdat$data_train$year_fac = factor(simdat$data_train$year)
+
+mod <- mvgam(y ~ - 1 + s(time, by = series, bs = 'cr', k = 20) + x,
+             trend_formula = ~ s(year_fac, bs = 're') - 1,
+             trend_model = AR(cor = TRUE),
+             family = gaussian(),
+             data = simdat$data_train,
+             silent = 2)
+
+tidy(mod, probs = c(0.2, 0.5, 0.8))
+}
+
+}
+\seealso{
+Other tidiers: 
+\code{\link{augment.mvgam}()}
+}
+\concept{tidiers}
diff --git a/tests/mvgam_examples.R b/tests/mvgam_examples.R
index 2cd17a47..e6da3afb 100644
--- a/tests/mvgam_examples.R
+++ b/tests/mvgam_examples.R
@@ -1,7 +1,8 @@
 # Small mvgam examples for testing post-fitting functions such as
 # predict, forecast, hindcast etc...
 testthat::skip_on_cran()
-set.seed(1234)
+SEED = 1234
+set.seed(SEED)
 library(mvgam)
 mvgam_examp_dat <- list(data_train =
                           structure(list(y = c(-0.317295790188251, 0.220334092025582,
@@ -151,7 +152,8 @@ mvgam_example1 <- mvgam(y ~ s(season, k = 5),
                         burnin = 300,
                         samples = 30,
                         chains = 1,
-                        backend = 'rstan')
+                        backend = 'rstan',
+                        seed = SEED)
 
 # Univariate process with trend_formula, trend_map and correlated process errors
 trend_map <- data.frame(series = unique(mvgam_examp_dat$data_train$series),
@@ -165,7 +167,8 @@ mvgam_example2 <- mvgam(y ~ 1,
                         burnin = 300,
                         samples = 30,
                         chains = 1,
-                        backend = 'rstan')
+                        backend = 'rstan',
+                        seed = SEED)
 
 # Multivariate process without trend_formula
 mvgam_example3 <- mvgam(y ~ s(season, k = 5),
@@ -175,7 +178,8 @@ mvgam_example3 <- mvgam(y ~ s(season, k = 5),
                         burnin = 300,
                         samples = 30,
                         chains = 1,
-                        backend = 'rstan')
+                        backend = 'rstan',
+                        seed = SEED)
 
 # Multivariate process with trend_formula and moving averages
 mvgam_example4 <- mvgam(y ~ series,
@@ -186,7 +190,8 @@ mvgam_example4 <- mvgam(y ~ series,
                         burnin = 300,
                         samples = 30,
                         chains = 1,
-                        backend = 'rstan')
+                        backend = 'rstan',
+                        seed = SEED)
 
 # GP dynamic factors (use list format to ensure it works in tests)
 list_data <- list()
@@ -203,7 +208,42 @@ mvgam_example5 <- mvgam(y ~ series + s(season, k = 5),
                         burnin = 300,
                         samples = 30,
                         chains = 1,
-                        backend = 'rstan')
+                        backend = 'rstan',
+                        seed = SEED)
+
+
+# Hierarchical dynamics example adapted from RW documentation example.
+# The difference is that this uses 4 species rather than 3.
+simdat1 <- sim_mvgam(trend_model = VAR(cor = TRUE),
+                     prop_trend = 0.95,
+                     n_series = 4,
+                     mu = c(1, 2, 3, 4))
+simdat2 <- sim_mvgam(trend_model = VAR(cor = TRUE),
+                     prop_trend = 0.95,
+                     n_series = 4,
+                     mu = c(1, 2, 3, 4))
+simdat3 <- sim_mvgam(trend_model = VAR(cor = TRUE),
+                     prop_trend = 0.95,
+                     n_series = 4,
+                     mu = c(1, 2, 3, 4))
+
+simdat_all <- rbind(simdat1$data_train %>%
+                      dplyr::mutate(region = 'qld'),
+                    simdat2$data_train %>%
+                      dplyr::mutate(region = 'nsw'),
+                    simdat3$data_train %>%
+                      dplyr::mutate(region = 'vic')) %>%
+  dplyr::mutate(species = gsub('series', 'species', series),
+                species = as.factor(species),
+                region = as.factor(region)) %>%
+  dplyr::arrange(series, time) %>%
+  dplyr::select(-series)
+
+mvgam_example6 <- mvgam(formula = y ~ species,
+                        trend_model = AR(gr = region, subgr = species),
+                        data = simdat_all,
+                        backend = 'rstan',
+                        seed = SEED)
 
 # Save examples as internal data
 usethis::use_data(
@@ -213,6 +253,7 @@ usethis::use_data(
   mvgam_example3,
   mvgam_example4,
   mvgam_example5,
+  mvgam_example6,
   internal = TRUE,
   overwrite = TRUE
 )
diff --git a/tests/testthat/_snaps/tidier_methods.md b/tests/testthat/_snaps/tidier_methods.md
new file mode 100644
index 00000000..3c7e5091
--- /dev/null
+++ b/tests/testthat/_snaps/tidier_methods.md
@@ -0,0 +1,360 @@
+# `tidy()` snapshot value of `mvgam_example1`
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["parameter", "type", "mean", "sd", "2.5%", "50%", "97.5%"]
+        },
+        "row.names": {
+          "type": "integer",
+          "attributes": {},
+          "value": [1, 2, 3, 4, 5, 6, 7]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["tbl_df", "tbl", "data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["sigma_obs[1]", "sigma_obs[2]", "sigma_obs[3]", "(Intercept)", "sigma[1]", "sigma[2]", "sigma[3]"]
+        },
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["observation_family_extra_param", "observation_family_extra_param", "observation_family_extra_param", "observation_beta", "trend_model_param", "trend_model_param", "trend_model_param"]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.5, 0.49, 0.62, 0.43, 0.15, 0.11, 0.041]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.079, 0.092, 0.09, 0.12, 0.061, 0.05, 0.02]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.35, 0.37, 0.48, 0.25, 0.056, 0.061, 0.023]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.5, 0.47, 0.6, 0.4, 0.17, 0.1, 0.034]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.68, 0.67, 0.8, 0.68, 0.23, 0.22, 0.091]
+        }
+      ]
+    }
+
+# `tidy()` snapshot value of `mvgam_example2`
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["parameter", "type", "mean", "sd", "2.5%", "50%", "97.5%"]
+        },
+        "row.names": {
+          "type": "integer",
+          "attributes": {},
+          "value": [1, 2, 3, 4, 5, 6, 7, 8, 9]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["tbl_df", "tbl", "data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["sigma_obs[1]", "sigma_obs[2]", "sigma_obs[3]", "(Intercept)", "Sigma[1,1]", "Sigma[2,1]", "Sigma[1,2]", "Sigma[2,2]", "(Intercept)_trend"]
+        },
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["observation_family_extra_param", "observation_family_extra_param", "observation_family_extra_param", "observation_beta", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_beta"]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.49, 0.44, 0.6, 0.099, 0.19, 0.083, 0.083, 0.16, 0.025]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.089, 0.079, 0.16, 1.7, 0.12, 0.097, 0.097, 0.099, 1.7]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.32, 0.3, 0.34, -2.1, 0.057, -0.068, -0.068, 0.012, -3.8]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.49, 0.45, 0.59, 0.0068, 0.15, 0.06, 0.06, 0.16, 0.15]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.67, 0.59, 0.92, 3.2, 0.5, 0.27, 0.27, 0.35, 2.3]
+        }
+      ]
+    }
+
+# `tidy()` snapshot value of `mvgam_example3`
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["parameter", "type", "mean", "sd", "2.5%", "50%", "97.5%"]
+        },
+        "row.names": {
+          "type": "integer",
+          "attributes": {},
+          "value": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["tbl_df", "tbl", "data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["sigma_obs[1]", "sigma_obs[2]", "sigma_obs[3]", "(Intercept)", "A[1,1]", "A[2,1]", "A[3,1]", "A[1,2]", "A[2,2]", "A[3,2]", "A[1,3]", "A[2,3]", "A[3,3]", "Sigma[1,1]", "Sigma[2,1]", "Sigma[3,1]", "Sigma[1,2]", "Sigma[2,2]", "Sigma[3,2]", "Sigma[1,3]", "Sigma[2,3]", "Sigma[3,3]"]
+        },
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["observation_family_extra_param", "observation_family_extra_param", "observation_family_extra_param", "observation_beta", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param"]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.54, 0.39, 0.52, 0.46, 0.035, 0.39, 0.41, -0.0021, 0.11, 0.061, -0.036, -0.18, 0.064, 0.0045, 0.0013, -0.0018, 0.0013, 0.12, 0.022, -0.0018, 0.022, 0.15]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.088, 0.14, 0.13, 0.058, 0.48, 2.2, 2.5, 0.19, 0.28, 1, 0.12, 0.57, 0.39, 0.0032, 0.0085, 0.015, 0.0085, 0.1, 0.027, 0.015, 0.027, 0.1]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.42, 0.18, 0.27, 0.36, -0.69, -3.6, -4, -0.32, -0.36, -1.4, -0.34, -1.4, -0.76, 0.0013, -0.016, -0.037, -0.016, 0.0053, -0.017, -0.037, -0.017, 0.0089]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.52, 0.41, 0.55, 0.46, -0.036, 0.48, 0.17, -0.0068, 0.12, 0.031, -0.0024, 0.0066, 0.13, 0.0037, 0.0014, 0.00075, 0.0014, 0.092, 0.026, 0.00075, 0.026, 0.16]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.72, 0.57, 0.71, 0.57, 0.88, 4.5, 5.4, 0.44, 0.68, 1.8, 0.084, 0.48, 0.66, 0.011, 0.017, 0.022, 0.017, 0.32, 0.072, 0.022, 0.072, 0.32]
+        }
+      ]
+    }
+
+# `tidy()` snapshot value of `mvgam_example4`
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["parameter", "type", "mean", "sd", "2.5%", "50%", "97.5%"]
+        },
+        "row.names": {
+          "type": "integer",
+          "attributes": {},
+          "value": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["tbl_df", "tbl", "data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["sigma_obs[1]", "sigma_obs[2]", "sigma_obs[3]", "(Intercept)", "seriesseries_2", "seriesseries_3", "A[1,1]", "A[2,1]", "A[3,1]", "A[1,2]", "A[2,2]", "A[3,2]", "A[1,3]", "A[2,3]", "A[3,3]", "theta[1,1]", "theta[2,1]", "theta[3,1]", "theta[1,2]", "theta[2,2]", "theta[3,2]", "theta[1,3]", "theta[2,3]", "theta[3,3]", "Sigma[1,1]", "Sigma[2,1]", "Sigma[3,1]", "Sigma[1,2]", "Sigma[2,2]", "Sigma[3,2]", "Sigma[1,3]", "Sigma[2,3]", "Sigma[3,3]", "(Intercept)_trend"]
+        },
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["observation_family_extra_param", "observation_family_extra_param", "observation_family_extra_param", "observation_beta", "observation_beta", "observation_beta", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_beta"]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.49, 0.34, 0.64, 0.31, -0.17, -0.22, 0.16, 1.6, 0.78, 0.16, 0.025, 0.1, 0.6, 1.1, -0.0088, -0.12, -0.8, -0.38, -0.024, 0.2, 0.11, 0.21, 0.25, 0.19, 0.015, 0.0035, -0.00015, 0.0035, 0.071, 0.0029, -0.00015, 0.0029, 0.0058, 0.24]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.071, 0.11, 0.13, 1.5, 0.3, 0.23, 0.48, 1.1, 0.74, 0.2, 0.38, 0.23, 1.2, 1.6, 0.36, 0.82, 2.7, 1.9, 0.31, 0.63, 0.29, 0.99, 2.2, 0.73, 0.021, 0.019, 0.004, 0.019, 0.05, 0.0093, 0.004, 0.0093, 0.0038, 1.4]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.38, 0.17, 0.46, -2, -0.77, -0.59, -0.63, 0.27, -0.028, -0.12, -0.65, -0.25, -1.1, -1.3, -0.56, -1.5, -8.5, -5.9, -0.72, -0.9, -0.4, -1.4, -4.4, -1.1, 0.00032, -0.02, -0.0089, -0.02, 0.012, -0.0097, -0.0089, -0.0097, 0.0016, -2.2]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.48, 0.34, 0.64, 0.05, -0.17, -0.21, 0.19, 1.1, 0.66, 0.088, 0.025, 0.059, 0.1, 0.8, -0.023, -0.21, -0.16, -0.032, 0.029, 0.33, 0.11, 0.14, -0.00025, 0.3, 0.0056, -0.00068, 0.00052, -0.00068, 0.059, 0.0015, 0.00052, 0.0015, 0.0057, 0.48]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.64, 0.53, 0.89, 2.8, 0.4, 0.21, 0.9, 4.2, 2.2, 0.58, 0.66, 0.48, 3.1, 4.2, 0.68, 1.4, 2.9, 1.2, 0.38, 1.3, 0.58, 2.1, 4.8, 1.5, 0.066, 0.047, 0.007, 0.047, 0.19, 0.025, 0.007, 0.025, 0.017, 2.4]
+        }
+      ]
+    }
+
+# `tidy()` snapshot value of `mvgam_example5`
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["parameter", "type", "mean", "sd", "2.5%", "50%", "97.5%"]
+        },
+        "row.names": {
+          "type": "integer",
+          "attributes": {},
+          "value": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["tbl_df", "tbl", "data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["sigma_obs[1]", "sigma_obs[2]", "sigma_obs[3]", "(Intercept)", "seriesseries_2", "seriesseries_3", "rho_gp[1]", "rho_gp[2]", "alpha_gp[1]", "alpha_gp[2]"]
+        },
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["observation_family_extra_param", "observation_family_extra_param", "observation_family_extra_param", "observation_beta", "observation_beta", "observation_beta", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param"]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.49, 0.48, 0.64, 0.56, -0.14, -0.28, 15, 4.7, 0.25, 0.25]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.09, 0.085, 0.082, 0.12, 0.23, 0.22, 16, 5, 0, 0]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.36, 0.36, 0.5, 0.31, -0.46, -0.65, 0.78, 1.6, 0.25, 0.25]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.47, 0.45, 0.62, 0.57, -0.21, -0.22, 6.4, 2.7, 0.25, 0.25]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.7, 0.67, 0.77, 0.76, 0.4, 0.016, 44, 19, 0.25, 0.25]
+        }
+      ]
+    }
+
+# `tidy()` snapshot value of `mvgam_example6`
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["parameter", "type", "mean", "sd", "2.5%", "50%", "97.5%"]
+        },
+        "row.names": {
+          "type": "integer",
+          "attributes": {},
+          "value": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["tbl_df", "tbl", "data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["(Intercept)", "speciesspecies_2", "speciesspecies_3", "speciesspecies_4", "ar1[1]", "ar1[2]", "ar1[3]", "ar1[4]", "ar1[5]", "ar1[6]", "ar1[7]", "ar1[8]", "ar1[9]", "ar1[10]", "ar1[11]", "ar1[12]", "alpha_cor", "Sigma_1[1,1]", "Sigma_1[2,1]", "Sigma_1[3,1]", "Sigma_1[4,1]", "Sigma_1[1,2]", "Sigma_1[2,2]", "Sigma_1[3,2]", "Sigma_1[4,2]", "Sigma_1[1,3]", "Sigma_1[2,3]", "Sigma_1[3,3]", "Sigma_1[4,3]", "Sigma_1[1,4]", "Sigma_1[2,4]", "Sigma_1[3,4]", "Sigma_1[4,4]", "Sigma_2[1,1]", "Sigma_2[2,1]", "Sigma_2[3,1]", "Sigma_2[4,1]", "Sigma_2[1,2]", "Sigma_2[2,2]", "Sigma_2[3,2]", "Sigma_2[4,2]", "Sigma_2[1,3]", "Sigma_2[2,3]", "Sigma_2[3,3]", "Sigma_2[4,3]", "Sigma_2[1,4]", "Sigma_2[2,4]", "Sigma_2[3,4]", "Sigma_2[4,4]", "Sigma_3[1,1]", "Sigma_3[2,1]", "Sigma_3[3,1]", "Sigma_3[4,1]", "Sigma_3[1,2]", "Sigma_3[2,2]", "Sigma_3[3,2]", "Sigma_3[4,2]", "Sigma_3[1,3]", "Sigma_3[2,3]", "Sigma_3[3,3]", "Sigma_3[4,3]", "Sigma_3[1,4]", "Sigma_3[2,4]", "Sigma_3[3,4]", "Sigma_3[4,4]"]
+        },
+        {
+          "type": "character",
+          "attributes": {},
+          "value": ["observation_beta", "observation_beta", "observation_beta", "observation_beta", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param", "trend_model_param"]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [1, 1.1, 2, 2.9, -0.34, 0.42, -0.084, 0.67, 0.13, 0.63, 0.012, 0.26, 0.51, 0.077, 0.053, 0.76, 0.11, 0.73, -0.17, 0.11, -0.38, -0.17, 0.98, 0.49, -0.017, 0.11, 0.49, 0.76, -0.062, -0.38, -0.017, -0.062, 0.59, 0.79, 0.081, 0.097, -0.074, 0.081, 0.52, 0.39, -0.024, 0.097, 0.39, 0.94, -0.36, -0.074, -0.024, -0.36, 0.95, 0.71, -0.019, 0.028, -0.083, -0.019, 0.81, -0.27, 0.54, 0.028, -0.27, 0.84, -0.13, -0.083, 0.54, -0.13, 0.67]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.081, 0.11, 0.09, 0.16, 0.16, 0.16, 0.17, 0.1, 0.21, 0.093, 0.11, 0.12, 0.2, 0.11, 0.14, 0.12, 0.07, 0.18, 0.13, 0.14, 0.1, 0.13, 0.21, 0.15, 0.092, 0.14, 0.15, 0.15, 0.096, 0.1, 0.092, 0.096, 0.1, 0.21, 0.097, 0.12, 0.11, 0.097, 0.12, 0.099, 0.089, 0.12, 0.099, 0.17, 0.11, 0.11, 0.089, 0.11, 0.16, 0.22, 0.1, 0.12, 0.11, 0.1, 0.17, 0.11, 0.13, 0.12, 0.11, 0.16, 0.1, 0.11, 0.13, 0.1, 0.13]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.87, 0.88, 1.8, 2.6, -0.64, 0.075, -0.41, 0.47, -0.29, 0.44, -0.22, 0.039, 0.074, -0.14, -0.22, 0.51, 0.017, 0.44, -0.44, -0.17, -0.6, -0.44, 0.65, 0.21, -0.2, -0.17, 0.21, 0.52, -0.24, -0.6, -0.2, -0.24, 0.42, 0.46, -0.1, -0.12, -0.3, -0.1, 0.34, 0.22, -0.2, -0.12, 0.22, 0.67, -0.6, -0.3, -0.2, -0.6, 0.69, 0.37, -0.24, -0.22, -0.31, -0.24, 0.54, -0.51, 0.32, -0.22, -0.51, 0.58, -0.34, -0.31, 0.32, -0.34, 0.45]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [1, 1.1, 2, 2.9, -0.35, 0.43, -0.076, 0.67, 0.13, 0.63, 0.016, 0.26, 0.52, 0.076, 0.05, 0.77, 0.091, 0.71, -0.17, 0.11, -0.37, -0.17, 0.96, 0.49, -0.017, 0.11, 0.49, 0.74, -0.064, -0.37, -0.017, -0.064, 0.58, 0.77, 0.079, 0.091, -0.072, 0.079, 0.51, 0.38, -0.023, 0.091, 0.38, 0.92, -0.35, -0.072, -0.023, -0.35, 0.93, 0.68, -0.016, 0.031, -0.078, -0.016, 0.8, -0.27, 0.53, 0.031, -0.27, 0.82, -0.12, -0.078, 0.53, -0.12, 0.66]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [1.2, 1.3, 2.2, 3.2, -0.032, 0.71, 0.24, 0.86, 0.54, 0.81, 0.24, 0.49, 0.84, 0.3, 0.33, 0.96, 0.27, 1.2, 0.084, 0.39, -0.19, 0.084, 1.5, 0.82, 0.16, 0.39, 0.82, 1.1, 0.13, -0.19, 0.16, 0.13, 0.83, 1.3, 0.27, 0.35, 0.15, 0.27, 0.78, 0.62, 0.15, 0.35, 0.62, 1.3, -0.16, 0.15, 0.15, -0.16, 1.3, 1.2, 0.18, 0.26, 0.12, 0.18, 1.2, -0.07, 0.81, 0.26, -0.07, 1.2, 0.061, 0.12, 0.81, 0.061, 0.96]
+        }
+      ]
+    }
+
diff --git a/tests/testthat/test-tidier_methods.R b/tests/testthat/test-tidier_methods.R
index 3aa6810c..12c629b2 100644
--- a/tests/testthat/test-tidier_methods.R
+++ b/tests/testthat/test-tidier_methods.R
@@ -1,5 +1,37 @@
 context("tidier methods")
 
+# `tidy()` tests
+test_that("`tidy()` snapshot value of `mvgam_example1`", {
+  local_edition(3)
+  expect_snapshot_value(tidy.mvgam(mvgam_example1), style = "json2")
+})
+
+test_that("`tidy()` snapshot value of `mvgam_example2`", {
+  local_edition(3)
+  expect_snapshot_value(tidy.mvgam(mvgam_example2), style = "json2")
+})
+
+test_that("`tidy()` snapshot value of `mvgam_example3`", {
+  local_edition(3)
+  expect_snapshot_value(tidy.mvgam(mvgam_example3), style = "json2")
+})
+
+test_that("`tidy()` snapshot value of `mvgam_example4`", {
+  local_edition(3)
+  expect_snapshot_value(tidy.mvgam(mvgam_example4), style = "json2")
+})
+
+test_that("`tidy()` snapshot value of `mvgam_example5`", {
+  local_edition(3)
+  expect_snapshot_value(tidy.mvgam(mvgam_example5), style = "json2")
+})
+
+test_that("`tidy()` snapshot value of `mvgam_example6`", {
+  local_edition(3)
+  expect_snapshot_value(tidy.mvgam(mvgam_example6), style = "json2")
+})
+
+# `augment()` tests
 test_that("augment doesn't error", {
   expect_no_error(augment(mvgam:::mvgam_example1))
   expect_no_error(augment(mvgam:::mvgam_example5))