Dynamic generalised additive models (DGAMs) for forecasting discrete ecological time series

Generalised additive models (GAMs) are increasingly popular tools for estimating smooth nonlinear relationships between predictors and response variables. GAMs are particularly relevant in ecology for representing hierarchical functions for discrete responses that encompass complex features including zero inflation, truncation and uneven sampling. However, GAMs are less useful for producing forecasts as their smooth functions provide unstable predictions outside the range of training data. We introduce dynamic generalised additive models (DGAMs), where the GAM linear predictor is jointly estimated with unobserved dynamic components to model time series that evolve as a function of nonlinear predictor associations and latent temporal processes. These models are especially useful for analysing multiple series, as they can estimate hierarchical smooth functions while learning complex temporal associations via dimension‐reduced latent factor processes. We implement our models in the mvgam R package, which estimates unobserved parameters for smoothing splines and latent temporal processes in a probabilistic framework. Using simulations, we illustrate how our models outperform competing formulations in realistic ecological forecasting tasks while identifying important smooth predictor functions. We use a real‐world case study to highlight some of mvgam's key features, which include functions for calculating correlations among series' latent trends, performing model selection using rolling window forecasts and posterior predictive checks, online data augmentation via a recursive particle filter and visualising probabilistic uncertainties for smooth functions and predictions. Dynamic GAMs (DGAMs) offer a solution to the challenge of forecasting discrete time series while estimating ecologically relevant nonlinear predictor associations. Our Bayesian latent factor approach will be particularly useful for exploring competing dynamic ecological models that encompass hierarchical smoothing structures while providing forecasts with robust uncertainties, tasks that are becoming increasingly important in applied ecology.


| INTRODUC TI ON
Rapidly changing climates and landscape modification are impacting species and ecosystems at all micro-and macroecological levels, incurring substantial economic and environmental costs (Kennedy et al., 2019;United Nations, 2015;World Health Organization, 2005).
There is broad consensus among scientists, parliamentarians and applied decision makers that anticipating probable future states is vital to mitigate the impacts of environmental change on ecosystem functionality and services Intergovernmental Panel on Climate Change, 2018;Schmidt et al., 2010).
Two challenges impede the improvement and adoption of common forecasting tools in ecology. First, ecosystems are driven by networks of interacting biotic and abiotic processes (Choler et al., 2001;Levin, 1998;Massoud et al., 2018). These dynamic natural processes are the products of multiple sources of variation including long-term trends, seasonal and other cyclic oscillations, environmental forcing, temporal dependence or species interactions (Auger-Méthé et al., 2021;Choler et al., 2001;Dietze, 2017). Second, ecological time series are often integer-valued variables, such as observations of species presence or abundance, that exhibit complex features including observation error, zero inflation, overdispersion, truncation at hard bounds, missing values and uneven sampling (Kowal & Canale, 2020;Lindén & Mäntyniemi, 2011;Simpson, 2018;Warton, 2018). Such discrete time series are far less supported in existing software than are real-valued series that can be readily modelled using assumptions of Gaussian error (Hyndman & Khandakar, 2008). Moreover, ecological observations are almost always multivariate when contextual information, such as data from environmental predictors or observations of non-target species, is considered. These features make it difficult to analyse ecological time series while sufficiently accounting for important systematic temporal components and multivariate dependencies (Auger-Méthé et al., 2021).
Time-series analyses are often concerned with decomposing temporal variation into components representing trend, seasonality and other cyclic changes. Generalised additive models (GAMs), which are increasingly used in ecology to identify nonlinear functional relationships (Guisan et al., 2002;Hughes et al., 2018;Pedersen et al., 2019;Simpson, 2018), offer a way to accomplish this decomposition. Outlined in detail previously (Hastie & Tibshirani, 1990;Wood, 2004), GAMs can briefly be described as modified generalised linear models (GLMs) in which the linear predictor includes a sum of smooth functions representing functional relationships between covariates and the response: where E(Y) is the conditional expec tation of a response assumed to be drawn from an exponential family distribution, 0 is the unknown intercept, the s i 's are a set of smooth functions over one or several predictor variables (the x's) and g is a monotonic link function. Each smooth function s i is composed of basis functions whose coefficients, which must be estimated, act as weights for the basis functions to control the function's shape. The total number of basis functions limits the potential complexity of the smooth function, with a larger set of basis functions allowing greater flexibility. In addition to their ability to represent complex and nonlinear ecological relationships, several other advantages of GAMs are that they can model a diversity of response families that accommodate ecological features (such as zero inflation) and that they can be formulated to include hierarchical smoothing for multivariate responses (Pedersen et al., 2019;Wood, 2017).
Given the set of basis coefficients that comprise each smooth function, a GAM can in principle be estimated as a GLM. However due to their incredible flexibility, GAMs will invariably overfit if left unconstrained (Hastie & Tibshirani, 1990;Marra & Wood, 2011;Wood, 2004). Penalised likelihood estimation avoids this overfitting by placing quadratic penalties on the basis coefficients (referred to as smoothing penalties), which penalise the function's 'wiggliness' and control the trade-off between fit and smoothness (Wood, 2004(Wood, , 2016. From a Bayesian perspective, another way to represent a smooth function is to draw the set of basis coefficients from a multivariate Gaussian distribution with the penalty acting on the prior precision. Larger penalties shrink the coefficient covariances, effectively forcing the function towards a straight line when the data do not justify a nonlinear relationship (Marra & Wood, 2011;Wood, 2016). GAMs are particularly sought after for modelling time series to identify nonlinear or time-varying covariate effects, perform smoothing of historical time series and uncover periods of rapid change, though strong temporal autocorrelation can make it challenging to estimate key parameters (Camara et al., 2021;Knape, 2016;Simpson, 2018;Spooner et al., 2018;Yang et al., 2012).
For many ecological studies that employ GAMs, a primary objective is predicting future states (Clark et al., 2020;Kaplan et al., 2016;Koolhof et al., 2021;Malick et al., 2020;Ward et al., 2014). However, a lingering issue in using GAMs for forecasting is the way in which smooth functions predict outside the range of training data. Many of the smooth functions used in ecological GAMs have zero second derivatives at the boundaries, meaning they will linearly extrapolate beyond the last observation (Elith et al., 2010;Zurell et al., 2012).
This projection of a straight line indefinitely into the future can produce unrealistic forecasts, particularly if the estimated function 'wiggles' (i.e. exhibits a pronounced change in the response-predictor relationship) near the boundary (Figure 1 top). There are technical solutions to help with this problem, for example, by extending the

