Skip to contents

Whole of ecosystem model

James Thorson

Source: vignettes/too_slow/whole_of_ecosystem.Rmd
whole_of_ecosystem.Rmd

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 an entire ecosystem, using inputs from a previous Rpath model from Whitehouse et al. 2021, and fitting it to biomass and catch time-series for pollock, cod, and arrowtooth.

Eastern Bering Sea

We first load various inputs that are used by Rpath:

data( whitehouse_2021 )
Ecopath = whitehouse_2021$Ecopath
Diet = whitehouse_2021$Diet
stanzas = whitehouse_2021$stanzas
stanza_groups = whitehouse_2021$stanza_groups

We also load time-series of surveys and catches that are used in the MICE model demonstration:

data( eastern_bering_sea )
Catch = eastern_bering_sea$Catch
catch = data.frame( "Year"=Catch$Year, "Mass"=Catch$Mass, "Taxon"=Catch$Taxon )
catch$Taxon = sapply( catch$Taxon, FUN=switch, "Arrowtooth" = "Arrowtooth flounder adult",
                                               "Cod" = "Pacific cod adult",
                                               "Pollock" = "Walleye pollock adult" )

#
Fish = eastern_bering_sea$Survey
biomass = data.frame(
  "Mass" = Fish$Mass,
  "Year" = Fish$Year,
  "Taxon" = sapply( Fish$Taxon, FUN=switch, "Arrowtooth" = "Arrowtooth flounder adult",
                                              "Cod" = "Pacific cod adult",
                                              "Pollock" = "Walleye pollock adult",
                                              NA)
)
biomass = na.omit(biomass)

We then reformat the various Rpath inputs. In particular, we eliminate the pelagic detritus pool (because EcoState currently does not allow more than one detritus functional group), and add benthic_detritus to the Ecopath inputs (because Ecopath solves for equilibrium values for detritus jointly with other groups):

#
which_ecopath = setdiff(1:nrow(Ecopath), 66)
taxa = c( Ecopath[which_ecopath,'Functional.group'], "benthic_detritus" )
type = sapply( taxa, FUN = switch, 
                     "benthic_detritus" = "detritus", 
                     "Large phytoplankton" = "auto", 
                     "Small phytoplankton" = "auto", 
                     "hetero" )
PB = c( Ecopath[which_ecopath,'P.B'], 0.5 )
QB = c( Ecopath[which_ecopath,'Q.B'], NA )
B = c( Ecopath[which_ecopath,'Biomass'], NA )
EE = c( Ecopath[which_ecopath,'EE'], 0.5 )
U = c( rep(0.2,length(which_ecopath)), 0 )

# Modify EE
EE = ifelse( taxa %in% c("benthic_detritus","Walleye pollock adult"), EE, NA )

# 
which_cols = which_ecopath + 1
which_rows = c( which_ecopath, 72 )
DC = cbind( Diet[which_rows, which_cols], "benthic_detritus"=0 )
DC = ifelse( is.na(DC), 0, as.matrix(DC) )
X = array(2, dim=rep(length(taxa),2) )

#
stgroups = stanza_groups[,'StanzaGroup'][stanzas[,'StGroupNum']]
Amax = ceiling((stanzas[,'Last']+1) / 12)

#
K = stanza_groups[,'VBGF_Ksp']
d = stanza_groups[,'VBGF_d']
Wmat = stanza_groups[,'Wmat']
SpawnX = rep(2, length(K))

# Add names
names(PB) = names(QB) = names(B) = names(EE) = names(U) = names(type) = taxa
dimnames(DC) = dimnames(X) = list( taxa, taxa )
names(Amax) = names(stgroups) = taxa[stanzas[,'GroupNum']]
names(K) = names(d) = names(Wmat) = names(SpawnX) = unique(stgroups)

# Eliminate juvenile values
which_juv = stanzas[which(!stanzas[,'Leading']),'GroupNum']
QB[which_juv] = B[which_juv] = NA

We then specify which parameters to estimate:

# Fit catchability coefficient for adult pollock
fit_Q = c( "Walleye pollock adult" )

# fit equilibrium biomass for adult arrowtooth and cod
fit_B = c( "Arrowtooth flounder adult", "Pacific cod adult" )

Finally, we will build inputs without running the model, modify the TMB map manually, and then estimate parameters

run = function( map, nlminb_loops, getsd ){
  out = ecostate(
    taxa = taxa,
    years = 1981:2020,
    catch = catch,
    biomass = biomass,
    PB = PB,
    QB = QB,
    DC = DC,
    B = B,
    EE = EE,
    X = X,
    type = type,
    U = U,
    fit_B = fit_B,
    fit_Q = fit_Q,
    control = ecostate_control( nlminb_loops = nlminb_loops,
                                getsd = getsd,
                                n_steps = 30,
                                profile = NULL,
                                random = NULL,
                                map = map,
                                start_tau = 0.1 ),
    settings = stanza_settings(
      taxa = taxa,
      stanza_groups = stgroups,
      K = K,
      Wmat = Wmat,
      d = d,
      Amax = Amax,
      SpawnX = SpawnX,
      STEPS_PER_YEAR = 12 )
  )
}
out0 = run( map = NULL,
           nlminb_loops = 0,
           getsd = FALSE )
map = out0$tmb_inputs$map
  map$logtau_i = factor(rep(NA,length(map$logtau_i)))
out = run( map = map,
           nlminb_loops = 1,
           getsd = TRUE )

We can compare equilibrium biomass between packages. These are identical for many taxa (because) we specify the same values for biomass for many taxa, and also differ for arrowtooth and cod because we are estimating equilibrium biomass for these two taxa:

# runtime
out$run_time
#> Time difference of 20.10183 mins

# Compile stuff
Borig = c( Ecopath[which_ecopath,'Biomass'], NA )
Bhat = print_ecopars(out, silent=TRUE)[[1]][,'B']
df = data.frame( taxa, 
                 "Whitehouse" = Borig, 
                 "EcoState" = Bhat, 
                 "Difference" = Bhat - Borig )
knitr::kable( df, digits=3 )
taxa Whitehouse EcoState Difference
Transient killer whales 0.000 0.000 0.000
Toothed whales 0.035 0.035 0.000
Gray whales 0.033 0.033 0.000
Baleen whales 0.508 0.508 0.000
Pinnipeds 0.016 0.016 0.000
Walrus and bearded seals 0.114 0.114 0.000
Northern fur seal juvenile 0.001 0.001 0.000
Northern fur seal adult 0.033 0.033 0.000
Ice seals 0.030 0.030 0.000
Other birds 0.001 0.001 0.000
Murres and puffins 0.009 0.009 0.000
Kittiwakes 0.001 0.001 0.000
Auklets 0.002 0.002 0.000
Fulmars 0.000 0.001 0.000
Albatross 0.000 0.000 0.000
Sharks 0.053 0.053 0.000
Walleye pollock juvenile 4.567 4.590 0.023
Walleye pollock adult 18.486 18.486 0.000
Pacific cod juvenile 0.187 0.164 -0.023
Pacific cod adult 2.465 2.183 -0.282
Herring 0.783 0.783 0.000
Arrowtooth flounder juvenile 0.008 0.013 0.004
Arrowtooth flounder adult 0.944 0.420 -0.525
Kamchatka flounder 0.056 0.056 0.000
Greenland turbot juvenile 0.007 0.007 0.000
Greenland turbot adult 0.349 0.349 0.000
Pacific halibut juvenile 0.002 0.002 0.000
Pacific halibut adult 0.222 0.222 0.000
Yellowfin sole 5.258 5.258 0.000
Flathead sole 1.341 1.341 0.000
Northern rock sole 3.383 3.383 0.000
Alaska plaice 1.068 1.068 0.000
Other flatfish 0.262 0.262 0.000
Skates 0.773 0.773 0.000
Sablefish 0.035 0.035 0.000
Eelpouts 2.372 2.372 0.000
Deep demersal fish 0.976 0.976 0.000
Pacific Ocean perch 0.167 0.167 0.000
Other rockfish 0.018 0.018 0.000
Northern rockfish 0.028 0.028 0.000
Shortraker rockfish 0.010 0.010 0.000
Rougheye rockfish 0.004 0.004 0.000
Atka mackerel 0.144 0.144 0.000
Shallow demersal fish 2.306 2.306 0.000
Large sculpins 0.540 0.540 0.000
Octopi 0.192 0.193 0.000
Squids 0.927 0.927 0.000
Salmon returning 0.164 0.164 0.000
Salmon smolts 0.014 0.014 0.000
Myctophids and Bathylagids 0.988 0.988 0.000
Capelin 1.245 1.245 0.000
Sandlance 2.528 2.528 0.000
Other forage fish 2.102 2.102 0.000
Tanner crab (C. bairdi ) 0.522 0.522 0.000
King crabs 0.238 0.238 0.000
Snow crab (C. opilio) 2.179 2.179 0.000
Pandalid shrimp 6.742 6.743 0.000
Benthic zooplankton 35.885 35.885 0.000
Motile epifauna 9.880 9.880 0.000
Structural epifauna 0.781 0.781 0.000
Infauna 87.230 87.230 0.000
Jellyfish 1.112 1.112 0.000
Pelagic zooplankton 2.643 2.643 0.000
Euphausiids 17.831 17.831 0.000
Copepods 26.857 26.857 0.000
Benthic microbes 22.061 22.061 0.000
Large phytoplankton 9.100 9.100 0.000
Small phytoplankton 39.504 39.504 0.000
benthic_detritus NA 12695.821 NA

Runtime for this vignette: 21 mins