Title:	A General Framework for Combining Ecosystem Models
Version:	1.1.2
Description:	Fit and sample from the ensemble model described in Spence et al (2018): "A general framework for combining ecosystem models"<doi:10.1111/faf.12310>.
License:	GPL (≥ 3)
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.1
Biarch:	true
URL:	https://github.com/CefasRepRes/EcoEnsemble
BugReports:	https://github.com/CefasRepRes/EcoEnsemble/issues
Depends:	R (≥ 3.5.0)
Imports:	methods, Rcpp, matrixcalc, RcppParallel (≥ 5.0.1), rstan (≥ 2.26.0), rstantools (≥ 2.1.1), dplyr, ggplot2, reshape2, tibble, cowplot
LinkingTo:	BH (≥ 1.66.0), Rcpp, RcppEigen (≥ 0.3.3.3.0), RcppParallel (≥ 5.0.1), rstan (≥ 2.26.0), StanHeaders (≥ 2.26.0)
SystemRequirements:	GNU make
Suggests:	rmarkdown, knitr, testthat (≥ 3.0.0), mgcv
Config/testthat/edition:	3
VignetteBuilder:	knitr
Collate:	'DiscrepancyPrior-IndLTPrior-class.R' 'DiscrepancyPrior-IndSTPrior-class.R' 'DiscrepancyPrior-ShaSTPrior-class.R' 'DiscrepancyPrior-TruthPrior-class.R' 'EcoEnsemble-package.R' 'EnsemblePrior-class.R' 'EnsembleData-class.R' 'EnsembleFit-class.R' 'EnsembleSample-class.R' 'Prior_sampling.R' 'RcppExports.R' 'data.R' 'fit_ensemble.R' 'fit_ensemble_withdrivers.R' 'generate_simulator_stan_data.R' 'generate_simulator_stan_data_withdrivers.R' 'get_stan_outputs.R' 'get_stan_outputs_withdrivers.R' 'plot_functions.R' 'plot_functions_withdrivers.R' 'prior_functions.R' 'stanmodels.R' 'validation_functions.R' 'validation_functions_withdrivers.R' 'validation_prior_functions.R'
NeedsCompilation:	yes
Packaged:	2025-03-18 16:48:05 UTC; MS23
Author:	Michael A. Spence [aut, cre], James A. Martindale [aut], Michael J. Thomson [aut], Thomas I. J. Bartos [aut]
Maintainer:	Michael A. Spence <michael.spence@cefas.gov.uk>
Repository:	CRAN
Date/Publication:	2025-03-18 18:20:02 UTC

A general framework for combining ecosystem models

Description

The EcoEnsemble package implements the framework for combining ecosystem models laid out in Spence et al (2018).

Details

The ensemble model can be implemented in three main stages:

1. Eliciting priors on discrepancy terms: This is done by using the EnsemblePrior constructor.

2. Fitting the ensemble model: Using fit_ensemble_model with simulator outputs, observations and prior information. The ensemble model can be fit, obtaining either the point estimate, which maximises the posterior density, or running Markov chain Monte Carlo to generate a sample from the posterior denisty of the ensemble model.

3. Sampling the latent variables from the fitted model: Using generate_sample with the fitted ensemble object, the discrepancy terms and the ensemble's best guess of the truth can be generated. Similarly to fit_ensemble_model, this can either be a point estimate or a full sample.

Author(s)

Maintainer: Michael A. Spence michael.spence@cefas.gov.uk (ORCID)

Authors:

• James A. Martindale (ORCID)

• Michael J. Thomson (ORCID)

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. https://mc-stan.org

Spence, M. A., J. L. Blanchard, A. G. Rossberg, M. R. Heath, J. J. Heymans, S. Mackinson, N. Serpetti, D. C. Speirs, R. B. Thorpe, and P. G. Blackwell. 2018. "A General Framework for Combining Ecosystem Models." Fish and Fisheries 19: 1013-42. https://onlinelibrary.wiley.com/doi/abs/10.1111/faf.12310

See Also

Useful links:

• https://github.com/CefasRepRes/EcoEnsemble

• Report bugs at https://github.com/CefasRepRes/EcoEnsemble/issues

Constructor for the EnsembleData class

Description

A constructor for the EnsembleData class. This is used to convert input data into the required form for fit_ensemble_model.

Usage

EnsembleData(observations, simulators, priors, drivers = FALSE, MMod)

Arguments

Value

An object of class EnsembleData

See Also

EnsembleData, EnsemblePrior, fit_ensemble_model

Examples

ensemble_data <- EnsembleData(observations = list(SSB_obs, Sigma_obs),
                             simulators = list(list(SSB_ewe, Sigma_ewe, "EwE"),
                                           list(SSB_fs,  Sigma_fs, "FishSUMS"),
                                           list(SSB_lm,  Sigma_lm, "LeMans"),
                                           list(SSB_miz, Sigma_miz, "mizer")),
                              priors = EnsemblePrior(4))

A class to hold the Ensemble data

Description

A class that holds the observation data, simulator outputs, and prior information to convert into the required form for fit_ensemble_model.

Slots

stan_input

A list of parameters in the correct form to fit the ensemble model in Stan.

observations

A list of length 2 containing observations and a covariance matrix. The first element is a data.frame or matrix with each column giving observations of each output of interest and each row a time. Rows should be named with the times and columns should be named the variables. The second element is the covariance matrix of the observations.

simulators

A list with length equal to the number of simulators. For each simulator, there is a list of 2 objects containing the simulator output and covariance matrix. The first element is a data.frame or matrix with each column giving a simulator outputs of interest and each row a time. Rows should be named with the times and columns should be named the variables. The second element is the covariance matrix of the simulator outputs.

priors

An EnsemblePrior object specifying the prior distributions for the ensemble.

See Also

EnsembleData, EnsemblePrior, fit_ensemble_model

The constructor for the EnsembleFit object

Description

The constructor for an EnsembleFit object. This function need not be called as an EnsembleFit object is constructed automatically by the fit_ensemble_model function. The samples slot contains the samples from the MCMC if a full sampling was completed, otherwise the point_estimate slot contains information about a point estimate.

Usage

EnsembleFit(ensemble_data, samples = NULL, point_estimate = NULL)

Arguments

Value

An object of class EnsembleFit

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. https://mc-stan.org

See Also

EnsembleFit, fit_ensemble_model,

A class to hold the samples or point estimates from the ensemble model.

Description

An EnsembleFit object is returned by the fit_ensemble_model function. The object contains a slot for the EnsembleData object originally used to fit the ensemble model. The samples slot contains the samples from the MCMC if a full sampling was completed, otherwise the point_estimate slot contains information about a point estimate.

Slots

ensemble_data

An EnsembleData object encapsulating the data used to fit the ensemble model.

samples

A stanfit object containing the samples drawn from the fitted model. The default value is NULL.

point_estimate

A list output of the optimised model. The default value is NULL.

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. https://mc-stan.org

See Also

EnsembleFit, fit_ensemble_model, generate_sample

A class to hold the priors for the ensemble model.

Description

An EnsemblePrior object encapsulates the prior information for the ensemble model.

Slots

d

A numeric specifying the number of variables of interest in the ensemble.

ind_st_params

A list containing a prior specification for the individual short-term discrepancies z_k^{(t)}. See details of the EnsemblePrior() constructor.

ind_lt_params

A list containing a prior specification for the individual long-term discrepancies \gamma_k. See details of the EnsemblePrior() constructor.

sha_st_params

A list containing a prior specification for the shared short-term discrepancies \eta^{(t)}. See details of the EnsemblePrior() constructor.

sha_lt_params

A numeric containing the standard deviations for the normal prior used on the shared short-term discrepancy \mu. If a single value is supplied, this is repeated for each variable

truth_params

A list containing a prior specification for the processes on the truth y^{(t)}. See details of the EnsemblePrior() constructor. The default value is TruthPrior(d).

priors_stan_input

A list containing the prior data in the correct form to fit the model in Stan. This information is automatically generated by the constructor.

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. https://mc-stan.org

Lewandowski, Daniel, Dorota Kurowicka, and Harry Joe. 2009. “Generating Random Correlation Matrices Based on Vines and Extended Onion Method.” Journal of Multivariate Analysis 100: 1989–2001.

A constructor for the EnsembleSample object

Description

A constructor for the EnsembleSample class. These objects are generated automatically using the generate_sample function.

Usage

EnsembleSample(ensemble_fit, mle, samples)

Arguments

Value

An object of class EnsembleSample

See Also

EnsembleSample, generate_sample

A class to hold samples of the ensemble model

Description

EnsembleSample objects are generated using the generate_sample function.

Slots

ensemble_fit

An EnsembleFit object containing the fitted ensemble model.

mle

An array of dimension T \times (M + 2)\times N_{sample} containing MLE point estimates from the ensemble_fit object, where T is the total time, M is the number of simulators and N_{sample} is the number of samples. For each time step, the tth element of the array is a matrix where each column is a sample and the rows are the variables:

\left( y^{(t)}, \eta^{(t)}, z_1^{(t)}, z_2^{(t)}, \ldots, z_M^{(t)}\right)'

where y^{(t)} is the ensemble model's prediction of the latent truth value at time t, \eta^{(t)} is the shared short-term discrepancy at time t, z_i^{(t)} is the individual short-term discrepancy of simulator i at time t.

samples

An array of dimension T \times (M + 2)\times N_{sample} containing samples from the ensemble_fit object, where T is the total time, M is the number of simulators and N_{sample} is the number of samples. For each time step, the tth element of the array is a matrix where each column is a sample and the rows are the variables:

\left( y^{(t)}, \eta^{(t)}, z_1^{(t)}, z_2^{(t)}, \ldots, z_M^{(t)}\right)'

where y^{(t)} is the ensemble model's prediction of the latent truth value at time t, \eta^{(t)} is the shared short-term discrepancy at time t, z_i^{(t)} is the individual short-term discrepancy of simulator i at time t.

See Also

EnsembleSample, generate_sample

Constructor for the EnsemblePrior class

Description

Constructors for the EnsemblePrior class and related classes. These functions are used to encode prior information for the ensemble model. The IndSTPrior, IndLTPrior, ShaSTPrior, and TruthPrior constructors encapsulate prior information.

Usage

IndLTPrior(
  parametrisation_form = "lkj",
  var_params = list(1, 1),
  cor_params = 1
)

IndSTPrior(
  parametrisation_form = "hierarchical",
  var_params = list(-3, 1, 8, 4),
  cor_params = list(0.1, 0.1, 0.1, 0.1),
  AR_params = c(2, 2)
)

ShaSTPrior(
  parametrisation_form = "lkj",
  var_params = list(1, 10),
  cor_params = 1,
  AR_params = c(2, 2)
)

TruthPrior(
  d,
  initial_mean = 0,
  initial_var = 100,
  rw_covariance = list(2 * d, diag(d))
)

EnsemblePrior(
  d,
  ind_st_params = IndSTPrior(),
  ind_lt_params = IndLTPrior(),
  sha_st_params = ShaSTPrior(),
  sha_lt_params = 5,
  truth_params = TruthPrior(d)
)

Arguments

Details

IndSTPrior and ShaSTPrior discrepancy prior parameter objects contain 4 slots corresponding to:

1. parametrisation_form - A character specifying how the priors are parametrised. Currently supported priors are 'lkj', 'inv_wishart', 'beta', 'hierarchical', or 'hierarchical_beta_conjugate' ('hierarchical' and 'hierarchical_beta_conjugate' are only supported for IndSTPrior objects).

2. var_params - The prior parameters for the discrepancy variances, either a list of length 2 or a numeric of length 4. See below.

3. cor_params - The correlation matrix parameters, either a list of length 2, a numeric of length 3 or a numeric of length 4. See below.

4. AR_params - Parameters for the autoregressive parameter as a numeric of length 2.

IndLTPrior discrepancy prior parameter objects contain the slots parametrisation_form, var_params, and cor_params.

There are currently five supported prior distributions on covariance matrices. As in Spence et. al. (2018), the individual and shared short-term discrepancy covariances, \Lambda_k and \Lambda_\eta, as well as the individual long-term discrepancy covariance, \Lambda_\gamma, are decomposed into a vector of variances and a correlation matrix

\Lambda = \sqrt{\mathrm{diag}(\pi)} P \sqrt{\mathrm{diag}(\pi)},

where \pi is the vector of variances for each variable of interest (VoI), and P is the correlation matrix.

Selecting 'lkj', 'inv_wishart', 'beta', 'hierarchical' or 'hierarchical_beta_conjugate' refers to setting LKJ, inverse Wishart, beta or hierarchical (with gamma-distributed hyperparameters or beta-conjugate-distributed hyperparameters) prior distributions on the covariance matrix respectively. The variance parameters should be passed through as the var_params slot of the object and the correlation parameters should be passed through as the cor_params. For 'lkj', 'inv_wishart', and 'beta' selections, variances are parameterised by gamma distributions, so the var_params slot should be a list of length two, where each element gives the shape and rate parameters for each VoI (either as a single value which is the same for each VoI or a numeric with the same length as the number of VoI). For example, setting var_params = list(c(5,6,7,8), c(4,3,2,1)) would correspond to a Gamma(5, 4) prior on the variance of the first VoI, a Gamma(6, 3) prior on the variance of the second VoI, etc... The correlations should be in the following form:

• If 'lkj' is selected, then cor_params should be a numeric \eta giving the LKJ shape parameter, such that the probability density is given by (Lewandowski et. al. 2009)

  f(\Sigma | \eta)\propto \mathrm{det} (\Sigma)^{\eta - 1}.

  Variances are parameterised by gamma distributions.

• If 'inv_wishart' is selected, then cor_params should be a list containing a scalar value \eta (giving the degrees of freedom) and a symmetric, positive definite matrix S (giving the scale matrix). The dimensions of S should be the same as the correlation matrix it produces (i.e d \times d where d is the number of VoI). The density of an inverse Wishart is given by

  f(W|\eta, S) = \frac{1}{2^{\eta d/2} \Gamma_N \left( \frac{\eta}{2} \right)} |S|^{\eta/2} |W|^{-(\eta + d + 1)/2} \exp \left(- \frac{1}{2} \mathrm{tr}\left(SW^{-1} \right) \right),

  where \Gamma_N is the multivariate gamma function and \mathrm{tr \left(X \right)} is the trace of X. Note that inverse Wishart distributions act over the space of all covariance matrices. When used for a correlation matrix, only the subset of valid covariance matrices that are also valid correlation matrices are considered. Variances are parameterised by gamma distributions.

• If 'beta' is selected, then cor_params should be a list containing two symmetric d\timesd matrices A and B giving the prior success parameters and prior failure parameters respectively. The correlation between the ith and jth VoI is \rho_{i, j} with

  \frac{1}{\pi} \tan^{-1} \frac{\rho_{i, j}}{\sqrt{1-\rho_{i, j}^2}} + \frac{1}{2} \sim \mathrm{beta}(A_{i, j}, B_{i, j}).

  Variances are parameterised by gamma distributions.

• If 'hierarchical' or 'hierarchical_beta_conjugate' is selected, then variances are parameterised by log-normal distributions:

  \log \pi_{k, i} \sim \mathrm{N}(\mu_i, \sigma^2_i)

  with priors

  \mu_i \sim \mathrm{N}(\alpha_\pi, \beta_\pi),

  \sigma^2_i \sim \mathrm{InvGamma}(\gamma_\pi, \delta_\pi).

  The var_params slot should then be a numeric of length 4, giving the \alpha_\pi, \beta_\pi, \gamma_\pi, \delta_\pi hyperparameters respectively. Correlations (\rho_{k, i, j} where \rho_{k, i, j} is the correlation between VoI i and j for the kth simulator) are parameterised by hierarchical beta distributions.

  \frac{\rho_{k, i, j} + 1}{2} \sim \mathrm{beta}(c_{k, i, j}, d_{k, i, j})

  with priors

  c_{k, i, j} \sim \mathrm{gamma}(\alpha_\rho, \beta_\rho),

  d_{k, i, j} \sim \mathrm{gamma}(\gamma_\rho, \delta_\rho).

  NOTE: These options is only supported for the individual short-term discrepancy terms.

• If 'hierarchical' is selected, then the cor_params slot should be a numeric of length 4 giving the \alpha_\rho, \beta_\rho, \gamma_\rho, \delta_\rho hyperparameters. respectively. NOTE: This option is only supported for the individual short-term discrepancy terms.

• If 'hierarchical_beta_conjugate' is selected, then the cor_params slot should be a numeric of length 3. Denoting the values by r,s,k, they map to the hyperparameters p, q, k of the beta conjugate distribution via k = k, p^{-1/k} = (1+e^{-s})(1+e^{-r}) and q^{1/k} = e^{-r}(1+e^{-s})^{-1}(1+e^{-r})^{-1}. The density of the beta conjugate distribution is defined up to a constant of proportionality by

  p(\alpha_\rho, \beta_\rho\,|\,p, q, k) \propto \frac{\Gamma(\alpha_\rho + \beta_\rho)^{k}p^{\alpha_\rho}q^{\beta_\rho}}{\Gamma(\alpha_\rho)^{k}\Gamma(\beta_\rho)^{k}}\,.

  NOTE: This option is only supported for the individual short-term discrepancy terms. Priors may also be specified for the autoregressive parameters for discrepancies modelled using autoregressive processes (i.e. for IndSTPrior and ShaSTPrior objects). These are parametrised via beta distributions such that the autoregressive parameter R \in (-1,1) satisfies

  \frac{R+1}{2} \sim \mathrm{Beta}(\alpha, \beta)

  .

In addition to priors on the discrepancy terms, it is also possible to add prior information on the truth. We require priors on the truth at t=0. By default, a N(0, 10) prior is used on the initial values., however this can be configured by the truth_params argument. The covariance matrix of the random walk of the truth \Lambda_y can be configured using an inverse-Wishart prior. The truth_params argument should be a TruthPrior object.

Value

EnsemblePrior returns an object of class EnsemblePrior. IndSTPrior returns an object of class IndSTPrior. IndLTPrior returns an object of class IndLTPrior. ShaSTPrior returns an object of class ShaSTPrior. TruthPrior returns an object of class TruthPrior.

References

Spence et. al. (2018). A general framework for combining ecosystem models. Fish and Fisheries, 19(6):1031-1042.

Examples

##### Different forms of the individual long term discrepancy priors
#LKJ(10) priors on correlation matrices and gamma(5, 3) priors on the variances
ist_lkj <- IndSTPrior("lkj", list(5, 3), 10)#

#Same as above but with an additional beta(2, 4) prior on
#the autoregressive parameter of the AR process.
ist_lkj <- IndSTPrior("lkj", list(5, 3), 10, AR_params = c(2, 4))

#Same as above but with different variance priors for 5 different variables of interest.
#This encodes that there is a gamma(1, 1) prior on the variance of the first variable,
#a gamma(23, 1) on the second variable etc...
ist_lkj <- IndSTPrior("lkj", list(c(1,23,24,6,87), c(1,1,1,1,5)), 10, AR_params = c(2, 4))

#Hierarchical priors with gamma(1,2) and gamma(10, 1) on the variance hyperparameters and
#gamma(3,4), gamma(5,6) on the correlation hyperparameters
ist_hie <- IndSTPrior("hierarchical", list(1,2,10,1), list(3,4,5,6))

#Hierarchical priors with gamma(1,2) and gamma(10, 1) on the variance hyperparameters and
#the beta conjugate prior with parameters (p = 0.75, q = 0.75, k = 0.2) on the
#correlation hyperparameters
ist_hie_beta_conj <- IndSTPrior("hierarchical_beta_conjugate",
list(1,2,10,1), list(0.75,0.75,0.2))

#Inverse Wishart correlation priors. Gamma(2, 1/3) priors are on the variances and
#inv-Wishart(5, diag(5)) on the correlation matrices.
ist_inW <- IndSTPrior("inv_wishart", list(2, 1/3),list(5, diag(5)))

##### TruthPrior
#Simple default truth prior with 7 variables of interest
truth_def <- TruthPrior(7)
# A more fine-tuned truth prior for an ensemble with 7 species.
truth_cus <- TruthPrior(7, initial_mean = 2, initial_var = 10, rw_covariance = list(10, diag(7)))

#The default priors for an ensemble with 8 variables of interest
priors <- EnsemblePrior(8)

#With 4 variables of interest.
priors <- EnsemblePrior(4)

#Defining custom priors for a model with 4 species.
num_species <- 5
priors <- EnsemblePrior(
  d = num_species,
  ind_st_params = IndSTPrior("lkj",  list(3, 2), 3, AR_params = c(2,4)),
  ind_lt_params = IndLTPrior(
    "beta",
    list(c(10,4,8, 7,6),c(2,3,1, 4,4)),
    list(matrix(5, num_species, num_species),
         matrix(0.5, num_species, num_species))
  ),
  sha_st_params = ShaSTPrior("inv_wishart",list(2, 1/3),list(5, diag(num_species))),
  sha_lt_params = 5,
  truth_params = TruthPrior(d = num_species, initial_mean = 5, initial_var = 10,
                            rw_covariance = list(10, diag(10)))
)

A class to hold the priors for the ensemble model.

Description

An IndLTPrior object encapsulates the prior information for the long-term discrepancies of the individual simulators within the ensemble model.

Details

Individual long-term discrepancies \gamma_k are drawn from a distribtion

\gamma_k \sim N(0, C_\gamma).

Accepted parametrisation forms for this discrepancy are lkj, beta, or inv_wishart. See details of the EnsemblePrior() constructor for more details.

Slots

parametrisation_form

The parametrisation by which the covariance matrix of the noise of the AR process is decomposed.

var_params

The parameters characterising the variance of the distribution of the individual long-term discrepancy.

cor_params

The parameters characterising the correlations of the distribution of the individual long-term discrepancy.

See Also

IndSTPrior, ShaSTPrior, TruthPrior,

A class to hold the priors for the ensemble model.

Description

An IndSTPrior object encapsulates the prior information for the short-term discrepancies of the individual simulators within the ensemble model.

Details

Individual short-term discrepancies z_k^{(t)} are modelled as an AR(1) process so that

z_k^{(t+1)} \sim N(R_k z_k^{(t)}, \Lambda_k).

Accepted parametrisation forms for this discrepancy are lkj, beta, inv_wishart, hierarchical or hierarchical_beta_conjugate. See details of the EnsemblePrior() constructor for more details.

Slots

AR_param

The parameters giving the beta parameters for the autoregressive parameter of the AR(1) process.

parametrisation_form

The parametrisation by which the covariance matrix of the noise of the AR process is decomposed.

var_params

The parameters characterising the variance of the AR process on the individual short-term discrepancy.

cor_params

The parameters characterising the correlations of the AR process on the individual short-term discrepancy. .

See Also

IndLTPrior, ShaSTPrior, TruthPrior,

Backwards Kalman filter

Description

Finds the most likely path through a dynamical linear model.

Usage

KalmanFilter_back(rhos, dee, R, Q, C, P, xhat, Time, y, obs)

Arguments

Details

For the model with the evolution process

x_{t+1}\sim{}N(\rho{}\cdot{}x_{t},Q)

and observation process

y_{t}\sim{}N(\rho{}(x_{t} + \delta),diag(R))

.

Using the sequential Kalman filter, the function gives the mostly path of x_{t} for all t.

Value

A matrix with dimensions nrow(time) and length(xhat) representing the most likely values of the latent variables.

References

Chui, C.K. & Chen, G. (2009) Kalman Filtering with Real-Time Applications. Springer, Berlin, Heidelberg, Fourth Edtion.

Kalman, R. E. (1960) A new approach to linear filtering and prediction problems. Trans. ASME, J. Basic Eng., 82, pp. 35-45.

Examples

fit <- fit_ensemble_model(observations = list(SSB_obs, Sigma_obs),
               simulators = list(list(SSB_ewe, Sigma_ewe, "EwE"),
                                 list(SSB_fs,  Sigma_fs, "FishSUMS"),
                                 list(SSB_lm,  Sigma_lm, "LeMans"),
                                 list(SSB_miz, Sigma_miz, "Mizer")),
               priors = EnsemblePrior(4,
               ind_st_params = IndSTPrior(parametrisation_form = "lkj",
               var_params= list(1,1), cor_params = 10, AR_params = c(2, 2))),
               full_sample = FALSE) #Only optimise in this case
transformed_data <- get_transformed_data(fit)
ex.fit <- fit@point_estimate$par
params <- get_parameters(ex.fit)
ret <- KalmanFilter_back(params$AR_params, params$lt_discrepancies,
                          transformed_data$all_eigenvalues_cov,params$SIGMA,
                          transformed_data$bigM, params$SIGMA_init, params$x_hat,
                          fit@ensemble_data@stan_input$time,transformed_data$new_data,
                          transformed_data$observation_available)

Ecopath with EcoSim SSB

Description

A dataset containing the predictions for spawning stock biomass of Norway Pout, Herring, Cod, and Sole in the North Sea between 1991-2050 under an MSY fishing scenario from EcoPath with EcoSim.

Usage

SSB_ewe

Format

A data.frame with 60 rows and 4 variables:

N.pout

Spawning stock biomass of Norway Pout, in log tonnes.

Herring

Spawning stock biomass of Herring, in log tonnes.

Cod

Spawning stock biomass of Cod, in log tonnes.

Sole

Spawning stock biomass of Sole, in log tonnes.

References

ICES (2016). Working Group on Multispecies Assessment Methods (WGSAM). Technical report, International Council for Exploration of the Seas.

FishSUMS SSB

Description

A ⁠data frame⁠ containing the predictions for spawning stock biomass of Norway Pout, Herring, and Cod in the North Sea between 1984-2050 under an MSY fishing scenario from FishSUMS. Note that FishSUMS does not produce outputs for Sole.

Usage

SSB_fs

Format

A data.frame with 67 rows and 3 variables:

N.pout

Spawning stock biomass of Norway Pout, in log tonnes.

Herring

Spawning stock biomass of Herring, in log tonnes.

Cod

Spawning stock biomass of Cod, in log tonnes.

References

Speirs, D., Greenstreet, S., and Heath, M. (2016). Modelling the effects of fishing on the North Sea fish community size composition. Ecological Modelling, 321, 35–45

LeMans SSB

Description

A data.frame containing the predictions for spawning stock biomass of Norway Pout, Herring, Cod, and Sole in the North Sea between 1986-2050 under an MSY fishing scenario from the LeMaRns package.

Usage

SSB_lm

Format

A data.frame with 65 rows and 4 variables:

N.pout

Spawning stock biomass of Norway Pout, in log tonnes.

Herring

Spawning stock biomass of Herring, in log tonnes.

Cod

Spawning stock biomass of Cod, in log tonnes.

Sole

Spawning stock biomass of Sole, in log tonnes.

References

Thorpe, R. B., Le Quesne, W. J. F., Luxford, F., Collie, J. S., and Jennings, S. (2015). Evaluation and management implications of uncertainty in a multispecies size-structured model of population and community responses to fishing. Methods in Ecology and Evolution, 6(1), 49–58.

mizer SSB

Description

A data.frame containing the predictions for spawning stock biomass of Norway Pout, Herring, Cod, and Sole in the North Sea between 1984-2050 under an MSY fishing scenario from mizer.

Usage

SSB_miz

Format

A data frame with 67 rows and 4 variables:

N.pout

Spawning stock biomass of Norway Pout, in log tonnes.

Herring

Spawning stock biomass of Herring, in log tonnes.

Cod

Spawning stock biomass of Cod, in log tonnes.

Sole

Spawning stock biomass of Sole, in log tonnes.

References

Blanchard, J. L., Andersen, K. H., Scott, F., Hintzen, N. T., Piet, G., and Jennings, S. (2014). Evaluating targets and trade-offs among fisheries and conservation objectives using a multispecies size spectrum model. Journal of Applied Ecology, 51(3), 612–622

Stock assessment SSB

Description

A data.frame containing the single species stock assessment estimates of spawning stock biomass of Norway Pout, Herring, Cod, and Sole in the North Sea between 1984-2017.

Usage

SSB_obs

Format

A data frame with 34 rows and 4 variables:

N.pout

Spawning stock biomass of Norway Pout, in log tonnes.

Herring

Spawning stock biomass of Herring, in log tonnes.

Cod

Spawning stock biomass of Cod, in log tonnes.

Sole

Spawning stock biomass of Sole, in log tonnes.

References

Herring Assessment Working Group for the Area South of 62 N (HAWG)(2020).Technical report, ICES Scientific Reports. ACOM:07. 960 pp, ICES, Copenhagen.

Report of the Working Group on the Assessment of Demersal Stocks in the North Sea and Skagerrak (2020). Technical report, ICES Scientific Reports. ACOM:22.pp, ICES, Copenhagen.

A class to hold the priors for the ensemble model.

Description

An ShaSTPrior object encapsulates the prior information for the short-term discrepancies of the shared discrepancy of the ensemble model.

Details

Shared short-term discrepancies \mathbf{\eta}^{(t)} are modelled as an AR(1) process so that

\mathbf{\eta}^{(t+1)} \sim N(R_\eta \mathbf{\eta}, \Lambda_\eta).

Accepted parametrisation forms for this discrepancy are lkj, beta, or inv_wishart. See details of the EnsemblePrior() constructor for more details.

Slots

AR_param

The parameters giving the beta parameters for the main parameter of the AR(1) process.

parametrisation_form

The parametrisation by which the covariance matrix of the noise of the AR process is decomposed.

var_params

The parameters characterising the variance of the AR process on the shared short-term discrepancy.

cor_params

The parameters characterising the correlations of the AR process on the shared short-term discrepancy.

Ecopath with EcoSim Sigma

Description

A 4x4 covariance matrix quantifying the parameter uncertainty of Ecopath with EcoSim

Usage

Sigma_ewe

Format

A 4x4 matrix.

References

Mackinson, S., Platts, M., Garcia, C., and Lynam, C. (2018). Evaluating the fishery and ecological consequences of the proposed North Sea multi-annual plan. PLOS ONE, 13(1), 1–23.

FishSUMS Sigma

Description

A 3x3 covariance matrix quantifying the parameter uncertainty of FishSUMS

Usage

Sigma_fs

Format

A 3x3 matrix.

References

Spence, M. A., Blanchard, J. L., Rossberg, A. G., Heath, M. R., Heymans, J. J., Mackinson, S., Serpetti, N., Speirs, D. C., Thorpe, R. B., and Blackwell, P. G. (2018). A general framework for combining ecosystem models. Fish and Fisheries, 19(6), 1031–1042.

LeMans Sigma

Description

LeMans Sigma

Usage

Sigma_lm

Format

A 4x4 matrix. A 4x4 covariance matrix quantifying the parameter uncertainty of LeMans

References