| Univariate models for a single ecological time series
A Bayesian framework to model fitting and parameter estimation involves defining a joint probability distribution over all observable and unobservable quantities in a statistical model that aligns with expert beliefs about the data generating process . A DGAM is naturally viewed from a Bayesian perspective, where prior beliefs about the nonlinearity of a function can be elicited to inform the complexity and penalisation of the smooth (Miller, 2019;Pedersen et al., 2019;Wood, 2013) while accounting for possible unobserved temporal dependence in line with the expectation that time series evolves as serially autocorrelated dynamic processes (Hyndman & Athanasopoulos, 2018). In its basic form, the DGAM for a discrete integer-valued time series is written as: where E Y t is the conditional expectation of the response at time t and z t is the (latent) dynamic process estimate at time t. Readers familiar with state-space models will recognise the benefits of separating the temporal and observation processes (Auger-Méthé (3) z t = 0 + z t−1 + e t F I G U R E 1 Estimated trends and forecasts from two GAMs applied to a discrete time series. In the top panel, a thin plate regression spline with a penalised second derivative is used for the trend, leading to a smooth function (top left) and linear extrapolation when forecasting (top right). In the bottom panel, the trend penalty is placed on the first derivative, resulting in flat extrapolation when forecasting. Trend shading shows 95% confidence intervals, while forecast shading shows empirical quantiles. Both models were fitted to a simulated seasonal discrete time series in R using the mgcv package with the general formula: Y ~ s(year, bs = 'tp') + s(season, bs = 'cc') + ti(season, year), family = nb()).  , 2021;Heilman et al., 2022), but it is worth clarifying these advantages explicitly. First, estimating the trend as a dynamic random variable avoids problems that can occur in competing autoregressive observation models where measurement error or outliers can have large influences on estimated AR parameters and cause highly unstable forecasts (see an example in Appendix S1). Second, it is far easier to handle missing or irregularly sampled observations using latent processes. Because the z t are unobserved latent variables, they will continue to evolve, even when an observation Y t is missing, via dynamic equations that conveniently provide recursive expressions for h-step ahead prediction, historical filtering and updating of forecasts (Durbin & Koopman, 2012). In contrast, a missing observation in an AR3 observation model will result in NAs for four rows of the design matrix (one missing Y t and three missing AR predictors) that can make parameter estimation difficult for software that automatically excludes rows with missing values (such as commonly used linear modelling packages in R). Other advantages of a state-space form are that trend dynamics provide a probabilistic model for the temporal evolution of a process, which can often be more useful than a smoothed trend (such as a penalised spline) by facilitating simulation and comparison with other processes, allowing new observations to be assimilated to adapt a forecast distribution via recursive Kalman or particle filtering (Massoud et al., 2018) and providing a means for multiple observation processes to depend on shared latent processes (Ward et al., 2021).
In its simplest form, temporal dependence can be modelled as a random walk with possible drift, where 0 in Equation (3) is the drift parameter and the residual error e t is drawn from a zerocentred Gaussian distribution with a fixed (time-invariant) standard deviation. This can easily be expanded to include autoregressive (AR) processes. For example, the following specifies a latent AR2 model: The assumption of a fixed standard deviation for the process error could potentially be a limitation if the series of interest displays nonconstant volatility with perturbations that may be evidence of responses to 'shocks'.
The time series literature is rich with model specifications for accommodating dynamic distributional models, including stochastic volatility, GARCH or Lévy processes (Bartumeus, 2007;Carrasco & Chen, 2002). In sharp contrast, temporal dependence could also be modelled via a latent Gaussian process (or other stochastic process), which provides a nonparametric probability distribution over functions. Gaussian processes are particularly suitable for ecological time series where we often expect dynamics to evolve as a smooth function and we wish to estimate the co-

