Calculating lengths at ‘RMIS region - Brood year - Age’ level

length_ref_val <- read_xlsx("GAM_reference_values.xlsx", sheet = 1)
length_ref_val <- length_ref_val %>% distinct(region,age, .keep_all = TRUE)
head(length_ref_val)
## # A tibble: 6 × 3
##   region   age ref_mean
##   <chr>  <dbl>    <dbl>
## 1 CECR       1     611.
## 2 CECR       2     708.
## 3 CECR       3     841.
## 4 CECR       4     919.
## 5 CECR       5     906.
## 6 FRTH       1     622.
gam_predictions  <- read_xlsx("GAM_predictions.xlsx")
head(gam_predictions)
## # A tibble: 6 × 4
##   RMIS_region ocen_age brood_year rel_pred
##   <chr>          <dbl>      <dbl>    <dbl>
## 1 CECR               1       1977    -2.33
## 2 CECR               1       1978    -2.47
## 3 CECR               1       1979    -2.61
## 4 CECR               1       1980    -2.75
## 5 CECR               1       1981    -2.88
## 6 CECR               1       1982    -3.01
size_at_age <- left_join(gam_predictions,length_ref_val, by = c("RMIS_region" = "region", "ocen_age" = "age"))
size_at_age$pred.avg.length <- (size_at_age$rel_pred*(size_at_age$ref_mean/100)) + size_at_age$ref_mean
head(size_at_age)
## # A tibble: 6 × 6
##   RMIS_region ocen_age brood_year rel_pred ref_mean pred.avg.length
##   <chr>          <dbl>      <dbl>    <dbl>    <dbl>           <dbl>
## 1 CECR               1       1977    -2.33     611.            597.
## 2 CECR               1       1978    -2.47     611.            596.
## 3 CECR               1       1979    -2.61     611.            595.
## 4 CECR               1       1980    -2.75     611.            594.
## 5 CECR               1       1981    -2.88     611.            593.
## 6 CECR               1       1982    -3.01     611.            592.

Size at age and region

RMIS_region brood_year 1 2 3 4 5
CECR 2003 603.1848 706.7292 851.5607 943.8654 938.957672
CECR 2004 605.2615 707.0123 849.2048 937.7984 929.173088
CECR 2005 607.3994 707.2954 846.4283 930.9042 918.120132
CECR 2006 609.5984 707.5785 843.2310 923.4585 906.161196
CECR 2007 611.9195 707.8615 839.7814 915.4612 893.74927
CECR 2008 614.3628 708.1446 835.9952 907.0962 NA
CECR 2009 616.8671 708.4277 832.1248 898.4554 NA
CECR 2010 619.4326 708.7108 828.2545 889.6308 NA
FRTH 1977 562.7028 664.9061 837.6596 976.9968 NA
FRTH 1978 564.8189 666.8125 840.1071 975.3354 NA
FRTH 1979 566.8728 668.7189 842.5546 973.7718 NA

In this case, the NAs belonging to the CECR 2008-2010 and FRTH 1977 - 1979 at age 5 are due to the ‘rel_pred’ and ‘reference size’ value being missing at this age.

Let’s check where the NAs are

print(NA_size_at_age)
## # A tibble: 22 × 2
## # Groups:   RMIS_region [21]
##    RMIS_region ocen_age
##    <chr>          <dbl>
##  1 CECR               5
##  2 FRTH               5
##  3 GRAY               1
##  4 GRAY               5
##  5 GST                5
##  6 HOOD               5
##  7 JNST               5
##  8 JUAN               5
##  9 KLTR               5
## 10 LOCR               5
## # ℹ 12 more rows

All regions have at least one missing value at ‘age-5’, and only the ‘GRAY’ region has missing values at ages 1 and 5. Since the FRAM abundance database does not contain information for age 1, the missing values at ‘age-1’ will not be relevant in the next calculations.

The ‘region-brood year-age-4’ predicted average lengths were used to fill the NAs at ‘age-5’ respectively. However, this was not the best way to do so, as you will see later, therefore I needed to fill the gaps on a different way that will be explained below, thus keep reading :)

We are assuming that that ocean-5s are the same size as ocean-4 fish. The image from the supplement of our 2018 Fish and Fisheries paper shows that while ocean-5s are larger the difference is relatively small and has decreased over time to a point where ocean-4s and ocean-5s are pretty much the same length – these are coast-wide averages and differences will vary by population, but we could use it to justify such an assumption. We would introduce a small bias in those lengths (<4% percent based on where we do have data), but that might not matter much considering the low proportions of ocean-5s (I’d be surprised if any of the populations have more than 1-2% ocean-5s returning) (Jan)

