diff --git a/DESCRIPTION b/DESCRIPTION
index 3cb84a50..c96db513 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: mvgam
 Title: Multivariate (Dynamic) Generalized Additive Models
-Version: 1.0.5
-Date: 2023-06-03
+Version: 1.0.6
+Date: 2023-10-06
 Authors@R:
     person("Nicholas J", "Clark", , "nicholas.j.clark1214@gmail.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-7131-3301"))
diff --git a/NAMESPACE b/NAMESPACE
index 3c8ad740..d5a027d1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -78,6 +78,7 @@ export(plot_mvgam_series)
 export(plot_mvgam_smooth)
 export(plot_mvgam_trend)
 export(plot_mvgam_uncertainty)
+export(plot_predictions.mvgam)
 export(ppc)
 export(roll_eval_mvgam)
 export(score)
@@ -131,6 +132,7 @@ importFrom(magrittr,"%>%")
 importFrom(marginaleffects,get_coef)
 importFrom(marginaleffects,get_predict)
 importFrom(marginaleffects,get_vcov)
+importFrom(marginaleffects,plot_predictions)
 importFrom(marginaleffects,set_coef)
 importFrom(methods,new)
 importFrom(mgcv,bam)
@@ -150,6 +152,7 @@ importFrom(posterior,as_draws_rvars)
 importFrom(posterior,rhat)
 importFrom(posterior,variables)
 importFrom(rlang,missing_arg)
+importFrom(rlang,parse_expr)
 importFrom(rstantools,posterior_epred)
 importFrom(rstantools,posterior_linpred)
 importFrom(rstantools,posterior_predict)
@@ -169,7 +172,9 @@ importFrom(stats,dlnorm)
 importFrom(stats,dnbinom)
 importFrom(stats,dnorm)
 importFrom(stats,dpois)
+importFrom(stats,drop.terms)
 importFrom(stats,ecdf)
+importFrom(stats,fitted)
 importFrom(stats,formula)
 importFrom(stats,frequency)
 importFrom(stats,gaussian)
diff --git a/R/compute_edf.R b/R/compute_edf.R
index 1bd7ae35..b802bc55 100644
--- a/R/compute_edf.R
+++ b/R/compute_edf.R
@@ -1,4 +1,5 @@
 #' Compute approximate EDFs of smooths
+#' @importFrom stats fitted
 #'@noRd
 compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names){
 
@@ -91,7 +92,11 @@ compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names){
         (mgcv_model$off[i] + ncol(mgcv_model$S[[i]]) - 1)
       XWXS[ind, ind] <- XWXS[ind, ind] + mgcv_model$S[[i]] * lambda[i]
     }
-    edf <- diag(solve(XWXS, XWX))
+    suppressWarnings(edf <- try(diag(solve(XWXS, XWX)), silent = TRUE))
+    if(inherits(edf, 'try-error')){
+      edf <- rep(1, length(coef(mgcv_model)))
+      names(edf) <- names(coef(mgcv_model))
+    }
     mgcv_model$edf <- edf
   }
 
diff --git a/R/dynamic.R b/R/dynamic.R
index cd681b75..e0196341 100644
--- a/R/dynamic.R
+++ b/R/dynamic.R
@@ -145,6 +145,24 @@ interpret_mvgam = function(formula, N){
 
   attr(newformula, '.Environment') <- attr(formula, '.Environment')
 
+  # Check if any terms use the gp wrapper, as mvgam cannot handle
+  # multivariate GPs yet
+  response <- terms.formula(newformula)[[2]]
+  tf <- terms.formula(newformula, specials = c("gp"))
+  which_gp <- attr(tf,"specials")$gp
+  if(length(which_gp) != 0L){
+    gp_details <- vector(length = length(which_gp),
+                          mode = 'list')
+      for(i in seq_along(which_gp)){
+        gp_details[[i]] <- eval(parse(text = rownames(attr(tf,
+                                                            "factors"))[which_gp[i]]))
+      }
+    if(any(unlist(lapply(purrr::map(gp_details, 'term'), length)) > 1)){
+      stop('mvgam cannot yet handle multidimensional gps',
+           call. = FALSE)
+    }
+  }
+
   # Check if any terms use the dynamic wrapper
   response <- terms.formula(newformula)[[2]]
   tf <- terms.formula(newformula, specials = c("dynamic"))
diff --git a/R/get_linear_predictors.R b/R/get_linear_predictors.R
index 248946bc..22d4bbf2 100644
--- a/R/get_linear_predictors.R
+++ b/R/get_linear_predictors.R
@@ -26,6 +26,47 @@ obs_Xp_matrix = function(newdata, mgcv_model){
                                     type = 'lpmatrix'))
   }
 
+  # Check for any gp() terms and update the design matrix
+  # accordingly
+  if(!is.null(attr(mgcv_model, 'gp_att_table'))){
+    # Compute the eigenfunctions from the supplied attribute table,
+    # and add them to the Xp matrix
+
+    # Extract GP attributes
+    gp_att_table <- attr(mgcv_model, 'gp_att_table')
+    gp_covariates <- unlist(purrr::map(gp_att_table, 'covariate'))
+    by <- unlist(purrr::map(gp_att_table, 'by'))
+    level <- unlist(purrr::map(gp_att_table, 'level'))
+    k <- unlist(purrr::map(gp_att_table, 'k'))
+    scale <- unlist(purrr::map(gp_att_table, 'scale'))
+    mean <- unlist(purrr::map(gp_att_table, 'mean'))
+    max_dist <- unlist(purrr::map(gp_att_table, 'max_dist'))
+    boundary <- unlist(purrr::map(gp_att_table, 'boundary'))
+    L <- unlist(purrr::map(gp_att_table, 'L'))
+
+    # Compute eigenfunctions
+    test_eigenfunctions <- lapply(seq_along(gp_covariates), function(x){
+      prep_eigenfunctions(data = newdata,
+                          covariate = gp_covariates[x],
+                          by = by[x],
+                          level = level[x],
+                          k = k[x],
+                          boundary = boundary[x],
+                          L = L[x],
+                          mean = mean[x],
+                          scale = scale[x],
+                          max_dist = max_dist[x])
+    })
+
+    # Find indices to replace in the design matrix and replace with
+    # the computed eigenfunctions
+    starts <- purrr::map(gp_att_table, 'first_coef')
+    ends <- purrr::map(gp_att_table, 'last_coef')
+    for(i in seq_along(starts)){
+      Xp[,c(starts[[i]]:ends[[i]])] <- test_eigenfunctions[[i]]
+    }
+  }
+
   return(Xp)
 }
 
@@ -46,31 +87,6 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all',
   trend_test$series <- trend_indicators
   trend_test$y <- NULL
 
-  # Because these are set up inherently as dynamic factor models,
-  # we ALWAYS need to forecast the full set of trends, regardless of
-  # which series (or set of series) is being forecast
-  # data.frame(series = trend_test$series,
-  #            time = trend_test$time,
-  #            row_num = 1:length(trend_test$time)) %>%
-  #   dplyr::group_by(series, time) %>%
-  #   dplyr::slice_head(n = 1) %>%
-  #   dplyr::ungroup() %>%
-  #   dplyr::arrange(time, series) %>%
-  #   dplyr::pull(row_num) -> inds_keep
-  #
-  # if(inherits(newdata, 'list')){
-  #   trend_test <- lapply(trend_test, function(x){
-  #     if(is.matrix(x)){
-  #       matrix(x[inds_keep,], ncol = NCOL(x))
-  #     } else {
-  #       x[inds_keep]
-  #     }
-  #
-  #   })
-  # } else {
-  #   trend_test <- trend_test[inds_keep, ]
-  # }
-
   suppressWarnings(Xp_trend  <- try(predict(mgcv_model,
                                             newdata = trend_test,
                                             type = 'lpmatrix'),
@@ -95,6 +111,48 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all',
                                           newdata = testdat,
                                           type = 'lpmatrix'))
   }
+
+  # Check for any gp() terms and update the design matrix
+  # accordingly
+  if(!is.null(attr(mgcv_model, 'gp_att_table'))){
+    # Compute the eigenfunctions from the supplied attribute table,
+    # and add them to the Xp matrix
+
+    # Extract GP attributes
+    gp_att_table <- attr(mgcv_model, 'gp_att_table')
+    gp_covariates <- unlist(purrr::map(gp_att_table, 'covariate'))
+    by <- unlist(purrr::map(gp_att_table, 'by'))
+    level <- unlist(purrr::map(gp_att_table, 'level'))
+    k <- unlist(purrr::map(gp_att_table, 'k'))
+    scale <- unlist(purrr::map(gp_att_table, 'scale'))
+    mean <- unlist(purrr::map(gp_att_table, 'mean'))
+    max_dist <- unlist(purrr::map(gp_att_table, 'max_dist'))
+    boundary <- unlist(purrr::map(gp_att_table, 'boundary'))
+    L <- unlist(purrr::map(gp_att_table, 'L'))
+
+    # Compute eigenfunctions
+    test_eigenfunctions <- lapply(seq_along(gp_covariates), function(x){
+      prep_eigenfunctions(data = newdata,
+                          covariate = gp_covariates[x],
+                          by = by[x],
+                          level = level[x],
+                          k = k[x],
+                          boundary = boundary[x],
+                          L = L[x],
+                          mean = mean[x],
+                          scale = scale[x],
+                          max_dist = max_dist[x])
+    })
+
+    # Find indices to replace in the design matrix and replace with
+    # the computed eigenfunctions
+    starts <- purrr::map(gp_att_table, 'first_coef')
+    ends <- purrr::map(gp_att_table, 'last_coef')
+    for(i in seq_along(starts)){
+      Xp_trend[,c(starts[[i]]:ends[[i]])] <- test_eigenfunctions[[i]]
+    }
+  }
+
   return(Xp_trend)
 }
 
diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R
index d1b06ea1..3118dc1e 100644
--- a/R/get_mvgam_priors.R
+++ b/R/get_mvgam_priors.R
@@ -223,6 +223,31 @@ get_mvgam_priors = function(formula,
   # Ensure series and time variables are present
   data_train <- validate_series_time(data_train, name = 'data')
 
+  # Check for gp and mo terms in the formula
+  orig_formula <- gp_terms <- mo_terms <- NULL
+  if(any(grepl('gp(', attr(terms(formula), 'term.labels'), fixed = TRUE))){
+
+    # Check that there are no multidimensional gp terms
+    formula <- interpret_mvgam(formula, N = max(data_train$time))
+    orig_formula <- formula
+
+    # Keep intercept?
+    keep_intercept <- attr(terms(formula), 'intercept') == 1
+
+    # Indices of gp() terms in formula
+    gp_terms <- which_are_gp(formula)
+
+    # Extract GP attributes
+    gp_details <- get_gp_attributes(formula)
+
+    # Replace with s() terms so the correct terms are included
+    # in the model.frame
+    formula <- gp_to_s(formula)
+    if(!keep_intercept){
+      formula <- update(formula, trend_y  ~ . -1)
+    }
+  }
+
   # Check for missing rhs in formula
   drop_obs_intercept <- FALSE
   if(length(attr(terms(formula), 'term.labels')) == 0 &
@@ -240,7 +265,9 @@ get_mvgam_priors = function(formula,
            call. = FALSE)
     }
   }
-  orig_formula <- formula
+  if(is.null(orig_formula)){
+    orig_formula <- formula
+  }
 
   # Validate observation formula
   formula <- interpret_mvgam(formula, N = max(data_train$time))
