diff --git a/DESCRIPTION b/DESCRIPTION
index 343968ef..9385b638 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,16 +1,16 @@
 Package: mvgam
 Title: Multivariate (Dynamic) Generalized Additive Models
 Version: 1.1.3
-Date: 2024-07-02
+Date: 2024-09-03
 Authors@R: person("Nicholas J", "Clark", , "nicholas.j.clark1214@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-7131-3301"))
 Description: Fit Bayesian Dynamic Generalized Additive Models to sets of time series. Users can build dynamic nonlinear State-Space models that can incorporate semiparametric effects in observation and process components, using a wide range of observation families. Estimation is performed using Markov Chain Monte Carlo with Hamiltonian Monte Carlo in the software 'Stan'. References: Clark & Wells (2022) <doi:10.1111/2041-210X.13974>.
 URL: https://github.com/nicholasjclark/mvgam, https://nicholasjclark.github.io/mvgam/
 BugReports: https://github.com/nicholasjclark/mvgam/issues
 License: MIT + file LICENSE
 Depends: 
-    R (>= 3.6.0),
-    brms (>= 2.21.0)
+    R (>= 3.6.0)
 Imports: 
+    brms (>= 2.21.0),
     methods,
     mgcv (>= 1.8-13),
     insight (>= 0.19.1),
diff --git a/NAMESPACE b/NAMESPACE
index 58ed9bd6..ea9ce1a2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -27,6 +27,7 @@ S3method(get_data,mvgam_prefit)
 S3method(get_predict,mvgam)
 S3method(get_vcov,mvgam)
 S3method(hindcast,mvgam)
+S3method(irf,mvgam)
 S3method(lfo_cv,mvgam)
 S3method(logLik,mvgam)
 S3method(log_posterior,mvgam)
@@ -41,6 +42,7 @@ S3method(pairs,mvgam)
 S3method(plot,mvgam)
 S3method(plot,mvgam_conditional_effects)
 S3method(plot,mvgam_forecast)
+S3method(plot,mvgam_irf)
 S3method(plot,mvgam_lfo)
 S3method(posterior_epred,mvgam)
 S3method(posterior_linpred,mvgam)
@@ -73,9 +75,14 @@ export(PW)
 export(RW)
 export(VAR)
 export(add_residuals)
+export(bernoulli)
+export(beta_binomial)
 export(betar)
 export(code)
 export(compare_mvgams)
+export(comparisons)
+export(conditional_effects)
+export(datagrid)
 export(drawDotmvgam)
 export(dynamic)
 export(ensemble)
@@ -85,14 +92,22 @@ export(eval_smoothDotmodDotsmooth)
 export(eval_smoothDotmoiDotsmooth)
 export(forecast)
 export(get_mvgam_priors)
+export(gp)
 export(hindcast)
+export(hypotheses)
+export(irf)
 export(lfo_cv)
+export(lognormal)
+export(loo)
+export(loo_compare)
 export(lv_correlations)
+export(mcmc_plot)
 export(mvgam)
 export(nb)
 export(neff_ratio)
 export(nmix)
 export(nuts_params)
+export(plot_comparisons)
 export(plot_mvgam_factors)
 export(plot_mvgam_fc)
 export(plot_mvgam_pterms)
@@ -102,13 +117,30 @@ export(plot_mvgam_series)
 export(plot_mvgam_smooth)
 export(plot_mvgam_trend)
 export(plot_mvgam_uncertainty)
+export(plot_predictions)
+export(plot_slopes)
+export(posterior_epred)
+export(posterior_linpred)
+export(posterior_predict)
 export(pp_check)
 export(ppc)
+export(predictions)
+export(prior)
+export(prior_)
+export(prior_string)
+export(rhat)
 export(roll_eval_mvgam)
+export(s)
 export(score)
 export(series_to_mvgam)
+export(set_prior)
 export(sim_mvgam)
+export(slopes)
+export(stancode)
+export(standata)
+export(student)
 export(student_t)
+export(t2)
 export(te)
 export(ti)
 export(tweedie)
@@ -122,22 +154,27 @@ importFrom(bayesplot,neff_ratio)
 importFrom(bayesplot,nuts_params)
 importFrom(bayesplot,pp_check)
 importFrom(brms,bernoulli)
+importFrom(brms,beta_binomial)
 importFrom(brms,conditional_effects)
 importFrom(brms,dbeta_binomial)
 importFrom(brms,do_call)
 importFrom(brms,dstudent_t)
 importFrom(brms,get_prior)
+importFrom(brms,gp)
 importFrom(brms,logm1)
 importFrom(brms,lognormal)
 importFrom(brms,mcmc_plot)
 importFrom(brms,ndraws)
 importFrom(brms,pbeta_binomial)
+importFrom(brms,prior)
+importFrom(brms,prior_)
 importFrom(brms,prior_string)
 importFrom(brms,pstudent_t)
 importFrom(brms,qstudent_t)
 importFrom(brms,rbeta_binomial)
 importFrom(brms,read_csv_as_stanfit)
 importFrom(brms,rstudent_t)
+importFrom(brms,set_prior)
 importFrom(brms,stancode)
 importFrom(brms,standata)
 importFrom(brms,student)
@@ -171,11 +208,18 @@ importFrom(loo,is.loo)
 importFrom(loo,loo)
 importFrom(loo,loo_compare)
 importFrom(magrittr,"%>%")
+importFrom(marginaleffects,comparisons)
+importFrom(marginaleffects,datagrid)
 importFrom(marginaleffects,get_coef)
 importFrom(marginaleffects,get_predict)
 importFrom(marginaleffects,get_vcov)
+importFrom(marginaleffects,hypotheses)
+importFrom(marginaleffects,plot_comparisons)
 importFrom(marginaleffects,plot_predictions)
+importFrom(marginaleffects,plot_slopes)
+importFrom(marginaleffects,predictions)
 importFrom(marginaleffects,set_coef)
+importFrom(marginaleffects,slopes)
 importFrom(methods,cbind2)
 importFrom(mgcv,Predict.matrix)
 importFrom(mgcv,Rrank)
@@ -186,8 +230,12 @@ importFrom(mgcv,get.var)
 importFrom(mgcv,initial.sp)
 importFrom(mgcv,interpret.gam)
 importFrom(mgcv,nb)
+importFrom(mgcv,s)
 importFrom(mgcv,smooth.construct)
 importFrom(mgcv,smoothCon)
+importFrom(mgcv,t2)
+importFrom(mgcv,te)
+importFrom(mgcv,ti)
 importFrom(parallel,clusterEvalQ)
 importFrom(parallel,clusterExport)
 importFrom(parallel,setDefaultCluster)
diff --git a/NEWS.md b/NEWS.md
index cb8e9e17..67f5ab72 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,7 @@
 * Added `standata.mvgam_prefit`, `stancode.mvgam` and `stancode.mvgam_prefit` methods for better alignment with 'brms' workflows
 * Added 'gratia' to *Enhancements* to allow popular methods such as `draw()` to be used for 'mvgam' models if 'gratia' is already installed
 * Added an `ensemble.mvgam_forecast` method to generate evenly weighted combinations of probabilistic forecast distributions
+* Added an `irf.mvgam` method to compute Generalized and Orthogonalized Impulse Response Functions (IRFs) from models fit with Vector Autoregressive dynamics
 
 ## Deprecations
 * The `drift` argument has been deprecated. It is now recommended for users to include parametric fixed effects of "time" in their respective GAM formulae to capture any expected drift effects
diff --git a/R/backends.R b/R/backends.R
index c1637a9e..a672f1b0 100644
--- a/R/backends.R
+++ b/R/backends.R
@@ -197,7 +197,8 @@
                                                          'sigma_obs_vec',
                                                          'shape_vec',
                                                          'phi_inv',
-                                                         'lv_coefs_raw'),
+                                                         'lv_coefs_raw',
+                                                         'L_Sigma'),
                                              algorithm = algorithm)
   }
 
diff --git a/R/conditional_effects.R b/R/conditional_effects.R
index e38b9b73..d0e15cc9 100644
--- a/R/conditional_effects.R
+++ b/R/conditional_effects.R
@@ -2,9 +2,7 @@
 #'
 #' Display conditional effects of one or more numeric and/or categorical
 #' predictors in `mvgam` models, including two-way interaction effects.
-#' @importFrom brms conditional_effects
 #' @importFrom ggplot2 scale_colour_discrete scale_fill_discrete theme_classic
-#' @importFrom marginaleffects plot_predictions
 #' @importFrom graphics plot
 #' @importFrom grDevices devAskNewPage
 #' @inheritParams brms::conditional_effects.brmsfit
@@ -219,6 +217,10 @@ conditional_effects.mvgam = function(x,
   return(out)
 }
 
+#' @export
+#' @importFrom brms conditional_effects
+brms::conditional_effects
+
 #' @rdname conditional_effects.mvgam
 #' @export
 plot.mvgam_conditional_effects = function(x,
diff --git a/R/families.R b/R/families.R
index 1795a9d8..a1d895c0 100644
--- a/R/families.R
+++ b/R/families.R
@@ -2,7 +2,7 @@
 #' @importFrom stats make.link dgamma pgamma rgamma qnorm plnorm runif pbeta dlnorm dpois
 #' @importFrom stats pnorm ppois plogis gaussian poisson Gamma dnbinom rnbinom dnorm dbeta
 #' @importFrom stats binomial rbinom pbinom dbinom qbinom qlogis
-#' @importFrom brms lognormal student bernoulli rstudent_t qstudent_t dstudent_t pstudent_t dbeta_binomial rbeta_binomial pbeta_binomial
+#' @importFrom brms lognormal student beta_binomial bernoulli rstudent_t qstudent_t dstudent_t pstudent_t dbeta_binomial rbeta_binomial pbeta_binomial
 #' @importFrom mgcv betar nb
 #' @param link a specification for the family link function. At present these cannot
 #' be changed
@@ -80,6 +80,30 @@ nb = function(...){
   mgcv::nb(...)
 }
 
+#' @rdname mvgam_families
+#' @export
+lognormal = function(...){
+  brms::lognormal(...)
+}
+
+#' @rdname mvgam_families
+#' @export
+student = function(...){
+  brms::student(...)
+}
+
+#' @rdname mvgam_families
+#' @export
+bernoulli = function(...){
+  brms::bernoulli(...)
+}
+
+#' @rdname mvgam_families
+#' @export
+beta_binomial = function(...){
+  brms::beta_binomial(...)
+}
+
 #' @rdname mvgam_families
 #' @examples
 #' \donttest{
diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R
index 266bd974..4d213c7e 100644
--- a/R/get_mvgam_priors.R
+++ b/R/get_mvgam_priors.R
@@ -1254,6 +1254,22 @@ get_mvgam_priors = function(formula,
   return(out)
 }
 
+#' @export
+#' @importFrom brms prior
+brms::prior
+
+#' @export
+#' @importFrom brms prior_
+brms::prior_
+
+#' @export
+#' @importFrom brms set_prior
+brms::set_prior
+
+#' @export
+#' @importFrom brms prior_string
+brms::prior_string
+
 #' Use informative scale and intercept priors following brms example
 #' @importFrom stats mad qcauchy setNames
 #' @importFrom brms logm1 prior_string get_prior
diff --git a/R/irf.mvgam.R b/R/irf.mvgam.R
new file mode 100644
index 00000000..5335f368
--- /dev/null
+++ b/R/irf.mvgam.R
@@ -0,0 +1,229 @@
+#' Calculate latent VAR impulse response functions
+#'
+#' Compute Generalized or Orthogonalized Impulse Response Functions (IRFs) from
+#' \code{mvgam} models with Vector Autoregressive dynamics
+#'
+#'@name irf.mvgam
+#'@param object \code{list} object of class \code{mvgam} resulting from a call to [mvgam()]
+#'that used a Vector Autoregressive latent process model (either as `VAR(cor = FALSE)` or
+#'`VAR(cor = TRUE)`)
+#'@param h Positive \code{integer} specifying the forecast horizon over which to calculate
+#'the IRF
+#'@param cumulative \code{Logical} flag indicating whether the IRF should be cumulative
+#'@param orthogonal \code{Logical} flag indicating whether orthogonalized IRFs should be
+#'calculated. Note that the order of the variables matters when calculating these
+#'@param ... ignored
+#'@details Generalized or Orthogonalized Impulse Response Functions can be computed
+#'using the posterior estimates of Vector Autoregressive parameters. This function
+#'generates a positive "shock" for a target process at time `t = 0` and then
+#'calculates how  each of the remaining processes in the latent VAR are expected
+#'to respond over the forecast horizon `h`. The function computes IRFs for all
+#'processes in the object and returns them in an array that can be plotted using
+#'the S3 `plot` function
+#'@return An object of class \code{mvgam_irf} containing the posterior IRFs. This
+#'object can be used with the supplied S3 functions \code{plot}
+#'@author Nicholas J Clark
+#'@seealso \code{\link{VAR}}, \code{\link{plot.mvgam_irf}}
+#' @examples
+#' \donttest{
+#' # Simulate some time series that follow a latent VAR(1) process
+#' simdat <- sim_mvgam(family = gaussian(),
+#'                     n_series = 4,
+#'                     trend_model = VAR(cor = TRUE),
+#'                     prop_trend = 1)
+#' plot_mvgam_series(data = simdat$data_train, series = 'all')
+#'
+#' # Fit a model that uses a latent VAR(1)
+#' mod <- mvgam(y ~ -1,
+#'              trend_formula = ~ 1,
+#'              trend_model = VAR(cor = TRUE),
+#'              family = gaussian(),
+#'              data = simdat$data_train,
+#'              silent = 2)
+#'
+#' # Calulate Generalized IRFs for each series
+#' irfs <- irf(mod, h = 12, cumulative = FALSE)
+#'
+#' # Plot them
+#' plot(irfs, series = 1)
+#' plot(irfs, series = 2)
+#' plot(irfs, series = 3)
+#' }
+#'@export
+irf <- function(object, ...){
+  UseMethod("irf", object)
+}
+
+#'@rdname irf.mvgam
+#'@method irf mvgam
+#'@export
+irf.mvgam <- function(object,
+                     h = 1,
+                     cumulative = FALSE,
+                     orthogonal = FALSE,
+                     ...){
+  validate_pos_integer(h)
+  trend_model <- attr(object$model_data, 'trend_model')
+  if(!trend_model %in% c('VAR', 'VARcor', 'VAR1', 'VAR1cor')){
+    stop('Only VAR(1) models currently supported for calculating IRFs',
+         call. = FALSE)
+  }
+  beta_vars <- mcmc_chains(object$model_output, 'A')
+  sigmas <- mcmc_chains(object$model_output, 'Sigma')
+  n_series <- object$n_lv
+  all_irfs <- lapply(seq_len(NROW(beta_vars)), function(draw){
+
+    # Get necessary VAR parameters into a simple list format
+    x <- list(K = n_series,
+              A = matrix(beta_vars[draw,],
+                         nrow = n_series,
+                         ncol = n_series,
+                         byrow = TRUE),
+              Sigma = matrix(sigmas[draw,],
+                             nrow = n_series,
+                             ncol = n_series,
+                             byrow = TRUE),
+              p = 1)
+
+    # Calculate the IRF
+    gen_irf(x, h = h, cumulative = cumulative, orthogonal = orthogonal)
+  })
+  class(all_irfs) <- 'mvgam_irf'
+  return(all_irfs)
+}
+
+#### Functions to compute Generalized Impulse Response functions
+# Much of this code is modified from R code generously provided by Clinton Watkins:
+# https://www.clintonwatkins.com/posts/2021-generalised-impulse-response-function-R/ ####
+
+#' Calculate impulse response functions
+#' @noRd
+gen_irf = function(x, h = 6, cumulative = TRUE, orthogonal = FALSE){
+
+  impulse <- paste0('process_', 1:x$K)
+
+  # Create arrays to hold calculations
+  IRF_o <- array(data = 0, dim = c(h, x$K, x$K),
+                 dimnames = list(NULL, impulse, impulse))
+  IRF_g <- array(data = 0, dim = c(h, x$K, x$K),
+                 dimnames = list(NULL, impulse, impulse))
+  IRF_g1 <- array(data = 0, dim = c(h, x$K, x$K))
+
+
+  # Estimation of orthogonalised or generalised IRFs
+  if(orthogonal){
+    var_ma <- var_psi(x, h)
+  } else {
+    var_ma <- var_phi(x, h)
+  }
+
+  sigma_u <- x$Sigma
+  P <- t(chol(sigma_u))
+  sig_jj <- diag(sigma_u)
+
+  for (jj in 1:x$K){
+    indx_ <- matrix(0, x$K, 1)
+    indx_[jj, 1] <- 1
+
+    for (kk in 1:h){
+      IRF_o[kk, ,jj] <- var_ma[, , kk] %*% P %*% indx_  # Peseran-Shin eqn 7 (OIRF)
+      IRF_g1[kk, ,jj] <- var_ma[, , kk] %*% sigma_u %*% indx_
+      IRF_g[kk, ,jj] <- sig_jj[jj] ^ (-0.5) * IRF_g1[kk, ,jj]  # Peseran-Shin eqn 10 (GIRF)
+
+    }
+  }
+
+  if(orthogonal == TRUE){
+    irf <- IRF_o
+  } else if(orthogonal == FALSE) {
+    irf <- IRF_g
+  } else {
+    stop("\nError! Orthogonalised or generalised IRF?\n")
+  }
+
+  idx <- length(impulse)
+  irs <- list()
+  for (ii in 1:idx) {
+    irs[[ii]] <- matrix(irf[1:(h), impulse, impulse[ii]], nrow = h)
+    colnames(irs[[ii]]) <- impulse
+    if (cumulative) {
+      if (length(impulse) > 1)
+        irs[[ii]] <- apply(irs[[ii]], 2, cumsum)
+      if (length(impulse) == 1) {
+        tmp <- matrix(cumsum(irs[[ii]]))
+        colnames(tmp) <- impulse
+        irs[[ii]] <- tmp
+      }
+    }
+  }
+  names(irs) <- impulse
+  result <- irs
+  return(irs)
+}
+
+#' Convert a VAR A matrix to its moving average representation
+#' @noRd
+var_phi = function(x, h = 10){
+  h <- abs(as.integer(h))
+  K <- x$K
+  p <- x$p
+  A <- as.array(x$A)
+  if(h >= p){
+    As <- array(0, dim = c(K, K, h + 1))
+    for(i in (p + 1):(h + 1)){
+      As[, , i] <- matrix(0, nrow = K, ncol = K)
+    }
+  } else {
+    As <- array(0, dim = c(K, K, p))
+  }
+  As[, , 1] <- A
+  Phi <- array(0, dim=c(K, K, h + 1))
+  Phi[, ,1] <- diag(K)
+  Phi[, , 2] <- Phi[, , 1] %*% As[, , 1]
+  if (h > 1) {
+    for (i in 3:(h + 1)) {
+      tmp1 <- Phi[, , 1] %*% As[, , i-1]
+      tmp2 <- matrix(0, nrow = K, ncol = K)
+      idx <- (i - 2):1
+      for (j in 1:(i - 2)) {
+        tmp2 <- tmp2 + Phi[, , j+1] %*% As[, , idx[j]]
+      }
+      Phi[, , i] <- tmp1 + tmp2
+    }
+  }
+  return(Phi)
+}
+
+#' Convert a VAR A matrix to its orthogonalised moving average representation
+#' @noRd
+var_psi = function(x, h=10){
+  h <- abs(as.integer(h))
+  Phi <- var_phi(x, h = h)
+  Psi <- array(0, dim = dim(Phi))
+  sigma_u <- x$Sigma
+  P <- t(chol(sigma_u))
+  dim3 <- dim(Phi)[3]
+  for(i in 1:dim3){
+    Psi[, , i] <- Phi[, , i] %*% P
+  }
+  return(Psi)
+}
+
+#' Forecast error decomposition
+#' @noRd
+var_fecov = function(x, h) {
+  sigma_u <- x$Sigma
+  sigma_yh <- array(NA, dim = c(x$K, x$K, h))
+  sigma_yh[, , 1] <- sigma_u
+  Phi <- var_phi(x, h = h)
+  if (h > 1) {
+    for (i in 2:h) {
+      temp <- matrix(0, nrow = x$K, ncol = x$K)
+      for (j in 2:i) {
+        temp <- temp + Phi[, , j] %*% sigma_u %*% t(Phi[, , j])
+      }
+      sigma_yh[, , i] <- temp + sigma_yh[, , 1]
+    }
+  }
+  return(sigma_yh)
+}
diff --git a/R/loo.mvgam.R b/R/loo.mvgam.R
index 799f3106..9e0a85be 100644
--- a/R/loo.mvgam.R
+++ b/R/loo.mvgam.R
@@ -144,6 +144,14 @@ loo_compare.mvgam <- function(x, ...,
   loo_compare(loos)
 }
 
+#' @export
+#' @importFrom loo loo
+loo::loo
+
+#' @export
+#' @importFrom loo loo_compare
+loo::loo_compare
+
 #'@noRd
 split_mod_dots = function (x, ..., model_names = NULL, other = TRUE) {
 
diff --git a/R/marginaleffects.mvgam.R b/R/marginaleffects.mvgam.R
index 8b7abc4d..abd8f154 100644
--- a/R/marginaleffects.mvgam.R
+++ b/R/marginaleffects.mvgam.R
@@ -21,6 +21,38 @@
 #' @author Nicholas J Clark
 NULL
 
+#' @export
+#' @importFrom marginaleffects predictions
+marginaleffects::predictions
+
+#' @export
+#' @importFrom marginaleffects plot_predictions
+marginaleffects::plot_predictions
+
+#' @export
+#' @importFrom marginaleffects slopes
+marginaleffects::slopes
+
+#' @export
+#' @importFrom marginaleffects plot_slopes
+marginaleffects::plot_slopes
+
+#' @export
+#' @importFrom marginaleffects comparisons
+marginaleffects::comparisons
+
+#' @export
+#' @importFrom marginaleffects plot_comparisons
+marginaleffects::plot_comparisons
+
+#' @export
+#' @importFrom marginaleffects datagrid
+marginaleffects::datagrid
+
+#' @export
+#' @importFrom marginaleffects hypotheses
+marginaleffects::hypotheses
+
 #' Functions needed for working with marginaleffects
 #' @rdname mvgam_marginaleffects
 #' @export
diff --git a/R/mcmc_plot.mvgam.R b/R/mcmc_plot.mvgam.R
index 66432bad..1631917e 100644
--- a/R/mcmc_plot.mvgam.R
+++ b/R/mcmc_plot.mvgam.R
@@ -2,7 +2,6 @@
 #'
 #' Convenient way to call MCMC plotting functions
 #' implemented in the \pkg{bayesplot} package
-#' @importFrom brms mcmc_plot
 #' @importFrom bayesplot color_scheme_set color_scheme_get
 #' @inheritParams brms::mcmc_plot
 #' @inheritParams as.data.frame.mvgam
@@ -119,3 +118,7 @@ mcmc_plot.mvgam = function(object,
   # Return the plot
   return(out_plot)
 }
+
+#' @export
+#' @importFrom brms mcmc_plot
+brms::mcmc_plot
diff --git a/R/mvgam_diagnostics.R b/R/mvgam_diagnostics.R
index 969b3b08..1ade4ac8 100644
--- a/R/mvgam_diagnostics.R
+++ b/R/mvgam_diagnostics.R
@@ -52,6 +52,7 @@ log_posterior.mvgam <- function(object, ...) {
 
 #' @rdname mvgam_diagnostics
 #' @importFrom posterior rhat
+#' @export rhat
 #' @export
 rhat.mvgam <- function(x, pars = NULL, ...) {
   # bayesplot uses outdated rhat code from rstan
@@ -100,31 +101,3 @@ neff_ratio.mvgam <- function(object, pars = NULL, ...) {
   names(ess) <- tmp$variable
   ess / brms::ndraws(draws)
 }
-
-#' @rdname mvgam_diagnostics
-#' @importFrom bayesplot neff_ratio
-#' @export neff_ratio
-#' @export
-neff_ratio.mvgam <- function(object, pars = NULL, ...) {
-  insight::check_if_installed("matrixStats",
-                              reason = 'to calculate effective sample sizes')
-  # bayesplot uses outdated ess code from rstan
-  # bayesplot::neff_ratio(object$fit, pars = pars, ...)
-  if(is.null(pars)){
-    vars_extract <- unlist(purrr::map(variables(object), 'orig_name'))
-    vars_extract <- vars_extract[-grep('ypred', vars_extract)]
-    draws <- as_draws_array(object,
-                            variable = vars_extract,
-                            use_alias = FALSE)
-  } else {
-    draws <- as_draws_array(object,
-                            variable = pars)
-  }
-  tmp <- posterior::summarise_draws(
-    draws, ess_bulk = posterior::ess_bulk, ess_tail = posterior::ess_tail
-  )
-  # min of ess_bulk and ess_tail mimics definition of posterior::rhat.default
-  ess <- matrixStats::rowMins(cbind(tmp$ess_bulk, tmp$ess_tail))
-  names(ess) <- tmp$variable
-  ess / ndraws(draws)
-}
diff --git a/R/mvgam_formulae.R b/R/mvgam_formulae.R
index 00429e09..11651b90 100644
--- a/R/mvgam_formulae.R
+++ b/R/mvgam_formulae.R
@@ -36,7 +36,28 @@
 #' \code{\link[mgcv]{jagam}},
 #' \code{\link[mgcv]{gam}},
 #' \code{\link[mgcv]{s}},
+#' \code{\link[brms]{gp}},
 #' \code{\link[stats]{formula}}
 #' @author Nicholas J Clark
 #' @name mvgam_formulae
 NULL
+
+#' @export
+#' @importFrom brms gp
+brms::gp
+
+#' @export
+#' @importFrom mgcv s
+mgcv::s
+
+#' @export
+#' @importFrom mgcv te
+mgcv::te
+
+#' @export
+#' @importFrom mgcv ti
+mgcv::ti
+
+#' @export
+#' @importFrom mgcv t2
+mgcv::t2
diff --git a/R/mvgam_irf-class.R b/R/mvgam_irf-class.R
new file mode 100644
index 00000000..f904559f
--- /dev/null
+++ b/R/mvgam_irf-class.R
@@ -0,0 +1,116 @@
+#' `mvgam_irf` object description
+#'
+#' A \code{mvgam_irf} object returned by function \code{\link{irf}}.
+#' Run `methods(class = "mvgam_irf")` to see an overview of available methods.
+#' @details A `mvgam_irf` object contains a list of posterior IRFs, each stored as
+#' its own list
+#' @seealso [mvgam], [VAR]
+#' @author Nicholas J Clark
+#' @name mvgam_irf-class
+NULL
+
+#'Plot impulse responses from an mvgam_irf object
+
+#'This function takes an \code{mvgam_irf} object and produces plots of Impulse Response Functions
+#'
+#'@param x \code{list} object of class \code{mvgam_irf}. See [irf()]
+#'@param series \code{integer} specifying which process series should be given the shock
+#'@param ... ignored
+#'@return A base R plot or set of plots
+#'@author Nicholas J Clark
+#'@export
+plot.mvgam_irf = function(x, series = 1, ...){
+
+  all_irfs <- x
+  n_processes <- dim(all_irfs[[1]][[1]])[2]
+  h <- dim(all_irfs[[1]][[1]])[1]
+
+  # Extract IRFs for the specific series
+  impulse_responses <- lapply(seq_along(all_irfs), function(j){
+    all_irfs[[j]][series]
+  })
+
+  # Plot imp responses for all series apart from the impacted one
+  c_light <- c("#DCBCBC")
+  c_light_highlight <- c("#C79999")
+  c_mid <- c("#B97C7C")
+  c_mid_highlight <- c("#A25050")
+  c_dark <- c("#8F2727")
+  c_dark_highlight <- c("#7C0000")
+
+
+  n_plots <- n_processes
+  pages <- 1
+
+  if (n_plots > 4) pages <- 2
+  if (n_plots > 8) pages <- 3
+  if (n_plots > 12) pages <- 4
+  if (pages != 0)  {
+    ppp <- n_plots %/% pages
+
+    if (n_plots %% pages != 0) {
+      ppp <- ppp + 1
+      while (ppp * (pages - 1) >= n_plots) pages <- pages - 1
+    }
+
+    # Configure layout matrix
+    c <- r <- trunc(sqrt(ppp))
+    if (c<1) r <- c <- 1
+    if (c*r < ppp) c <- c + 1
+    if (c*r < ppp) r <- r + 1
+
+    .pardefault <- par(no.readonly=T)
+    on.exit(par(.pardefault))
+    oldpar <- par(mfrow = c(r,c))
+
+  } else { ppp <- 1; oldpar <- par()}
+
+  for(x in 1:n_plots){
+    if(x == series){
+      plot(x = 1:h, y = 1:h, type = "n", bty = 'L',
+           ylim = c(-2, 2),
+           yaxt = 'n',
+           xaxt = 'n',
+           ylab = '',
+           xlab = '')
+      abline(h = 0, lwd = 2.5)
+      box(bty = 'l', lwd = 2)
+
+    } else {
+      responses <- do.call(rbind, lapply(seq_along(impulse_responses), function(j){
+        impulse_responses[[j]][[1]][,x]
+      }))
+
+      probs = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
+      cred <- sapply(1:NCOL(responses),
+                     function(n) quantile(responses[,n],
+                                          probs = probs, na.rm = TRUE))
+      pred_vals <- 1:h
+      plot(1, type = "n", bty = 'L',
+           xlab = '',
+           xaxt = 'n',
+           ylab = '',
+           xlim = c(1, h),
+           ylim = c(-1.1 * max(abs(cred)), 1.1 * max(abs(cred))))
+      polygon(c(pred_vals, rev(pred_vals)), c(cred[1,], rev(cred[9,])),
+              col = c_light, border = NA)
+      polygon(c(pred_vals, rev(pred_vals)), c(cred[2,], rev(cred[8,])),
+              col = c_light_highlight, border = NA)
+      polygon(c(pred_vals, rev(pred_vals)), c(cred[3,], rev(cred[7,])),
+              col = c_mid, border = NA)
+      polygon(c(pred_vals, rev(pred_vals)), c(cred[4,], rev(cred[6,])),
+              col = c_mid_highlight, border = NA)
+      lines(pred_vals, cred[5,], col = c_dark, lwd = 2.5)
+      abline(h = 0, lwd = 3, col = 'white')
+      abline(h = 0, lwd = 2.5)
+      box(bty = 'L', lwd = 2)
+    }
+
+    if(x != series){
+      axis(1, cex.axis = 1, tck = -0.05)
+      title(paste0('process_', series, ' -> process_', x), line = 0)
+    }
+
+  }
+  layout(1)
+}
diff --git a/R/posterior_epred.mvgam.R b/R/posterior_epred.mvgam.R
index 56ea36ac..c314c6a8 100644
--- a/R/posterior_epred.mvgam.R
+++ b/R/posterior_epred.mvgam.R
@@ -10,7 +10,6 @@
 #' incorporated in the draws computed by \code{posterior_epred} while the
 #' residual error is ignored there. However, the estimated means of both methods
 #' averaged across draws should be very similar.
-#' @importFrom rstantools posterior_epred
 #' @inheritParams predict.mvgam
 #' @param ndraws Positive `integer` indicating how many posterior draws should be used.
 #' If `NULL` (the default) all draws are used.
@@ -37,7 +36,6 @@
 #' where \code{n_samples} is the number of posterior samples from the fitted object
 #' and \code{n_obs} is the number of observations in \code{newdata}
 #' @seealso \code{\link{hindcast.mvgam}} \code{\link{posterior_linpred.mvgam}} \code{\link{posterior_predict.mvgam}}
-#' @aliases posterior_epred
 #' @examples
 #' \donttest{
 #' # Simulate some data and fit a model
@@ -85,7 +83,6 @@ posterior_epred.mvgam = function(object,
 #' Compute posterior draws of the linear predictor, that is draws before
 #' applying any link functions or other transformations. Can be performed for
 #' the data used to fit the model (posterior predictive checks) or for new data.
-#' @importFrom rstantools posterior_linpred
 #' @inheritParams posterior_epred.mvgam
 #' @param transform Logical; if \code{FALSE}
 #'  (the default), draws of the linear predictor are returned.
@@ -166,7 +163,6 @@ posterior_linpred.mvgam = function(object,
 #' \code{\link{posterior_epred.mvgam}}. This is because the residual error
 #' is incorporated in \code{posterior_predict}. However, the estimated means of
 #' both methods averaged across draws should be very similar.
-#' @importFrom rstantools posterior_predict
 #' @inheritParams predict.mvgam
 #' @inheritParams posterior_epred.mvgam
 #' @param process_error Logical. If \code{TRUE} and \code{newdata} is supplied,
@@ -229,6 +225,19 @@ posterior_predict.mvgam = function(object,
   return(out)
 }
 
+#' @export
+#' @importFrom rstantools posterior_predict
+rstantools::posterior_predict
+
+#' @export
+#' @importFrom rstantools posterior_epred
+rstantools::posterior_epred
+
+#' @export
+#' @importFrom rstantools posterior_linpred
+rstantools::posterior_linpred
+
+
 #' Expected Values of the Posterior Predictive Distribution
 #'
 #' This method extracts posterior estimates of the fitted values
diff --git a/R/ppc.mvgam.R b/R/ppc.mvgam.R
index 7e4b9d8c..d2c2ce23 100644
--- a/R/ppc.mvgam.R
+++ b/R/ppc.mvgam.R
@@ -772,7 +772,7 @@ ppc.mvgam = function(object, newdata, data_test, series = 1, type = 'hist',
 #' predictions are generated for the original observations used for the model fit.
 #' @param ... Further arguments passed to \code{\link{predict.mvgam}}
 #'   as well as to the PPC function specified in \code{type}.
-#' @inheritParams prepare_predictions.brmsfit
+#' @inheritParams brms::prepare_predictions.brmsfit
 #'
 #' @return A ggplot object that can be further
 #'  customized using the \pkg{ggplot2} package.
diff --git a/R/stan_utils.R b/R/stan_utils.R
index d5b9238c..f57acbb3 100644
--- a/R/stan_utils.R
+++ b/R/stan_utils.R
@@ -41,6 +41,9 @@ print.mvgammodel = function(x, ...){
 
 #' @export
 #' @importFrom brms stancode
+brms::stancode
+
+#' @export
 #' @param ... ignored
 #' @rdname code
 stancode.mvgam_prefit = function(object, ...){
@@ -57,6 +60,9 @@ stancode.mvgam = function(object, ...){
 
 #' @export
 #' @importFrom brms standata
+brms::standata
+
+#' @export
 #' @param ... ignored
 #' @rdname code
 standata.mvgam_prefit = function(object, ...){
diff --git a/R/ti.R b/R/ti.R
deleted file mode 100644
index 086cb643..00000000
--- a/R/ti.R
+++ /dev/null
@@ -1,48 +0,0 @@
-# This file contains functions dealing with the extended
-# formula syntax to specify smooth terms via mgcv
-
-#' Defining smooths in \pkg{mvgam} formulae
-#'
-#' Functions used in definition of smooth terms within model formulae.
-#' The functions do not evaluate a (spline) smooth - they exist purely
-#' to help set up mvgam models using spline based smooths.
-#'
-#' @param ... Arguments passed to \code{\link[mgcv:ti]{mgcv::ti}} or
-#'  \code{\link[mgcv:ti]{mgcv::te}}
-#'
-#' @details The functions defined here are just simple wrappers of the respective
-#'   functions of the \pkg{mgcv} package. When using them, please cite the
-#'   appropriate references obtained via \code{citation("mgcv")}.
-#'
-#' @seealso \code{\link[mgcv:ti]{mgcv::ti}}, \code{\link[mgcv:ti]{mgcv::te}}
-#'
-#'
-#' @examples
-#' \donttest{
-#' # Simulate some data
-#' dat <- mgcv::gamSim(1, n = 200, scale = 2)
-#'
-#' # Fit univariate smooths for all predictors
-#' fit1 <- mvgam(y ~ s(x0) + s(x1) + s(x2) + s(x3),
-#'               data = dat, chains = 2, family = gaussian())
-#' summary(fit1)
-#' conditional_effects(fit1)
-#'
-#' # Fit a more complicated smooth model
-#' fit2 <- mvgam(y ~ te(x0, x1) + s(x2, by = x3),
-#'               data = dat, chains = 2, family = gaussian())
-#' summary(fit2)
-#' conditional_effects(fit2)
-#' }
-#'
-#' @rdname ti
-#' @export
-ti <- function(...) {
-  mgcv::ti(...)
-}
-
-#' @rdname ti
-#' @export
-te <- function(...) {
-  mgcv::te(...)
-}
diff --git a/inst/doc/data_in_mvgam.Rmd b/inst/doc/data_in_mvgam.Rmd
index a704cb38..b8240edd 100644
--- a/inst/doc/data_in_mvgam.Rmd
+++ b/inst/doc/data_in_mvgam.Rmd
@@ -24,7 +24,7 @@ knitr::opts_chunk$set(
 ```{r setup, include=FALSE}
 knitr::opts_chunk$set(
   echo = TRUE,   
-  dpi = 100,
+  dpi = 150,
   fig.asp = 0.8,
   fig.width = 6,
   out.width = "60%",
diff --git a/inst/doc/mvgam_overview.Rmd b/inst/doc/mvgam_overview.Rmd
index c3e2602f..508c63d3 100644
--- a/inst/doc/mvgam_overview.Rmd
+++ b/inst/doc/mvgam_overview.Rmd
@@ -23,7 +23,7 @@ knitr::opts_chunk$set(
 ```{r setup, include=FALSE}
 knitr::opts_chunk$set(
   echo = TRUE,   
-  dpi = 100,
+  dpi = 150,
   fig.asp = 0.8,
   fig.width = 6,
   out.width = "60%",
@@ -146,7 +146,7 @@ z_{i,t} & \sim \text{Normal}(ar1_i^{distance} * z_{i,t-1}, \sigma_i) \end{align*
 Where $distance$ is a vector of non-negative measurements of the time differences between successive observations. See the **Examples** section in `?CAR` for an illustration of how to set these models up. 
 
 ## Regression formulae
-`mvgam` supports an observation model regression formula, built off the `mgcv` package, as well as an optional process model regression formula. The formulae supplied to \code{\link{mvgam}} are exactly like those supplied to `glm()` except that smooth terms, `s()`,
+`mvgam` supports an observation model regression formula, built off the `mvgcv` package, as well as an optional process model regression formula. The formulae supplied to \code{\link{mvgam}} are exactly like those supplied to `glm()` except that smooth terms, `s()`,
 `te()`, `ti()` and `t2()`, time-varying effects using `dynamic()`, monotonically increasing (using `s(x, bs = 'moi')`) or decreasing splines (using `s(x, bs = 'mod')`; see `?smooth.construct.moi.smooth.spec` for details), as well as Gaussian Process functions using `gp()`, can be added to the right hand side (and `.` is not supported in `mvgam` formulae). See `?mvgam_formulae` for more guidance.
   
 For setting up State-Space models, the optional process model formula can be used (see [the State-Space model vignette](https://nicholasjclark.github.io/mvgam/articles/trend_formulas.html) and [the shared latent states vignette](https://nicholasjclark.github.io/mvgam/articles/trend_formulas.html) for guidance on using trend formulae).
@@ -223,7 +223,15 @@ You can also summarize multiple variables, which is helpful to search for data r
 summary(model_data)
 ```
 
-We have some `NA`s in our response variable `count`. These observations will generally be thrown out by most modelling packages in \R. But as you will see when we work through the tutorials, `mvgam` keeps these in the data so that predictions can be automatically returned for the full dataset. The time series and some of its descriptive features can be plotted using `plot_mvgam_series()`:
+We have some `NA`s in our response variable `count`. Let's visualize the data as a heatmap to get a sense of where these are distributed (`NA`s are shown as red bars in the below plot)
+```{r}
+image(is.na(t(model_data %>%
+                dplyr::arrange(dplyr::desc(time)))), axes = F,
+      col = c('grey80', 'darkred'))
+axis(3, at = seq(0,1, len = NCOL(model_data)), labels = colnames(model_data))
+```
+
+These observations will generally be thrown out by most modelling packages in \R. But as you will see when we work through the tutorials, `mvgam` keeps these in the data so that predictions can be automatically returned for the full dataset. The time series and some of its descriptive features can be plotted using `plot_mvgam_series()`:
 ```{r}
 plot_mvgam_series(data = model_data, series = 1, y = 'count')
 ```
@@ -306,6 +314,7 @@ mcmc_plot(object = model1,
 We can also use the wide range of posterior checking functions available in `bayesplot` (see `?mvgam::ppc_check.mvgam` for details):
 ```{r}
 pp_check(object = model1)
+pp_check(model1, type = "rootogram")
 ```
 
 There is clearly some variation in these yearly intercept estimates. But how do these translate into time-varying predictions? To understand this, we can plot posterior hindcasts from this model for the training period using `plot.mvgam` with `type = 'forecast'`
@@ -325,6 +334,13 @@ hc <- hindcast(model1, type = 'link')
 range(hc$hindcasts$PP)
 ```
 
+Objects of class `mvgam_forecast` have an associated plot function as well:
+```{r Plot hindcasts on the linear predictor scale}
+plot(hc)
+```
+
+This plot can look a bit confusing as it seems like there is linear interpolation from the end of one year to the start of the next. But this is just due to the way the lines are automatically connected in base \R plots  
+  
 In any regression analysis, a key question is whether the residuals show any patterns that can be indicative of un-modelled sources of variation. For GLMs, we can use a modified residual called the [Dunn-Smyth, or randomized quantile, residual](https://www.jstor.org/stable/1390802){target="_blank"}. Inspect Dunn-Smyth residuals from the model using `plot.mvgam` with `type = 'residuals'`
 ```{r Plot posterior residuals}
 plot(model1, type = 'residuals')
@@ -349,12 +365,22 @@ model1b <- mvgam(count ~ s(year_fac, bs = 're') - 1,
 
 ```{r eval=FALSE}
 model1b <- mvgam(count ~ s(year_fac, bs = 're') - 1,
-                 family = poisson(),
-                 data = data_train,
-                 newdata = data_test)
+                family = poisson(),
+                data = data_train,
+                newdata = data_test)
 ```
 
-We can view the test data in the forecast plot to see that the predictions do not capture the temporal variation in the test set
+Repeating the plots above gives insight into how the model's hierarchical prior formulation provides all the structure needed to sample values for un-modelled years
+```{r}
+plot(model1b, type = 're')
+```
+
+
+```{r}
+plot(model1b, type = 'forecast')
+```
+
+We can also view the test data in the forecast plot to see that the predictions do not capture the temporal variation in the test set
 ```{r Plotting predictions against test data}
 plot(model1b, type = 'forecast', newdata = data_test)
 ```
@@ -402,7 +428,12 @@ Rather than printing the summary each time, we can also quickly look at the post
 coef(model2)
 ```
 
-Look at the estimated effect of `ndvi` using using a histogram. This can be done by first extracting the posterior coefficients:
+Look at the estimated effect of `ndvi` using `plot.mvgam` with `type = 'pterms'`
+```{r Plot NDVI effect}
+plot(model2, type = 'pterms')
+```
+
+This plot indicates a positive linear effect of `ndvi` on `log(counts)`. But it may be easier to visualise using a histogram, especially for parametric (linear) effects. This can be done by first extracting the posterior coefficients as we did in the first example:
 ```{r}
 beta_post <- as.data.frame(model2, variable = 'betas')
 dplyr::glimpse(beta_post)
@@ -424,9 +455,19 @@ abline(v = 0, lwd = 2.5)
 ```
 
 ### `marginaleffects` support
-Given our model used a nonlinear link function (log link in this example), it can still be difficult to fully understand what relationship our model is estimating between a predictor and the response. Fortunately, the `marginaleffects` package makes this relatively straightforward. Objects of class `mvgam` can be used with `marginaleffects` to inspect contrasts, scenario-based predictions, conditional and marginal effects, all on the outcome scale. Like `brms`, `mvgam` has the simple `conditional_effects` function to make quick and informative plots for main effects, which rely on `marginaleffects` support. This will likely be your go-to function for quickly understanding patterns from fitted `mvgam` models
+Given our model used a nonlinear link function (log link in this example), it can still be difficult to fully understand what relationship our model is estimating between a predictor and the response. Fortunately, the `marginaleffects` package makes this relatively straightforward. Objects of class `mvgam` can be used with `marginaleffects` to inspect contrasts, scenario-based predictions, conditional and marginal effects, all on the outcome scale. Here we will use the `plot_predictions` function from `marginaleffects` to inspect the conditional effect of `ndvi` (use `?plot_predictions` for guidance on how to modify these plots):
 ```{r warning=FALSE}
-conditional_effects(model2)
+plot_predictions(model2, 
+                 condition = "ndvi",
+                 # include the observed count values
+                 # as points, and show rugs for the observed
+                 # ndvi and count values on the axes
+                 points = 0.5, rug = TRUE)
+```
+
+Now it is easier to get a sense of the nonlinear but positive relationship estimated between `ndvi` and `count`. Like `brms`, `mvgam` has the simple `conditional_effects` function to make quick and informative plots for main effects. This will likely be your go-to function for quickly understanding patterns from fitted `mvgam` models
+```{r warning=FALSE}
+plot(conditional_effects(model2), ask = FALSE)
 ```
 
 ## Adding predictors as smooths
@@ -463,9 +504,33 @@ Where the smooth function $f_{time}$ is built by summing across a set of weighte
 summary(model3)
 ```
 
-The summary above now contains posterior estimates for the smoothing parameters as well as the basis coefficients for the nonlinear effect of `time`. We can visualize `conditional_effects` as before:
+The summary above now contains posterior estimates for the smoothing parameters as well as the basis coefficients for the nonlinear effect of `time`. We can visualize the conditional `time` effect using the `plot` function with `type = 'smooths'`:
+```{r}
+plot(model3, type = 'smooths')
+```
+
+By default this plots shows posterior empirical quantiles, but it can also be helpful to view some realizations of the underlying function (here, each line is a different potential curve drawn from the posterior of all possible curves):
+```{r}
+plot(model3, type = 'smooths', realisations = TRUE,
+     n_realisations = 30)
+```
+
+### Derivatives of smooths
+A useful question when modelling using GAMs is to identify where the function is changing most rapidly. To address this, we can plot estimated 1st derivatives of the spline:
+```{r Plot smooth term derivatives, warning = FALSE, fig.asp = 1}
+plot(model3, type = 'smooths', derivatives = TRUE)
+```
+
+Here, values above `>0` indicate the function was increasing at that time point, while values `<0` indicate the function was declining. The most rapid declines appear to have been happening around timepoints 50 and again toward the end of the training period, for example.
+  
+Use `conditional_effects` again for useful plots on the outcome scale:
 ```{r warning=FALSE}
-conditional_effects(model3, type = 'link')
+plot(conditional_effects(model3), ask = FALSE)
+```
+
+Or on the link scale:
+```{r warning=FALSE}
+plot(conditional_effects(model3, type = 'link'), ask = FALSE)
 ```
 
 Inspect the underlying `Stan` code to gain some idea of how the spline is being penalized:
@@ -526,6 +591,13 @@ Here the term $z_t$ captures autocorrelated latent residuals, which are modelled
 summary(model4)
 ```
 
+View conditional smooths for the `ndvi` effect:
+```{r warning=FALSE, message=FALSE}
+plot_predictions(model4, 
+                 condition = "ndvi",
+                 points = 0.5, rug = TRUE)
+```
+
 View posterior hindcasts / forecasts and compare against the out of sample test data
 ```{r}
 plot(model4, type = 'forecast', newdata = data_test)
diff --git a/inst/doc/nmixtures.Rmd b/inst/doc/nmixtures.Rmd
index f7585ab1..62d4a08e 100644
--- a/inst/doc/nmixtures.Rmd
+++ b/inst/doc/nmixtures.Rmd
@@ -23,7 +23,7 @@ knitr::opts_chunk$set(
 ```{r setup, include=FALSE}
 knitr::opts_chunk$set(
   echo = TRUE,
-  dpi = 100,
+  dpi = 150,
   fig.asp = 0.8,
   fig.width = 6,
   out.width = "60%",
diff --git a/man/irf.mvgam.Rd b/man/irf.mvgam.Rd
new file mode 100644
index 00000000..14436390
--- /dev/null
+++ b/man/irf.mvgam.Rd
@@ -0,0 +1,75 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/irf.mvgam.R
+\name{irf.mvgam}
+\alias{irf.mvgam}
+\alias{irf}
+\title{Calculate latent VAR impulse response functions}
+\usage{
+irf(object, ...)
+
+\method{irf}{mvgam}(object, h = 1, cumulative = FALSE, orthogonal = FALSE, ...)
+}
+\arguments{
+\item{object}{\code{list} object of class \code{mvgam} resulting from a call to \code{\link[=mvgam]{mvgam()}}
+that used a Vector Autoregressive latent process model (either as \code{VAR(cor = FALSE)} or
+\code{VAR(cor = TRUE)})}
+
+\item{...}{ignored}
+
+\item{h}{Positive \code{integer} specifying the forecast horizon over which to calculate
+the IRF}
+
+\item{cumulative}{\code{Logical} flag indicating whether the IRF should be cumulative}
+
+\item{orthogonal}{\code{Logical} flag indicating whether orthogonalized IRFs should be
+calculated. Note that the order of the variables matters when calculating these}
+}
+\value{
+An object of class \code{mvgam_irf} containing the posterior IRFs. This
+object can be used with the supplied S3 functions \code{plot}
+}
+\description{
+Compute Generalized or Orthogonalized Impulse Response Functions (IRFs) from
+\code{mvgam} models with Vector Autoregressive dynamics
+}
+\details{
+Generalized or Orthogonalized Impulse Response Functions can be computed
+using the posterior estimates of Vector Autoregressive parameters. This function
+generates a positive "shock" for a target process at time \code{t = 0} and then
+calculates how  each of the remaining processes in the latent VAR are expected
+to respond over the forecast horizon \code{h}. The function computes IRFs for all
+processes in the object and returns them in an array that can be plotted using
+the S3 \code{plot} function
+}
+\examples{
+\donttest{
+# Simulate some time series that follow a latent VAR(1) process
+simdat <- sim_mvgam(family = gaussian(),
+                    n_series = 4,
+                    trend_model = VAR(cor = TRUE),
+                    prop_trend = 1)
+plot_mvgam_series(data = simdat$data_train, series = 'all')
+
+# Fit a model that uses a latent VAR(1)
+mod <- mvgam(y ~ -1,
+             trend_formula = ~ 1,
+             trend_model = VAR(cor = TRUE),
+             family = gaussian(),
+             data = simdat$data_train,
+             silent = 2)
+
+# Calulate Generalized IRFs for each series
+irfs <- irf(mod, h = 12, cumulative = FALSE)
+
+# Plot them
+plot(irfs, series = 1)
+plot(irfs, series = 2)
+plot(irfs, series = 3)
+}
+}
+\seealso{
+\code{\link{VAR}}, \code{\link{plot.mvgam_irf}}
+}
+\author{
+Nicholas J Clark
+}
diff --git a/man/mvgam_diagnostics.Rd b/man/mvgam_diagnostics.Rd
index e483b900..7fe9c1d7 100644
--- a/man/mvgam_diagnostics.Rd
+++ b/man/mvgam_diagnostics.Rd
@@ -17,8 +17,6 @@
 
 \method{rhat}{mvgam}(x, pars = NULL, ...)
 
-\method{neff_ratio}{mvgam}(object, pars = NULL, ...)
-
 \method{neff_ratio}{mvgam}(object, pars = NULL, ...)
 }
 \arguments{
diff --git a/man/mvgam_families.Rd b/man/mvgam_families.Rd
index 6d993699..1e0a439f 100644
--- a/man/mvgam_families.Rd
+++ b/man/mvgam_families.Rd
@@ -6,6 +6,10 @@
 \alias{student_t}
 \alias{betar}
 \alias{nb}
+\alias{lognormal}
+\alias{student}
+\alias{bernoulli}
+\alias{beta_binomial}
 \alias{nmix}
 \title{Supported mvgam families}
 \usage{
@@ -17,6 +21,14 @@ betar(...)
 
 nb(...)
 
+lognormal(...)
+
+student(...)
+
+bernoulli(...)
+
+beta_binomial(...)
+
 nmix(link = "log")
 }
 \arguments{
diff --git a/man/mvgam_formulae.Rd b/man/mvgam_formulae.Rd
index 5926d28b..1b27c3a6 100644
--- a/man/mvgam_formulae.Rd
+++ b/man/mvgam_formulae.Rd
@@ -46,6 +46,7 @@ extensive documentation for the \code{mgcv} package.
 \code{\link[mgcv]{jagam}},
 \code{\link[mgcv]{gam}},
 \code{\link[mgcv]{s}},
+\code{\link[brms]{gp}},
 \code{\link[stats]{formula}}
 }
 \author{
diff --git a/man/mvgam_irf-class.Rd b/man/mvgam_irf-class.Rd
new file mode 100644
index 00000000..e70f46fd
--- /dev/null
+++ b/man/mvgam_irf-class.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mvgam_irf-class.R
+\name{mvgam_irf-class}
+\alias{mvgam_irf-class}
+\title{\code{mvgam_irf} object description}
+\description{
+A \code{mvgam_irf} object returned by function \code{\link{irf}}.
+Run \code{methods(class = "mvgam_irf")} to see an overview of available methods.
+}
+\details{
+A \code{mvgam_irf} object contains a list of posterior IRFs, each stored as
+its own list
+}
+\seealso{
+\link{mvgam}, \link{VAR}
+}
+\author{
+Nicholas J Clark
+}
diff --git a/man/mvgam_marginaleffects.Rd b/man/mvgam_marginaleffects.Rd
index 16e6648e..5564cfdb 100644
--- a/man/mvgam_marginaleffects.Rd
+++ b/man/mvgam_marginaleffects.Rd
@@ -86,8 +86,6 @@ arguments.}
 \item \code{newdata = datagrid(cyl = c(4, 6))}: \code{cyl} variable equal to 4 and 6 and other regressors fixed at their means or modes.
 \item See the Examples section and the \code{\link[marginaleffects:datagrid]{datagrid()}} documentation.
 }
-\item \code{\link[=subset]{subset()}} call with a single argument to select a subset of the dataset used to fit the model, ex: \code{newdata = subset(treatment == 1)}
-\item \code{\link[dplyr:filter]{dplyr::filter()}} call with a single argument to select a subset of the dataset used to fit the model, ex: \code{newdata = filter(treatment == 1)}
 \item string:
 \itemize{
 \item "mean": Marginal Effects at the Mean. Slopes when each predictor is held at its mean or mode.
diff --git a/man/plot.mvgam_irf.Rd b/man/plot.mvgam_irf.Rd
new file mode 100644
index 00000000..2eaaffb1
--- /dev/null
+++ b/man/plot.mvgam_irf.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mvgam_irf-class.R
+\name{plot.mvgam_irf}
+\alias{plot.mvgam_irf}
+\title{Plot impulse responses from an mvgam_irf object
+This function takes an \code{mvgam_irf} object and produces plots of Impulse Response Functions}
+\usage{
+\method{plot}{mvgam_irf}(x, series = 1, ...)
+}
+\arguments{
+\item{x}{\code{list} object of class \code{mvgam_irf}. See \code{\link[=irf]{irf()}}}
+
+\item{series}{\code{integer} specifying which process series should be given the shock}
+
+\item{...}{ignored}
+}
+\value{
+A base R plot or set of plots
+}
+\description{
+Plot impulse responses from an mvgam_irf object
+This function takes an \code{mvgam_irf} object and produces plots of Impulse Response Functions
+}
+\author{
+Nicholas J Clark
+}
diff --git a/man/posterior_epred.mvgam.Rd b/man/posterior_epred.mvgam.Rd
index cbf354ce..4c8f7fee 100644
--- a/man/posterior_epred.mvgam.Rd
+++ b/man/posterior_epred.mvgam.Rd
@@ -2,7 +2,6 @@
 % Please edit documentation in R/posterior_epred.mvgam.R
 \name{posterior_epred.mvgam}
 \alias{posterior_epred.mvgam}
-\alias{posterior_epred}
 \title{Draws from the Expected Value of the Posterior Predictive Distribution}
 \usage{
 \method{posterior_epred}{mvgam}(
diff --git a/man/reexports.Rd b/man/reexports.Rd
new file mode 100644
index 00000000..87ce027e
--- /dev/null
+++ b/man/reexports.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/conditional_effects.R, R/get_mvgam_priors.R,
+%   R/loo.mvgam.R, R/marginaleffects.mvgam.R, R/mcmc_plot.mvgam.R,
+%   R/mvgam_formulae.R, R/posterior_epred.mvgam.R, R/stan_utils.R
+\docType{import}
+\name{reexports}
+\alias{reexports}
+\alias{conditional_effects}
+\alias{prior}
+\alias{prior_}
+\alias{set_prior}
+\alias{prior_string}
+\alias{loo}
+\alias{loo_compare}
+\alias{predictions}
+\alias{plot_predictions}
+\alias{slopes}
+\alias{plot_slopes}
+\alias{comparisons}
+\alias{plot_comparisons}
+\alias{datagrid}
+\alias{hypotheses}
+\alias{mcmc_plot}
+\alias{gp}
+\alias{s}
+\alias{te}
+\alias{ti}
+\alias{t2}
+\alias{posterior_predict}
+\alias{posterior_epred}
+\alias{posterior_linpred}
+\alias{stancode}
+\alias{standata}
+\title{Objects exported from other packages}
+\keyword{internal}
+\description{
+These objects are imported from other packages. Follow the links
+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{loo}{\code{\link[loo]{loo}}, \code{\link[loo]{loo_compare}}}
+
+  \item{marginaleffects}{\code{\link[marginaleffects]{comparisons}}, \code{\link[marginaleffects]{datagrid}}, \code{\link[marginaleffects]{hypotheses}}, \code{\link[marginaleffects]{plot_comparisons}}, \code{\link[marginaleffects]{plot_predictions}}, \code{\link[marginaleffects]{plot_slopes}}, \code{\link[marginaleffects]{predictions}}, \code{\link[marginaleffects]{slopes}}}
+
+  \item{mgcv}{\code{\link[mgcv]{s}}, \code{\link[mgcv]{t2}}, \code{\link[mgcv]{te}}, \code{\link[mgcv:te]{ti}}}
+
+  \item{rstantools}{\code{\link[rstantools]{posterior_epred}}, \code{\link[rstantools]{posterior_linpred}}, \code{\link[rstantools]{posterior_predict}}}
+}}
+
diff --git a/man/ti.Rd b/man/ti.Rd
deleted file mode 100644
index 6631b3e0..00000000
--- a/man/ti.Rd
+++ /dev/null
@@ -1,47 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/ti.R
-\name{ti}
-\alias{ti}
-\alias{te}
-\title{Defining smooths in \pkg{mvgam} formulae}
-\usage{
-ti(...)
-
-te(...)
-}
-\arguments{
-\item{...}{Arguments passed to \code{\link[mgcv:ti]{mgcv::ti}} or
-\code{\link[mgcv:ti]{mgcv::te}}}
-}
-\description{
-Functions used in definition of smooth terms within model formulae.
-The functions do not evaluate a (spline) smooth - they exist purely
-to help set up mvgam models using spline based smooths.
-}
-\details{
-The functions defined here are just simple wrappers of the respective
-functions of the \pkg{mgcv} package. When using them, please cite the
-appropriate references obtained via \code{citation("mgcv")}.
-}
-\examples{
-\donttest{
-# Simulate some data
-dat <- mgcv::gamSim(1, n = 200, scale = 2)
-
-# Fit univariate smooths for all predictors
-fit1 <- mvgam(y ~ s(x0) + s(x1) + s(x2) + s(x3),
-              data = dat, chains = 2, family = gaussian())
-summary(fit1)
-conditional_effects(fit1)
-
-# Fit a more complicated smooth model
-fit2 <- mvgam(y ~ te(x0, x1) + s(x2, by = x3),
-              data = dat, chains = 2, family = gaussian())
-summary(fit2)
-conditional_effects(fit2)
-}
-
-}
-\seealso{
-\code{\link[mgcv:ti]{mgcv::ti}}, \code{\link[mgcv:ti]{mgcv::te}}
-}
diff --git a/src/mvgam.dll b/src/mvgam.dll
index 119225b3..6180c635 100644
Binary files a/src/mvgam.dll and b/src/mvgam.dll differ
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
index 3046445b..843b188b 100644
Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