Annual predictions for the size-at-age for each ocean age, Ohlberger2018

full_size_at_age <- size_at_age%>%
  full_join(check.grid,size_at_age, by = c("RMIS_region","brood_year","ocen_age"))%>%
  group_by(RMIS_region,brood_year)%>%
  mutate(pred.avg.length.fixed = if_else(is.na(pred.avg.length) & ocen_age == 5, pred.avg.length[ocen_age == 4], pred.avg.length, missing = pred.avg.length))%>%
   ungroup()

Lipid ranking and Length-Lipid calculation

I manually organized this database in Excel. O’neill et al. 2014, calculated the lipid tier for each stock at the FRAM stock level, and Ohlberger et al. 2018 estimated the salmon stock lengths at the RMIS level. Hence, I matched FRAM stocks to RMIS regions by overlaying the FRAM stock origin river with the RMIS atlas.

lipid_ranking_stock <- read_xlsx("SRKW_prey_kcal.xlsx", sheet = 9) #Tab 'Lipid calls clean'
lipid_ranking_stock <- lipid_ranking_stock%>%distinct(FRAM.long.names,RMIS.Region, .keep_all = TRUE) #Collapsing by FRAM Stock and RMIS regions since the tier does not change among age classes.

lipid_ranking_params <-read_xlsx("SRKW_prey_kcal.xlsx", sheet = 12) #Parameters to calculate the lipid content relative to lipid ranking and length
lipid_ranking_stock[,2:7]%>%head(10)%>%kable(format = "html")%>%kable_styling()
Lipid.calls.names FRAM.names FRAM.long.names RMIS.Region Age Lipid Ranking
Nk/Sm Fall Nooksack/Samish Fall Marked Nooksack/Samish Fall NOWA 3 low
NF Nook Spr NF Nooksack Spr Marked NF Nooksack Spr NOWA 3 medium
SF Nook Spr SF Nooksack Spr Marked SF Nooksack Spr NOWA 3 medium
Skag Su/Fl Fing Skagit Summer/Fall Fing Marked Skagit Summer/Fall Fing SKAG 3 low
Skag Su/Fl Year Skagit Summer/Fall Year Marked Skagit Summer/Fall Year SKAG 3 low
Skag Sprng Year Skagit Spring Year Marked Skagit Spring Year SKAG 3 medium
Snoh Fall Fing Snohomish Fall Fing Marked Snohomish Fall Fing NPS 3 low
Snoh Fall Year Snohomish Fall Year Marked Snohomish Fall Year NPS 3 low
Stil Fall Fing Stillaguamish Fall Fing Marked Stillaguamish Fall Fing NPS 3 low
Tula Fall Fing Tulalip Fall Fing Marked Tulalip Fall Fing NPS 3 low

Then, full_size_at_age is linked to lipid_ranking_stock by RMIS region.

lipids <- left_join(full_size_at_age,lipid_ranking_stock, by = c("RMIS.Region"))
rmis.region age brood_year rel_pred ref_mean pred_length fram.long.names lipid_ranking
CECR 1 1977 -2.33 610.82 596.5879 Marked CR Oregon Hatchery Tule low
CECR 2 1977 -1.17 707.72 699.4397 Marked CR Oregon Hatchery Tule low
CECR 3 1977 1.71 841.38 855.7676 Marked CR Oregon Hatchery Tule low
CECR 4 1977 6.45 919.23 978.5203 Marked CR Oregon Hatchery Tule low
CECR 5 1977 NA NA 978.5203 Marked CR Oregon Hatchery Tule low
CECR 1 1978 -2.47 610.82 595.7327 Marked CR Oregon Hatchery Tule low

Next, the lipid content is calculated according to the lipid ranking, as follows:

For ‘high´lipid ranking \[\text{kcal}^{-1} = 1.8034e^{-05}*(\text{length})^{3.0796}\] For ’medium’ lipid ranking \[\text{kcal}^{-1} = 1.1051e^{-05}*(\text{length})^{3.122}\] For ‘low’ lipid ranking \[\text{kcal}^{-1} = 7.2074e^{-06}*(\text{length})^{3.143}\]

The results are shown in the next table:

rmis.region age brood_year rel_pred ref_mean pred_length fram.long.names lipid_ranking lipid_content
CECR 1 1977 -2.33 610.82 596.5879 Marked CR Oregon Hatchery Tule low 3817.011
CECR 2 1977 -1.17 707.72 699.4397 Marked CR Oregon Hatchery Tule low 6292.579
CECR 3 1977 1.71 841.38 855.7676 Marked CR Oregon Hatchery Tule low 11862.413
CECR 4 1977 6.45 919.23 978.5203 Marked CR Oregon Hatchery Tule low 18077.555
CECR 5 1977 NA NA 978.5203 Marked CR Oregon Hatchery Tule low 18077.555
CECR 1 1978 -2.47 610.82 595.7327 Marked CR Oregon Hatchery Tule low 3799.841

Calculating Cohorts

The ‘Cohort’ year is necessary since the lipid content database needs to be linked to the FRAM abundance database. This way, it is possible to calculate the total lipids per stock/ocean age/year.

The assumption here is that the Cohort year is calculated as:

\[\text{Cohort} = \text{Brood year} + \text{Freshwater residence} + \text{Ocean age}\]

The freshwater age is estimated as follows: fish returning in spring (March-June) are 2 years old in freshwater (FW age-2), whereas fish returning in summer and fall (July-November) are 1 year old in freshwater (FW age-1).

stock_fw_age <- read_xlsx("FRAM stock_season run and fw age.xlsx", sheet = 1)
colnames(stock_fw_age) <- tolower(names(stock_fw_age))

head(stock_fw_age)
## # A tibble: 6 × 3
##   fram.long.names                    season.run freshwater.age
##   <chr>                              <chr>               <dbl>
## 1 Marked Central Valley Fall         Fall                    1
## 2 Marked Columbia R Upriver Bright   Fall                    1
## 3 Marked Columbia R Upriver Summer   Summer                  1
## 4 Marked Cowlitz River Spring        Spring                  2
## 5 Marked CR Bonneville Pool Hatchery Fall                    1
## 6 Marked CR Oregon Hatchery Tule     Fall                    1
lipids <- left_join(lipids,stock_fw_age, by = "fram.long.names")
lipids[,c(1:3,7,9:11)]%>%head()%>%kable(format = "html")%>%kable_styling()
rmis.region age brood_year fram.long.names lipid_content season.run freshwater.age
CECR 1 1977 Marked CR Oregon Hatchery Tule 3817.011 Fall 1
CECR 2 1977 Marked CR Oregon Hatchery Tule 6292.579 Fall 1
CECR 3 1977 Marked CR Oregon Hatchery Tule 11862.413 Fall 1
CECR 4 1977 Marked CR Oregon Hatchery Tule 18077.555 Fall 1
CECR 5 1977 Marked CR Oregon Hatchery Tule 18077.555 Fall 1
CECR 1 1978 Marked CR Oregon Hatchery Tule 3799.841 Fall 1 
lipids <- lipids%>%mutate(cohort = brood_year+freshwater.age+age)

Filling missing values after 2016

stock_abundance <- read_csv("Cohort_Stock_Shelton_seasons.csv")
head(stock_abundance[,c(2,4,8,12,14,20)])
## # A tibble: 6 × 6
##   RunID   Age StartCohort Year.run Year.Season.run StockLongName                
##   <dbl> <dbl>       <dbl>    <dbl> <chr>           <chr>                        
## 1    31     2     107436.     1991 1991.Spr        UnMarked Nooksack/Samish Fall
## 2    31     2      98215.     1991 1991.Sum        UnMarked Nooksack/Samish Fall
## 3    31     2     145806.     1991 1991.Fall       UnMarked Nooksack/Samish Fall
## 4    31     3      66961.     1991 1991.Spr        UnMarked Nooksack/Samish Fall
## 5    31     3      60356.     1991 1991.Sum        UnMarked Nooksack/Samish Fall
## 6    31     3      82951.     1991 1991.Fall       UnMarked Nooksack/Samish Fall

This plot shows the abundance data available for each FRAM stock by age (columns) and cohort (rows). As we can see, there are some gaps of data which were not filled since we assumed no fish were released.

The following figure displays the abundance over time for each FRAM stock in Spring and age (line color)

The ‘lipids’ database only have data for ‘Marked’ stocks, therefore, to ensure accurate abundance and energy content estimates, we assume that the marked and unmarked fish from the same tributaries have the same lipid content.

## # A tibble: 10 × 12
##    rmis.region   age brood_year rel_pred ref_mean pred_length fram.long.names   
##    <chr>       <dbl>      <dbl>    <dbl>    <dbl>       <dbl> <chr>             
##  1 CECR            1       1977    -2.33     611.        597. Marked CR Oregon …
##  2 CECR            2       1977    -1.17     708.        699. Marked CR Oregon …
##  3 CECR            3       1977     1.71     841.        856. Marked CR Oregon …
##  4 CECR            4       1977     6.45     919.        979. Marked CR Oregon …
##  5 CECR            5       1977    NA         NA         979. Marked CR Oregon …
##  6 CECR            1       1977    -2.33     611.        597. UnMarked CR Orego…
##  7 CECR            2       1977    -1.17     708.        699. UnMarked CR Orego…
##  8 CECR            3       1977     1.71     841.        856. UnMarked CR Orego…
##  9 CECR            4       1977     6.45     919.        979. UnMarked CR Orego…
## 10 CECR            5       1977    NA         NA         979. UnMarked CR Orego…
## # ℹ 5 more variables: lipid_ranking <chr>, lipid_content <dbl>,
## #   season.run <chr>, freshwater.age <dbl>, cohort <dbl>

The next plot shows the data available of lipids for each FRAM stock by age (columns) and the calculated cohorts (rows). The ages are lagged throughout time due to the ‘cohort’ calculation in the previous steps.

To fill the blank spaces, it was assumed that the lengths remained the same the last years for each stock/age.

The next plot shows the data available (blue) after filling the gaps from the last years. Notice that we are missing values before 1983, this is not important because the abundance data starts at 1988, thus, these values will be filtered out eventually.

When we plot the lipid content for each FRAM stock by age we got the following: ___

Up to this part, I realized that for age 5 the ‘lipid_content’ suddenly goes up which is odd. This happened due to added lipid content from age 4 to NAs at age 5.

I tackled this issue based on two conditions:

1) If there was data available for age 5, I kept the ‘reference means’ and the ‘rel_pred’ constant for the previous years relative to the last year with information.

2) For the years for age 5 without the ‘reference means’ and “rel_pred”, the values were taken relative to age 4.

Then, both length and lipids were recalculated.

x <- lipids.full%>%
  filter(cohort >= 1988)%>%
  group_by(fram.long.names, age)%>%
  arrange(fram.long.names,age)%>%
  fill(rel_pred,ref_mean,lipid_ranking,rmis.region,freshwater.age, .direction = 'down')%>%
  mutate(pred_length_fixed = (rel_pred*(ref_mean/100)) + ref_mean)%>%#Recalculate the pred_length.
  ungroup()%>%
  group_by(fram.long.names,cohort)%>%
  arrange(fram.long.names)%>%
  mutate(pred_length_fixed2 = vctrs::vec_fill_missing(pred_length_fixed, direction = 'down'), #vctrs::vec_fill_missing does the same as fill but within mutate()
         lipid_content_fixed = case_when(lipid_ranking == "high" ~ lipid_ranking_params$a[1]*pred_length_fixed2^(lipid_ranking_params$b[1]),
                                   lipid_ranking == "medium" ~ lipid_ranking_params$a[2]*pred_length_fixed2^(lipid_ranking_params$b[2]),
                                   lipid_ranking == "low" ~ lipid_ranking_params$a[3]*pred_length_fixed2^(lipid_ranking_params$b[3])))%>%
  select(fram.long.names,rmis.region,cohort,brood_year,age,freshwater.age,rel_pred,ref_mean,pred_length = pred_length_fixed2,lipid_ranking,lipid_content = lipid_content_fixed)%>%
  ungroup()



y <- x%>%
  group_by(fram.long.names)%>%
  arrange(fram.long.names,cohort)%>%
  fill(rel_pred)%>%
  mutate(ref_mean_fill = vctrs::vec_fill_missing(ref_mean, direction = 'down'),
         lipid_content_theoretical = case_when(lipid_ranking == "high" ~ lipid_ranking_params$a[1]*ref_mean_fill^(lipid_ranking_params$b[1]),
                                   lipid_ranking == "medium" ~ lipid_ranking_params$a[2]*ref_mean_fill^(lipid_ranking_params$b[2]),
                                   lipid_ranking == "low" ~ lipid_ranking_params$a[3]*ref_mean_fill^(lipid_ranking_params$b[3])))%>%
  select(fram.long.names,rmis.region,cohort,brood_year,age,freshwater.age,rel_pred,ref_mean = ref_mean_fill,pred_length,lipid_ranking,lipid_content,lipid_content_theoretical)%>%
  ungroup()

lipids.full <-  y
rm(x)

This way the lipid content are constant for the last years, as shown in the next plot:

____

Wrapping up

lipid_content refers to lipids calculated from predicted lengths, which were previously calculated using Ohlberger’s reference values and rates of change.

lipid_content_theoretical refers to lipids calculated from the reference mean lengths.

lipid_content_constant refers to lipids calculated assuming that the lengths of the fish have not changed from 1988.

The following plots show the total lipid values for age 4 (solid line) and 5 (dashed line) under two assumptions: predicted lengths from Jan’s model (blue) and a constant length from 1988 (red).

## Joining with `by = join_by(fram.long.names, rmis.region, age)`

Predict (Blue), and Constant (Red) lipid content over time for age 4 - 5 

lipids.full.all.ages <- lipids.full%>%
  group_by(fram.long.names,cohort)%>%
  summarise(lipid_content = sum(lipid_content),
            lipid_content_constant = sum(lipid_content_constant))%>%
  ungroup()
## `summarise()` has grouped output by 'fram.long.names'. You can override using
## the `.groups` argument.

Predict (Blue), and Constant (Red) lipid content over time for ALL ages

Interestingly, some stocks exhibited a much larger increase in size and consequently in lipids compared to the reference and constant values (e.g. Fraser river stocks)

Calculating the TOTAL energy content by Stock!

Now the lipids.fulldatabase is linked to the FRAM’s stock_abundance database.

lipids_abundance <- left_join(lipids.full, stock_abundance, by = c("fram.long.names" = "StockLongName","cohort"="Year.run", "age"="Age"))

The next heat map shows the data available in the FRAM stock abundance database. In this case, the data was broken down for each stock by ‘Year.Season.run’ (columns) and age (rows).

The total energy content is now calculated by multiplying the energy content of each FRAM stock by the abundance of salmon (i.e. StartCohort variable) in each cohort, age and season.

x <- lipids_abundance%>%
  filter(age != 1)%>% #Remove age 1 because there are not abundance values
  mutate(total_lipid_theoretical = lipid_content_theoretical*StartCohort,
         total_lipid_content = lipid_content*StartCohort,
         total_lipid_constant = lipid_content_constant*StartCohort)%>%
  relocate(total_lipid_theoretical,total_lipid_content,total_lipid_constant, .after = lipid_content_constant)

lipids_abundance <- x
rm(x)

The next plots show the difference between total_lipid_content (i.e. lipids calculated from predicted lengths) and total_lipid_constant (i.e. lipids calculated from constant length) by each FRAM stock for age 4 (red) and 5 only (blue) in the Spring run. Solid line and dashed line depicts the total_lipid_contentand total_lipid_constant respectively.

lipids_total_abundance_spr <- lipids_abundance%>%
  filter(Season.run =='.Spr')%>%
  group_by(cohort)%>%
  summarise(total_lipid_content = sum(total_lipid_content, na.rm = TRUE),
            total_lipid_constant = sum(total_lipid_constant, na.rm = TRUE))

w <- ggplot(lipids_total_abundance_spr, aes(x = cohort))+
  geom_line(aes(y = total_lipid_content/100000000000), color = 'blue')+
  geom_line(aes(y = total_lipid_constant/100000000000), color = 'red')+
  geom_vline(xintercept = 2007)+
  labs(title = 'Lipid content over time - Spring abundance - all stocks - all ages. Blue = Predicted, Red = constant', x = 'Year run', y='Lipid content (x10^11)')+
  #facet_wrap(~Season.run)+
  theme_minimal()

w

Finally, the Shelton stocks are linked to the FRAM stocks. The sum of lipids is calculated by Shelton Stock/Cohort/Season/Age level, and the total energy content is shown in the following plots:

The following plot shows the Lipid content over time (Only Spring season) - Dashed lines: lipids calculated using constant length, Solid lines: lipids calculated using predicted length”, x= “Year run”, y = “Total lipid content”.

Lipid content over time - all ages - Spring season - Dashed lines: lipids calculated using constant length, Solid lines: lipids calculated using predicted length