| Dynamic factor DGAMs for analysing multiple ecological time series
Here, we describe how a DGAM can be modified into a joint multivariate statistical model for collections of time series with potentially common dynamics. Dynamic factor models that account for relationships in time-series data are closely aligned with static latent factor models, which are used in quantitative ecology to jointly model abundances of multiple species by estimating shared responses to unmeasured ecological drivers (Ovaskainen et al., 2017;Thorson et al., 2016;Ward et al., 2021;Warton et al., 2015). A latent factor model is a function of unmeasured random predictors (factors) that induce correlations between multiple responses via factor loadings while exercising dimension reduction. Often, species do demonstrate correlated responses to environmental gradients, meaning that a smaller set of factors (i.e. a low-dimensional representation) than the total number of possible species-predictor relationships can adequately capture the main axes of covariation (Letten et al., 2015;Warton et al., 2015). A dynamic factor model assumes the factors evolve as time series. The strength of this approach is that a small number of common factors can often model the temporal behaviours of a much larger set of series. This dimension reduction simplifies the estimation and forecasting tasks, as only the smaller set of factors and the series' specific factor loadings need to be estimated to generate forecasts (De Stefani et al., 2019).
In a dynamic factor DGAM, each series' latent trend is composed of a linear combination of these common factors: where E Y j,t is the expected response for series j at time t, the z m,t 's are estimates for the M factors at time t and the j 's are factor loadings. As in the univariate case, the factors can evolve either as random walks with drift or as autoregressive processes up to order 3.
A challenge with any factor model is the need to determine the number of factors M (Bhattacharya & Dunson, 2011;Fox et al., 2009;Thorson et al., 2016;Tobler et al., 2019). Setting M too small prevents temporal dependencies from being adequately modelled, leading to poor convergence and difficulty estimating smooth parameters. By contrast, setting M too large leads to unnecessary computation. The problem can be approached by formulating a prior distribution that enforces exponentially increasing penalties on the factor variances to allow any un-needed factors to evolve as flat lines. Following Welty et al. (2009) andWells et al. (2016), one such prior assumes that factors up to a certain threshold number π have precisions of similar magnitudes, after which they increase exponentially (leading to variances that shrink towards zero). Along with π, two other hyperparameters can be estimated to control the baseline penalty and the rate at which penalties exponentially increase, respectively, allowing the data to inform the selection of dynamic factors. We caution, however,  (Wood, 2016). The model is modified to include dynamic components (either as random walk, AR trends up to order 3 or Gaussian processes) and to update any prior distributions specified by the user, while all data reformatting necessary for modelling is done automatically. Employing either the JAGS software through the R interface rjags (Plummer, 2003)  mvgam extends functions available in existing software in several ways. First, while fully Bayesian GAMs can be estimated using a variety of software including brms (Bürkner, 2017), BayesX (Brezger et al., 2005) and bamlss (Umlauf et al., 2018), mvgam is the only software we are aware of that can simultaneously estimate any smooth function available in mgcv together with latent dynamic trends (bamlss and BayesX can estimate a diversity of smooth functions, but to our knowledge, dynamic latent processes cannot be estimated; brms offers more flexibility for time series and can accommodate dynamic latent processes, including AR and ARMA processes of order 1, but we are not aware of extensions to dynamic factors). Second, our software can employ Hamiltonian Monte Carlo using Stan for much more efficient and unbiased MCMC sampling compared to Gibbs samplers (BayesX uses its own custom Gibbs samplers, while bamlss does not employ full MCMC). Perhaps the most important advantage of Hamiltonian Monte Carlo is the powerful diagnostics it provides for detecting posterior degeneracies, which can help uncover model inadequacies or incompatibilities between model and observed data (Betancourt, 2017). Finally, our package is designed for analysing and forecasting sets of discrete time series, and as such the additional utilities we offer for working with time series (including options to compare models using rolling forecast evaluation as well as routines to assimilate new observations 'online' for automatic forecast updating; Appendix S1) make our software attractive for a range of applied forecasting tasks.
It is notable that our design permits any formula allowed in mgcv to be used for the GAM component of the linear predictor, providing a user-friendly way to explore dynamic ecological models that encompass nonlinear smooth functions. Other advantages of our framework are (1) missing values are allowed for the responses; (2) upper bounds can be used via truncated likelihoods; (3) smooth distributed lag covariate functions can be estimated alongside latent temporal components to form complex dynamic nonlinear models (Gasparrini, 2011); and (4) dynamic components can easily be forecasted via their autoregressive equations (for random walk and AR trends) or via their estimated covariance functions (for Gaussian process trends), providing robust probabilistic uncertainties.
While the mvgam package does not currently support stochastic volatility or moving average trends, these processes could be added by the user at any time (the package can be used to generate all model files, data objects and initial values, so that a model can be easily modified for conditioning outside of mvgam, i.e. with rstan, rjags or other interfaces directly).

