Georges Bank is a highly productive sub-marine plateau located directly east of Cape Cod, Massachusetts. The bank has a long history of commercial fishing and is currently the focus of an example Fisheries Ecosystem Plan (eFEP) for the New England Fisheries Management Council. This repository will house all of the code that goes into creating a mass balance food web model version of Georges Bank. The model will be built using Rpath, the R implementaion of the mass balance algorithms.
The geographic extent of the model will be the Georges Bank Ecological Production Unit developed by the Ecosystem Assessment and Dynamics Branch (EDAB) of the Northeast Fisheries Science Center (NEFSC).
The base structure of the model is the EMAX model of Link et al. (2006). However, it will be much more disaggregated. Here are the model nodes:
| Group | type |
|---|---|
| Seabirds | Bird |
| Seals | Mammal |
| BalWhale | Whale |
| ToothWhale | Whale |
| HMS | HMS |
| Sharks | Shark |
| AtlHerring | Forage |
| AtlMackerel | Forage |
| Butterfish | Forage |
| SmPelagics | Forage |
| Mesopelagics | Forage |
| OtherPelagics | Forage |
| Cod | DemRound |
| Haddock | DemRound |
| Goosefish | DemRound |
| OffHake | DemRound |
| SilverHake | DemRound |
| RedHake | DemRound |
| WhiteHake | DemRound |
| Redfish | DemRound |
| Pollock | DemRound |
| OceanPout | DemRound |
| BlackSeaBass | DemRound |
| Bluefish | DemRound |
| Scup | DemRound |
| OtherDemersals | DemRound |
| SouthernDemersals | DemRound |
| Fourspot | DemFlat |
| SummerFlounder | DemFlat |
| AmPlaice | DemFlat |
| Windowpane | DemFlat |
| WinterFlounder | DemFlat |
| WitchFlounder | DemFlat |
| YTFlounder | DemFlat |
| OtherFlatfish | DemFlat |
| SmFlatfishes | DemFlat |
| SpinyDogfish | DemRound |
| SmoothDogfish | DemRound |
| Barndoor | Dem |
| WinterSkate | Dem |
| LittleSkate | Dem |
| OtherSkates | Dem |
| Illex | PelInvert |
| Loligo | PelInvert |
| OtherCephalopods | PelInvert |
| AmLobster | BenInvert |
| Macrobenthos | BenInvert |
| Megabenthos | BenInvert |
| AtlScallop | BenInvert |
| Clams | BenInvert |
| OtherShrimps | PelInvert |
| Krill | PelInvert |
| Micronekton | PelInvert |
| GelZooplankton | Zoop |
| Mesozooplankton | Zoop |
| Microzooplankton | Zoop |
| Phytoplankton | Phyto |
Several other groups showed up in the biomass data pull describe below but were not present for the whole time series. Therefore they were omitted.
| Group | Years |
|---|---|
| AmShad | 3 |
| AtlHalibut | 2 |
| NShrimp | 1 |
| RedCrab | 2 |
| RiverHerring | 2 |
| StripedBass | 2 |
| Tilefish | 1 |
| Group | Type |
|---|---|
| Detritus | Detrital Group |
| Discards | Detrital Group |
| DredgeScallop | Fleet |
| DredgeClam | Fleet |
| Gillnet | Fleet |
| Longline | Fleet |
| Seine | Fleet |
| PotTrap | Fleet |
| OttertrawlSmall | Fleet |
| OttertrawlLarge | Fleet |
| Midwater | Fleet |
| OtherFisheries | Fleet |
For non-bivalve species that are encountered in our surveys I am using Bigelow Fall data 2012 to 2016 - without conversion factors. Fall 2017 survey was conducted on the Pisces with no Mid coverage which will be an issue for future work…also spring/fall 2018 compromised due to ship issues.
I am using the scallop survey and clam survey for Atlantic Scallops and Ocean Quhogs/Surfclams respectively. Note that the clam survey samples regions on a rotating basis. See below table:
| Year | ntows |
|---|---|
| 2012 | 0 |
| 2013 | 115 |
| 2014 | 29 |
| 2015 | 8 |
| 2016 | 166 |
Biomass for surveyed species is calculated as the mean from 2013 to 2015. Clams will be an average of 2013 and 2016 as those two years are more representive than the others.
2012 through 2016 are used to see if there is a significant trend for a BA term. Clams were not assessed for a trend because there are only two data points.
Species not encountered by our survey are using EMAX values (Link et al. 2006). I also assumed no BA term.
Three groups showed significant trends. All three were positive which was a surprise. The three groups were Haddock, Southern Demersals, and Atlantic Mackerel. Haddock and Southern Demersal make sense but not Mackerel. An increasing biomass would be counter to the resent assessmets. So in talking with Sarah, I’m going to hold off on BA terms unless model doesn’t balance right.
Need to double check “Freshwater” designation
Pulling data from the stockefficiency tables for 2013 - 2015. Landings will be an average of those years.
Need to deal with skates which are landed as skates(ns). Human consumption mainly winter or thorny while bait mainly little skate (44th SAW). Set all utilcd of 0 to winter skates and utilcd of 7 to little skates. End fate of observer data not known but only 7.6% of skates are listed as skates(ns). The rest are identified to species. Therefore I’m just allowing the “OtherSkates” node absorb the unknown skates.
Dividing otter trawl landings by large and small mesh. The regulation for July 2015 was 6.5 inches. Approximately 18% of otter trawl records (~11% by weight) do not have a mesh associated with them.
Total landings by gear type are divided by the area of the NAFO Statistical Areas that have been selected as Georges Bank. This area is much larger than the survey foot print.
Figuring out which otter trawl trips were small mesh was challenging. Originally I was using the otter trawl gear table (OBOTG). Unfortunately it was not to the haul level and some trips switched gear during the trip. Some even switched between small mesh and large mesh. After extensive searching and frustration I founs the OBDBS table OBOTGH which contains a link to both the haul (OBHAU) and gear tables (OBOTG). The species table (OBSPP) which contains the hail weights for both caught and discarded species contains the primary key for hauls (LINK3). OBOTGH has LINK3 and LINK4 (Gear key) which can be used then to assign mesh size. Not all catches had an associated gear entry but 95.8% did (which is better than the landings tables!). For the purposes of calculating a DK ratio I am not worrying about the unassigned 5%.
Had to manually assign Rpath groups for the incidental takes in the observer database. There were 16 sea turtles in the >500K records. I decided this was not enough to warrent their inclusion. Also, incidental takes are measured as number of individuals not weight. Estimate length is recorded but unfortunately nearly 89% of records do not have on. I did add animal condition to drop incidental takes were the animal was still alive (~11% of records).
I applied a average weight based on the values in Trites and Pauly (1998). I averaged male and female mean weights as I did not know what sex many of the animals were.
First step is assigning prey to Rpath nodes. I am using Jason’s SASprey12B file which has various categories for prey. There are 1456 different prey items including blown and empty. One to one relationships are easy to assign. Then I used the MODCAT to deal with chuncks of prey species at once. Other categories such as ANALCAT, ANALCOM, and Collcom have been useful to assign blocks of species at once.
The largest issue is what to do with prey such as bony fish uncl and hake unclassified. Obviously don’t want to get rid of them or assign to OtherDemersals.
To ensure that I had enough stomach sample for each predator I used the entire time series of food habits.
| Group | Stomachs Full | Stomachs 12-16 |
|---|---|---|
| AmPlaice | 439 | 83 |
| AmShad | 119 | 29 |
| AtlHalibut | 33 | 2 |
| AtlHerring | 1734 | 203 |
| AtlMackerel | 1214 | 195 |
| Barndoor | 1619 | 807 |
| BlackSeaBass | 46 | 24 |
| Bluefish | 417 | 23 |
| Butterfish | 660 | 222 |
| Cod | 6538 | 291 |
| Fourspot | 2222 | 350 |
| Freshwater | 29 | 5 |
| Goosefish | 682 | 140 |
| Haddock | 4446 | 995 |
| Illex | 684 | 0 |
| LargePelagics | 2 | 0 |
| LittleSkate | 10209 | 1481 |
| Loligo | 475 | 0 |
| Mesopelagics | 10 | 0 |
| OceanPout | 906 | 200 |
| OffHake | 24 | 3 |
| OtherDemersals | 9223 | 1182 |
| OtherPelagics | 2 | 0 |
| OtherSkates | 671 | 138 |
| Pollock | 602 | 22 |
| Redfish | 45 | 14 |
| RedHake | 3319 | 457 |
| RiverHerring | 91 | 20 |
| Scup | 174 | 47 |
| SilverHake | 5016 | 617 |
| SmFlatfishes | 302 | 69 |
| SmoothDogfish | 392 | 72 |
| SpinyDogfish | 8118 | 467 |
| StripedBass | 62 | 0 |
| SummerFlounder | 321 | 61 |
| Weakfish | 2 | 0 |
| WhiteHake | 912 | 73 |
| Windowpane | 3050 | 301 |
| WinterFlounder | 2199 | 356 |
| WinterSkate | 9009 | 734 |
| WitchFlounder | 338 | 78 |
| YTFlounder | 2440 | 299 |
Percent weight of prey was calculated using the cluster design explained in Nelson (2014):
\[ \hat{r} = \frac{\sum_{i = 1}^nM_i\hat{\mu_i}}{\sum_{i = 1}^nM_i}\] where \(\hat{r}\) is the mean attribute of interest (in our case weight per stomach), \(M_i\) is the total number of fish in each cluster (Station/Rpath group), and \(\hat{\mu_i}\) the mean attribute in the cluster calculated as:
\[ \hat{\mu_i} = \frac{\sum_{j = 1}^{M_i}{y_{ij}}}{M_i} \] where \(y_{ij}\) is attribute of fish j in cluster i. After calculating the mean weight per stomach of each prey item I converted to proportions as:
\[ \%prey_i = \frac{\hat{r_i}}{\sum_{i=1}^n\hat{r_i}} \]
Needed to add in diet from groups not encoutered in our survey. Will use EMAX as a starting point and proportion aggregate groups to new Rpath groups using the proportion of biomass. Several groups need to be merged from EMAX to Rpath (Macrobenthos and Mesozooplankton). I weighted the DCs by the EMAX groups respective biomasses.
Also noticed that the food habits data from the survey does not include discards or detritus. The EMAX groups did have some nodes feeding on them. I will try and balance without it and see where I get.
Note Calculated the clustering wrong which resulted in fish dominating diets rather than lower trophic levels. Fixed.
While more indepth methods exist (see Aydin et al. 200x), for this model I pulled the numbers from FishBase (citation?). Actually the FishBase values are not complete and I’m not happy with some of them (PB the same for herring and cod!). Instead I am going to use Buchheister et al. (2017). This is a shelf-wide disaggregated Ecopath model based off EMAX. Several groups have multiple stanzas so I will biomass weight their values. Their biological parameters were based on stock assessments.
There were still many groups that were not broken out by Buchheister et al. (2017). I used fishbase.org to pull their parameters where available. For PB I used \(1/Longevitywild\). QB was an average of the QBs on fishbase from different regions.
After mapping all Buchheister groups and individual new groups, SmFlatfish and OtherFlatfish still did not have parameters. I used SmPelagics for SmFlatfish due to their size and position in the food web. I use OtherDemersals for OtherFlatfish. Those parameters were an average of the three demersal fish groups from EMAX and should be agood approximation.
Rediscovered the R scripts for getting information out of Ecobase. I’m going to get ballpark estimates from across other models.
Initial balance saw 45/60 groups with EEs > 1. Going to use PreBal to identify issues.
Coding up the PreBal routines from Link (2010). Some of the checks are easy to make general. Some of them require knowing the classification of the species node (i.e. demersal flatfish or small pelagic). Will hard code for this project and revisit at a later date.
Biomass range - 6 orders of magnitude GOOD
Slope of biomass - -0.70 BAD
Decreased the biomass of many lower trophic levels by an order of magnitude and increased the biomass of everything else by a factor of 4. This gives a slope of -32% which I will accredit to GB being a highly productive but overfished system.
Predators to Prey
| Ratio | Value | Outcome |
|---|---|---|
| Sm pelagics:Zooplankton | 0.69 | GOOD |
| Zooplankton:Phytoplankton | 0.94 | OK |
| Sm Pelagics:Phytoplankton | 0.66 | GOOD |
| Demersals:Benthic Inverts | 0.14 | GOOD |
| Sharks/HMS:Sm Pelagics | 0.002 | GOOD |
| Marine Mammals/Birds:Sm Pelagics | 0.04 | GOOD |
| Whales:Zooplankton | 0.02 | GOOD |
Energy Pathways
| Ratio | Value | Outcome |
|---|---|---|
| Demersal:Pelagic | 3.23 | OK |
| Flatfish:Roundfish | 1.68 | GOOD |
| Sm Pelagics:All fish | 0.24 | GOOD |
| HMS:All fish | 5.3e-4 | GOOD |
| Sharks:All fish | 1.9e-5 | GOOD |
| Demersals:All fish | 0.76 | GOOD |
| Macroinverts:All inverts | 0.87 | OK |
| Gelatinous ZP:All inverts | 0.01 | GOOD |
| Shrimp/micronekton:All inverts | 0.02 | GOOD |
| Other Zooplankton:All inverts | 0.06 | GOOD |
| Zooplankton:Benthos | 0.06 | OK |
| Benthivores:Piscivores | 4.31 | OK |
| Benthivores:Planktivores | 2.63 | OK |
| Planktivores:Piscivores | 1.64 | GOOD |
| <TL3:TL4+ | 6.86 | OK |
| TL | mean.PB | GB.initialPB |
|---|---|---|
| 1-2 | 70.910 | 91.25 |
| 2-3 | 27.470 | 20.78 |
| 3-4 | 1.560 | 1.28 |
| 4+ | 0.605 | 0.27 |
ABC is approximate Bayesian computation - originally developed in genetics to estimate the posterior distribution when calculating the likelihood is impractical.
General reminder of Bayes Theorem:
\[ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} \] where \(p(\theta|D)\) is the posterior distribution, \(p(D|\theta)\) the likelihood, \(p(\theta)\) the prior distribution, and \(p(D)\) the evidence or marginal likelihood.
All ABCs approximate the likelihood through simulations where outcomes are compared to observations. The core is the ABC rejection algorithm where parameters are first sampled from the prior distribution. Given a sampled parameter point \(\hat{\theta}\), a data set \(\hat{D}\) is then simulated under the statistical model M specified by \(\hat{\theta}\). If the generated \(\hat{D}\) is too different from the observed data D, the sampled parameter value is discarded. The distance measure \(\rho(\hat{D}, D)\) determines the level of discrepancy between and D based on a given metric (e.g. Euclidean distance).
\[ \rho(\hat{D}, D) \leq \epsilon \] The probability of generating a \(\hat{D}\) that is within tolerance of D decreases as the dimensionality of the data increases. To avoid this a set of lower-dimensional summary statistics, S(D) are used that captures the informatin in D. The acceptance criteria then becomes:
\[ \rho(\hat{S(D)}, S(D)) \leq \epsilon \]
Package in R for doing approximate Bayesian computation: abc (Csilléry, François, and Blum 2012)
Buchheister, Andre, Thomas J Miller, Edward D Houde, and David A Loewensteiner. 2017. “Technical Documentation of the Northwest Atlantic Continental Shelf (NWACS) Ecosystem Model. Report to the Lenfest Ocean Program, Washington, D.C.” University of Maryland Center for Environmental Sciences Report TS -694 -17.
Csilléry, Katalin, Olivier François, and Michael G. B. Blum. 2012. “Abc: An R Package for Approximate Bayesian Computation (ABC).” Methods in Ecology and Evolution 3 (3): 475–79. https://doi.org/10.1111/j.2041-210X.2011.00179.x.
Link, Jason S. 2010. “Adding Rigor to Ecological Network Models by Evaluating a Set of Pre-Balance Diagnostics: A Plea for PREBAL.” Ecological Modelling 221 (12): 1580–91. https://doi.org/10.1016/j.ecolmodel.2010.03.012.
Link, Jason S., Carolyn A. Griswold, Elizabeth T. Methratta, and Jessie Gunnard. 2006. “Documentation for the Energy Modeling and Analysis eXercise (EMAX).” US Department of Commerce.
Nelson, Gary A. 2014. “Cluster Sampling: A Pervasive, yet Little Recognized Survey Design in Fisheries Research.” Transactions of the American Fisheries Society 143 (4): 926–38. https://doi.org/10.1080/00028487.2014.901252.
Trites, Andrew W, and Daniel Pauly. 1998. “Estimating Mean Body Masses of Marine Mammals from Maximum Body Lengths” 76: 11.