diff --git a/R/globals.R b/R/globals.R
index bba0c982..7960d0f1 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -8,4 +8,5 @@ utils::globalVariables(c("y", "year", "smooth_vals", "smooth_num",
                          "mod_call", "particles", "obs", "mgcv_model",
                          "param_name", "outcome", "mgcv_plottable",
                          "term", "data_test", "object", "row_num", "trends_test",
-                         "trend", "trend_series", "trend_y", ".", "gam"))
+                         "trend", "trend_series", "trend_y", ".", "gam",
+                         "group", "mod", "row_id"))
diff --git a/R/gp.R b/R/gp.R
new file mode 100644
index 00000000..edff9fc9
--- /dev/null
+++ b/R/gp.R
@@ -0,0 +1,683 @@
+#' Make gp() attributes table and necessary stan lines
+#' @noRd
+make_gp_additions = function(gp_details, data,
+                             newdata,
+                             model_data, mgcv_model,
+                             gp_terms){
+  # Need to expand combination of GPs if any of the by variables
+  # is a factor; mvgam will drop unused levels automatically
+  by <- gp_details$by
+  gp_details$row_id <- 1:NROW(gp_details)
+  gp_covariates <- gp_details$gp_covariates
+  gp_details_orig <- gp_details
+  if(any(!is.na(by))){
+    for(i in 1:length(by)){
+      if(!is.na(by[i])){
+        if(is.factor(data[[by[i]]])){
+          nlevels <- length(levels(droplevels(data[[by[i]]])))
+          new_details <- do.call(rbind, lapply(1:nlevels, function(x){
+            gp_details_orig[i,]
+          }))
+          new_details$level <- levels(droplevels(data[[by[i]]]))
+
+          rows_drop <- which(gp_details$gp_covariates == gp_covariates[i] &
+                               gp_details$by == by[i])
+
+          gp_details <- gp_details[-rows_drop,]
+          gp_details <- rbind(gp_details, new_details)
+        }
+      }
+    }
+  }
+
+  # Preserve ordering of terms
+  gp_details %>%
+    dplyr::arrange(row_id, level) %>%
+    dplyr::select(-row_id) -> gp_details
+
+  # Prepare the GP objects and Stan data lines;
+  # Data updates need to happen BEFORE calling vectorise
+  # model_file modifications can happen AFTER drop_obs_intercept
+  # to ensure everything behaves predictably
+  gp_covariates <- gp_details$gp_covariates
+  scale <- gp_details$scale
+  boundary <- gp_details$boundary
+  by <- gp_details$by
+  level <- gp_details$level
+
+  # Prep the covariate for GP modelling
+  gp_data <- lapply(seq_along(gp_covariates), function(x){
+
+    # Find the correct k to ensure that the total number of coefficients
+    # when using gp(k = k) is the same as when using s(k = k + 1)
+    smooth_terms <- unlist(purrr::map(mgcv_model$smooth, 'term'))
+    smooth_bys <- unlist(purrr::map(mgcv_model$smooth, 'by'))
+    if(any(smooth_bys == 'NA')){
+      smooth_bys[smooth_bys == 'NA'] <- NA
+    }
+    term_k <- mgcv_model$smooth[[min(which(smooth_bys %in% by[x] &
+                               smooth_terms == gp_covariates[x]))]]$df
+
+    prep_gp_covariate(data = data,
+                      covariate = gp_covariates[x],
+                      by = by[x],
+                      level = level[x],
+                      scale = scale[x],
+                      k = term_k,
+                      boundary = boundary[x])
+  })
+
+  # Consolidate Stan data objects and add to model_data
+  gp_stan_data <- do.call(c, purrr::map(gp_data, 'data_append'))
+  model_data <- append(model_data, gp_stan_data)
+
+  # Consolidate attribute tables
+  gp_att_table <- purrr::map(gp_data, 'att_table')
+
+  # Create updated design matrix by replacing the s() basis functions with
+  # the gp() eigenfunctions
+  coefs_replace <- list()
+    for(x in gp_terms){
+      label <- attr(terms(formula(mgcv_model)), 'term.labels')[x]
+      s_attributes <- eval(rlang::parse_expr(label))
+      if(s_attributes$by != 'NA'){
+        coef_name <- paste0('s(', s_attributes$term, '):', s_attributes$by)
+      } else {
+        coef_name <- paste0('s(', s_attributes$term, ')')
+      }
+      which_replace <- grep(coef_name, names(coef(mgcv_model)), fixed = TRUE)
+      names(mgcv_model$coefficients)[which_replace] <-
+        gsub('s(', 'gp(', names(mgcv_model$coefficients)[which_replace],
+             fixed = TRUE)
+      coefs_replace[[x]] <- which_replace
+    }
+
+  newX <- model_data$X
+
+  # Training data eigenfunctions
+  eigenfuncs <- do.call(cbind, purrr::map(gp_data, 'eigenfunctions'))
+
+  # Testing data eigenfunctions
+  if(!is.null(newdata)){
+    gp_covariates <- unlist(purrr::map(gp_att_table, 'covariate'))
+    by <- unlist(purrr::map(gp_att_table, 'by'))
+    level <- unlist(purrr::map(gp_att_table, 'level'))
+    k <- unlist(purrr::map(gp_att_table, 'k'))
+    scale <- unlist(purrr::map(gp_att_table, 'scale'))
+    mean <- unlist(purrr::map(gp_att_table, 'mean'))
+    max_dist <- unlist(purrr::map(gp_att_table, 'max_dist'))
+    boundary <- unlist(purrr::map(gp_att_table, 'boundary'))
+    L <- unlist(purrr::map(gp_att_table, 'L'))
+    test_eigenfunctions <- lapply(seq_along(gp_covariates), function(x){
+      prep_eigenfunctions(data = newdata,
+                          covariate = gp_covariates[x],
+                          by = by[x],
+                          level = level[x],
+                          k = k[x],
+                          boundary = boundary[x],
+                          L = L[x],
+                          mean = mean[x],
+                          scale = scale[x],
+                          max_dist = max_dist[x])
+    })
+
+    eigenfuncs <- rbind(eigenfuncs,
+                        do.call(cbind, test_eigenfunctions))
+  }
+
+  newX[unlist(coefs_replace), ] <- t(eigenfuncs)
+  model_data$X <- newX
+
+  # Consolidate Stan data lines
+  gp_stan_lines <- paste0(purrr::map(gp_data, 'data_lines'), collapse = '')
+
+  # Add coefficient indices to attribute table and to Stan data
+  for(covariate in seq_along(gp_att_table)){
+    coef_indices <- grep(gp_att_table[[covariate]]$name,
+                         names(coef(mgcv_model)), fixed = TRUE)
+    gp_att_table[[covariate]]$first_coef <- min(coef_indices)
+    gp_att_table[[covariate]]$last_coef <- max(coef_indices)
+
+    gp_names <- gp_att_table[[covariate]]$name
+    gp_names <- gsub('(', '_', gp_names, fixed = TRUE)
+    gp_names <- gsub(')', '_', gp_names, fixed = TRUE)
+    gp_names <- gsub(':', 'by', gp_names, fixed = TRUE)
+
+    gp_stan_lines <- paste0(gp_stan_lines,
+                            paste0('array[',  gp_att_table[[covariate]]$k,
+                                   '] int b_idx_',
+                                   gp_names,
+                                   '; // gp basis coefficient indices\n'))
+    gp_idx_data <- list(coef_indices)
+    names(gp_idx_data) <- paste0('b_idx_',
+                                 gp_names)
+    model_data <- append(model_data, gp_idx_data)
+  }
+
+  # Add the GP attribute table to the mgcv_model
+  attr(mgcv_model, 'gp_att_table') <- gp_att_table
+
+  return(list(model_data = model_data,
+              mgcv_model = mgcv_model,
+              gp_stan_lines = gp_stan_lines,
+              gp_att_table = gp_att_table))
+}
+
+#' Which terms are gp() terms?
+#' @noRd
+which_are_gp = function(formula){
+  tf <- terms.formula(formula, specials = c("gp"))
+  if(is.null(rlang::f_lhs(formula))){
+    out <- attr(tf,"specials")$gp
+  } else {
+    out <- attr(tf,"specials")$gp - 1
+  }
+  return(out)
+}
+
+#' Convert gp() terms to s() terms for initial model construction#'
+#' @importFrom stats drop.terms
+#' @noRd
+gp_to_s <- function(formula){
+  # Extract details of gp() terms
+  gp_details <- get_gp_attributes(formula)
+
+  # Drop these terms from the formula
+  which_gp <- which_are_gp(formula)
+  response <- rlang::f_lhs(formula)
+
+  suppressWarnings(tt <- try(drop.terms(terms(formula),
+                                        which_gp,
+                                        keep.response = TRUE),
+                             silent = TRUE))
+  if(inherits(tt, 'try-error')){
+    newformula <- as.formula(paste(response, '~ 1'))
+  } else {
+    tt <- drop.terms(terms(formula), which_gp, keep.response = TRUE)
+    newformula <- reformulate(attr(tt, "term.labels"), rlang::f_lhs(formula))
+  }
+
+  # Now replace the gp() terms with s() for constructing the initial model
+  s_terms <- vector()
+  for(i in 1:NROW(gp_details)){
+    if(!is.na(gp_details$by[i])){
+      s_terms[i] <- paste0('s(',
+                           gp_details$gp_covariates[i],
+                           ', by = ',
+                           gp_details$by[i],
+                           ', k = ',
+                           gp_details$k[i] + 1, ')')
+    } else {
+      s_terms[i] <- paste0('s(',
+                           gp_details$gp_covariates[i],
+                           ', k = ',
+                           gp_details$k[i] + 1, ')')
+    }
+  }
+
+  if(length(attr(terms(newformula), 'term.labels')) == 0){
+    rhs <- '1'
+  } else {
+    rhs <- attr(terms(newformula), 'term.labels')
+  }
+
+  newformula <- as.formula(paste(response, '~',
+                                 paste(paste(rhs,
+                                             collapse = '+'), '+',
+                                       paste(s_terms, collapse = '+'))))
+  attr(newformula, '.Environment') <- attr(formula, '.Environment')
+  return(newformula)
+}
+
+
+#' Store attributes of the gp terms
+#' @importFrom rlang parse_expr
+#' @noRd
+get_gp_attributes = function(formula){
+  gp_terms <- rownames(attr(terms(formula), 'factors'))[
+    grep('gp', rownames(attr(terms(formula), 'factors')))]
+  gp_attributes <- lapply(seq_along(gp_terms), function(x){
+    eval(rlang::parse_expr(gp_terms[x]))
+  })
+
+  # Extract information necessary to construct the GP terms
+  gp_covariates <- unlist(purrr::map(gp_attributes, 'term'))
+  k <- unlist(purrr::map(gp_attributes, 'k'))
+  if(any(is.na(k))){
+    k[is.na(k)] <- 10
+  }
+  scale <- unlist(purrr::map(gp_attributes, 'scale'))
+  boundary <- unlist(purrr::map(gp_attributes, 'c'))
+  if(any(is.na(boundary))){
+    boundary[is.na(boundary)] <- 5.0/4
+  }
+  by <- unlist(purrr::map(gp_attributes, 'by'))
+  if(any(by == 'NA')){
+    by[by == 'NA'] <- NA
+  }
+
+  # Return as a data.frame
+  return(data.frame(gp_covariates,
+                    k,
+                    scale,
+                    boundary,
+                    by,
+                    level = NA))
+}
+
+
+#' Propagate a Hilbert basis GP
+#' Evaluate Laplacian eigenfunction for a given GP basis function
+#' @noRd
+phi = function(boundary, m, centred_covariate) {
+  1 / sqrt(boundary) * sin((m * pi)/(2 * boundary) *
+                             (centred_covariate + boundary))
+}
+
+#' Evaluate eigenvalues for a given GP basis function
+#' @noRd
+lambda = function(boundary, m) {
+  ((m * pi)/(2 * boundary))^2
+}
+
+#' Spectral density squared exponential Gaussian Process kernel
+#' @noRd
+spd = function(alpha_gp, rho_gp, eigenvalues) {
+  (alpha_gp^2) * sqrt(2 * pi) * rho_gp *
+    exp(-0.5 * (rho_gp^2) * (eigenvalues^2))
+}
+
+#' @noRd
+sim_hilbert_gp = function(alpha_gp,
+                          rho_gp,
+                          b_gp,
+                          last_trends,
+                          fc_times,
+                          train_times,
+                          mean_train_times){
+
+  num_gp_basis <- length(b_gp)
+
+  # Get vector of eigenvalues of covariance matrix
+  eigenvalues <- vector()
+  for(m in 1:num_gp_basis){
+    eigenvalues[m] <- lambda(boundary = (5.0/4) *
+                               (max(train_times) - min(train_times)),
+                             m = m)
+  }
+
+  # Get vector of eigenfunctions
+  eigenfunctions <- matrix(NA, nrow = length(fc_times),
+                           ncol = num_gp_basis)
+  for(m in 1:num_gp_basis){
+    eigenfunctions[, m] <- phi(boundary = (5.0/4) *
+                                 (max(train_times) - min(train_times)),
+                               m = m,
+                               centred_covariate = fc_times - mean_train_times)
+  }
+
+  # Compute diagonal of covariance matrix
+  diag_SPD <- sqrt(spd(alpha_gp = alpha_gp,
+                       rho_gp = rho_gp,
+                       sqrt(eigenvalues)))
+
+  # Compute GP trend forecast
+  as.vector((diag_SPD * b_gp) %*% t(eigenfunctions))
+}
+
+#' Mean-center and scale the particular covariate of interest
+#' so that the maximum Euclidean distance between any two points is 1
+#' @noRd
+scale_cov <- function(data, covariate, by, level, scale = TRUE,
+                      mean, max_dist){
+  if(!is.na(by)){
+    if(!is.na(level)){
+      Xgp <- data[[covariate]][data[[by]] == level]
+    }
+  } else {
+    Xgp <- data[[covariate]]
+  }
+
+  if(is.na(mean)){
+    Xgp_mean <- mean(Xgp, na.rm = TRUE)
+  } else {
+    Xgp_mean <- mean
+  }
+
+  if(is.na(max_dist)){
+    Xgp_max_dist <- (abs(max(Xgp, na.rm = TRUE) -
+                           min(Xgp, na.rm = TRUE)))
+  } else {
+    Xgp_max_dist <- max_dist
+  }
+
+  if(scale){
+    # Mean center and divide by max euclidean distance
+    (Xgp - Xgp_mean) / Xgp_max_dist
+
+  } else {
+    # Just mean center
+    Xgp - Xgp_mean
+  }
+}
+
+#' prep GP eigenfunctions
+#' @noRd
+prep_eigenfunctions = function(data,
+                               covariate,
+                               by = NA,
+                               level = NA,
+                               k,
+                               boundary,
+                               mean = NA,
+                               max_dist = NA,
+                               scale = TRUE,
+                               L){
+
+  # Extract and scale covariate (scale set to FALSE if this is a prediction
+  # step so that we can scale by the original training covariate values supplied
+  # in mean and max_dist)
+  covariate_cent <- scale_cov(data = data,
+                              covariate = covariate,
+                              by = by,
+                              level = level,
+                              mean = mean,
+                              max_dist = max_dist,
+                              scale = scale)
+
+  # Construct matrix of eigenfunctions
+  eigenfunctions <- matrix(NA, nrow = length(covariate_cent),
+                           ncol = k)
+  if(missing(L)){
+    L <- brms:::choose_L(covariate_cent, boundary)
+  }
+
+  for(m in 1:k){
+    # eigenfunctions[, m] <- phi(boundary = boundary *
+    #                                      (max(covariate_cent) -
+    #                                         min(covariate_cent)),
+    #                                    m = m,
+    #                                    centred_covariate = covariate_cent)
+
+    eigenfunctions[, m] <- brms:::eigen_fun_cov_exp_quad(x = matrix(covariate_cent),
+                                                         m = m,
+                                                         L = L)
+  }
+
+  # Multiply eigenfunctions by the 'by' variable if one is supplied
+  if(!is.na(by)){
+    if(!is.na(level)){
+      # no multiplying needed as this is a factor by variable,
+      # but we need to pad the eigenfunctions with zeros
+      # for the observations where the by is a different level
+      full_eigens <- matrix(0, nrow = length(data[[by]]),
+                            ncol = NCOL(eigenfunctions))
+      full_eigens[(1:length(data[[by]]))[
+        data[[by]] == level],] <- eigenfunctions
+      eigenfunctions <- full_eigens
+    } else {
+      eigenfunctions <- eigenfunctions * data[[by]]
+    }
+  }
+  eigenfunctions
+}
+
+#' Prep Hilbert Basis GP covariates
+#' @noRd
+prep_gp_covariate = function(data,
+                             covariate,
+                             by = NA,
+                             level = NA,
+                             scale = TRUE,
+                             boundary = 5.0/4,
+                             k = 20){
+
+  covariate_cent <- scale_cov(data = data,
+                              covariate = covariate,
+                              scale = scale,
+                              by = by,
+                              mean = NA,
+                              max_dist = NA,
+                              level = level)
+
+  if(!is.na(by)){
+    if(!is.na(level)){
+      Xgp <- data[[covariate]][data[[by]] == level]
+    }
+  } else {
+    Xgp <- data[[covariate]]
+  }
+
+  covariate_mean <- mean(Xgp, na.rm = TRUE)
+  covariate_max_dist <- ifelse(scale,
+                               abs(max(Xgp,
+                                  na.rm = TRUE) -
+                                 min(Xgp,
+                                     na.rm = TRUE)),
+                               1)
+  # Check k
+  if(k > length(unique(covariate_cent))){
+    warning('argument "k" > number of unique covariate values;\ndropping to the maximum allowed "k"',
+            call. = FALSE)
+    k <- length(unique(covariate_cent))
+  }
+
+  # Construct vector of eigenvalues for GP covariance matrix; the
+  # same eigenvalues are always used in prediction, so we only need to
+  # create them when prepping the data. They will need to be included in
+  # the Stan data list
+  L <- brms:::choose_L(covariate_cent, boundary)
+  eigenvalues <- vector()
+  for(m in 1:k){
+    eigenvalues[m] <- sqrt(lambda(boundary = L,
+                                     m = m))
+  }
+
+  # Construct matrix of eigenfunctions; this will change depending on the values
+  # of the covariate, so it needs to be computed and included as data but also needs
+  # to be computed to make predictions
+  eigenfunctions <- prep_eigenfunctions(data = data,
+                                        covariate = covariate,
+                                        by = by,
+                                        level = level,
+                                        k = k,
+                                        boundary = boundary,
+                                        mean = NA,
+                                        max_dist = covariate_max_dist,
+                                        scale = scale)
+
+  # Make attributes table
+  byname <- ifelse(is.na(by), '', paste0(':', by))
+  covariate_name <- paste0('gp(', covariate, ')', byname)
+  if(!is.na(level)){
+    covariate_name <- paste0(covariate_name, level)
+  }
+  att_table <- list(effect = 'gp',
+                    name = covariate_name,
+                    covariate = covariate,
+                    by = by,
+                    level = level,
+                    k = k,
+                    boundary = boundary,
+                    L = L,
+                    scale = scale,
+                    mean = covariate_mean,
+                    max_dist = covariate_max_dist,
+                    eigenvalues = eigenvalues)
+
+  # Items to add to Stan data
+  # Number of basis functions
+  covariate_name <- gsub('(', '_', covariate_name, fixed = TRUE)
+  covariate_name <- gsub(')', '_', covariate_name, fixed = TRUE)
+  covariate_name <- gsub(':', 'by', covariate_name, fixed = TRUE)
+  data_lines <- paste0('int<lower=1> k_', covariate_name, '; // basis functions for approximate gp\n')
+  append_dat <- list(k = k)
+  names(append_dat) <- paste0('k_', covariate_name, '')
+
+  # Approximate GP eigenvalues
+  data_lines <- paste0(data_lines, paste0(
+    'vector[',
+    'k_', covariate_name,
+    '] l_', covariate_name, '; // approximate gp eigenvalues\n'),
+    collapse = '\n')
+  append_dat2 <- list(slambda = eigenvalues)
+  names(append_dat2) <- paste0('l_', covariate_name, '')
+  append_dat <- append(append_dat, append_dat2)
+
+  # Return necessary objects in a list
+  list(att_table = att_table,
+       data_lines = data_lines,
+       data_append = append_dat,
+       eigenfunctions = eigenfunctions)
+}
+
+add_gp_model_file = function(model_file, model_data, mgcv_model, gp_additions){
+
+  # Add data lines
+  model_file[grep('int<lower=0> ytimes[n, n_series];',
+                  model_file, fixed = TRUE)] <-
+    paste0(model_file[grep('int<lower=0> ytimes[n, n_series];',
+                           model_file, fixed = TRUE)],
+           '\n',
+           gp_additions$gp_stan_lines)
+  model_file <- readLines(textConnection(model_file), n = -1)
+
+  # Replace the multi_normal_prec lines with spd_cov_exp_quad
+  gp_names <- unlist(purrr::map(attr(mgcv_model, 'gp_att_table'), 'name'))
+  gp_names_clean <- gsub('(', '_', gp_names, fixed = TRUE)
+  gp_names_clean <- gsub(')', '_', gp_names_clean, fixed = TRUE)
+  gp_names_clean <- gsub(':', 'by', gp_names_clean, fixed = TRUE)
+  s_to_remove <- list()
+  for(i in seq_along(gp_names)){
+    s_name <- gsub('gp(', 's(', gp_names[i], fixed = TRUE)
+    to_replace <- grep(paste0('// prior for ', s_name, '...'),
+                       model_file, fixed = TRUE) + 1
+    pattern <- "S\\s*(.*?)\\s*\\["
+    result <- regmatches(model_file[to_replace],
+                         regexec(pattern, model_file[to_replace]))[[1]]
+    s_to_remove[[i]] <- unique(unlist(regmatches(result,
+                                                 gregexpr("[[:digit:]]+", result))))
+
+    model_file[grep(paste0('// prior for ', s_name, '...'),
+         model_file, fixed = TRUE)] <-
+      gsub('s(', 'gp(', model_file[grep(paste0('// prior for ', s_name, '...'),
+                                        model_file, fixed = TRUE)], fixed = TRUE)
+
+    model_file[to_replace] <-
+      paste0('z_',
+             gp_names_clean[i],
+             ' ~ std_normal();\n',
+             'alpha_',
+             gp_names_clean[i],
+             ' ~ normal(0, 0.5);\n',
+             'rho_',
+             gp_names_clean[i],
+             ' ~ inv_gamma(1.65, 5.97);\n',
+             'b_raw[b_idx_',
+             gp_names_clean[i],
+             '] ~ std_normal();\n')
+  }
+  b_line <- max(grep('b[', model_file, fixed = TRUE))
+  b_edits <- paste0('b[b_idx_',
+                    gp_names_clean,
+                    '] = sqrt(spd_cov_exp_quad(l_',
+                    gp_names_clean,
+                    ', alpha_',
+                    gp_names_clean,
+                    ', rho_',
+                    gp_names_clean,
+                    ')) .* z_',
+                    gp_names_clean,
+                    ';',
+                    collapse = '\n')
+  model_file[b_line] <- paste0(model_file[b_line],
+                               '\n',
+                               b_edits)
+  model_file <- readLines(textConnection(model_file), n = -1)
+
+  # Remove un-needed penalty matrices from the model file and the
+  # model data
+  for(i in seq_along(unique(unlist(s_to_remove)))){
+    model_data[[paste0('S',
+                       unique(unlist(s_to_remove))[i])]] <- NULL
+    model_file <- model_file[-grep(paste0('mgcv smooth penalty matrix S',
+                                          unique(unlist(s_to_remove))[i]),
+                                   model_file, fixed = TRUE)]
+  }
+
+  # Add alpha, rho and z lines in parameters and model blocks
+  alpha_names <- paste(paste0('real<lower=0> alpha_', gp_names_clean,
+                              ';'),
+                       collapse = '\n')
+  rho_names <- paste(paste0('real<lower=0> rho_', gp_names_clean,
+                            ';'),
+                     collapse = '\n')
+  z_names <- paste(paste0('vector[k_',
+                          gp_names_clean,
+                          '] z_',
+                          gp_names_clean,
+                          ';'),
+                   collapse = '\n')
+  model_file[grep("vector[num_basis] b_raw;", model_file, fixed = TRUE)] <-
+    paste0("vector[num_basis] b_raw;\n\n",
+           '// gp term sd parameters\n',
+           alpha_names,
+           '\n\n// gp term length scale parameters\n',
+           rho_names,
+           '\n\n// gp term latent variables\n',
+           z_names,
+           '\n')
+  model_file <- readLines(textConnection(model_file), n = -1)
+
+  # Add spd_cov_exp_quad function from brms code
+  if(any(grepl('functions {', model_file, fixed = TRUE))){
+    model_file[grep('functions {', model_file, fixed = TRUE)] <-
+      paste0('functions {\n',
+             '/* Spectral density function of a Gaussian process\n',
+             '* with squared exponential covariance kernel\n',
+             '* Args:\n',
+             '*   l_gp: numeric eigenvalues of an SPD GP\n',
+             '*   alpha_gp: marginal SD parameter\n',
+             '*   rho_gp: length-scale parameter\n',
+             '* Returns:\n',
+             '*   numeric values of the GP function evaluated at l_gp\n',
+             '*/\n',
+             'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n',
+             'int NB = size(l_gp);\n',
+             'vector[NB] out;\n',
+             'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n',
+             'real neg_half_lscale2 = -0.5 * square(rho_gp);\n',
+             'for (m in 1:NB) {\n',
+             'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n',
+             '}\n',
+             'return out;\n',
+             '}\n')
+  } else {
+    model_file[grep('Stan model code', model_file)] <-
+      paste0('// Stan model code generated by package mvgam\n',
+             'functions {\n',
+             '/* Spectral density function of a Gaussian process\n',
+             '* with squared exponential covariance kernel\n',
+             '* Args:\n',
+             '*   l_gp: numeric eigenvalues of an SPD GP\n',
+             '*   alpha_gp: marginal SD parameter\n',
+             '*   rho_gp: length-scale parameter\n',
+             '* Returns:\n',
+             '*   numeric values of the GP function evaluated at l_gp\n',
+             '*/\n',
+             'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n',
+             'int NB = size(l_gp);\n',
+             'vector[NB] out;\n',
+             'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n',
+             'real neg_half_lscale2 = -0.5 * square(rho_gp);\n',
+             'for (m in 1:NB) {\n',
+             'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n',
+             '}\n',
+             'return out;\n',
+             '}\n}\n')
+  }
+  model_file <- readLines(textConnection(model_file), n = -1)
+
+  return(list(model_file = model_file,
+              model_data = model_data))
+}
diff --git a/R/marginaleffects.mvgam.R b/R/marginaleffects.mvgam.R
index 061f79e7..2fb734da 100644
--- a/R/marginaleffects.mvgam.R
+++ b/R/marginaleffects.mvgam.R
@@ -17,6 +17,66 @@
 #' @author Nicholas J Clark
 NULL
 
+#' Effect plot as implemented in \pkg{marginaleffects}
+#'
+#' Convenient way to call marginal or conditional effect plotting functions
+#' implemented in the \pkg{marginaleffects} package
+#' @importFrom marginaleffects plot_predictions
+#' @importFrom bayesplot color_scheme_set color_scheme_get
+#' @inheritParams marginaleffects::plot_predictions
+#' @return A \code{\link[ggplot2:ggplot]{ggplot}} object
+#' that can be further customized using the \pkg{ggplot2} package,
+#' or a `data.frame` (if `draw=FALSE`)
+#'
+#' @export
+plot_predictions.mvgam = function(model,
+                                  condition = NULL,
+                                  by = NULL,
+                                  newdata = NULL,
+                                  type = NULL,
+                                  vcov = NULL,
+                                  conf_level = 0.95,
+                                  wts = NULL,
+                                  transform = NULL,
+                                  points = 0,
+                                  rug = FALSE,
+                                  gray = FALSE,
+                                  draw = TRUE,
+                                  ...){
+  # Set red colour scheme
+  orig_col <- .Options$ggplot2.discrete.colour
+  orig_fill <- .Options$ggplot2.discrete.fill
+  options(ggplot2.discrete.colour = c("#B97C7C",
+                                      "#A25050",
+                                      "#8F2727",
+                                      "darkred",
+                                      "#630000",
+                                      "#300000",
+                                      "#170000"),
+          ggplot2.discrete.fill = c("#B97C7C",
+                                    "#A25050",
+                                    "#8F2727",
+                                    "darkred",
+                                    "#630000",
+                                    "#300000",
+                                    "#170000"))
+  plot_predictions(model = model,
+                   condition = condition,
+                   by = by,
+                   newdata = newdata,
+                   type = type,
+                   vcov = vcov,
+                   conf_level = conf_level,
+                   wts = wts,
+                   transform = transform,
+                   points = points,
+                   rug = rug,
+                   gray = gray,
+                   draw = draw,
+                   ...)
+
+}
+
 #' Functions needed for working with marginaleffects
 #' @rdname mvgam_marginaleffects
 #' @export