Thorpe, R. B., Le Quesne, W. J. F., Luxford, F., Collie, J. S., and Jennings, S. (2015). Evaluation and management implications of uncertainty in a multispecies size-structured model of population and community responses to fishing. Methods in Ecology and Evolution, 6(1), 49–58.

mizer Sigma

Description

mizer Sigma

Usage

Sigma_miz

Format

A 4x4 matrix. A 4x4 covariance matrix quantifying the parameter uncertainty of mizer

References

Spence, M. A., Blackwell, P. G., and Blanchard, J. L. (2016). Parameter uncertainty of a dynamic multispecies size spectrum model. Canadian Journal of Fisheries and Aquatic Sciences, 73(4), 589–59

Stock assessment Sigma

Description

A 4x4 covariance matrix quantifying the covariances of the stock assessment estimates of biomass.

Usage

Sigma_obs

Format

A 4x4 matrix.

References

Herring Assessment Working Group for the Area South of 62 N (HAWG).Technical report, ICES Scientific Reports. ACOM:07. 960 pp, ICES, Copenhagen.

Report of the Working Group on the Assessment of Demersal Stocks inthe North Sea and Skagerrak. Technical report, ICES Scientific Reports. ACOM:22.pp, ICES, Copenhagen.

A class to hold the priors for the truth model in the ensemble framework

Description

An TruthPrior object encapsulates the prior information for the short-term discrepancies of the shared discrepancy of the ensemble model.

Details

The truth \mathbf{y}^{(t)} is modelled as a random walk such that

\mathbf{y}^{(t+1)} \sim N(\mathbf{y}^{(t)}, \Lambda_y).

The covariance matrix \Lambda_y is parameterised by an inverse Wishart distribution (contained in the rw_covariance slot) and the initial value is modelled as drawn from a normal distribution.

Slots

d

A numeric giving the number of variables of interest in the ensemble model.

initial_mean

A numeric giving the standard deviation of the normal prior on the initial mean value of the random walk. This is the same standard deviation for each variable of interest.

initial_var

A list of length 2 containing the shape and scale parameters (respectively) for the gamma priors on the variance of the initial value of the truth.

rw_covariance

A list of length 2 containing the inverse-Wishart parameters for the covariance of the random walk of the truth.

Fits the ensemble model

Description

fit_ensemble_model runs an MCMC of the ensemble model. This process can take a long time depending on the size of the datasets.

Usage

fit_ensemble_model(
  observations,
  simulators,
  priors,
  full_sample = TRUE,
  control = list(adapt_delta = 0.95),
  drivers = FALSE,
  MMod,
  ...
)

Arguments

Value

An EnsembleFit object.

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. https://mc-stan.org

See Also

EnsembleFit, EnsembleSample

Examples

fit <- fit_ensemble_model(observations = list(SSB_obs, Sigma_obs),
               simulators = list(list(SSB_ewe, Sigma_ewe, "EwE"),
                                 list(SSB_fs,  Sigma_fs, "FishSUMS"),
                                 list(SSB_lm,  Sigma_lm, "LeMans"),
                                 list(SSB_miz, Sigma_miz, "Mizer")),
               priors = EnsemblePrior(4,
               ind_st_params = IndSTPrior(parametrisation_form = "lkj",
               var_params= list(1,1), cor_params = 10, AR_params = c(2, 2))),
               full_sample = FALSE) #Only optimise in this case

Generate samples from a fitted ensemble model.

Description

Methods to generates samples of the latent variables from a fitted ensemble model.

Usage

generate_sample(fit, num_samples = 1)

get_transformed_data(fit)

get_parameters(ex.fit, x = 1)

get_mle(x = 1, ex.fit, transformed_data, time)

gen_sample(x = 1, ex.fit, transformed_data, time)

get_transformed_data_dri(fit)

Arguments

Details

The samples are created using the methods described in Strickland et. al. (2009) and Durbin and Koopman (2002).

Value

generate_sample gives a list of length 2, the first element being the MLE of latent variables and the second element being a set of samples of the latent variables.

• If fit is a sampling of the ensemble model parameters, then:

  • mle is a time\times (3M + 2) \times num_samples array, where M is the number of simulators and num_samples is the number of samples from the ensemble model, giving the MLE of the latent variables for each available sample from the ensemble model.

  • sample is a time\times (3M + 2) \times num_samples array, giving a sample of the latent variables for each available sample of the ensemble model.

• If fit is a point estimate of the ensemble model parameters, then:

  • mle is a time\times (3M + 2) \times 1 array giving the MLE of the latent variables for the point estimate of the ensemble model.

  • sample is a time\times (3M + 2) \times num_samples array, giving num_samples samples of the latent variables for the single point estimate of the ensemble model.

get_transformed_data gives a list of transformed input data.

get_parameters gives a list of ensemble parameters from the requested sample.

get_mle If fit is a sampling of the ensemble model parameters, then this is a time\times (3M + 2) \times num_samples array. If fit is a point estimate of the ensemble model parameters, then this is a time\times (3M + 2) \times 1 array giving the MLE of the latent variables for the point estimate of the ensemble model.

gen_sample If fit is a sampling of the ensemble model parameters, then this is a time\times (3M + 2) \times num_samples array, giving a sample of the latent variables for each available sample of the ensemble model. If fit is a point estimate of the ensemble model parameters, then this is a time\times (3M + 2) \times num_samples array.

References

J. Durbin, S. J. Koopman (2002) A simple and efficient simulation smoother for state space time series analysis Biometrika, Volume 89, Issue 3, August 2002, Pages 603–616,

Chris M.Strickland, Ian. W.Turner, Robert Denhamb, Kerrie L.Mengersena. Efficient Bayesian estimation of multivariate state space models Computational Statistics & Data Analysis Volume 53, Issue 12, 1 October 2009, Pages 4116-4125

Examples

fit <- fit_ensemble_model(observations = list(SSB_obs, Sigma_obs),
               simulators = list(list(SSB_ewe, Sigma_ewe, "EwE"),
                                 list(SSB_fs,  Sigma_fs, "FishSUMS"),
                                 list(SSB_lm,  Sigma_lm, "LeMans"),
                                 list(SSB_miz, Sigma_miz, "Mizer")),
               priors = EnsemblePrior(4,
               ind_st_params = IndSTPrior(parametrisation_form = "lkj",
               var_params= list(1,1), cor_params = 10, AR_params = c(2, 2))),
               full_sample = FALSE) #Only optimise in this case
samples <- generate_sample(fit, num_samples = 2000)

Return the compiled ensemble model Stan object.

Description

Gets the unfit, compiled stanmodel object encoding the ensemble model. This allows for manual fitting of the ensemble model directly using rstan::sampling.

Usage

get_mcmc_ensemble_model(priors, likelihood = TRUE, drivers = FALSE)

Arguments

Value

The stanmodel object encoding the ensemble model.

Examples

priors <- EnsemblePrior(4)
mod <- get_mcmc_ensemble_model(priors)

ensemble_data <- EnsembleData(observations = list(SSB_obs, Sigma_obs),
                             simulators = list(list(SSB_ewe, Sigma_ewe, "EwE"),
                                           list(SSB_fs,  Sigma_fs, "FishSUMS"),
                                           list(SSB_lm,  Sigma_lm, "LeMans"),
                                           list(SSB_miz, Sigma_miz, "mizer")),
                              priors = priors)

out <- rstan::sampling(mod, ensemble_data@stan_input, chains = 1)

Plot the ensemble output

Description

Plots the latent variables predicted by the ensemble model, along with simulator outputs and observations.

Usage

## S3 method for class 'EnsembleSample'
plot(x, variable = NULL, quantiles = c(0.05, 0.95), ...)

Arguments

Value

The ggplot object.

Examples

priors <- EnsemblePrior(4)
prior_density <- prior_ensemble_model(priors, M = 4)
samples <- sample_prior(observations = list(SSB_obs, Sigma_obs),
             simulators = list(list(SSB_miz, Sigma_miz),
                               list(SSB_ewe, Sigma_ewe),
                               list(SSB_fs, Sigma_fs),
                               list(SSB_lm, Sigma_lm)),
             priors = priors,
             sam_priors = prior_density)
plot(samples) #Plot the prior predictive density.
plot(samples, variable="Herring")

Generate samples of parameters from prior distribution

Description

Methods to generates samples of the parameters from the prior distribution of the ensemble model.

Usage

prior_ensemble_model(
  priors,
  M = 1,
  MM = NULL,
  full_sample = TRUE,
  control = list(adapt_delta = 0.95),
  ...
)

Arguments

Value

A list containing two items named samples and point_estimate. If full_sample==TRUE, samples is a stanfit and point_estimate is a NULL object, else samples is a NULL and point_estimate is a list object. It is possible to generate a point estimate for the prior if the individual short-term discrepancy prior follows a hierarchical parameterisation.

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. https://mc-stan.org

See Also

EnsembleFit

Examples

priors <- EnsemblePrior(4)
prior_density <- prior_ensemble_model(priors, M = 4)

Generate samples of latent variables from prior predictive distribution

Description

Methods to generates samples of the latent variables from the prior predictive distribution of the ensemble model.

Usage

sample_prior(
  observations,
  simulators,
  priors,
  sam_priors,
  num_samples = 1,
  full_sample = TRUE,
  drivers = FALSE,
  ...
)

Arguments

Details

The samples are created using the methods described in Strickland et. al. (2009) and Durbin and Koopman (2002).

Value

An EnsembleSample object.

References

J. Durbin, S. J. Koopman (2002) A simple and efficient simulation smoother for state space time series analysis Biometrika, Volume 89, Issue 3, August 2002, Pages 603-616,

Chris M.Strickland, Ian. W.Turner, RobertDenhamb, Kerrie L.Mengersena. Efficient Bayesian estimation of multivariate state space models Computational Statistics & Data Analysis Volume 53, Issue 12, 1 October 2009, Pages 4116-4125

See Also

EnsembleFit, EnsembleSample, generate_sample, prior_ensemble_model

Examples

priors <- EnsemblePrior(4)
prior_density <- prior_ensemble_model(priors, M = 4)
samples <- sample_prior(observations = list(SSB_obs, Sigma_obs),
             simulators = list(list(SSB_miz, Sigma_miz),
                               list(SSB_ewe, Sigma_ewe),
                               list(SSB_fs, Sigma_fs),
                               list(SSB_lm, Sigma_lm)),
             priors = priors,
             sam_priors = prior_density)
plot(samples) #Plot the prior predictive density.