| S IMUL ATIONS
We used simulations to examine the performance of our DGAM formulation. Briefly, we simulated multiseries datasets with 72 time points (6 years of data for monthly series) consisting of Negative Binomial observations (size parameter = 5) for sets of series whose log-linear predictors included a hierarchical seasonal pattern (where each series' seasonal pattern was created by drawing from a global seasonal pattern with common Gaussian noise; see function sim_ mvgam() in the mvgam package for R code to produce simulations) and uncorrelated latent trends. Temporal dependences followed independent random walk processes. We investigated model sensitivity to missingness (proportion missing = 0, 10 or 50%), dimensionality (number of series = 2, 4 or 12) and the magnitude of the temporal component relative to seasonality (0.3 for moderate dynamics or 0.7 for strong dynamics; see Figure S1 for an example of two series with the same seasonality but different strengths of trend).
Each simulated dataset was fit with the same set of four models.
First, we fit a hierarchical GAM using mgcv that included a random intercept per series (s(series, bs = 're')), a cyclic smooth function for global seasonality (s(season, m = 2, k = 8, bs = 'cc')), local smooth functions for series-specific deviations from global seasonality (s(season, series, m = 1, k = 4, bs = 'fs')), a smooth function for a global trend (s(year, k = 4)) and local smooth functions for series-specific deviations from the global trend (s(year, series, m = 1, k = 4, bs = 'fs')). Our next model was a GAM (also fitted with mgcv) that used a stochastic trend via an autoregressive observation model. This model used same hierarchical seasonality smooths functions as the GAM above but replaced the trend smooths with an AR1 parametric term for the effect of log(y t−1 ), with separate AR1 terms estimated for each series. We chose to model the AR1 term on the log scale as this reduces sensitivity of the AR parameter estimates to outliers (see Appendix S1 for an investigation of the forecasting behaviours of autoregressive observation models for discrete time series). Note, however, that because each missing observation results in additional missing rows in the design matrix (due to missing values in AR predictors), we were unable to fit this model for the simulations where 50% of observations were missing. We next asked whether a dynamic factor process could capture the multiseries temporal dynamics by fitting a dynamic factor DGAM (with M = half the number of series) with identical random effect and seasonal smooth functions but no yearly smooth function. Finally, we fit a 'null' dynamic factor DGAM that only estimated random intercepts but no seasonal smooth function. Negative binomial distributions were specified for each model and AR1 models were used for modelling the DGAM dynamic factor processes. Each combination of missingness, dimensionality and strength of dynamics was used to generate five replicate datasets, yielding a total of 60 simulations. For mgcv models, estimation of smoothing penalties was performed using restricted maximum likelihood (method = 'REML'). Gaussian priors were specified for AR parameters ( ) (mean = 0; variance = 0.1) in the mvgam implementation. Following Wood (2016), zero-centred multivariate Gaussian priors were used for each smooth's ß parameters and exponential priors were used for the smoothing penalties.  (Gelman & Rubin, 1992) and by visual inspection of posterior chains.
The relative performances of each model were explored using out of sample forecasts. We trained models on the first 5 years of data (60 observations) and generated forecasts for the remaining year (12 observations). Probabilistic forecast performance was evaluated using a discrete version of the Rank Probability Score (DRPS; Gneiting & Raftery, 2007) and coverage of 90% prediction intervals.
Forecasts with lower DRPS and coverage closer to 0.9 were considered more accurate.

| C A S E S TUDY: FOREC A S TING TI CK ABUNDAN CE S
Amblyomma americanum and Ixodes scapularis are two widespread species of hard ticks capable of transmitting a diversity of parasites to animals and humans, many of which are zoonotic (Rochlin & Toledo, 2020). Due to the medical and ecological importance of these species, a common goal is to understand factors that in- Temperatures between −5°C and 5°C can affects various components of tick physiological diapause and host-seeking behaviours (Clark, 1995). We included a cumulative growing degree day (cum_gdd) variable using temperature records for each site's nearest weather station from NOAA's Daily Global Historical Climatology Network daily database as a covariate. The predictor was calculated as the total number of days up to the start of the tick season (1st June) in which the mean of the day's maximum and minimum temperatures was above 0°C. We fit species-specific DGAMs to 4 years of data (2015-2018) for 17 A. americanum plots (nested in seven NEON sites) and for eight I. scapularis plots (nested in three sites) using the most recent release of the NEON tick drag sampling product (National Ecological Observatory Network, 2022). Counts of ticks were aggregated at the temporal resolution of epidemiological week, a standardised method of counting weeks developed by the US Centers for Disease Control and Prevention to facilitate direct comparisons across years. Time points during winter (epidemiological weeks 1-14 and 41-53) had entirely missing observations as no sampling occurred during this period, but we kept these in the model as missing data. For each species, we fit four models representing different hypothetical dynamics, though we caution that our goal here was not to carry out a rigorous analysis but to highlight how DGAMs could be used to facilitate model selection and scrutiny: • Null: There is no seasonality, rather the latent factors/random site-level effects of cum_gdd fully influence the dynamics for the plot-level series. We hypothesised that the site-specific partial effects of cum_gdd could be mildly nonlinear, so we set k = 5 for this Comparisons of 90% interval coverages strongly favoured the two DGAMs (Figure 3). Intervals for the DGAMs frequently included 25%-35% more of the out of sample observations than did the intervals for the two GAMs. There was little distinction between the two DGAMs, even as the number of series and the strength of the underlying dynamics increased (Figure 3). Results were similar when inspecting 90% interval coverage as a function of missingness, with the DGAMs strongly outperforming the GAMs ( Figure S3).

| DGAM and NEON tick abundance forecasts
Our results suggested that Hyp3, which captured hierarchical seasonality by allowing individual plot-level seasonal patterns to deviate from a global seasonality function, was the best-performing model when forecasting I. scapularis nymphal abundance across NEON sites, while the null model that did not include seasonality was the worst performing ( Figure 4). Nominal coverages of 90% intervals were accurate for the three seasonal models (ranging from 87% to 88%), while the intervals for the null model were generally wider than they needed to be (97% coverage; Figure 4).
However, there was variation across plots in terms of forecast performance, suggesting that an ensemble forecast (which combines forecasts from multiple models) could improve performance ( Figure S4). Inspection of probability integral transform (PIT) histograms, which should be uniform if predictions are evenly distributed about the truth (Simonis et al., 2021), revealed that all models apart from the null tended to underpredict to some degree (left-skewed PIT histograms; Figure S5). Figure 5 shows example mvgam visualisations for a single plot, including estimated smooth functions, forecasts and dynamic trend estimates (along with their probabilistic uncertainties). When conditioning on seasonality and the trend, I. scapularis abundances demonstrated no apparent association with variation in cumulative growing degree days ( Figure 5). Inspection of the latent dynamic components for the three seasonal models revealed positive within-site correlations for sites SCBI and SERC ( Figure S6). Example mvgam visualisations of posterior checks for training (retrodictive) and forecast periods (predictive), useful for checking if a model is capable of simulating time series that resemble key aspects of the observed data without notable discrepancies, are shown in Figure S7. Examples highlighting how smooth function and trend realisations can be plotted, which can improve model interpretation over quantile or density plots, are shown in Figure S8.
In agreement with the I. scapularis models, A. americanum abundance was also best predicted by the Hyp3 model. Example visualisations of estimated plot-level seasonal functions are shown in Figure 6, while a visualisation of estimated random effect intercept distributions is shown in Figure S9. Our model estimated that tick abundances in some plots (i.e. SERC_001) tended to show earlier peaks around epidemiological week 24, while abundance in other plots (i.e. TALL_001) followed a broader curve with a peak around epidemiological week 30 ( Figure 6).

| Quantifying uncertainty contributions among mvgam model components
In addition to plotting smooth functions and forecasts, mvgam offers utilities to compute relative contributions of the latent dynamic and GAM components to forecast uncertainty. This process of partitioning uncertainty is an important step in analysing a model's forecasts to diagnose the main drivers of prediction uncertainty and prioritise aspects of models or data that require further investigation (Dietze, 2017;Heilman et al., 2022). Comparisons of uncertainty contributions for four of the A. americanum forecasts indicate that both components contribute to forecast uncertainty, but to varying degrees over time and across plots ( Figure 7). However, across all plots, dynamic trend uncertainty tended to increase over time, becoming relatively more important during the peak tick season (3-22 weeks ahead).

| DISCUSS ION
We have introduced an R package for fitting Bayesian Dynamic GAMs (DGAMs) that incorporate the flexibility of the widely popular penalised smoothing functions in mgcv with latent dynamic components for analysing and forecasting discrete time series. Keys to mvgam's performance are its ability to cope with substantial missing data, scale to large collections of discrete time series and provide robust uncertainty quantification. In recent years, there has been increased interest in using time-series models for uncertainty interval estimation as opposed to point predictions, a trend that lends well to Bayesian inference Makridakis et al., 2020). This is particularly relevant for ecological forecasts, F I G U R E 2 Log(discrete rank probability score) (DRPS) performance for out of sample forecasts from competing models fitted to sets of simulated discrete time series. Panels depict models fitted with different levels of data missingness (proportion of observations set to NA) and temporal dynamics strength. The seasonal GAM was fitted using R package mgcv, while the seasonal and nonseasonal DGAMs were fitted using the mvgam package (using the Hamiltonian Monte Carlo software Stan). Lower scores indicate better model performance.  where point estimates are less important for making informed decisions than are conditional probability statements Dietze, 2017;Dietze et al., 2018;White et al., 2019 processes and other models of stochastic processes, to increase the diversity of models that can be interrogated using mvgam • The addition of Markov-switching processes to allow dynamic factor loadings to be drawn from different sets of correlation 'regimes', allowing correlation structures to change over time in a principled way (Fox et al., 2010) • The incorporation of covariates into the latent temporal models (i.e. as dynamic linear models) to explicitly address broader hypotheses about the factors that influence temporal dynamics (Heilman et al., 2022)

| Challenges in estimating DGAM parameters
The joint estimation of smoothing parameters, basis coefficients, latent trend variances or overdispersion parameters is not without its challenges (Wood, 2016). Posterior geometries for such high-dimensional models can become complex enough that traditional MCMC samplers based on Random Walk proposals (e.g. Gibbs samplers) will not be able to sample the parameter space without reverting to painfully small step sizes that result in high posterior autocorrelation and very slow exploration (Betancourt, 2017). Maximum likelihood and related estimators will not likely produce better uncertainty quantification, as verifying how and when posterior geometries can be accurately approximated under an asymptotic regime is a huge and elusive challenge. mvgam's exploitation of Hamiltonian Monte Carlo is a major advantage for tackling DGAM parameter estimation, but we stress that there is no one-size-fits-all default F I G U R E 7 Output from the plot_ mvgam_uncertainty function in mvgam showing relative contributions of the dynamic temporal (grey) and GAM (red) components to forecast uncertainty for four Amblyomma americanum plots estimated from a dynamic GAM with hierarchical seasonality. Forecast horizons were varied over a '1-year' horizon (52 weeks matching data availability).   (Gelman et al., 2020). One illuminating situation that we have encountered is the difficulty in jointly estimating a latent trend and overdispersion parameters such as in the Negative Binomial or Tweedie distributions. This is because both processes (overdispersion and autocorrelation) may be able to explain the dispersion around the mean, particularly when using Random Walk or AR trends that can jump around easily. Users will need to use theory and judgement to decide how to tackle these challenges, for example, by assuming there is overdispersion in the observation process (with consultation from appropriate references; i.e. Bliss & Fisher, 1953, Lindén & Mäntyniemi, 2011 but that the trend is smooth, in which case a latent Gaussian process with suitable prior on the length scale would be appropriate. Smoothing splines are also challenging in a way because they do not readily facilitate principled prior modelling, where expert elicitation could help to constrain prior function shapes towards those that are compatible with domain expertise as part of a Bayesian workflow (Betancourt, 2021;Gelman et al., 2020). Users are recommended to refer to the wealth of material relating to the mgcv package for choosing a smoothing basis and basis dimension that are compatible with expected function shapes (Wood, 2004(Wood, , 2013(Wood, , 2017.

| CON CLUS ION
The R package mvgam provides a user-friendly tool for researchers and practitioners interested in fitting DGAMs to analyse and forecast ecological time series. The problems associated with smooth spline extrapolation are not limited to ecology however, as the need to forecast sets of discrete nonlinear time series is a common challenge in areas as diverse as speech recognition, tourism demand, natural language processing and finance (Hyndman & Athanasopoulos, 2018;Makridakis et al., 2018). Beyond the examples showcased here, the package can be especially useful to identify avenues for model improvement via its ability to assimilate new observations online to update forecast distributions (showcased in Appendix S1). With growing interest in both the application of hierarchical GAMs to ecological problems and the need to use iterative forecasts to make ecology a more predictive discipline, mvgam can become a vital addition to the applied ecologist's analytical toolbox.

AUTH O R CO NTR I B UTI O N S
Nicholas J. Clark involved in conceptualization, data curation, project administration, software, visualisation, writing-original draft.
Nicholas J. Clark and Konstans Wells involved in formal analysis, methodology, validation, writing-review & editing.

ACK N OWLED G EM ENTS
NOAA temperature data were supplied by Daniel Ruiz-Carrascal as part of the 2021 NEON Ecological Forecasting Challenge. This research was funded by an ARC DECRA fellowship to N. Clark (DE210101439).

CO N FLI C T O F I NTE R E S T
None of the authors have a conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.

DATA AVA I L A B I L I T Y S TAT E M E N T
The manuscript uses data that are archived by the National Ecological Observatory Network (https://data.neons cience.org/).
The data have been downloaded and converted into a usable format for modelling, and this version of the data is available, along with an archived version of the mvgam R package and all R scripts used to produce analyses in this manuscript (and in the Appendices), at a permanently archived Zenodo repository (https://doi.org/10.5281/ zenodo.6918047; Clark & Wells, 2022).