diff --git a/R/mvgam.R b/R/mvgam.R
index 7e898317..3f5669b3 100644
--- a/R/mvgam.R
+++ b/R/mvgam.R
@@ -616,6 +616,31 @@ mvgam = function(formula,
   validate_pos_integer(samples)
   validate_pos_integer(thin)
 
+  # Check for gp and mo terms in the formula
+  orig_formula <- gp_terms <- mo_terms <- NULL
+  if(any(grepl('gp(', attr(terms(formula), 'term.labels'), fixed = TRUE))){
+
+    # Check that there are no multidimensional gp terms
+    formula <- interpret_mvgam(formula, N = max(data_train$time))
+    orig_formula <- formula
+
+    # Keep intercept?
+    keep_intercept <- attr(terms(formula), 'intercept') == 1
+
+    # Indices of gp() terms in formula
+    gp_terms <- which_are_gp(formula)
+
+    # Extract GP attributes
+    gp_details <- get_gp_attributes(formula)
+
+    # Replace with s() terms so the correct terms are included
+    # in the model.frame
+    formula <- gp_to_s(formula)
+    if(!keep_intercept){
+      formula <- update(formula, trend_y  ~ . -1)
+    }
+  }
+
   # Check for missing rhs in formula
   drop_obs_intercept <- FALSE
   if(length(attr(terms(formula), 'term.labels')) == 0 &
@@ -633,7 +658,9 @@ mvgam = function(formula,
            call. = FALSE)
     }
   }
-  orig_formula <- formula
+  if(is.null(orig_formula)){
+    orig_formula <- formula
+  }
 
   # Check for brmspriors
   if(!missing("data")){
@@ -1475,6 +1502,27 @@ mvgam = function(formula,
     inits <- family_inits(family = family_char, trend_model,
                           smooths_included, model_data)
 
+    # Include any GP term updates
+    if(!is.null(gp_terms)){
+      gp_additions <- make_gp_additions(gp_details = gp_details,
+                                        data = data_train,
+                                        newdata = data_test,
+                                        model_data = stan_objects$model_data,
+                                        mgcv_model = ss_gam,
+                                        gp_terms = gp_terms)
+      stan_objects$model_data <- gp_additions$model_data
+      ss_gam <- gp_additions$mgcv_model
+
+      gp_names <- unlist(purrr::map(gp_additions$gp_att_table, 'name'))
+      gp_names <- gsub('gp(', 's(', gp_names, fixed = TRUE)
+      rhos_change <- list()
+      for(i in seq_along(gp_names)){
+        rhos_change[[i]] <- grep(gp_names[i], rho_names, fixed = TRUE)
+      }
+      rho_names[c(unique(unlist(rhos_change)))] <- gsub('s(', 'gp(',
+                                                      rho_names[c(unique(unlist(rhos_change)))],
+                                                      fixed = TRUE)
+    }
 
     # Vectorise likelihoods
     vectorised <- vectorise_stan_lik(model_file = stan_objects$stan_file,
@@ -1585,6 +1633,33 @@ mvgam = function(formula,
       vectorised$model_file <- readLines(textConnection(vectorised$model_file), n = -1)
     }
 
+    # Remaining model file updates for any GP terms
+    if(!is.null(gp_terms)){
+     final_gp_updates <- add_gp_model_file(model_file = vectorised$model_file,
+                                           model_data = vectorised$model_data,
+                                           mgcv_model = ss_gam,
+                                           gp_additions = gp_additions)
+     vectorised$model_file <- final_gp_updates$model_file
+     vectorised$model_data <- final_gp_updates$model_data
+    }
+
+    # Update monitor pars for any GP terms
+    if(any(grepl('real<lower=0> alpha_gp',
+                vectorised$model_file, fixed = TRUE)) &
+       !lfo){
+      alpha_params <- trimws(gsub(';', '', gsub('real<lower=0> ',
+                                                '',
+                                                grep('real<lower=0> alpha_gp',
+                                                     vectorised$model_file, fixed = TRUE,
+                                                     value = TRUE), fixed = TRUE)))
+      rho_params <- trimws(gsub(';', '', gsub('real<lower=0> ',
+                                              '',
+                                              grep('real<lower=0> rho_gp',
+                                                   vectorised$model_file, fixed = TRUE,
+                                                   value = TRUE), fixed = TRUE)))
+      param <- c(param, alpha_params, rho_params)
+    }
+
     # Tidy the representation
     clean_up <- vector()
     for(x in 1:length(vectorised$model_file)){
@@ -1674,7 +1749,7 @@ mvgam = function(formula,
   # outside of mvgam ####
   if(!run_model){
     unlink('base_gam.txt')
-      output <- structure(list(call = formula,
+      output <- structure(list(call = orig_formula,
                                trend_call = if(!missing(trend_formula)){
                                  trend_formula
                                } else {
@@ -2133,7 +2208,7 @@ mvgam = function(formula,
   }
 
   #### Return the output as class mvgam ####
-  output <- structure(list(call = formula,
+  output <- structure(list(call = orig_formula,
                            trend_call = if(!missing(trend_formula)){
                              trend_formula
                            } else {
diff --git a/R/stan_utils.R b/R/stan_utils.R
index 052b9e58..8b78c303 100644
--- a/R/stan_utils.R
+++ b/R/stan_utils.R
@@ -2570,6 +2570,9 @@ add_trend_predictors = function(trend_formula,
   b_trend_lines <- gsub('\\bb\\b', 'b_trend', b_trend_lines)
   b_trend_lines <- gsub('raw', 'raw_trend', b_trend_lines)
   b_trend_lines <- gsub('num_basis', 'num_basis_trend', b_trend_lines)
+  b_trend_lines <- gsub('idx', 'trend_idx', b_trend_lines)
+  b_trend_lines <- gsub('l_gp', 'trend_l_gp', b_trend_lines)
+  b_trend_lines <- gsub('k_gp', 'trend_k_gp', b_trend_lines)
   model_file[grep("// derived latent states", model_file, fixed = TRUE)] <-
     paste0('// process model basis coefficients\n',
            paste(b_trend_lines, collapse = '\n'),
@@ -2589,10 +2592,11 @@ add_trend_predictors = function(trend_formula,
 
   # Add any multinormal smooth lines
   if(any(grepl('multi_normal_prec', trend_model_file)) |
-     any(grepl('// priors for smoothing parameters', trend_model_file))){
+     any(grepl('// priors for smoothing parameters', trend_model_file)) |
+     any(grepl('// prior for gp', trend_model_file))){
     trend_smooths_included <- TRUE
 
-    # Replace any noncontiguous indices from trend model so names aren't
+    # Replace any indices from trend model so names aren't
     # conflicting with any possible indices in the observation model
     if(any(grepl('idx', trend_model_file))){
       trend_model_file <- gsub('idx', 'trend_idx', trend_model_file)
@@ -2600,7 +2604,8 @@ add_trend_predictors = function(trend_formula,
       names(idx_data) <- gsub('idx', 'trend_idx', names(idx_data))
       model_data <- append(model_data, idx_data)
 
-      idx_lines <- grep('int trend_idx', trend_model_file)
+      idx_lines <- c(grep('int trend_idx', trend_model_file),
+                     grep('// gp basis coefficient indices', trend_model_file))
       model_file[min(grep('data {', model_file, fixed = TRUE))] <-
         paste0('data {\n',
                paste(trend_model_file[idx_lines],
@@ -2608,6 +2613,107 @@ add_trend_predictors = function(trend_formula,
       model_file <- readLines(textConnection(model_file), n = -1)
     }
 
+    # Check for gp() terms
+    if(any(grepl('l_gp', trend_model_file))){
+
+      # Add spd_cov_exp_quad function from brms code
+      if(any(grepl('functions {', model_file, fixed = TRUE))){
+        model_file[grep('functions {', model_file, fixed = TRUE)] <-
+          paste0('functions {\n',
+                 '/* Spectral density function of a Gaussian process\n',
+                 '* with squared exponential covariance kernel\n',
+                 '* Args:\n',
+                 '*   l_gp: numeric eigenvalues of an SPD GP\n',
+                 '*   alpha_gp: marginal SD parameter\n',
+                 '*   rho_gp: length-scale parameter\n',
+                 '* Returns:\n',
+                 '*   numeric values of the GP function evaluated at l_gp\n',
+                 '*/\n',
+                 'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n',
+                 'int NB = size(l_gp);\n',
+                 'vector[NB] out;\n',
+                 'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n',
+                 'real neg_half_lscale2 = -0.5 * square(rho_gp);\n',
+                 'for (m in 1:NB) {\n',
+                 'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n',
+                 '}\n',
+                 'return out;\n',
+                 '}\n')
+      } else {
+        model_file[grep('Stan model code', model_file)] <-
+          paste0('// Stan model code generated by package mvgam\n',
+                 'functions {\n',
+                 '/* Spectral density function of a Gaussian process\n',
+                 '* with squared exponential covariance kernel\n',
+                 '* Args:\n',
+                 '*   l_gp: numeric eigenvalues of an SPD GP\n',
+                 '*   alpha_gp: marginal SD parameter\n',
+                 '*   rho_gp: length-scale parameter\n',
+                 '* Returns:\n',
+                 '*   numeric values of the GP function evaluated at l_gp\n',
+                 '*/\n',
+                 'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n',
+                 'int NB = size(l_gp);\n',
+                 'vector[NB] out;\n',
+                 'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n',
+                 'real neg_half_lscale2 = -0.5 * square(rho_gp);\n',
+                 'for (m in 1:NB) {\n',
+                 'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n',
+                 '}\n',
+                 'return out;\n',
+                 '}\n}\n')
+      }
+      model_file <- readLines(textConnection(model_file), n = -1)
+
+      trend_model_file <- gsub('l_gp', 'trend_l_gp', trend_model_file)
+      trend_model_file <- gsub('k_gp', 'trend_k_gp', trend_model_file)
+      idx_data <- trend_mvgam$model_data[grep('l_gp', names(trend_mvgam$model_data))]
+      names(idx_data) <- gsub('l_gp', 'trend_l_gp', names(idx_data))
+      model_data <- append(model_data, idx_data)
+
+      l_lines <- grep('// approximate gp eigenvalues', trend_model_file, fixed = TRUE)
+      model_file[min(grep('data {', model_file, fixed = TRUE))] <-
+        paste0('data {\n',
+               paste(trend_model_file[l_lines],
+                     collapse = '\n'))
+      model_file <- readLines(textConnection(model_file), n = -1)
+    }
+
+    if(any(grepl('k_gp', trend_model_file))){
+      idx_data <- trend_mvgam$model_data[grep('k_gp', names(trend_mvgam$model_data))]
+      names(idx_data) <- gsub('k_gp', 'trend_k_gp', names(idx_data))
+      model_data <- append(model_data, idx_data)
+
+      k_lines <- grep('// basis functions for approximate gp', trend_model_file, fixed = TRUE)
+      model_file[min(grep('data {', model_file, fixed = TRUE))] <-
+        paste0('data {\n',
+               paste(trend_model_file[k_lines],
+                     collapse = '\n'))
+      model_file <- readLines(textConnection(model_file), n = -1)
+
+      # Update the parameters block with gp params
+      start <- grep("// gp term sd parameters", trend_model_file,
+                    fixed = TRUE)
+      end <- grep("// gp term latent variables", trend_model_file,
+                  fixed = TRUE) + 1
+      last <- end
+      for(i in end:(end+50)){
+        if(grepl('vector[trend_k_gp', trend_model_file[i],
+                 fixed = TRUE)){
+          last <- i
+        } else {
+          break
+        }
+      }
+      gp_params <- paste(trend_model_file[start:last],
+                         collapse = '\n')
+
+      model_file[min(grep('parameters {', model_file, fixed = TRUE))] <-
+        paste0('parameters {\n',
+               gp_params)
+      model_file <- readLines(textConnection(model_file), n = -1)
+    }
+
     if(any(grepl("int<lower=0> n_sp; // number of smoothing parameters",
                  model_file, fixed = TRUE))){
       model_file[grep("int<lower=0> n_sp; // number of smoothing parameters",
@@ -2630,6 +2736,12 @@ add_trend_predictors = function(trend_formula,
                                trend_model_file[grep('normal(0, lambda',
                                                      trend_model_file, fixed = TRUE)-1])
     }
+
+    if(any(grepl('// prior for gp', trend_model_file))){
+      spline_coef_headers <- c(spline_coef_headers,
+                               trend_model_file[grep('// prior for gp',
+                                                     trend_model_file, fixed = TRUE)])
+    }
     spline_coef_headers <- gsub('...', '_trend...', spline_coef_headers,
                                 fixed = TRUE)
 
@@ -2646,6 +2758,16 @@ add_trend_predictors = function(trend_formula,
       }
     }
 
+    if(any(grepl('// prior for gp', trend_model_file))){
+      starts <- grep('// prior for gp', trend_model_file, fixed = TRUE) + 1
+      ends <- grep('// prior for gp', trend_model_file, fixed = TRUE) + 4
+      for(i in seq_along(starts)){
+        spline_coef_lines <- c(spline_coef_lines,
+                               paste(trend_model_file[starts[i]:ends[i]],
+                                     collapse = '\n'))
+      }
+    }
+
     spline_coef_lines <- gsub('_raw', '_raw_trend', spline_coef_lines)
     spline_coef_lines <- gsub('lambda', 'lambda_trend', spline_coef_lines)
     spline_coef_lines <- gsub('zero', 'zero_trend', spline_coef_lines)
@@ -2757,7 +2879,6 @@ add_trend_predictors = function(trend_formula,
                "// posterior predictions")
     }
 
-
     model_file <- readLines(textConnection(model_file), n = -1)
 
   }
diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R
index 1070b735..4c843f01 100644
--- a/R/summary.mvgam.R
+++ b/R/summary.mvgam.R
@@ -180,84 +180,82 @@ if(object$family == 'Gamma'){
     }
   }
 
-smooth_labs <- do.call(rbind, lapply(seq_along(object$mgcv_model$smooth), function(x){
-  data.frame(label = object$mgcv_model$smooth[[x]]$label,
-             term = paste(object$mgcv_model$smooth[[x]]$term, collapse = ','),
-             class = class(object$mgcv_model$smooth[[x]])[1])
-}))
-
-if(any(!is.na(object$sp_names)) & !all(smooth_labs$class == 'random.effect')){
-  if(!is.null(object$trend_call)){
-    cat("\nGAM smoothing parameter (rho) estimates:\n")
-  } else {
-    cat("\nGAM observation smoothing parameter (rho) estimates:\n")
-  }
-  rho_coefs <- mcmc_summary(object$model_output, 'rho',
-                            digits = digits,
-                            variational = object$algorithm %in% c('fullrank',
-                                                                  'meanfield'))[,c(3:7)]
-  rownames(rho_coefs) <- paste0(object$sp_names, '_rho')
+if(all(is.na(object$sp_names))){
+
+} else {
+  if(any(unlist(purrr::map(object$mgcv_model$smooth, inherits, 'random.effect')))){
+    re_labs <- unlist(lapply(purrr::map(object$mgcv_model$smooth, 'term'),
+                             paste, collapse = ','))[
+                               unlist(purrr::map(object$mgcv_model$smooth, inherits, 'random.effect'))]
+    re_labs <- gsub('series', 'trend', re_labs)
+    re_sds <- mcmc_summary(object$model_output, 'sigma_raw',
+                           ISB = TRUE, digits = digits,
+                           variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
 
-  # Don't print random effect lambdas as they follow the prior distribution
-  if(any(smooth_labs$class == 'random.effect')){
-    re_smooths <- smooth_labs %>%
-      dplyr::filter(class == 'random.effect') %>%
-      dplyr::pull(label)
+    re_mus <- mcmc_summary(object$model_output, 'mu_raw',
+                           ISB = TRUE, digits = digits,
+                           variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
 
-  print(rho_coefs[!rownames(rho_coefs) %in% re_smooths,])
+    rownames(re_sds) <- paste0('sd(',re_labs,')')
+    rownames(re_mus) <- paste0('mean(',re_labs,')')
 
-  } else {
-    print(rho_coefs)
+    if(!is.null(object$trend_call)){
+      cat("\nGAM observation model group-level estimates:\n")
+    } else {
+      cat("\nGAM group-level estimates:\n")
+    }
+    print(rbind(re_mus, re_sds))
   }
 }
 
-if(any(!is.na(object$sp_names))){
-  cat("\nApproximate significance of GAM observation smooths:\n")
-  suppressWarnings(printCoefmat(summary(object$mgcv_model)$s.table[, c(1,3,4), drop = FALSE],
-                                digits = min(3, digits + 1),
-                                signif.stars = getOption("show.signif.stars"),
-                                has.Pvalue = TRUE, na.print = "NA",
-                                cs.ind = 1))
-}
-
-if(any(smooth_labs$class == 'random.effect')){
-  re_smooths <- smooth_labs %>%
-    dplyr::filter(class == 'random.effect') %>%
-    dplyr::pull(term)
-
-  if(object$fit_engine == 'jags'){
-    re_sds <- mcmc_summary(object$model_output,
-                           paste0('sigma_raw',
-                                  seq_along(re_smooths)),
-                           digits = digits,
-                           variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+if(!is.null(attr(object$mgcv_model, 'gp_att_table'))){
+  gp_names <- unlist(purrr::map(attr(object$mgcv_model, 'gp_att_table'), 'name'))
+  alpha_params <- gsub(':', 'by', gsub(')', '_',
+                       gsub('(', '_', paste0('alpha_', gp_names),
+                            fixed = TRUE), fixed = TRUE))
+  alpha_summary <- mcmc_summary(object$model_output, alpha_params,
+                         ISB = TRUE, digits = digits,
+                         variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+
+  rownames(alpha_summary) <- paste0('alpha_', gp_names)
+
+  rho_params <- gsub(':', 'by', gsub(')', '_',
+                                       gsub('(', '_', paste0('rho_', gp_names),
+                                            fixed = TRUE), fixed = TRUE))
+  rho_summary <- mcmc_summary(object$model_output, rho_params,
+                                ISB = TRUE, digits = digits,
+                                variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+  rownames(rho_summary) <- paste0('rho_', gp_names)
 
-    re_mus <- mcmc_summary(object$model_output,
-                           paste0('mu_raw',
-                                  seq_along(re_smooths)),
-                           digits = digits,
-                           variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+  if(!is.null(object$trend_call)){
+    cat("\nGAM observation model gp term marginal deviation (alpha) and length scale (rho) estimates:\n")
   } else {
-    re_sds <- mcmc_summary(object$model_output, 'sigma_raw',
-                           ISB = TRUE,
-                           digits = digits,
-                           variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
-
-    re_mus <- mcmc_summary(object$model_output, 'mu_raw',
-                           ISB = TRUE,
-                           digits = digits,
-                           variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+    cat("\nGAM gp term marginal deviation (alpha) and length scale (rho) estimates:\n")
   }
+  print(rbind(alpha_summary, rho_summary))
+}
 
-  rownames(re_sds) <- paste0('sd(',re_smooths,')')
-  rownames(re_mus) <- paste0('mean(',re_smooths,')')
+if(any(!is.na(object$sp_names))){
+  gam_sig_table <- summary(object$mgcv_model)$s.table[, c(1,3,4), drop = FALSE]
+  if(!is.null(attr(object$mgcv_model, 'gp_att_table'))){
+    gp_names <- unlist(purrr::map(attr(object$mgcv_model, 'gp_att_table'), 'name'))
+    if(all(rownames(gam_sig_table) %in% gsub('gp(', 's(', gp_names, fixed = TRUE))){
 
-  if(!is.null(object$trend_call)){
-    cat("\nGAM observation model group-level estimates:\n")
-    print(rbind(re_mus, re_sds))
-  } else {
-    cat("\nGAM group-level estimates:\n")
-    print(rbind(re_mus, re_sds))
+    } else {
+      gam_sig_table <- gam_sig_table[!rownames(gam_sig_table) %in%
+                                       gsub('gp(', 's(', gp_names, fixed = TRUE),,drop = FALSE]
+
+      if(!is.null(object$trend_call)){
+        cat("\nApproximate significance of GAM observation smooths:\n")
+      } else {
+        cat("\nApproximate significance of GAM observation smooths:\n")
+      }
+      suppressWarnings(printCoefmat(gam_sig_table,
+                                    digits = min(3, digits + 1),
+                                    signif.stars = getOption("show.signif.stars"),
+                                    has.Pvalue = TRUE, na.print = "NA",
+                                    cs.ind = 1))
+      }
   }
 }
 
@@ -576,54 +574,68 @@ if(!is.null(object$trend_call)){
   if(all(is.na(object$trend_sp_names))){
 
   } else {
-      to_print <- vector(length = length(object$trend_mgcv_model$smooth))
-      for(i in 1:length(object$trend_mgcv_model$smooth)){
-        if(inherits(object$trend_mgcv_model$smooth[[i]], 'random.effect')){
-          to_print[i] <- FALSE
-        } else {
-          to_print[i] <- TRUE
-        }
-      }
+    if(any(unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect')))){
+      re_labs <- unlist(lapply(purrr::map(object$trend_mgcv_model$smooth, 'term'),
+                               paste, collapse = ','))[
+                                 unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect'))]
+      re_labs <- gsub('series', 'trend', re_labs)
+      re_sds <- mcmc_summary(object$model_output, 'sigma_raw_trend',
+                             ISB = TRUE, digits = digits,
+                             variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+
+      re_mus <- mcmc_summary(object$model_output, 'mu_raw_trend',
+                             ISB = TRUE, digits = digits,
+                             variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+
+      rownames(re_sds) <- paste0('sd(',re_labs,')_trend')
+      rownames(re_mus) <- paste0('mean(',re_labs,')_trend')
+
+      cat("\nGAM process model group-level estimates:\n")
+      print(rbind(re_mus, re_sds))
+    }
+  }
+
+  if(!is.null(attr(object$trend_mgcv_model, 'gp_att_table'))){
+    gp_names <- unlist(purrr::map(attr(object$trend_mgcv_model, 'gp_att_table'), 'name'))
+    alpha_params <- gsub(':', 'by', gsub(')', '_',
+                                         gsub('(', '_', paste0('alpha_', gp_names),
+                                              fixed = TRUE), fixed = TRUE))
+    alpha_summary <- mcmc_summary(object$model_output, alpha_params,
+                                  ISB = TRUE, digits = digits,
+                                  variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
 
-      if(all(to_print == FALSE)){
+    rownames(alpha_summary) <- gsub('series', 'trend', paste0('alpha_', gp_names))
+
+    rho_params <- gsub(':', 'by', gsub(')', '_',
+                                       gsub('(', '_', paste0('rho_', gp_names),
+                                            fixed = TRUE), fixed = TRUE))
+    rho_summary <- mcmc_summary(object$model_output, rho_params,
+                                ISB = TRUE, digits = digits,
+                                variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
+    rownames(rho_summary) <- gsub('series', 'trend', paste0('rho_', gp_names))
+
+    cat("\nGAM process model gp term marginal deviation (alpha) and length scale (rho) estimates:\n")
+    print(rbind(alpha_summary, rho_summary))
+  }
+
+  if(any(!is.na(object$trend_sp_names))){
+    gam_sig_table <- summary(object$trend_mgcv_model)$s.table[, c(1,3,4), drop = FALSE]
+    if(!is.null(attr(object$trend_mgcv_model, 'gp_att_table'))){
+      gp_names <- unlist(purrr::map(attr(object$trend_mgcv_model, 'gp_att_table'), 'name'))
+      if(all(rownames(gam_sig_table) %in% gsub('gp(', 's(', gp_names, fixed = TRUE))){
 
       } else {
-        cat("\nGAM process smoothing parameter (rho) estimates:\n")
-        rho_coefs <- mcmc_summary(object$model_output, 'rho_trend',
-                                  digits = digits,
-                                  variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
-        rownames(rho_coefs) <- paste0(object$trend_sp_names, '_rho_trend')
-        print(rho_coefs[to_print,])
+        gam_sig_table <- gam_sig_table[!rownames(gam_sig_table) %in%
+                                         gsub('gp(', 's(', gp_names, fixed = TRUE),,drop = FALSE]
 
         cat("\nApproximate significance of GAM process smooths:\n")
-        suppressWarnings(printCoefmat(summary(object$trend_mgcv_model)$s.table[,c(1,3,4), drop = FALSE],
+        suppressWarnings(printCoefmat(gam_sig_table,
                                       digits = min(3, digits + 1),
                                       signif.stars = getOption("show.signif.stars"),
                                       has.Pvalue = TRUE, na.print = "NA",
                                       cs.ind = 1))
-
-      }
-
-      if(any(to_print == FALSE)){
-        re_labs <- unlist(lapply(purrr::map(object$trend_mgcv_model$smooth, 'term'),
-                                 paste, collapse = ','))[
-          unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect'))]
-        re_labs <- gsub('series', 'trend', re_labs)
-        re_sds <- mcmc_summary(object$model_output, 'sigma_raw_trend',
-                               ISB = TRUE, digits = digits,
-                               variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
-
-        re_mus <- mcmc_summary(object$model_output, 'mu_raw_trend',
-                               ISB = TRUE, digits = digits,
-                               variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)]
-
-        rownames(re_sds) <- paste0('sd(',re_labs,')_trend')
-        rownames(re_mus) <- paste0('mean(',re_labs,')_trend')
-
-        cat("\nGAM process model group-level estimates:\n")
-        print(rbind(re_mus, re_sds))
-
       }
+    }
   }
 }
 
@@ -649,8 +661,6 @@ if(object$fit_engine == 'jags'){
   } else {
     cat('\nRhat looks reasonable for all parameters\n')
   }
-
-
 }
 
 }
diff --git a/R/trends.R b/R/trends.R
index ba5e7aff..32377390 100644
--- a/R/trends.R
+++ b/R/trends.R
@@ -71,66 +71,6 @@ sim_gp = function(last_trends, h, rho_gp, alpha_gp){
                    h))
 }
 
-#' Propagate a Hilbert basis GP
-#' Evaluate Laplacian eigenfunction for a given GP basis function
-#' @noRd
-phi = function(boundary, m, centred_covariate) {
-  1 / sqrt(boundary) * sin((m * pi)/(2 * boundary) *
-                             (centred_covariate + boundary))
-}
-
-#' Evaluate eigenvalues for a given GP basis function
-#' @noRd
-lambda = function(boundary, m) {
-  ((m * pi)/(2 * boundary))^2
-}
-
-#' Spectral density squared exponential Gaussian Process kernel
-#' @noRd
-spd = function(alpha_gp, rho_gp, eigenvalues) {
-  (alpha_gp^2) * sqrt(2 * pi) * rho_gp *
-    exp(-0.5 * (rho_gp^2) * (eigenvalues^2))
-}
-
-#' @noRd
-sim_hilbert_gp = function(alpha_gp,
-                          rho_gp,
-                          b_gp,
-                          last_trends,
-                          fc_times,
-                          train_times,
-                          mean_train_times){
-
-  num_gp_basis <- length(b_gp)
-
-  # Get vector of eigenvalues of covariance matrix
-  eigenvalues <- vector()
-  for(m in 1:num_gp_basis){
-    eigenvalues[m] <- lambda(boundary = (5.0/4) *
-                               (max(train_times) - min(train_times)),
-                             m = m)
-  }
-
-  # Get vector of eigenfunctions
-  eigenfunctions <- matrix(NA, nrow = length(fc_times),
-                           ncol = num_gp_basis)
-  for(m in 1:num_gp_basis){
-    eigenfunctions[, m] <- phi(boundary = (5.0/4) *
-                                 (max(train_times) - min(train_times)),
-                               m = m,
-                               centred_covariate = fc_times - mean_train_times)
-  }
-
-  # Compute diagonal of covariance matrix
-  diag_SPD <- sqrt(spd(alpha_gp = alpha_gp,
-                       rho_gp = rho_gp,
-                       sqrt(eigenvalues)))
-
-  # Compute GP trend forecast
-  as.vector((diag_SPD * b_gp) %*% t(eigenfunctions))
-}
-
-
 #' AR3  simulation function
 #' @param last_trends Vector of trend estimates leading up to the current timepoint
 #' @param h \code{integer} specifying the forecast horizon
diff --git a/R/update_priors.R b/R/update_priors.R
index 5bb2745d..30845d86 100644
--- a/R/update_priors.R
+++ b/R/update_priors.R
@@ -39,8 +39,8 @@ update_priors = function(model_file,
                     model_file, fixed = TRUE))){
 
         # Updating parametric effects
-        if(any(grepl(paste0(priors$param_name[i], '...'), model_file, fixed = TRUE))){
-          header_line <- grep(paste0(priors$param_name[i], '...'), model_file, fixed = TRUE)
+        if(any(grepl(paste0('// prior for ', priors$param_name[i], '...'), model_file, fixed = TRUE))){
+          header_line <- grep(paste0('// prior for ', priors$param_name[i], '...'), model_file, fixed = TRUE)
           newprior <- paste(trimws(strsplit(priors$prior[i], "[~]")[[1]][2]))
           model_file[header_line + 1] <-
             paste(trimws(strsplit(model_file[header_line + 1], "[~]")[[1]][1]), '~',
diff --git a/man/plot_predictions.mvgam.Rd b/man/plot_predictions.mvgam.Rd
new file mode 100644
index 00000000..77f37263
--- /dev/null
+++ b/man/plot_predictions.mvgam.Rd
@@ -0,0 +1,107 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/marginaleffects.mvgam.R
+\name{plot_predictions.mvgam}
+\alias{plot_predictions.mvgam}
+\title{Effect plot as implemented in \pkg{marginaleffects}}
+\usage{
+plot_predictions.mvgam(
+  model,
+  condition = NULL,
+  by = NULL,
+  newdata = NULL,
+  type = NULL,
+  vcov = NULL,
+  conf_level = 0.95,
+  wts = NULL,
+  transform = NULL,
+  points = 0,
+  rug = FALSE,
+  gray = FALSE,
+  draw = TRUE,
+  ...
+)
+}
+\arguments{
+\item{model}{Model object}
+
+\item{condition}{Conditional predictions
+\itemize{
+\item Character vector (max length 3): Names of the predictors to display.
+\item Named list (max length 3): List names correspond to predictors. List elements can be:
+\itemize{
+\item Numeric vector
+\item Function which returns a numeric vector or a set of unique categorical values
+\item Shortcut strings for common reference values: "minmax", "quartile", "threenum"
+}
+\item 1: x-axis. 2: color/shape. 3: facets.
+\item Numeric variables in positions 2 and 3 are summarized by Tukey's five numbers \code{?stats::fivenum}
+}}
+
+\item{by}{Marginal predictions
+\itemize{
+\item Character vector (max length 3): Names of the categorical predictors to marginalize across.
+\item 1: x-axis. 2: color. 3: facets.
+}}
+
+\item{newdata}{When \code{newdata} is \code{NULL}, the grid is determined by the \code{condition} argument. When \code{newdata} is not \code{NULL}, the argument behaves in the same way as in the \code{predictions()} function.}
+
+\item{type}{string indicates the type (scale) of the predictions used to
+compute contrasts or slopes. This can differ based on the model
+type, but will typically be a string such as: "response", "link", "probs",
+or "zero". When an unsupported string is entered, the model-specific list of
+acceptable values is returned in an error message. When \code{type} is \code{NULL}, the
+default value is used. This default is the first model-related row in
+the \code{marginaleffects:::type_dictionary} dataframe.}
+
+\item{vcov}{Type of uncertainty estimates to report (e.g., for robust standard errors). Acceptable values:
+\itemize{
+\item FALSE: Do not compute standard errors. This can speed up computation considerably.
+\item TRUE: Unit-level standard errors using the default \code{vcov(model)} variance-covariance matrix.
+\item String which indicates the kind of uncertainty estimates to return.
+\itemize{
+\item Heteroskedasticity-consistent: \code{"HC"}, \code{"HC0"}, \code{"HC1"}, \code{"HC2"}, \code{"HC3"}, \code{"HC4"}, \code{"HC4m"}, \code{"HC5"}. See \code{?sandwich::vcovHC}
+\item Heteroskedasticity and autocorrelation consistent: \code{"HAC"}
+\item Mixed-Models degrees of freedom: "satterthwaite", "kenward-roger"
+\item Other: \code{"NeweyWest"}, \code{"KernHAC"}, \code{"OPG"}. See the \code{sandwich} package documentation.
+}
+\item One-sided formula which indicates the name of cluster variables (e.g., \code{~unit_id}). This formula is passed to the \code{cluster} argument of the \code{sandwich::vcovCL} function.
+\item Square covariance matrix
+\item Function which returns a covariance matrix (e.g., \code{stats::vcov(model)})
+}}
+
+\item{conf_level}{numeric value between 0 and 1. Confidence level to use to build a confidence interval.}
+
+\item{wts}{string or numeric: weights to use when computing average contrasts or slopes. These weights only affect the averaging in \verb{avg_*()} or with the \code{by} argument, and not the unit-level estimates themselves. Internally, estimates and weights are passed to the \code{weighted.mean()} function.
+\itemize{
+\item string: column name of the weights variable in \code{newdata}. When supplying a column name to \code{wts}, it is recommended to supply the original data (including the weights variable) explicitly to \code{newdata}.
+\item numeric: vector of length equal to the number of rows in the original data or in \code{newdata} (if supplied).
+}}
+
+\item{transform}{A function applied to unit-level adjusted predictions and confidence intervals just before the function returns results. For bayesian models, this function is applied to individual draws from the posterior distribution, before computing summaries.}
+
+\item{points}{Number between 0 and 1 which controls the transparency of raw data points. 0 (default) does not display any points.}
+
+\item{rug}{TRUE displays tick marks on the axes to mark the distribution of raw data.}
+
+\item{gray}{FALSE grayscale or color plot}
+
+\item{draw}{\code{TRUE} returns a \code{ggplot2} plot. \code{FALSE} returns a \code{data.frame} of the underlying data.}
+
+\item{...}{Additional arguments are passed to the \code{predict()} method
+supplied by the modeling package.These arguments are particularly useful
+for mixed-effects or bayesian models (see the online vignettes on the
+\code{marginaleffects} website). Available arguments can vary from model to
+model, depending on the range of supported arguments by each modeling
+package. See the "Model-Specific Arguments" section of the
+\code{?marginaleffects} documentation for a non-exhaustive list of available
+arguments.}
+}
+\value{
+A \code{\link[ggplot2:ggplot]{ggplot}} object
+that can be further customized using the \pkg{ggplot2} package,
+or a \code{data.frame} (if \code{draw=FALSE})
+}
+\description{
+Convenient way to call marginal or conditional effect plotting functions
+implemented in the \pkg{marginaleffects} package
+}
diff --git a/src/mvgam.dll b/src/mvgam.dll
index b52e44d4..764c3c83 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 fdc1786a..6eb07ea5 100644
Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