diff --git a/NAMESPACE b/NAMESPACE
index e51682f5..7a0f6529 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -158,5 +158,7 @@ importFrom(stats,ts)
 importFrom(stats,update)
 importFrom(stats,update.formula)
 importFrom(stats,var)
+importFrom(utils,getFromNamespace)
 importFrom(utils,head)
+importFrom(utils,lsf.str)
 importFrom(utils,tail)
diff --git a/R/dynamic.R b/R/dynamic.R
index d868dc2f..8ebdba14 100644
--- a/R/dynamic.R
+++ b/R/dynamic.R
@@ -38,7 +38,7 @@ dynamic = function(variable, rho = 5, stationary = TRUE){
 
   # Check rho
   if(rho <= 0){
-    stop('argument "rho" in dynamic() must be a positive value',
+    stop('Argument "rho" in dynamic() must be a positive value',
          call. = FALSE)
   }
 
@@ -112,7 +112,7 @@ interpret_mvgam = function(formula, N){
     dyn_to_gp = function(term, N){
 
       if(term$rho > N - 1){
-        stop('argument "rho" in dynamic() cannot be larger than (max(time) - 1)',
+        stop('Argument "rho" in dynamic() cannot be larger than (max(time) - 1)',
              call. = FALSE)
       }
 
diff --git a/R/formula.mvgam.R b/R/formula.mvgam.R
index 31032148..45ce5a71 100644
--- a/R/formula.mvgam.R
+++ b/R/formula.mvgam.R
@@ -1,7 +1,7 @@
 #'Extract model.frame from a fitted mvgam object
 #'
 #'
-#'@inheritParams stats::formula
+#'@param x `mvgam` object
 #'@param trend_effects \code{logical}, return the model.frame from the
 #'observation model (if \code{FALSE}) or from the underlying process
 #'model (if\code{TRUE})
diff --git a/R/globals.R b/R/globals.R
index c624ea58..38cbceb8 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -7,4 +7,5 @@ utils::globalVariables(c("y", "year", "smooth_vals", "smooth_num",
                          "assimilated", "eval_horizon", "label",
                          "mod_call", "particles", "obs", "mgcv_model",
                          "param_name", "outcome", "mgcv_plottable",
-                         "term", "data_test", "object", "row_num", "trends_test"))
+                         "term", "data_test", "object", "row_num", "trends_test",
+                         "trend"))
diff --git a/R/marginaleffects.mvgam.R b/R/marginaleffects.mvgam.R
index fdc4fb05..cedf8b37 100644
--- a/R/marginaleffects.mvgam.R
+++ b/R/marginaleffects.mvgam.R
@@ -2,12 +2,15 @@
 #' @importFrom stats coef model.frame
 #' @importFrom insight find_predictors get_data
 #' @importFrom marginaleffects get_coef set_coef get_vcov get_predict
+#' @importFrom utils getFromNamespace
 #' @inheritParams marginaleffects::get_coef
 #' @inheritParams marginaleffects::set_coef
 #' @inheritParams marginaleffects::get_vcov
 #' @inheritParams marginaleffects::get_predict
 #' @inheritParams insight::get_data
 #' @inheritParams insight::find_predictors
+#' @param trend_effects `logical`, extract from the process model component
+#' (only applicable if a `trend_formula` was specified in the model)
 #' @param process_error `logical`. If `TRUE`, uncertainty in the latent
 #' process (or trend) model is incorporated in predictions
 #' @param n_cores `Integer` specifying number of cores to use for
@@ -102,36 +105,55 @@ get_data.mvgam = function (x, source = "environment", verbose = TRUE, ...) {
     # there won't be an easy way to match them up (for example if multiple
     # series depend on a shared latent trend)
     if(!is.null(x$trend_call)){
+
+      # Original series, time and outcomes
+      orig_dat <- data.frame(series = x$obs_data$series,
+                             time = x$obs_data$time,
+                             y = x$obs_data$y)
+
       # Add indicators of trend names as factor levels using the trend_map
-      trend_indicators <- vector(length = length(x$obs_data$time))
-      for(i in 1:length(x$obs_data$time)){
+      trend_indicators <- vector(length = length(orig_dat$time))
+      for(i in 1:length(orig_dat$time)){
         trend_indicators[i] <- x$trend_map$trend[which(x$trend_map$series ==
-                                                         x$obs_data$series[i])]
+                                                         orig_dat$series[i])]
       }
       trend_indicators <- as.factor(paste0('trend', trend_indicators))
 
-      # Only keep one time observation per trend
-      data.frame(series = trend_indicators,
-                 time = x$obs_data$time,
-                 y = x$obs_data$y,
-                 row_num = 1:length(x$obs_data$time)) %>%
-        dplyr::group_by(series, time) %>%
+      # Trend-level data, before any slicing that took place
+      trend_level_data <- data.frame(trend_series = trend_indicators,
+                                     series = orig_dat$series,
+                                     time = orig_dat$time,
+                                     y = orig_dat$y,
+                                     row_num = 1:length(x$obs_data$time))
+
+      # # We only kept one time observation per trend
+      trend_level_data %>%
+        dplyr::group_by(trend_series, time) %>%
         dplyr::slice_head(n = 1) %>%
         dplyr::pull(row_num) -> idx
+
+      # Extract model.frame for trend_level effects and add the
+      # trend indicators
       mf_data <- model.frame(x, trend_effects = TRUE)
-      mf_obs <- model.frame(x, trend_effects = FALSE)[idx, , drop = FALSE]
-      mf_data <- cbind(mf_obs, mf_data)
+      mf_data$trend_series <- trend_level_data$trend_series[idx]
+      mf_data$time <- trend_level_data$time[idx]
+
+      # Now join with the original data so the original observations can
+      # be included
+      trend_level_data %>%
+        dplyr::left_join(mf_data, by = c('trend_series', 'time')) %>%
+        dplyr::select(-trend_series, -row_num, -trend_y) -> mf_data
+
+      # Extract any predictors from the observation level model and
+      # bind to the trend level model.frame
+      mf_obs <- model.frame(x, trend_effects = FALSE)
+      mf_data <- cbind(mf_obs, mf_data) %>%
+        subset(., select = which(!duplicated(names(.))))
 
       # Now get the observed response, in case there are any
       # NAs there that need to be updated
-      data.frame(series = trend_indicators,
-                 time = x$obs_data$time,
-                 y = x$obs_data$y,
-                 row_num = 1:length(x$obs_data$time)) %>%
-        dplyr::group_by(series, time) %>%
-        dplyr::slice_head(n = 1) %>%
-        dplyr::pull(y) -> obs_response
-      mf_data[,resp] <- obs_response
+      mf_data[,resp] <- x$obs_data$y
+
     } else {
       mf_data <- model.frame(x, trend_effects = FALSE)
       mf_data[,resp] <- x$obs_data[[resp]]
@@ -140,7 +162,9 @@ get_data.mvgam = function (x, source = "environment", verbose = TRUE, ...) {
   }, error = function(x) {
     NULL
   })
-  insight:::.prepare_get_data(x, mf, effects = "all", verbose = verbose)
+
+  prep_data <- utils::getFromNamespace(".prepare_get_data", "insight")
+  prep_data(x, mf, effects = "all", verbose = verbose)
 }
 
 #' @rdname mvgam_marginaleffects
diff --git a/R/predict.mvgam.R b/R/predict.mvgam.R
index 9110adfe..57780c76 100644
--- a/R/predict.mvgam.R
+++ b/R/predict.mvgam.R
@@ -1,5 +1,6 @@
 #'Predict from the GAM component of an mvgam model
 #'@importFrom parallel clusterExport stopCluster setDefaultCluster
+#'@importFrom utils lsf.str
 #'@importFrom stats predict
 #'@param object \code{list} object returned from \code{mvgam}
 #'@param newdata Optional \code{dataframe} or \code{list} of test data containing the
diff --git a/R/update.mvgam.R b/R/update.mvgam.R
index eeeecc4e..677d3c8c 100644
--- a/R/update.mvgam.R
+++ b/R/update.mvgam.R
@@ -111,7 +111,7 @@ update.mvgam = function(object, formula,
     trend_model <- object$trend_model
 
     if(trend_model == 'VAR1'){
-      if(any(is.na(mvgam:::mcmc_summary(object$model_output, 'Sigma')[,6]))){
+      if(any(is.na(mcmc_summary(object$model_output, 'Sigma')[,6]))){
         trend_model <- 'VAR1'
       } else {
         trend_model <- 'VAR1cor'
diff --git a/R/validations.R b/R/validations.R
index 5a6fd0fd..d6f0aac1 100644
--- a/R/validations.R
+++ b/R/validations.R
@@ -123,7 +123,7 @@ validate_proportional = function(x){
   s <- substitute(x)
   x <- base::suppressWarnings(as.numeric(x))
   if (length(x) != 1L || anyNA(x)) {
-    stop("argument '", s, "' must be a single numeric value",
+    stop("Argument '", s, "' must be a single numeric value",
          call. = FALSE)
   }
 
@@ -138,16 +138,16 @@ validate_pos_integer = function(x){
   s <- substitute(x)
   x <- base::suppressWarnings(as.numeric(x))
   if (length(x) != 1L || anyNA(x)) {
-    stop("argument '", s, "' must be a single numeric value",
+    stop("Argument '", s, "' must be a single numeric value",
          call. = FALSE)
   }
 
   if(sign(x) != 1){
-    stop("argument '", s, "' must be a positive integer",
+    stop("Argument '", s, "' must be a positive integer",
          call. = FALSE)
   } else {
     if(x%%1 != 0){
-      stop("argument '", s, "' must be a positive integer",
+      stop("Argument '", s, "' must be a positive integer",
            call. = FALSE)
     }
   }
@@ -173,13 +173,13 @@ validate_trendmap = function(trend_map,
 
   # trend_map must have an entry for each unique time series
   if(!all(sort(trend_map$series) == sort(unique(data_train$series)))){
-    stop('argument "trend_map" must have an entry for every unique time series in "data"',
+    stop('Argument "trend_map" must have an entry for every unique time series in "data"',
          call. = FALSE)
   }
 
   # trend_map must not specify a greater number of trends than there are series
   if(max(trend_map$trend) > length(unique(data_train$series))){
-    stop('argument "trend_map" specifies more latent trends than there are series in "data"',
+    stop('Argument "trend_map" specifies more latent trends than there are series in "data"',
          call. = FALSE)
   }
 
@@ -189,7 +189,7 @@ validate_trendmap = function(trend_map,
   }
 
   if(!all(drop_zero(sort(unique(trend_map$trend))) == seq(1:max(trend_map$trend)))){
-    stop('argument "trend_map" must link at least one series to each latent trend',
+    stop('Argument "trend_map" must link at least one series to each latent trend',
          call. = FALSE)
   }
 
diff --git a/man/formula.mvgam.Rd b/man/formula.mvgam.Rd
index 23951702..630b22a3 100644
--- a/man/formula.mvgam.Rd
+++ b/man/formula.mvgam.Rd
@@ -7,6 +7,8 @@
 \method{formula}{mvgam}(x, trend_effects = FALSE, ...)
 }
 \arguments{
+\item{x}{\code{mvgam} object}
+
 \item{trend_effects}{\code{logical}, return the model.frame from the
 observation model (if \code{FALSE}) or from the underlying process
 model (if\code{TRUE})}
diff --git a/man/mvgam_marginaleffects.Rd b/man/mvgam_marginaleffects.Rd
index b41fa5c3..2c7ba12e 100644
--- a/man/mvgam_marginaleffects.Rd
+++ b/man/mvgam_marginaleffects.Rd
@@ -40,6 +40,9 @@
 \arguments{
 \item{model}{Model object}
 
+\item{trend_effects}{\code{logical}, extract from the process model component
+(only applicable if a \code{trend_formula} was specified in the model)}
+
 \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
diff --git a/tests/testthat/test-dynamic.R b/tests/testthat/test-dynamic.R
index 0c9669e0..1a0df406 100644
--- a/tests/testthat/test-dynamic.R
+++ b/tests/testthat/test-dynamic.R
@@ -28,7 +28,7 @@ test_that("rho argument must be positive numeric", {
                             data = data,
                             family = gaussian(),
                             run_model = FALSE),
-               'argument "rho" in dynamic() must be a positive value',
+               'Argument "rho" in dynamic() must be a positive value',
                fixed = TRUE)
 })
 
@@ -40,11 +40,11 @@ test_that("rho argument cannot be larger than N - 1", {
                             data = data,
                             family = gaussian(),
                             run_model = FALSE),
-               'argument "rho" in dynamic() cannot be larger than (max(time) - 1)',
+               'Argument "rho" in dynamic() cannot be larger than (max(time) - 1)',
                fixed = TRUE)
 
   expect_error(mvgam:::interpret_mvgam(formula = y ~ dynamic(covariate, rho = 120),
                                        N = 100),
-               'argument "rho" in dynamic() cannot be larger than (max(time) - 1)',
+               'Argument "rho" in dynamic() cannot be larger than (max(time) - 1)',
                fixed = TRUE)
 })
diff --git a/tests/testthat/test-mvgam.R b/tests/testthat/test-mvgam.R
index 25de983e..d12e667d 100644
--- a/tests/testthat/test-mvgam.R
+++ b/tests/testthat/test-mvgam.R
@@ -45,7 +45,7 @@ test_that("rho argument must be positive numeric", {
                             data = data,
                             family = gaussian(),
                             run_model = FALSE),
-               'argument "rho" in dynamic() must be a positive value',
+               'Argument "rho" in dynamic() must be a positive value',
                fixed = TRUE)
 })
 
diff --git a/tests/testthat/test-sim_mvgam.R b/tests/testthat/test-sim_mvgam.R
index 162402db..8877c9e8 100644
--- a/tests/testthat/test-sim_mvgam.R
+++ b/tests/testthat/test-sim_mvgam.R
@@ -18,7 +18,7 @@ test_that("trend_rel must be a valid proportion", {
   expect_error(sim_mvgam(family = gaussian(),
                          trend_model = 'AR2',
                          trend_rel = -0.1),
-               'Argument "trend_rel" must be a proportion ranging from 0 to 1, inclusive')
+               "Argument 'trend_rel' must be a proportion ranging from 0 to 1, inclusive")
 })
 
 test_that("n_lv must be a positive integer", {
@@ -26,5 +26,6 @@ test_that("n_lv must be a positive integer", {
                          trend_model = 'AR2',
                          trend_rel = 0.4,
                          n_lv = 0.5),
-               'Argument "n_lv" must be a positive integer')
+               "Argument 'n_lv' must be a positive integer")
 })
+