Model of Intermediate Complexity"

In ecostate: State-Space Mass-Balance Model for Marine Ecosystems

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
start_time = Sys.time()
# Install locally
#  devtools::install_local( R'(C:\Users\James.Thorson\Desktop\Git\ecostate)', force=TRUE )
#  devtools::install_github( 'James-Thorson-NOAA/ecostate', force=TRUE )
# Build
# setwd(R'(C:\Users\James.Thorson\Desktop\Git\ecostate)'); 
# devtools::build_rmd("vignettes/model_of_intermediate_complexity.Rmd"); rmarkdown::render( "vignettes/model_of_intermediate_complexity.Rmd", rmarkdown::pdf_document())

library(ecostate)

ecostate is an R package for fitting the mass-balance dynamics specified by EcoSim as a state-space model. We here demonstrate how it can be fitted to a real-world data set as a "Model of Intermediate Complexity" while including 10 functional groups and 1 detritus pool, using data across four trophic levels and representing both pelagic and demersal energy pathways.

Eastern Bering Sea

We first load the Survey, Catch, PB, and QB values, and define other biological inputs:

# load data
data(eastern_bering_sea)

# Reformat inputs
years = 1982:2021 # Catch only goes through 2021, and starting pre-data in 1982 doesn't play well with fit_B0
taxa = c( "Pollock", "Cod", "Arrowtooth", "Copepod", "Other_zoop", "Chloro", "NFS", "Krill", "Benthic_invert", "Benthos", "Detritus" )

# Define types
type_i = sapply( taxa, FUN=switch, "Detritus" = "detritus",
                                   "Chloro" = "auto",
                                   "hetero" )

# Starting values
U_i = EE_i = B_i = array( NA, dim=length(taxa), 
                    dimnames=list(names(eastern_bering_sea$P_over_B)))
B_i[c("Cod", "Arrowtooth", "NFS")] = c(1, 0.5, 0.02)
EE_i[] = 1
U_i[] = 0.2

# Define default vulnerability, except for primary producers
X_ij = array( 2, dim=c(length(taxa),length(taxa)) )
dimnames(X_ij) = list(names(B_i),names(B_i))
X_ij[,'Chloro'] = 91

We then fit the function call. This is very slow:

# Define parameters to estimate
fit_Q = c("Pollock", "Copepod", "Chloro", "Other_zoop", "Krill")
fit_B0 = c("Pollock", "Cod", "Arrowtooth", "NFS")
fit_B = c("Cod", "Arrowtooth", "NFS")  

# Define process errors to estimate
# Only estimating Pollock to speed up demonstration
fit_eps = "Pollock"

# Which taxa to include
taxa_to_include = c( "NFS", "Pollock", "Copepod", "Chloro", "Krill" )
# To run full model use:
# taxa_to_include = taxa

# Define priors
log_prior = function(p){
  # Prior on process-error log-SD to stabilize model
  logp = sum(dnorm( p$logtau_i, mean=log(0.2), sd=1, log=TRUE ), na.rm=TRUE)
}

# Run model
out = ecostate( taxa = taxa_to_include,
                years = years,
                catch = eastern_bering_sea$Catch,
                biomass = eastern_bering_sea$Survey,
                PB = eastern_bering_sea$P_over_B,
                QB = eastern_bering_sea$Q_over_B,
                DC = eastern_bering_sea$Diet_proportions,
                B = B_i,
                EE = EE_i,
                U = U_i,
                type = type_i,
                X = X_ij,
                fit_B = fit_B,
                fit_Q = fit_Q,
                fit_eps = fit_eps,
                fit_B0 = fit_B0,
                log_prior = log_prior,
                control = ecostate_control( n_steps = 20, # using 15 by default
                                            start_tau = 0.01,
                                            tmbad.sparse_hessian_compress = 0 ))

# print output
out

We can then plot the estimated foodweb:

# Plot foodweb at equilibrium
# using pelagic producers as x-axis and trophic level as y-axis
plot_foodweb( out$rep$out_initial$Qe_ij,  
              xtracer_i = ifelse(taxa_to_include=="Krill",1,0),
              B_i = out$rep$out_initial$B_i,
              type_i = type_i[taxa_to_include] )

run_time = Sys.time() - start_time

Runtime for this vignette: r paste( round(unclass(run_time),2), attr(run_time, "units") )

ecostate documentation built on April 3, 2025, 5:25 p.m.