diff --git a/R/globals.R b/R/globals.R
index e90a1c69..7c47f8b3 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -19,4 +19,5 @@ utils::globalVariables(c("y", "year", "smooth_vals", "smooth_num",
                          "time_lag", "dis_time", "maxt", "orig_rows",
                          "matches", "time.", "file_name", ".data",
                          "horizon", "target", "Series", "evd", "mean_evd",
-                         "total_evd"))
+                         "total_evd", "smooth_label", "by_variable",
+                         "gr", "tot_subgrs", "subgr"))
diff --git a/R/jsdgam.R b/R/jsdgam.R
index 10b0ea14..b3e8ee06 100644
--- a/R/jsdgam.R
+++ b/R/jsdgam.R
@@ -8,6 +8,7 @@
 #'\pkg{mvgam} can handle)
 #'
 #'@inheritParams mvgam
+#'@inheritParams ZMVN
 #'@param factor_formula A \code{character} string specifying the linear predictor
 #'effects for the latent factors. These are exactly like the formula
 #'for a GLM except that smooth terms, `s()`, `te()`, `ti()`, `t2()`, as well as time-varying
@@ -45,7 +46,7 @@
 #'pooled and any latent, unmodelled residual associations can also be learned. In \pkg{mvgam}, both of
 #'these effects can be modelled with the full power of latent factor Hierarchical GAMs, providing unmatched
 #'flexibility to model full communities of species
-#'@seealso \code{\link[mgcv]{mvgam}}
+#'@seealso [mvgam]
 #'@return A \code{list} object of classes \code{jsdgam} and \code{mvgam} containing model output,
 #'the text representation of the model file,
 #'the mgcv model output (for easily generating simulations at
diff --git a/man/jsdgam.Rd b/man/jsdgam.Rd
index 30787021..fe726e2e 100644
--- a/man/jsdgam.Rd
+++ b/man/jsdgam.Rd
@@ -86,6 +86,27 @@ of trials
 \item\code{beta_binomial()} as for \code{binomial()} but allows for overdispersion}
 Default is \code{poisson()}. See \code{\link{mvgam_families}} for more details}
 
+\item{unit}{The unquoted name of the variable that represents the unit of analysis in \code{data} over
+which latent residuals should be correlated. This variable should be either a
+\code{numeric} or \code{integer} variable in the supplied \code{data}.
+Defaults to \code{time} to be consistent with other functionalities
+in \pkg{mvgam}, though note that the data need not be time series in this case. See examples below
+for further details and explanations}
+
+\item{subgr}{A subgrouping \code{factor} variable specifying which element in \code{data} represents the
+different observational units. Defaults to \code{series} to be consistent with other functionalities
+in \pkg{mvgam}, though note that the data need not be time series in this case.
+But note that models that use the hierarchical correlations (by supplying a value for \code{gr})
+\emph{should not} include a \code{series} element in \code{data}. Rather, this element will be created internally based
+on the supplied variables for \code{gr} and \code{subgr}. For example, if you are modelling
+counts for a group of species (labelled as \code{species} in the data) across sampling sites
+(labelled as \code{site} in the data) in three
+different geographical regions (labelled as \code{region}), and you would like the residuals to be correlated
+within regions, then you should specify
+\code{unit = site},  \code{gr = region}, and \code{subgr = species}. Internally, \code{mvgam()} will appropriately order
+the data by \code{unit} (in this case, by \code{site}) and create
+the \code{series} element for the data using something like: \code{series = as.factor(paste0(group, '_', subgroup))}}
+
 \item{share_obs_params}{\code{logical}. If \code{TRUE} and the \code{family}
 has additional family-specific observation parameters (e.g. variance components in
 \code{student_t()} or \code{gaussian()}, or dispersion parameters in \code{nb()} or \code{betar()}),
@@ -352,7 +373,7 @@ loo(mod)
 }
 }
 \seealso{
-\code{\link[mgcv]{mvgam}}
+\link{mvgam}
 }
 \author{
 Nicholas J Clark
diff --git a/src/mvgam.dll b/src/mvgam.dll
index d78c66bd..c69c4397 100644
Binary files a/src/mvgam.dll and b/src/mvgam.dll differ