Commit 57d9216b authored by Poulet Camille's avatar Poulet Camille
Browse files

Comments on Growth Calibration

parent d43fcf3e
......@@ -18,6 +18,7 @@ rm(list = ls())
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggrepel)
library(readr)
library(XML)
library(forcats)
......@@ -46,11 +47,27 @@ source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
```
## Temperature effect experienced by fish for growth
During their life at sea, shads travel together in the marine environment, and population mixing occurs in different places according to season, suggesting that all populations relied/experienced various temperature to growth.
After hatching, juveniles spent their first summer in their natal river before emigrated to coastal marine waters. Juveniles reach their ‘inshore basin’ located at the main entrance of the river when fall begins.
At the end of fall, juveniles move seaward and reach the closest ‘wintering offshore basins’ to spend the winter months.
In early spring, juveniles and sub-adults then move to one of the three ‘summering offshore basin’ to grow according to the migration routes defined by Dadswell et al. (1987) and come back each year until they reach their size at maturity. After spending 4-5 years at sea, ripe individuals started their spawning migration and go back to the ‘inshore basin’ near to their natal river, and enter a ‘river basin’ to spawn.
To account for differences in temperatures experienced by fish to growth, we use the seasonal migration timing describe above.
## Temperature in offshore
For each river basin composing the physical environment of the model, we compute the temperature encountered y fish for each season, according to the connections described above.
```{r Build the connections file with temperature experienced by fish according to their offshore basins, echo = FALSE, warning=FALSE}
#Overview of connections existing between summering and wintering offshore basins
connections %>%
distinct(wintering_offshore_basin_name,summering_offshore_basin_name)
distinct(wintering_offshore_basin_name,summering_offshore_basin_name) %>%
flextable() %>%
autofit()
#Compute an average temperature experienced by juveniles during their life at sea according to connections files i.e. wintering and summering areas.
......@@ -91,22 +108,24 @@ tempEffectInOffshore <- bind_rows(wint,sum) %>%
#mutate(temperatureEffect = temperatureEffect(Tmean, growPar$tempMinGrow, growPar$tempOptGrow, growPar$tempMaxGrow))
#delete intermediate files
remove(sum,wint)
#see the tempertaure range for each season
tempEffectInOffshore %>%
#ggplot(aes(temperature_O, inshore_basin_name), col = season)+
ggplot()+
geom_line(aes(temperature_O, season))
#geom_label(aes(label =season, col = season), size = 3)
```
## Temperature in river and inshore basin
```{r Build the temperature expericenced by fish during their residency in river and inshore basins, echo = FALSE, warning=FALSE}
#TODO
#Contruire un tableau, du même style que celui-ci avec les temperatures du natal basins
After hatching, juveniles spent their first summer in river, and at the end of summer, migrate to the inshore basin to spend the fall.
Then, for each river basin, we pick the corresponding water temperature according to season i.e. the river temperature in summer, and the inshore temperature in fall.
#season = c("spring","summer","fall","winter")
#zone = c("river_tempertaure","inshore_temperature")
```{r Build the temperature expericenced by fish during their residency in river and inshore basins, echo = FALSE, warning=FALSE}
tempEfffectInRiver <- tempInriver %>%
tempEfffectInRI <- tempInriver %>%
pivot_longer(cols = contains('temperature'), names_to = "season", values_to = "river_temperature") %>%
mutate(season = str_remove(season, "_river_temperature")) %>%
inner_join(tempInInshore %>%
......@@ -116,15 +135,15 @@ tempEfffectInRiver <- tempInriver %>%
filter(between(year, 1981, 2010)) %>%
pivot_longer(c("river_temperature", "inshore_temperature"), names_to = "zone", values_to = "temperature")
#Optimzation
tempRiver <- tempEfffectInRiver %>%
#TODO: code a optimiser
tempRiver <- tempEfffectInRI %>%
filter(season %in% c("spring","summer") & zone == "river_temperature")
tempInshore <- tempEfffectInRiver %>%
tempInshore <- tempEfffectInRI %>%
filter(season %in% c("fall","winter") & zone == "inshore_temperature")
tempEffectInFreshWat <- bind_rows(tempRiver,tempInshore) %>%
tempEffectInRI <- bind_rows(tempRiver,tempInshore) %>%
arrange(basin_name,year,season) %>%
select(basin_name,year, season, temperature) %>%
group_by(basin_name,season) %>%
......@@ -138,8 +157,19 @@ tempEffectInFreshWat <- bind_rows(tempRiver,tempInshore) %>%
```
# Growth curves
In the GR3D model, the growth parameters (Table 1) give the following curves when random and temperature effects are not considered.
```{r growth parameters include in the GR3D model, include = TRUE, fig.cap = "Growth parameters include into the GR3D model t"}
growPar %>%
pivot_longer(tempMinGrow:lFirstMaturityForFemale, names_to = "parameters", values_to = "nominal values") %>%
flextable() %>%
autofit()
```
The growth parameters used in GR3D give the following curves when random and temperature effects are not considered.
```{r, include = TRUE, fig.cap = "Growth curve without temperature effect"}
......@@ -158,9 +188,37 @@ dfGrowth %>% filter(basin_name == 'reference') %>%
```
To optimize the growth parameters included in GR3D for the Amercian shad, we used values gave by Stich et al, 2020.
## Optimization of growth parameters
To optimize the growth parameters included in GR3D for the American shad, we used growth parameters values provided by Stich et al, 2020 (ASFMC, 2020)(Table 2)
```{r, fig.cap = "Growth parameters from Stich et al (2020) and reported in the last ASFCM report on the Amercian shad (2020)"}
Stich2020_sel <- read_csv("../NEA_calibration_offline/Stich_Table 9.csv") %>%
filter(parameter == "K" | parameter == "Linf"|parameter == "t0") %>%
dplyr::rename("catchment"="cachtment") %>%
select(parameter,catchment,mean) %>%
pivot_wider(names_from = parameter, values_from = mean) %>%
Stich2020_sel %>%
flextable() %>%
autofit()
```
Stich et al (2020) used the parameter t0, instead of the length at hatching (L0), as in commonly used.
To optimization on growth, required having values for length at hatching.
Then, we computed the length at hatching derived from Stich et al (2020) based on the growth coefficient, K, and the age at hatching available using the following equation:
$$
L_0 = L_{inf} . (1-{\exp^{t_0.k})}
$$
where L_0 is the length at hatching, t0 is the age at hatching, and k is the growth coefficient. t_0 and k wer both derived from Stich et al (2020).
```{r Growth from Stich et al, 2020, echo = FALSE, warning=FALSE}
```{r Growth parameters from Stich et al, 2020, echo = FALSE, warning=FALSE}
#------------------------------------------#
#Compute L0 for Stich et al (2020)
#-----------------------------------------#
#function to compute L0 based on growth paramters from Stich et al.
#K, Linf and t0 were provided for the three metapopulations.
......@@ -173,16 +231,30 @@ vonBertalanffyLengthAtHatch = function(Linf, K, t0){
return(L0)
}
#vonBertalanffyTry(454,0.49,-0.21)
#parameters from Stich
Stich2020_sel <- read_csv("../NEA_calibration_offline/Stich_Table 9.csv") %>%
filter(parameter == "K" | parameter == "Linf"|parameter == "t0") %>%
dplyr::rename("catchment"="cachtment") %>%
select(parameter,catchment,mean) %>%
pivot_wider(names_from = parameter, values_from = mean) %>%
#Results
Stich2020_sel<-Stich2020_sel %>%
mutate(L0_theo = vonBertalanffyLengthAtHatch(Linf, K, t0))
#create a data frame and compute L at age with values from Stich
#Table
Stich2020_sel %>%
flextable() %>%
autofit()
```
```{r Growth curves from Stich et al, 2020, echo = FALSE, warning=FALSE}
#------------------------------------------#
#Compute growth curves from Stich et al (2020)
#------------------------------------------#
#create a data frame (growthInOffshore) with 4 variables:
#basin_name = name of the river_basin included in GR3D
#age: vector ranging from 0 to 10
#season: 4 values, spring, summer, fall, winter
#latitude : lat_outlet of the river basin, used to define the 3 metapopulations.
#L = compute L at age with values from Stich
growthInOffshore <- expand_grid(basin_name = nea_riverBasinFeatures$basin_name , age = seq(0,9.75,.25)) %>%
mutate(season = case_when (
......@@ -190,7 +262,7 @@ growthInOffshore <- expand_grid(basin_name = nea_riverBasinFeatures$basin_name ,
(age - floor(age)) == 0.25 ~ 'summer',
(age - floor(age)) == 0.50 ~ 'fall',
(age - floor(age)) == 0.75 ~ 'winter'
)) %>% inner_join(nea_riverBasinFeatures %>% select(basin_name, lat_outlet), by ='basin_name') %>%
)) %>% inner_join(nea_riverBasinFeatures %>% select(basin_name, lat_outlet), by ='basin_name') %>% #add latitude to identify the 3 regions
select(basin_name, lat_outlet, age, season)%>%
mutate(metapop = case_when
(lat_outlet < 33.8 ~ "semelparous",
......@@ -198,40 +270,51 @@ growthInOffshore <- expand_grid(basin_name = nea_riverBasinFeatures$basin_name ,
lat_outlet > 41.28793 ~ "northern iteroparous")) %>%
inner_join(Stich2020_sel, by = c("metapop" = "catchment")) %>%
group_by(basin_name,metapop) %>%
#compute the length at each age based on the growth parameters from Stich et al (2020) provided above.
mutate(LStich = vonBertalanffyGrowth(age, L0_theo, Linf, K) / 10) %>% # in cm
ungroup()
#------------------------------------------#
#Add temperature effect for growth curves (2020)
#------------------------------------------#
#Add the temperature time-series for offshore, river, and inshore basins
growthInBasin <-growthInOffshore %>%
inner_join(tempEffectInOffshore,
by = c("basin_name" = "inshore_basin_name", "season")) %>%
right_join(tempEffectInFreshWat,
right_join(tempEffectInRI,
by = c("basin_name","season"))
#for each season, pick the temperature depending on the location of fish
growthInBasin <- growthInBasin %>%
mutate(temperature_RIO = ifelse(season == "spring" & age == 0, temperature_RI,
ifelse(season == "summer" & age == 0.25, temperature_RI,
ifelse(season == "fall" & age == 0.50, temperature_RI, temperature_O)))) %>%
arrange(basin_name,age,season) %>%
select(basin_name, offshore_basin_name, metapop,age, season, temperature_O, temperature_RIO, LStich)
select(basin_name, offshore_basin_name, metapop,age, season, temperature_O, temperature_RI, temperature_RIO, LStich)
#group_by(metapop,age,season) %>%
#summarize(LStich = unique(LStich))
```
Stich et al (2020) work on 11 rivers, and 1 river ( North (MA)) is not referenced in our database.
Optimization was run only for rivers in which presences were recorded, and for which growth parameters were available.
```{r join between Stich and GR3D database}
# NEW =====
# add presence and reorder level for metapop and season
growthInBasin <- growthInBasin %>%
inner_join(nea_presence %>% select(-basin_id), by ='basin_name') %>%
mutate(metapop = factor(metapop, levels = c('northern iteroparous', 'southern iteroparous', 'semelparous')),
season = factor(season, levels = c('winter', 'spring', 'summer', 'fall')))
```
There is 1 river ( North (MA)) in Stich et al paper that is not referenced in our database.
```{r join between Stich and GR3D database}
# NEW =====
#rename the basin_name to fit with our list of river_basin name
Stich2020_sel <- Stich2020_sel %>% mutate(basin_name = catchment) %>%
mutate(basin_name = replace(basin_name, basin_name == "Tar-Pamlico", "Pamlico-Tar"),
basin_name = replace(basin_name, basin_name == "Albemarle", "Roanoke"),
basin_name = replace(basin_name, basin_name == "Upper Chesapeake Bay", "Susquehanna"))
basin_name = replace(basin_name, basin_name == "Upper Chesapeake Bay", "Susquehanna"))
Stich2020_sel %>% filter(!catchment %in% c('coastwide', 'northern iteroparous', 'semelparous', 'southern iteroparous' )) %>%
select(catchment, basin_name) %>%
......@@ -241,10 +324,11 @@ Stich2020_sel %>% filter(!catchment %in% c('coastwide', 'northern iteroparous',
autofit()
```
Temperature ranged between -0.89 and 28.15°C. If we considered only basins in which presences are recorded, the thermal range is narrow, ranging from 1.66 to 27.91.
Recorded presence of shad inside these thermal ranges implicitly suggested that temperature between 1.66 and 27.91 are suitable and preferentially used by shads.
As such, we could use this thermal limits as starting values for growth optimisation.
```{r range of temperature}
```{r range of temperature, include = TRUE}
# NEW =====
growthInBasin %>%
group_by(obs_1900_1950) %>%
......@@ -264,7 +348,7 @@ growthInBasin %>%
#define parameters for both sex based on value from the XML file
growParModifed <- growPar %>%
growParUnisex <- growPar %>%
transmute(tempMinGrow,tempMaxGrow,tempOptGrow, sigmaDeltaLVonBert,lengthAtHatching,
kOpt = mean(c(kOptForFemale,kOptForMale)),
Linf = mean(c(linfVonBertForFemale,linfVonBertForMale)))
......@@ -280,7 +364,7 @@ growParModifed <- growPar %>%
```{r}
# new ====
save(growParModifed,growthInBasin, nea_presence, Stich2020_sel, file = 'SOS.rdata' )
save(growParUnisex,growthInBasin, nea_presence, Stich2020_sel, file = 'SOS.rdata' )
rm(list = ls())
load('SOS.rdata')
......@@ -599,11 +683,12 @@ write_csv(growParOptim, "growParOptim.csv")
# NEW
# ===============================================================================
#parGrowthUsed = c(res_B$par) %>% enframe() %>% pivot_wider()
parGrowthUsed = c(res_C$par, fixedPar_C) %>% enframe() %>% pivot_wider()
growParOptim = c(res_C$par, fixedPar_C) %>% enframe() %>% pivot_wider()
# ================================================================================
dataVerif <- computeGrowAllBasins(data = dataCalibration,
growPar = parGrowthUsed %>% mutate(sigmaDeltaLVonBert = 0) )
growPar = growParOptim %>%
mutate(sigmaDeltaLVonBert = 0))
# graph by meta pop
dataVerif %>% group_by(metapop, age) %>%
......@@ -619,22 +704,19 @@ save(dataVerif, growth_optimal,res, res_B, res_C, file = "growth_optimal.RData")
dataTemperatureMetapop <- dataVerif %>%
group_by(metapop, season) %>%
summarise(temperature = mean(temperature), .groups = 'drop') %>%
mutate(tempEffect = temperatureEffect(temperature, parGrowthUsed$tempMinGrow, parGrowthUsed$tempOptGrow, parGrowthUsed$tempMaxGrow ),
mutate(tempEffect = temperatureEffect(temperature, growParOptim$tempMinGrow, growParOptim$tempOptGrow, growParOptim$tempMaxGrow ),
metapop = factor(metapop, levels = c('northern iteroparous', 'southern iteroparous', 'semelparous')),
season = factor(season, levels = c('winter', 'spring', 'summer', 'fall')))
data.frame(temperature = seq(0,40,0.1)) %>%
mutate(tempEffect = temperatureEffect(temperature, parGrowthUsed$tempMinGrow, parGrowthUsed$tempOptGrow, parGrowthUsed$tempMaxGrow )) %>%
mutate(tempEffect = temperatureEffect(temperature, growParOptim$tempMinGrow, growParOptim$tempOptGrow, growParOptim$tempMaxGrow )) %>%
ggplot() + ylim(0,1)+
geom_line(aes(x = temperature, y = tempEffect)) +
geom_point(data = dataTemperatureMetapop, aes(x = temperature, y = tempEffect, color = metapop, shape = season ), size=3)
```
```{r draw growth curves}
#load("growth_optimal.RData")
......@@ -680,17 +762,71 @@ df %>%
```
#Growth Per Gender
Growth curves were defined for sex combined according to Stich et al (2020). To account variations in growth coefficients and size at maturity between males and females, we need distinct growth curve for each sex.
Then, we look for the kOpt for males and females so that both curves corresponded on average to the growth curve (sex-combined) obtained through optimisation for each metapopulation.
According to the literature, males and females have differences of 0.5 in age at maturity (TODO: a vérifier)
```{r compute growth per gender}
objFUN = function(par, ages, parGrow){
meanCurve = vonBertalanffyGrowth(ages, parGrow$L0, parGrow$Linf, parGrow$K)
femaleCurve = vonBertalanffyGrowth(ages, parGrow$L0, parGrow$Linf, par['Kfemale'])
maleCurve = vonBertalanffyGrowth(ages, parGrow$L0, parGrow$Linf, par['Kmale'])
SCE1 = sum(((maleCurve + femaleCurve)/2 - meanCurve)^2) / length(ages)
femaleAgeAtMaturity = vonBertalanffyInverse(parGrow$LatMaturity, parGrow$L0, parGrow$Linf, par['Kfemale'])
maleAgeAtMaturity = vonBertalanffyInverse(parGrow$LatMaturity, parGrow$L0, parGrow$Linf, par['Kmale'])
SCE2 = ((femaleAgeAtMaturity - maleAgeAtMaturity) - parGrow$diffMeanAgeAtMaturity)^2
return(unname(SCE1 + SCE2))
}
parGrow = data.frame(Linf = 45.4,
K = 0.49,
L0 = 4.4,
LatMaturity = 40 ,
diffMeanAgeAtMaturity = 1.5)
parGrow = parGrowthUsed %>%
select(lengthAtHatching,Linf,kOpt) %>%
rename(c("L0" = "lengthAtHatching",
"K" ="kOpt")) %>%
mutate(Lmaturity = 40,
diffMeanAgeAtMaturity = 0.5)
ages = seq_len(12)
objFUN(par = c(Kfemale = 0.45, Kmale = .52),
ages = ages,
parGrow = parGrow)
res = optim(par = c(Kfemale = 0.45, Kmale = .52),
fn = objFUN,
ages = ages,
parGrow = parGrow)
data.frame(age = seq(0,15, .1)) %>% mutate(Lmale = vonBertalanffyGrowth(age, parGrow$L0, parGrow$Linf, res$par['Kmale']),
Lfemale = vonBertalanffyGrowth(age, parGrow$L0, parGrow$Linf, res$par['Kfemale']),
Lboth = vonBertalanffyGrowth(age, parGrow$L0, parGrow$Linf, parGrow$K)
) %>%
pivot_longer(cols = starts_with('L'), names_to = 'gender', values_to = 'L') %>%
mutate(gender = str_replace(gender, 'L', '')) %>%
ggplot(aes(x = age, y = L, color = gender)) +
geom_line() +
geom_abline(slope = 0, intercept = parGrow$LatMaturity)
data.frame(age = seq(0,15, .1)) %>% mutate(LStich = vonBertalanffyGrowth(age, L0 = 4.44, Linf = 45.4 , K = 0.49)) %>%
ggplot(aes( x = age, y = LStich)) +
geom_path() +
geom_abline(slope = 0, intercept = c(40.3, 42))
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment