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

length_ref_val <- read_xlsx('/Users/domingpa/Documents/Github Repositories/Chinook_energetics/R/output/GAM_reference_values.xlsx')
head(length_ref_val)
## # A tibble: 6 × 3
##   region   age ref_mean
##   <chr>  <dbl>    <dbl>
## 1 CECR       1     612.
## 2 CECR       2     707.
## 3 CECR       3     839.
## 4 CECR       4     918.
## 5 CECR       5     906.
## 6 FRTH       1     622.
gam_predictions  <- read_xlsx('/Users/domingpa/Documents/Github Repositories/Chinook_energetics/R/output/GAM_predictions.xlsx')
head(gam_predictions)
## # A tibble: 6 × 7
##   region   age  year   fit se.fit ref_mean rel_pred
##   <chr>  <dbl> <dbl> <dbl>  <dbl>    <dbl>    <dbl>
## 1 CECR       1  1980  600.  10.1      612.    -1.97
## 2 CECR       1  1981  599.   8.89     612.    -2.15
## 3 CECR       1  1982  598.   7.85     612.    -2.34
## 4 CECR       1  1983  597.   6.99     612.    -2.52
## 5 CECR       1  1984  595.   6.31     612.    -2.71
## 6 CECR       1  1985  594.   5.78     612.    -2.91
colnames(gam_predictions)[2:3] <- c("ocean_age","brood_year")

Size at age and region

region brood_year 1 2 3 4 5
CECR 2008 614.2402 706.8939 833.1406 906.277004042424 NA
CECR 2009 615.5813 706.5605 828.6538 898.200344675202 NA
CECR 2010 616.9442 706.2197 824.0969 890.041501740687 NA
FRTH 1980 567.1555 673.5965 845.8200 NA NA
FRTH 1981 569.8844 674.9115 847.4117 NA NA
FRTH 1982 572.5491 676.2243 849.0288 NA NA
FRTH 1983 575.1021 677.5355 850.6648 NA NA
FRTH 1984 577.4986 678.8474 852.3196 NA NA
FRTH 1985 579.7072 680.1636 853.9999 NA NA
FRTH 1986 581.7178 681.4891 855.7181 NA NA
FRTH 1987 583.5443 682.8299 857.4908 NA NA

In this case, the NAs belonging to the CECR 2008-2010 and FRTH 1980 - 1987 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: 23 × 2
## # Groups:   region [21]
##    region ocean_age
##    <chr>      <dbl>
##  1 CECR           5
##  2 FRTH           4
##  3 FRTH           5
##  4 GRAY           1
##  5 GRAY           5
##  6 GST            5
##  7 HOOD           5
##  8 JNST           5
##  9 JUAN           5
## 10 KLTR           5
## # ℹ 13 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.

But before filling the NAs, let’s see how the data looks so far.

To fill the NAs, I tackled this issue based on two assumptions:

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

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

Additionally, JUAN-age 5 is removed since the size changes were estimated from few samples, resulting in a dramatic size decrease. The age-5 sizes were calculated based nder assumption # 2.

gam_predictions <- gam_predictions%>%
  filter(!(region == 'JUAN' & ocean_age == 5))

full_gam_predictions <- check.grid%>%
  left_join(gam_predictions, by = c("region","brood_year","ocean_age"))%>%
  group_by(region, ocean_age)%>%
  arrange(region,ocean_age)%>%
  fill(fit, se.fit, ref_mean, rel_pred, .direction = 'down')%>% #Fill NAs under assumption 1
  ungroup()%>%
  group_by(region, brood_year)%>%
  arrange(region,brood_year)%>%
  fill(fit, se.fit, ref_mean, rel_pred, .direction = 'down')%>% #Fill NAs under assumption 2
  ungroup()

avg_size_at_year <- full_gam_predictions%>%
  group_by(ocean_age,brood_year)%>%
  summarise(avg.size = mean(fit, na.rm = TRUE))%>%
  ungroup()

Chinook salmon size over brood year - Each line is a region and color indicates ocean age

Chinook salmon size changes divided by ocean age - Color indicates region and black line is the mean size trend for each age

Size changes over brood year for each region - Color indicates ocean age

Add Lipid tier to FRAM abundance database

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.

The following stocks did not have lipid tier and were assigned to the ‘medium’ category:

“Unmarked Nooksack Spr Natural” “Marked Nooksack Spr Hatchery,” and “Unmarked Nooksack Spr Hatchery,” which belong to the NOWA RMIS region.

Additionally, the stocks “Unmarked Mid Oregon Coast Fall” and “Marked Mid Oregon Coast Fall” were also assigned to the ‘medium’ lipid tier and are associated with the NOOR RMIS region.

Calculating Cohorts

The ‘Cohort’ year is necessary since the size estimates 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 (or Run 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).

We could derive the Brood year as:

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

Now the abundance and the size predictions databases are linked by ‘RMIS region’, ‘brood year’ and ‘ocean age’

temp <- stock_abundance%>%
  select(StockID, Age, StartCohort,Shelton.TimeStep,Year.run, Season.run, StockLongName,StockName = FRAM.names.fix, RMIS.Region.fix, Lipid.Ranking.fix,Brood.year)%>%
  left_join(full_gam_predictions, by = c("RMIS.Region.fix"='region','Brood.year'='brood_year','Age'='ocean_age'))

SPS has abundance data but not size estimates. Therefore, it is assumed that the lengths of SPS salmon are similar and follow trends comparable to those from nearby areas, such as HOOD.

check1 <- data.frame(region = c(unique(full_gam_predictions$region)))
check2 <- data.frame(region = c(unique(stock_abundance$RMIS.Region.fix)))

anti_join(check2,check1, by = 'region')
##   region
## 1    SPS
sps_values <- full_gam_predictions %>%
  filter(region == "HOOD")%>%
  mutate(region = "SPS")

stock_abundance <- temp %>%
  left_join(sps_values, by = c('RMIS.Region.fix'='region',"Brood.year"='brood_year', "Age"='ocean_age'), suffix = c("", ".sps")) %>% # Merge SPS's values
  mutate(
    fit = if_else(RMIS.Region.fix == "SPS" & is.na(fit), fit.sps, fit),
    se.fit = if_else(RMIS.Region.fix == "SPS" & is.na(se.fit), se.fit.sps, se.fit),
    ref_mean = if_else(RMIS.Region.fix == "SPS" & is.na(ref_mean), ref_mean.sps, ref_mean),
    rel_pred = if_else(RMIS.Region.fix == "SPS" & is.na(rel_pred), rel_pred.sps, rel_pred)
  ) %>%
  select(-ends_with(".sps")) # Remove temporary columns

rm(temp)

Since the abundance data span up to 2020, the size estimates from 2011 to 2017 (Brood year) are assumed to remain constant after 2010 which is the last year with size estimations. An example of this is shown in the plot below:

stock_abundance_full <-stock_abundance%>%
  arrange(StockLongName,Season.run,Age)%>%
  group_by(StockLongName,Season.run,Age)%>%
  fill(fit, se.fit, ref_mean, rel_pred, .direction = 'down')%>%
  ungroup()%>%
  arrange(StockLongName, Year.run, Shelton.TimeStep, Age)

ggplot(stock_abundance_full[stock_abundance_full$Season.run == '.Spr',], aes(x = Year.run, y = fit, color = as.factor(Age)))+
  geom_line()+
  facet_wrap(~StockLongName, ncol = 4)+
  theme_minimal()

Calculate lipid content

Three lipid indexes were calculated:

lipid_content_f refers to lipids calculated from predicted lengths (fit),

lipid_content_t refers to lipids calculated from the reference mean lengths (ref_mean).

lipid_content_c refers to lipids calculated assuming that the lengths of the fish have not changed relative to the first brood year (contant_length).

The indexes were calculated according to the lipid ranking and the length-lipid relationship parameters found in O’neill et al., 2014:

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}\]

Then, the TOTAL lipid content (predicted, theoretical and constant) was calculated by multiplying the lipid context index by the abundance.

constantlengthref <- stock_abundance_full %>%
  select(StockLongName,Year.run, Age, fit)%>%
  rename(constant_length = fit)%>%
  distinct(StockLongName, Year.run,Age, .keep_all = TRUE)%>%
  group_by(StockLongName,Age)%>%
  filter(Year.run == min(Year.run))%>%
  ungroup()%>%
  select(-Year.run)

stock_abundance_full <- stock_abundance_full%>%
  left_join(constantlengthref, by = c('StockLongName','Age'))%>%
  mutate(lipid_content_f = case_when(Lipid.Ranking.fix == "high" ~ lipid_ranking_params$a[1]*fit^(lipid_ranking_params$b[1]),
                                   Lipid.Ranking.fix == "medium" ~ lipid_ranking_params$a[2]*fit^(lipid_ranking_params$b[2]),
                                   Lipid.Ranking.fix == "low" ~ lipid_ranking_params$a[3]*fit^(lipid_ranking_params$b[3])),
         
         lipid_content_t = case_when(Lipid.Ranking.fix == "high" ~ lipid_ranking_params$a[1]*ref_mean^(lipid_ranking_params$b[1]),
                                   Lipid.Ranking.fix == "medium" ~ lipid_ranking_params$a[2]*ref_mean^(lipid_ranking_params$b[2]),
                                   Lipid.Ranking.fix == "low" ~ lipid_ranking_params$a[3]*ref_mean^(lipid_ranking_params$b[3])),
         
         lipid_content_c = case_when(Lipid.Ranking.fix == "high" ~ lipid_ranking_params$a[1]*constant_length^(lipid_ranking_params$b[1]),
                                   Lipid.Ranking.fix == "medium" ~ lipid_ranking_params$a[2]*constant_length^(lipid_ranking_params$b[2]),
                                   Lipid.Ranking.fix == "low" ~ lipid_ranking_params$a[3]*constant_length^(lipid_ranking_params$b[3])),
         
         total_lipid_content_f = lipid_content_f*StartCohort,
         total_lipid_content_t = lipid_content_t*StartCohort,
         total_lipid_content_c = lipid_content_c*StartCohort
         )

The following plot shows the comparison between the lipid content calculated from the predicted ( solid line ) and the constant lengths ( dashed line ) over time. Some lines appear cut because the size estimates were linked based on the availability of abundance data over time.

The plot shows the Total_lipid_content predicted (Solid line) and constant (Dashed line) over time for the different FRAM stocks. Each stock lipid content includes the Marked and UnMarked salmon abundance.

The plot show the Total_lipid_content predicted (Solid line) and constant (Dashed line) over time for all stocks by season.

### Calculate Biomass

XXXXXX

#The parameters were taken from fishbase. Medium values refers to the general model, while high and low refer to upper and lower credible intervals.
#Bayesian length-weight: a=0.01047 (0.00488 - 0.02245), b=3.04 (2.87 - 3.21), in cm total length, based on LWR estimates for this Genus-body shape 

lw <-data.frame(ranking = factor(c("high","medium","low"), levels = c("high","medium","low")) , a = c(0.02245,0.01047,0.00488), b = c(3.21,3.04,2.87))

stock_abundance_full <- stock_abundance_full%>%
  mutate(fit.cm = fit/10, #We need to transform the lengths to cm
         ref_mean.cm = ref_mean/10,
         constant_length.cm = constant_length/10,
         
         weight_f = case_when(Lipid.Ranking.fix == "high" ~ lw$a[1]*fit.cm^(lw$b[1]),
                                           Lipid.Ranking.fix == "medium" ~ lw$a[2]*fit.cm^(lw$b[2]),
                                           Lipid.Ranking.fix == "low" ~ lw$a[3]*fit.cm^(lw$b[3])),
               
         weight_t = case_when(Lipid.Ranking.fix == "high" ~ lw$a[1]*ref_mean.cm^(lw$b[1]),
                                           Lipid.Ranking.fix == "medium" ~ lw$a[2]*ref_mean.cm^(lw$b[2]),
                                           Lipid.Ranking.fix == "low" ~ lw$a[3]*ref_mean.cm^(lw$b[3])),
         
         weight_c = case_when(Lipid.Ranking.fix == "high" ~ lw$a[1]*constant_length.cm^(lw$b[1]),
                                           Lipid.Ranking.fix == "medium" ~ lw$a[2]*constant_length.cm^(lw$b[2]),
                                           Lipid.Ranking.fix == "low" ~ lw$a[3]*constant_length.cm^(lw$b[3])),
               
               biomass_f.tons = (weight_f*StartCohort)/1000, #Transform grams to tons
               biomass_t.tons = (weight_t*StartCohort)/1000,
               biomass_c.tons = (weight_c*StartCohort)/1000
  )

stock_abundance_full$Lipid.Ranking.fix <- factor(stock_abundance_full$Lipid.Ranking.fix, levels = c("high","medium","low"))
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Total Biomass over time for each FRAM stock (only age 3 - 5) 

prop.sex <- c(0.01,0.5) # Age-specific sex ratios. 0.01 for age 2, and 0.5 for age >3. Taken from Hendrix, 2017. Model Description for the Sacramento River Winter-run Chinook Salmon Life Cycle Model 
prop.mat <- c(0.1,0.9,1) # Maturation rate probs for age 2, 0.9 for age 3, and 1 for age 4 and 5. taken from Hendrix, 2017

#Fecundity alometric params
a <- 0.000935 #taken from Jan's paper for Chinook salmon from Yukon River https://cdnsciencepub.com/doi/full/10.1139/cjfas-2020-0012#sec-4
b <- 2.36 #taken from Jan's paper for Chinook salmon from Yukon River

#Egg mass alometric params
a.mass <- 8.71e-12 #taken from Jan's paper for Chinook salmon from Yukon River 
b.mass <- 4.83 #taken from Jan's paper for Chinook salmon from Yukon River 

temp.fecundity <- stock_abundance_full%>%
  mutate(fecundity_f = case_when(Age == 2 ~ prop.sex[1] * prop.mat[1] * (a*fit^b),
                                 Age == 3 ~ prop.sex[1] * prop.mat[2] * (a*fit^b),
                                 Age >= 4 ~ prop.sex[2] * prop.mat[3] * (a*fit^b)),
         
         fecundity_c = case_when(Age == 2 ~ prop.sex[1] * prop.mat[1] * (a*constant_length^b),
                                 Age == 3 ~ prop.sex[1] * prop.mat[2] * (a*constant_length^b),
                                 Age >= 4 ~ prop.sex[2] * prop.mat[3] * (a*constant_length^b)),
         
         egg_mass_f = case_when(Age == 2 ~ prop.sex[1] * prop.mat[1] * (a.mass*fit^b.mass),
                                 Age == 3 ~ prop.sex[1] * prop.mat[2] * (a.mass*fit^b.mass),
                                 Age >= 4 ~ prop.sex[2] * prop.mat[3] * (a.mass*fit^b.mass)),
         
         egg_mass_c = case_when(Age == 2 ~ prop.sex[1] * prop.mat[1] * (a.mass*constant_length^b.mass),
                                 Age == 3 ~ prop.sex[1] * prop.mat[2] * (a.mass*constant_length^b.mass),
                                 Age >= 4 ~ prop.sex[2] * prop.mat[3] * (a.mass*constant_length^b.mass)),
         
         total_fecundity_f = fecundity_f*StartCohort,
         total_fecundity_c = fecundity_c*StartCohort,
         
         total_egg_mass_f = egg_mass_f*StartCohort,
         total_egg_mass_c = egg_mass_c*StartCohort)

temp.fecundity <- temp.fecundity%>%
  group_by(StockName, Season.run, Year.run, Age)%>%
  summarise(total_fecundity_f = sum(total_fecundity_f),
            total_fecundity_c = sum(total_fecundity_c),
            total_egg_mass_f = sum(total_egg_mass_f),
            total_egg_mass_c = sum(total_egg_mass_c))%>%
  ungroup()
  
ggplot(temp.fecundity[temp.fecundity$Season.run == '.Spr' & temp.fecundity$Age >2,], aes (x = Year.run, color = as.factor(Age)))+
  geom_line(aes( y = total_fecundity_f), linetype = 'solid', size = 0.7)+
  geom_line(aes( y = total_fecundity_c), linetype = 'dashed', size = 0.7)+
  labs(y = 'Fecundity (# eggs)', x = 'Year Run', title = 'Fecundity (age 3 - 5) - Dashed: Constant, Solid: Fitted')+
  theme_minimal()+
  facet_wrap(~StockName, ncol = 3, scales = 'free')

ggplot(temp.fecundity[temp.fecundity$Season.run == '.Spr' & temp.fecundity$Age >2,], aes (x = Year.run, color = as.factor(Age)))+
  geom_line(aes( y = total_egg_mass_f), linetype = 'solid', size = 0.7)+
  geom_line(aes( y = total_egg_mass_c), linetype = 'dashed', size = 0.7)+
  labs(y = 'Total egg mass (g)', x = 'Year Run', title = 'Egg mass (age 3 - 5) - Dashed: Constant, Solid: Fitted')+
  theme_minimal()+
  facet_wrap(~StockName, ncol = 3, scales = 'free')

Linking data to Shelton’s stock regions

Finally, the Shelton stocks are linked to the FRAM stocks. The sum of lipids and is calculated by Shelton Stock/Cohort/Season/Age level.

The following plot shows the Lipid content over time (Only Spring Timestep) - Dashed lines: lipids calculated using constant length, Solid lines: lipids calculated using predicted length”.

The relative change between the fitted and constant lipid content was calculated and shown in the following plots: