EcoEnsemble: A general framework for combining ecosystem models in R

Often there are several complex ecosystem models available to address a specific question. However, structural differences, systematic discrepancies and uncertainties mean that they typically produce different outputs. Rather than selecting a single ‘best’ model, it is desirable to combine them to give a coherent answer to the question at hand. Many methods of combining ecosystem models assume that one of the models is exactly correct, which is unlikely to be the case. Furthermore, models may not be fitted to the same data, have the same outputs, nor be run for the same time period, making many common methods difficult to implement. In this paper, we use a statistical model to describe the relationship between the ecosystem models, prior beliefs and observations to make coherent predictions of the true state of the ecosystem with robust quantification of uncertainty. We introduce EcoEnsemble, an R package that takes advantage of the statistical model's structure to efficiently fit the ensemble model, either sampling from the posterior distribution or maximising the posterior density. We demonstrate EcoEnsemble by investigating what would happen to four fish species in the North Sea under future management scenarios. Although developed for applications in ecology, EcoEnsemble can be used to combine any group of mechanistic models, for example in climate modelling, epidemiology or biology.

. Furthermore, by ignoring other available models, the amount of information utilised is limited, leading to an increase in uncertainty (Spence et al., 2018), often to a greater extent than formally recognised.
Instead of choosing one model, it is possible to combine them using an ensemble model. However, many such methods, such as weighting schemes (e.g. Bayesian model averaging, Banner & Higgs, 2017;Dormann et al., 2018), explicitly assume that one of the models correctly describe the truth (Chandler, 2013). Not only is this a strong assumption, but treating a suite of models as containing the true model tends to underestimate the uncertainty because alternative, possibly as yet undeveloped, models could give predictions outside of the current range (Chandler, 2013). Under this scheme, adding a new model to the suite could increase the uncertainty of the ensemble despite the additional information contained in the new model (Dormann et al., 2018).
Applying the above methods to complex ecosystem models is not straightforward. Often these models are fitted to different data (Ianelli et al., 2016), describe different processes and produce outputs that are not directly comparable (Spence et al., 2018). At the same time, complex ecosystem models often share many similarities. They can have similar or even identical processes and forcing inputs (e.g. Tittensor et al., 2018), or can be fitted using similar data, suggesting that the discrepancies in the models are not independent (Christiansen, 2021;Knutti, 2010;Rougier et al., 2013). Spence et al. (2018) developed an ensemble model adopting a Bayesian approach that treats the individual models as exchangeable and drawn from a population of possible models. They separated the discrepancy between the models and the true variables of interest (VoI) into discrepancy shared between all of the models and discrepancies that are specific to each model. By modelling in this way, the ensemble model exploits the strengths whilst discounting the weaknesses of each individual simulator (Chandler, 2013) while rigorously quantifying the uncertainty allowing for better management decisions (Spence, Griffiths, et al., 2021).
In this paper, we introduce EcoEnsemble, an R package that implements the ensemble model of Spence et al. (2018) for time-series outputs of complex models, where the models output linear combinations of the VoI. We introduce EcoEnsemble and the ensemble model in Section 2 and then demonstrate it on case study of fishing in the North Sea in Section 3.

| US ING ECOEN S EMB LE
EcoEnsemble is an R package for fitting and sampling from the ensemble model described in Spence et al. (2018) For each simulator, the 'best guess' is the VoI plus a discrepancy The discrepancy term is split between discrepancies that are shared by all of the simulators and discrepancies that are specific to the kth simulator. The shared discrepancy is split between the long-term shared discrepancy, , and the short-term shared discrepancy, (t) . Simulator k's individual discrepancy is spilt between and long-term individual discrepancy, k , and short-term individual discrepancy, z (t) k , that is The long-term individual discrepancies are normally distributed so that the long-term individual discrepancy of the kth simulator is The short-term discrepancy terms, (t) and z (t) k , follow stationary autoregressive models of order one, with autoregressive parameters R and R k , and covariance parameters Λ and Λ k , respectively.
A summary of the ensemble model is found in Table 1. For more details on the model, see Spence et al. (2018) or the 'EcoEnsemble' vignette in EcoEnsemble.

| Data requirements
In EcoEnsemble, the user requires: • observations of the VoI (ŷ (t) ) and a covariance matrix describing observational uncertainty (Σ y ).
• outputs from each of the simulators (x (t) k ) and covariance matrices describing the parameter uncertainty of each simulator (Σ k ).
• prior beliefs about the discrepancy terms.

| Observations
Observations of the VoI at some times are included as an input. The ensemble model says noisy observations of the VoI are (1) (3)

| Simulator outputs
Simulators should output linear combinations of the VoI. They do not need to produce outputs for the whole period, but a given simulator should output the same variables for each time that it does. The simulator outputs are noisy linear combinations of the 'best guess'. If the kth simulator is evaluated at time t, its output is where Σ k reflects the parameter uncertainty of the kth simulator and M k is an n k × d matrix relating the 'best guess' to the simulator output. If, for all simulators, VoI are either present or not, then the fit_ensemble_ model() function infers M k for k = 1, … , m. If other transformations of the VoI are outputted from the simulators, then the M k matrices should be specified by the user.
Parameter uncertainty for simulator outputs should also be quantified in the covariance matrices Σ k . This can be calculated, for example, by calibrating the simulator to empirical data (e.g. Spence,

| Prior beliefs
To include prior beliefs in EcoEnsemble, the user should place priors on m + 3 different covariance matrices: the covariance of the random walk of the VoI Σ y ; the covariance of the auto-regressive process on the shared short-term discrepancy, Λ ; the covariance of the long-term individual discrepancies, C ; and the covariances of the auto-regressive processes on the individual short-term discrepancies, Λ 1 , … , Λ m . The prior on the long-term shared discrepancy, , also needs to be specified.
For Σ y , the only choice for prior distribution is an inverse-Wishart distribution.
In EcoEnsemble, we follow the approach of Spence et al. (2018) and decompose each of the discrepancy covariance matrices and P is a correlation matrix.
In EcoEnsemble, for C and Λ , there are three available prior distributions for covariance matrices. In each case, the variance terms are parameterised by a gamma distribution, while the correlation matrices can be parameterised either by an LKJ distribution (Lewandowski et al., 2009), an inverse-Wishart distribution, or by beta distributions, using the method of concordance (Gokhale & Press, 1982;Zondervan-Zwijnenburg et al., 2017).
TA B L E 1 A summary of the variables in the ensemble model. Values of n k and M k are simulator specific.

| Fitting the model
The model described in Table 1 can be written as a dynamical linear model (see Supporting Information), meaning that the latent varia- for all t, can be integrated out of the likelihood. EcoEnsemble uses the sequential Kalman filter (Chui & Chen, 2009;Kalman, 1960) to fit the ensemble model (see the Supporting Information for more details). The highest posterior density can be found or a full sample can be generated using the No U-

| Model outputs
EcoEnsemble includes functionality for sampling (t) , as well as the most likely (t) , from an EnsembleFit object. For each set of parameters from a fitted model, we calculate the most likely latent variables and sample them using the algorithm developed by Strickland et al. (2009), based on Durbin and Koopman (2002).
generate_sample() creates an EnsembleSample object, which contains samples of (t) . Calling the plot() function with an EnsembleSample object plots the VoI, y (t) , the simulators outputs, x (t) k , and the observations, ŷ (t) , for all t and k (see Figure 2). samples <-generate_sample(fit) plot(samples)

| C A S E S TUDY
We demonstrate EcoEnsemble using an example from fisheries. We investigated the change in spawning stock biomass (SSB) of four species in the North Sea, Norway pout Trisopterus esmarkii, Atlantic herring Clupea harengus, Atlantic cod Gadus morhua and common sole Solea solea if fished at maximum sustainable yield (MSY), a management strategy that maximises the long-term yield for each species, from 2017 (Spence, Griffiths, et al., 2021). We used single-species stock assessments and four multispecies simulators to predict the SSB from 1984 to 2050 under this scenario. The data for this study are available in EcoEnsemble.

| Observations
Noisy observations of the SSB were taken from single-species stock assessments (ICES, 2020a(ICES, , 2020b. Uncertainty provided in stock assessments formed the diagonal elements of Σ y , with the off-diagonal elements being zero.

| Simulators
SSB was also predicted using four simulators: 2. LeMans is a discrete time length-based model that describes growth and predation (Thorpe et al., 2015). It outputs all four species from 1986 until 2050, meaning that n 2 = 4 and M 2 = I 4 . Its covariance was calculated in Thorpe et al. (2015).

EcoPath with EcoSim
3. mizer is a size-based model that describes ontogenetic feeding and growth, mortality and reproduction driven by size-dependent predation and maturation processes (Blanchard et al., 2014).
It outputs all four species from 1984 until 2050, meaning that n 3 = 4 and M 3 = I 4 , and its covariance was calculated in Spence et al. (2016). 4. FishSUMS is a discrete time length-based model that describes growth, density-dependent mortality and losses due to fishing and predation by explicitly modelled species, and seasonal reproduction (Speirs et al., 2016). It outputs Norway pout, Atlantic herring and Atlantic cod, that is it does not output sole, from 1984 until 2050, meaning that n 4 = 3 and that Its covariance was calculated in Spence et al. (2018).

| Prior information
We used the default priors in EcoEnsemble to fit the model. See Supporting Information and Exploring priors for more information. The result of the prior predictive of the ensemble model is shown in Figure 1 using the prior_ensemble_model() and sample_prior_mode() functions.

| Results and discussion
The ensemble model prediction, simulator outputs and stock assessment estimates are shown in Figure 2. The probability that the SSB of each species increases is given in Table 2.
The simulators appear to capture the general trend correctly for most species but often struggle to capture short-term dynam-

| CON CLUS IONS
Complex models are rapidly emerging as an invaluable tool for managing natural resources (Hyder et al., 2015). However, when models F I G U R E 1 A plot of the natural logarithm of the prior predictive spawning-stock biomass of each species considered in the study, as predicted by each simulator, single-species assessments and the ensemble model. The edges of the shaded area represent the 5% and 95% quantiles of the prior predictive SSB.

TA B L E 2
The ensemble model probability that the spawning stock biomass for species in 2050 is larger than the year in the first column. All authors contributed critically to the drafts and gave final approval for publication.

ACK N O WLE D G E M ENTS
This work was supported by the Centre for the Environment, Fisheries and Aquaculture Science Seedcorn project DP436 'Operationalising ensemble modelling for multispecies stock assessments'. The authors would also like to thank Robert Thorpe for comments on earlier versions of the paper.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.

PE E R R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-review/10.1111/ 2041-210X.14148.

DATA AVA I L A B I L I T Y S TAT E M E N T
EcoEnsemble can be downloaded from GitHub (https://github.com/ Cefas RepRe s/EcoEn semble) and the Comprehensive R Archive Network CRAN (https://cran.r-proje ct.org/web/packa ges/EcoEn sembl e/index.html) (Spence et al., 2023).

F I G U R E 2
A plot of the posterior predictive natural logarithm of the spawning-stock biomass of each species considered in the study, as predicted by each simulator, single-species assessments and the ensemble model. The edges of the shaded area represent the 5% and 95% quantiles of the posterior predictive SSB.