Commit 77d934d6 authored by Poulet Camille's avatar Poulet Camille
Browse files

Merge branch 'develop' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D into develop

parents de111f56 38575a86
......@@ -280,8 +280,8 @@
<mortalityRateInOffshore>0.4</mortalityRateInOffshore>
</species.Survive>
<!--<species.WriteEffectiveAndBiomassImportFluxes> <synchronisationMode>ASYNCHRONOUS</synchronisationMode> <exportSeason>SPRING</exportSeason>
<fileNameOutput>effectiveBiomassFluxesBeforeReproduction</fileNameOutput> </species.WriteEffectiveAndBiomassImportFluxes> -->
<!--<analysis.WriteEffectiveAndBiomassImportFluxes> <synchronisationMode>ASYNCHRONOUS</synchronisationMode> <exportSeason>SPRING</exportSeason>
<fileNameOutput>effectiveBiomassFluxesBeforeReproduction</fileNameOutput> </anaysis.WriteEffectiveAndBiomassImportFluxes> -->
<species.ReproduceAndSurviveAfterReproductionWithDiagnose>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
......@@ -328,15 +328,15 @@
</processesEachStep>
<processesAtEnd>
<species.WriteNutrientImportFluxes>
<analysis.WriteNutrientImportFluxes>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<fileName>nutrientImportfFluxes</fileName>
</species.WriteNutrientImportFluxes>
</analysis.WriteNutrientImportFluxes>
<species.WriteNutrientExportFluxes>
<analysis.WriteNutrientExportFluxes>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<fileName>nutrientExportFluxes</fileName>
</species.WriteNutrientExportFluxes>
</analysis.WriteNutrientExportFluxes>
</processesAtEnd>
</processes>
......
......@@ -283,6 +283,19 @@
<!--<species.WriteEffectiveAndBiomassImportFluxes> <synchronisationMode>ASYNCHRONOUS</synchronisationMode> <exportSeason>SPRING</exportSeason>
<fileNameOutput>effectiveBiomassFluxesBeforeReproduction</fileNameOutput> </species.WriteEffectiveAndBiomassImportFluxes> -->
<analysis.AnalyseSpawnerFeatures>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<analysisSeason>SPRING</analysisSeason>
<memorySize>30</memorySize>
</analysis.AnalyseSpawnerFeatures>
<analysis.AnalyseFishDistribution>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<analysisSeason>SPRING</analysisSeason>
<memorySize>30</memorySize>
<minimumRecruitsForPopulatedBasin>50</minimumRecruitsForPopulatedBasin>
</analysis.AnalyseFishDistribution>
<species.ReproduceAndSurviveAfterReproductionWithDiagnose>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<reproductionSeason>SPRING</reproductionSeason>
......@@ -305,6 +318,7 @@
<displayFluxesOnConsole>false</displayFluxesOnConsole>
</species.ReproduceAndSurviveAfterReproductionWithDiagnose>
<species.MigrateFromRiverToInshore>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<migrationSeasonToReachInshore>SPRING</migrationSeasonToReachInshore>
......@@ -328,13 +342,15 @@
</processesEachStep>
<processesAtEnd>
<species.WriteBiomassFluxes>
<analysis.AnalyseLikelihoodOfPresence>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<exportSeason>SPRING</exportSeason>
<fileNameOutput>biomassFluxes</fileNameOutput>
</species.WriteBiomassFluxes>
<presenceFileName>data/input/northeastamerica/nea_presence.csv</presenceFileName>
<period>obs_1900_1950</period>
<minimumRecruitsForPopulatedBasin>50</minimumRecruitsForPopulatedBasin>
<epsilon>0.001</epsilon>
</analysis.AnalyseLikelihoodOfPresence>
<species.IdentifyPopulation>
<analysis.IdentifyPopulation>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<consoleDisplay>false</consoleDisplay>
<years>
......@@ -342,7 +358,7 @@
</years>
<fluxesSeason>SPRING</fluxesSeason>
<fileNameOutput>effectiveFluxes</fileNameOutput>
</species.IdentifyPopulation>
</analysis.IdentifyPopulation>
</processesAtEnd>
</processes>
<useCemetery>false</useCemetery>
......
......@@ -61,7 +61,7 @@ vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sig
}
else {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
increment = exp(rnorm(1, mu, sigma))
increment = exp(rnorm(length(mu), mu, sigma))
}
L = L + increment * tempEffect
return(L)
......@@ -101,7 +101,7 @@ spawnerSurvivalPreReproduction <- function(Triver,Tmin, Topt, Tmax){
#River survival
survProbOptRiver = survivePar$survProbOptRiver
spRiver = survProbOptRiver * TemperatureEffect(Triver, Tmin, Topt, Tmax)
spRiver = survProbOptRiver * temperatureEffect(Triver, Tmin, Topt, Tmax)
#StockAfterSurv = S * SpRiver
......@@ -159,7 +159,7 @@ stockRecruitementRelationship <- function(Triver,Tmin, Topt, Tmax, survivalStock
alleeEffect = 1/(1 + exp(-log(19)*(survivalStock - n/ratioS95_S50*surfaceWatershed)/(n*surfaceWatershed - n/ratioS95_S50*surfaceWatershed)))
Rj = (alphaj * SurvivalStock * alleeEffect)/(betaj + survivalStock * alleeEffect)
Rj = (alphaj * survivalStock * alleeEffect)/(betaj + survivalStock * alleeEffect)
#Rj = ((aj * S) * p)/(Bj +S * p)
......
......@@ -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,36 +270,47 @@ 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()
growthInOffshore <-growthInOffshore %>%
#------------------------------------------#
#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"))
growthInOffshore <- growthInOffshore %>%
mutate(temperature = ifelse(season == "spring" & age == 0, temperature_RI,
#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, 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"),
......@@ -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')
......@@ -587,7 +671,7 @@ growParOptim <- data.frame(value = growth_optimal$par) %>%
rownames_to_column("parameter")
#Rdata
save(growth_optimal, res, res_B, res_C, file = "growth_optimal.RData")
#save(growth_optimal, res, res_B, res_C, file = "growth_optimal.RData")
#csv
write_csv(growParOptim, "growParOptim.csv")
......@@ -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) %>%
......@@ -613,34 +698,31 @@ dataVerif %>% group_by(metapop, age) %>%
geom_line(aes(y = LStich)) +
geom_line(aes(y = L))
save(dataVerif, growth_optimal,res, res_B, res_C, file = "growth_optimal.RData")
# thermal response curve with average seasonal temperature experienced by fish
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")
#load("growth_optimal.RData")
growthCurves = computeGrowOpt(growth_optimal$par, growthInOffshore)
save(growthCurves, file = "growthCurves.RData")
#growth curve according to metapopulations
growthCurves %>%
......@@ -651,7 +733,8 @@ growthCurves %>%
#scale_x_continuous(limits = c(0, 10)) +
geom_hline(yintercept = 45, color = "pink", linetype =1)+ #female
geom_hline(yintercept = 40, color = "lightblue", linetype = 1)+
facet_wrap(.~ metapop, ncol =3)
geom_vline(aes(age), xintercept = 8, linetype = 1)+
facet_wrap(.~ metapop, ncol =1)
#temperature effect
temp <- seq(0,45,0.1)
......@@ -679,41 +762,71 @@ df %>%
```
#Growth Per Gender
```{r growth try, echo = FALSE, warning=FALSE}
#k depending on water temperature in ofsshore basins
kTempEffect = function(kopt, temp, Tmin, Topt, Tmax){
kTemp <- kopt * temperatureEffect(temp,
Tmin,
Topt,
Tmax)
return(kTemp)}
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)
#Generic function
vonBertalanffyGrowthWithTempEffect = function(temp, age, L0, Linf, Kopt, Tmin, Topt, Tmax){
```{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)
koptTemp = kTempEffect(Kopt, temp, Tmin, Topt, Tmax)
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
t0 = log(1 - L0 / Linf) / koptTemp
return(Linf * (1 - exp(-koptTemp * (age - t0))))
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))
```
---
title: "Maturity exploration"
author: "Camille POULET"
date: "08/03/2021"
output: word_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
```
```{r library, include=FALSE}
library(dplyr)
library(tidyr)
library(ggplot2)
library(readr)
library(forcats)
library(knitr)
library(officedown)
library(flextable)
library(stringr)
```
```{r data, include=FALSE}
#load("growthCurves.RData")
load('SOS.rdata')
load("growth_optimal.RData")
source("../GR3D_Rdescription/GR3Dfunction.R")
```
# 1. Estimation of American shad maturity schedules (ASMFC, 2020)
Maturity schedules or ogives were calculated by determining the proportion of mature and immature fish as a function of age for both male and female based on the capture of fish of a species.
When age data were not available, length was used as a surrogate.
In either case, this methods required satisfying the assumption of equal mixing of mature and immature fish and equal recruitment to gear.For most anadromous fish, this condition was not satisfied, and the assumption was then made for 'in-river' captures when all encountered fish were reproductively mature (ASMFC, 2020)
To estimate the age at maturity and identify previous spawning activity, two methods have been used: scales and otholiths.
However, Maki et al. (2001) developped a new approah to estimate maturity probability using scale-derived age in conjonction with the scale-derived patterns of spawning (Maki et al. 2001)
Consistent with this method, the ASMFC used at sea mortality estimates to inform a maturity ogive, by assuming that all fish observed returning to river system were reproductively mature.
## 1.1. Maturity ogives
Otholiths data were used as age informations and scales as a binary variable discriminating first time spawners from multiple spawners.
Following this, a matrix of age-at-capture and age_at first spawn was build for both male and female from each of the 3 regional metapopulations (semelparous, sourthern iteroparous, northern iteroparous)
Assuming a proportional capture of fish across years, and that all fish are susceptible to capture, the proportion of fish mature at age a (p_a) would then calculated as follow: