Commit 1b6f6d7a authored by Poulet Camille's avatar Poulet Camille
Browse files

Clean Rmd maturation

No related merge requests found
Showing with 290 additions and 394 deletions
+290 -394
--- ---
date: "`r Sys.Date()`" date: "22/04/2021"
author: "Your Name" author: "Camille POULET & Patrick LAMBERT"
title: "officedown template" title: "officedown template"
output: output:
officedown::rdocx_document: officedown::rdocx_document:
mapstyles: mapstyles:
Normal: ['First Paragraph'] Normal: ['First Paragraph']
...@@ -12,13 +12,15 @@ output: ...@@ -12,13 +12,15 @@ output:
knitr::opts_chunk$set(echo = TRUE, fig.cap = TRUE) knitr::opts_chunk$set(echo = TRUE, fig.cap = TRUE)
library(officedown) library(officedown)
library(flextable)
library(officer) library(officer)
library(readr)
library(tidyverse) library(tidyverse)
library(tictoc) # see timeR library(tictoc) # see timeR
library(data.table) library(data.table)
library(readr) library(ggplot2)
library(ggrepel)
library(tibble)
``` ```
...@@ -26,305 +28,150 @@ library(readr) ...@@ -26,305 +28,150 @@ library(readr)
rm(list = ls()) rm(list = ls())
source("../GR3D_Rdescription/GR3Dfunction.R") source("../GR3D_Rdescription/GR3Dfunction.R")
source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R") source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
source("../GR3D_Rdescription/GR3D_NEA_env_data.R")
#load("growth_optimal.RData") #load("growth_optimal.RData")
load('SOS.rdata') load('SOS.rdata')
load('res_optim.rdata') load('res_optim.rdata')
``` ```
```{r load growth parameters and tempertaure pattern }
regional_metapop = c("semelparous", "southern iteroparous", "northern iteroparous") # Maturation
growParXML = growPar ##1.Observation of American shad maturity schedules (ASMFC, 2020)
#Growth parameters from Stich et al, 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.
#Temperature relative values from growth optimization 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)
growParStich <- Stich2020_sel %>% To estimate the age at maturity and identify previous spawning activity, two methods have been used: scales and otholiths.
filter(catchment %in% regional_metapop) %>% 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)
select(catchment, Linf, K, L0_theo) %>% 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.
rename("metapop" = "catchment") %>%
mutate(tempMinGrow = 1.6, tempOptGrow = 5, tempMaxGrow = 27.9,
Linf = Linf/10, lengthAtHatching = L0_theo/10, kOpt = K,
sigmaDeltaLVonBert = 0.4)
temperaturePattern <- growthInBasin %>% ## 1.1. Maturity ogives
group_by(metapop, age, season) %>%
summarise(temperature = mean(temperature_RIO), .groups = 'drop') %>%
#filter(metapop == "semelparous") %>%
# mutate(season = case_when(
# (age - floor(age)) == 0.00 ~ 'spring',
# (age - floor(age)) == 0.25 ~ 'summer',
# (age - floor(age)) == 0.50 ~ 'fall',
# (age - floor(age)) == 0.75 ~ 'winter'
# )) %>%
select(metapop,age, season,temperature)
``` 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)
```{r test } 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:
#Do not run
tic() $$
computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growParStich) %>% p_a = \frac{observed_a}{(observed_a + immature_a)}
filter(season == 'spring') $$
toc()
tic() where p_a is the proportion of fish mature at age a, observed_a are mature fish observed at age a, with a assumed to be between years 1 and 9, and immature a are the immature fish observed each age.
computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 100000, growParStich) %>% As it was considered that returning fish were exclusively mature, there is no observations of immature fish at age.
filter(season == 'spring') Then, the proportion of immature fish was estimated based on all the returns greater than age a and the estimated mortality (M) at sea as below:
toc()
computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) , $$
Nind = 10000, growPar) %>% immature_a = \sum_{j=a+1}^{9}\frac{(observed_j)}{\exp^{-M(j-a)}}
filter(season == 'spring') %>% $$
computeOgive(lengthAtMaturity = 35) ## 1.2.Natural mortality
``` The natural mortality at sea (M) was calculated based on maximum ages following the method described in Then et al. (2015) for each region. Estimates of natural mortality by region and sex were provided in Table 10, and Table 27 from the ASMFC.
```{r echo = FALSE, fig.cap= "Table 1. Natural mortality estimates by region and sex. Maximum age data for the northern iteroparous (NI) region were borrowed from the sourthern iteroparous (SI) region due to more prevalent and long term impacts of man-made barriers in the NI region", out.width='100%'}
```{r ogive variability} knitr::include_graphics("American_Shad_natural_mortality.PNG")
#Do not run
tic(msg = 'ogive variability')
replicate(10, computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
Nind = 100000, growPar) %>%
filter(season == 'spring') %>%
computeOgive(lengthAtMaturity = 35) %>% select(p) %>% unlist(use.names = FALSE))
toc()
``` ```
```{r growth curve envelope} ## 1.3 Results
# envelope of growth curves
computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growParStich) %>% By considering both sex and the mortality correction approach, female have a later maturity-at-age than males in each region, as commonly observed in other shad groups(Taverny, 1991).
group_by(age) %>% Both sexes begin to mature at 3 years,with males maturing on average in 4.64 years and females maturing on average in 5,28 years (ASMFC, 2020).
summarise(Lmin = min(L), For a deeper explanation on how the mean age for male and female was computed, see the xls file "MaturityOgives".
L025 = quantile(L, .25),
L25 = quantile(L, .25),
Lmed = median(L), ```{r American shad maturity at-age- observed by region and sex using the mortality correction approach (ASMFC, 2020), include = TRUE }
L75 = quantile(L, .75),
L975 = quantile(L, .975), #American shad maturity at-age- observed by region and sex using the mortality correction approach (ASMFC, 2020)
Lmax = max(L)) %>% ogivesObs = expand_grid(age = seq.int(9), metapop = regional_metapop) %>%
ggplot(aes(x = age)) + arrange(metapop, age) %>%
geom_ribbon(aes(ymin = Lmin, ymax = Lmax), fill = "gray30", alpha = .2) + group_by(metapop) %>%
geom_ribbon(aes(ymin = L025, ymax = L975), fill = "gray30", alpha = .2) + mutate(p_female = fifelse(
geom_ribbon(aes(ymin = L25, ymax = L75), fill = "gray30", alpha = .2) + metapop == "semelparous", c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1), fifelse (metapop == "southern iteroparous", c(0, 0, 0,0.04,0.27,0.64,0.81,0.9,1), c(0,0,0,0.04,0.69,0.69,0.9,1,1)))) %>%
geom_line(aes(y = Lmed)) + mutate(p_male = fifelse(
labs(x = 'Age (y)', y = 'Length (cm)') + metapop == "semelparous", c(0, 0, 0.07, 0.29, 0.60, 0.86, 1, 1, 1), fifelse (metapop == "southern iteroparous", c(0,0,0,0.12,0.5,0.79,0.8,0.84,1), c(0,0,0.03,0.21,0.82,1,1,1,1)))) %>%
geom_abline(slope = 0, intercept = c(40.3, 42)) ungroup()
# verify that all L < Linf
computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growParStich) %>%
filter(L>growParStich$Linf)
``` ```
```{r optimisation - OBSOLETE} #Growth
#obsolete (presque) au k ou! #1. Observation of growth curves (ASMFC, 2020)
ogiveObs = data.frame(age = seq.int(9), p = c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1))
lengthTrajectoriesInSpring = computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growPar) %>%
filter(season == 'spring')
# calibrate only lengthAtMaturity
objFn_A = function(par, ogiveObs, yearlyLengthTrajectories){
SSE = ogiveObs %>%
inner_join(computeOgive(yearlyLengthTrajectories,
lengthAtMaturity = par[1]),
by = 'age', suffix = c('.obs','.pred')) %>%
mutate(squareError = (p.obs - p.pred)^2) %>%
summarise(SSE = sum(squareError)) %>%
unlist(use.names = FALSE)
return(SSE)
# epsilon = 1e-6 # to avoid log(0)
# minusSumLogL = ogiveObs %>%
# inner_join(computeOgive(yearlyLengthTrajectories,
# lengthAtMaturity = par[1]),
# by = 'age', suffix = c('.obs','.pred')) %>%
# summarise(minusSumLogL = - sum(p.obs * log(p.pred + epsilon) + (1 - p.obs) * log(1 - p.pred + epsilon))) %>%
# unlist(use.names = FALSE)
# return(minusSumLogL)
}
```{r load growth parameters and temperature pattern }
objFn_A(par = c(lengthAtMaturity = 40), #vector including the 3 regional metapopulation
ogiveObs, regional_metapop = c("semelparous", "southern iteroparous", "northern iteroparous")
yearlyLengthTrajectories = lengthTrajectoriesInSpring)
res_A <- optim(par = c(lengthAtMaturity = 35), #Growth parameters from Stich et al, 2020
fn = objFn_A, growParStich <- Stich2020_sel %>%
ogiveObs = ogiveObs, filter(catchment %in% regional_metapop) %>%
yearlyLengthTrajectories = lengthTrajectoriesInSpring, select(catchment, Linf, K, L0_theo) %>%
method = 'Brent', rename("metapop" = "catchment") %>%
lower = 0, mutate(tempMinGrow = 1.6, tempOptGrow = 5, tempMaxGrow = 27.9,
upper = growPar$Linf) Linf = Linf/10, lengthAtHatching = L0_theo/10, kOpt = K,
sigmaDeltaLVonBert = 0.4)
```
computeOgive(lengthTrajectories = lengthTrajectoriesInSpring, #Application to GR3D-US
lengthAtMaturity = res_A$par) %>% #1.Growth
inner_join(ogiveObs, by = 'age', suffix = c(".pred", ".obs")) ##1.2. Defintion of growth curve
In GR3D, growth of individuals was modeled as a Von Bertalanffy growth function (Von Bertalanffy 1938).
To account for variability in individuals, the seasonal growth increment is computed based on the growth coefficient K that is linked to water temperature.
The effect of water temperature (T°) on the growth coefficient K was introduced through a dome-shape relationship similar to the one for reproduction.
# ----------------------------------------------------------------------------------------------------
# calibrate lengthAtMaturity and Kopt
objFn_B = function(par, growPar, ogiveObs, temperaturePattern, Nind = 100000, seasonSelect = 'spring', sigmaDeltaLVonBert = .2){
growPar$K = par[2]
growPar$sigmaDeltaLVonBert = sigmaDeltaLVonBert
ogivePredict <- computeMultipleLengthTrajectories(temperaturePattern , Nind = Nind, growPar) %>%
filter(season == seasonSelect) %>%
computeOgive(lengthAtMaturity = par[1])
SSE = ogiveObs %>%
inner_join(ogivePredict,
by = 'age', suffix = c('.obs','.pred')) %>%
mutate(squareError = (p.obs - p.pred)^2) %>%
summarise(SSE = sum(squareError)) %>%
unlist(use.names = FALSE)
return(SSE)
}
objFn_B(par = c(lengthAtMaturity = 35, Kopt = 0.49 ), $$
growPar = growPar,
ogiveObs = ogiveObs,
temperaturePattern = temperaturePattern,
Nind = 100000,
sigmaDeltaLVonBert = .2)
res_B <- optim(par = c(lengthAtMaturity = 35, Kopt = 0.49 ),
fn = objFn_B,
growPar = growPar,
ogiveObs = ogiveObs,
temperaturePattern = temperaturePattern,
Nind = 100000,
sigmaDeltaLVonBert = .2 ,
control = list(trace = 1,
reltol = 1e-5))
save(res_B, file = 'res_B.Rdata') mu= log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2
$$
So, growth is defined by 8 parameters, described in Table 1
computeMultipleLengthTrajectories(temperaturePattern, ```{r growth parameters include in the GR3D model, include = TRUE, fig.cap = "Table 1. Growth parameters include into the GR3D model "}
Nind = 100000,
growPar %>% mutate(K = res_B$par[2])) %>% #intial growth parameters for A.sapidissima (see initial parameter value from the Workshop)
filter(season == 'spring') %>% growParXML = growPar
computeOgive(lengthAtMaturity = res_B$par[1]) %>%
inner_join(ogiveObs, by = 'age', suffix = c('.pred', '.obs')) growParXML%>%
pivot_longer(tempMinGrow:lFirstMaturityForFemale, names_to = "parameters", values_to = "nominal values") %>%
flextable() %>%
autofit()
computeMultipleLengthTrajectories(temperaturePattern,
Nind = 10000,
growPar %>% mutate(K = res_B$par[2])) %>%
group_by(age) %>%
summarise(Lmin = min(L),
L025 = quantile(L, .25),
L25 = quantile(L, .25),
Lmed = median(L),
L75 = quantile(L, .75),
L975 = quantile(L, .975),
Lmax = max(L)) %>%
ggplot(aes(x = age)) +
geom_ribbon(aes(ymin = Lmin, ymax = Lmax), fill = "gray30", alpha = .2) +
geom_ribbon(aes(ymin = L025, ymax = L975), fill = "gray30", alpha = .2) +
geom_ribbon(aes(ymin = L25, ymax = L75), fill = "gray30", alpha = .2) +
geom_line(aes(y = Lmed)) +
labs(x = 'Age (y)', y = 'Length (cm)') +
geom_abline(slope =0, intercept = 44)
``` ```
```{r calibration K and lengthAtMaturity by gender with ogive, warning = "FALSE"} ##1.2. Temperature effect experienced by fish for growth
#par = res_D$par
#----------------------------------------------# 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.
#------------ Optim function ------------------#
#-----------------------------------------------#
objFn_D = function(par, fixedPar, ogivesObs, parGrowRef, temperaturePattern, 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.
Nind = 1000, seasonSelected = 'spring', At the end of fall, juveniles move seaward and reach the closest ‘wintering offshore basins’ to spend the winter months.
RNGseed =1) { 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.
fullPar <- enframe(c(par, fixedPar)) %>% pivot_wider() #full set of parameters
#extract parameters by sex
parFemale = fullPar %>% select(tempMinGrow, tempOptGrow, tempMaxGrow,
lengthAtHatching,
Linf = linfVonBert,
kOpt = kOptForFemale,
sigmaDeltaLVonBert,
lengthAtMaturity = lFirstMaturityForFemale)
parMale = fullPar %>% select(tempMinGrow, tempOptGrow, tempMaxGrow,
lengthAtHatching,
Linf = linfVonBert,
kOpt = kOptForMale,
sigmaDeltaLVonBert,
lengthAtMaturity = lFirstMaturityForMale)
# SSE for female ogive
ogivePredictFemale <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern ,
Nind = Nind,
growPar = parFemale,
RNGseed = RNGseed) %>%
filter(season == seasonSelected) %>%
computeOgive(lengthAtMaturity = parFemale$lengthAtMaturity) %>%
# add a penalty when p(max(age)) < 1 - avoid immature individuals at maximuma age
mutate(p = if_else(age == max(age) & p < 1, 10*(p-1)^2, p))
SSEfemale = ogivesObs %>% select(age, p = p_female) %>%
inner_join(ogivePredictFemale,
by = 'age', suffix = c('.obs','.pred')) %>%
mutate(squareError = (p.obs - p.pred)^2) %>%
summarise(SSE = sum(squareError)) %>%
unlist(use.names = FALSE)
# SSE for male ogive
ogivePredictMale <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern ,
Nind = Nind,
growPar = parMale,
RNGseed = RNGseed) %>%
filter(season == seasonSelected) %>%
computeOgive(lengthAtMaturity = parMale$lengthAtMaturity) %>%
# add a penalty when p(max(age)) < 1
mutate(p = if_else(age == max(age) & p < 1, 10*(p-1)^2, p))
SSEmale = ogivesObs %>% select(age, p = p_male) %>%
inner_join(ogivePredictMale,
by = 'age', suffix = c('.obs','.pred')) %>%
mutate(squareError = (p.obs - p.pred)^2) %>%
summarise(SSE = sum(squareError)) %>%
unlist(use.names = FALSE)
### SSE for growth curve
#growth curve for female by using temperature pattern with sigma = 0 (mean_curve)
femaleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern ,
Nind = 1,
growPar = parFemale %>% mutate(sigmaDeltaLVonBert=0) ,
RNGseed = RNGseed) %>%
filter(season == seasonSelected) %>%
select(age, Lfemale = L)#rename L by Lfemale
maleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern ,
Nind = 1,
growPar = parMale %>% mutate(sigmaDeltaLVonBert=0) ,
RNGseed = RNGseed) %>%
filter(season == seasonSelected) %>% #1 value for a season
select(age, Lmale = L)
#SCE1 for growth curve from Stich et al, and growth curve by sex
SCE1growth <-
femaleCurve %>%
inner_join(maleCurve, by = 'age') %>%
mutate(Lref = vonBertalanffyGrowth(age, parGrowRef$L0, parGrowRef$Linf, parGrowRef$K),
L2genders = (Lfemale + Lmale)/2,
#SE = ((Lfemale + Lmale)/2 - Lref)^2 ) %>%
SE = (((Lfemale + Lmale)/2 - Lref)/Lref)^2) %>% #TODO TESTER
summarise(SSE = sum(SE)) %>% unlist(use.names = FALSE)
SSEgrowAndMat = (SSEmale + SSEfemale)*10 + SCE1growth
return(SSEgrowAndMat)# 10 for weigthed the sumSquare compared to the sum squared for gives so that both sum squared varied in the same interval
}
To account for differences in temperatures experienced by fish to growth, we use the seasonal migration timing describe above.
```{r load temperature pattern }
#Temperature pattern relative to thermal conditions experienced by fish
temperaturePattern <- growthInBasin %>%
group_by(metapop, age, season) %>%
summarise(temperature = mean(temperature_RIO), .groups = 'drop') %>%
#filter(metapop == "semelparous") %>%
# mutate(season = case_when(
# (age - floor(age)) == 0.00 ~ 'spring',
# (age - floor(age)) == 0.25 ~ 'summer',
# (age - floor(age)) == 0.50 ~ 'fall',
# (age - floor(age)) == 0.75 ~ 'winter'
# )) %>%
select(metapop,age, season,temperature)
``` ```
```{r Camille - calibration by gender and metapop with ogive, warnin = "FALSE" } #3. Optimization of growth and maturity
```{r calibration by gender and metapop with ogive, warning = "FALSE" }
#----------------------------------------------# #----------------------------------------------#
#------------ Optim function ------------------# #------------ Optim function ------------------#
...@@ -455,60 +302,33 @@ SSEforAllMetapop = function(par, fixedPar, regional_metapop,ogivesObs, parGrowRe ...@@ -455,60 +302,33 @@ SSEforAllMetapop = function(par, fixedPar, regional_metapop,ogivesObs, parGrowRe
} }
#Verification # #Verification
tic() # tic()
SSEforAllMetapop(par = par_0, # SSEforAllMetapop(par = par_0,
fixedPar = fixedPar, # fixedPar = fixedPar,
regional_metapop = regional_metapop, # regional_metapop = regional_metapop,
ogivesObs = ogivesObs, # ogivesObs = ogivesObs,
parGrowRef = parGrowRef, # parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern, # temperaturePattern = temperaturePattern,
Nind = 1000, # Nind = 1000,
seasonSelected = "spring", # seasonSelected = "spring",
RNGseed = 1) # RNGseed = 1)
toc() # toc()
tic() # tic()
metapopulation = "semelparous" # SSEforOneMetapop(fullPar = enframe(c(par_0, fixedPar)) %>% pivot_wider(),
SSEforOneMetapop(fullPar = enframe(c(par_0, fixedPar)) %>% pivot_wider(), # ogivesObs = ogivesObs %>% filter(metapop == metapopulation),
ogivesObs = ogivesObs %>% filter(metapop == metapopulation), # parGrowRef = parGrowRef %>% filter(metapop == metapopulation),
parGrowRef = parGrowRef %>% filter(metapop == metapopulation), # temperaturePatternforOneMetapop = temperaturePattern %>% filter(metapop == metapopulation),
temperaturePatternforOneMetapop = temperaturePattern %>% filter(metapop == metapopulation), # Nind = 1000,
Nind = 1000, # seasonSelected = 'spring',
seasonSelected = 'spring', # RNGseed = 1)
RNGseed = 1) # toc()
toc()
# computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern = temperaturePattern %>% filter(metapop == "semelparous"),
# Nind = 10,
# growPar = parFemale,
# RNGseed = 1)
``` ```
```{r set parameters values for calibration, warning = "FALSE"} ```{r set parameters values for calibration, warning = "FALSE"}
#----------------------------------------------#
#------------ Starting parameters values ------------------#
#-----------------------------------------------#
#
# ogivesObs = data.frame(age = seq.int(9), p_female = c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1),
# p_male = c(0, 0, 0.07, 0.29, 0.60, 0.86, 1, 1, 1))
# for 3 metapopulations
ogivesObs = expand_grid(age = seq.int(9), metapop = regional_metapop) %>%
arrange(metapop, age) %>%
group_by(metapop) %>%
mutate(p_female = fifelse(
metapop == "semelparous", c(0, 0, 0.01, 0.09, 0.33, 0.63, 0.92, 1, 1), fifelse (metapop == "southern iteroparous", c(0, 0, 0,0.04,0.27,0.64,0.81,0.9,1), c(0,0,0,0.04,0.69,0.69,0.9,1,1)))) %>%
mutate(p_male = fifelse(
metapop == "semelparous", c(0, 0, 0.07, 0.29, 0.60, 0.86, 1, 1, 1), fifelse (metapop == "southern iteroparous", c(0,0,0,0.12,0.5,0.79,0.8,0.84,1), c(0,0,0.03,0.21,0.82,1,1,1,1)))) %>%
ungroup()
#----------Growth parameters------------# #----------Growth parameters------------#
#---------------------------------------- #----------------------------------------
...@@ -523,7 +343,7 @@ parGrowRef = growParStich %>% ...@@ -523,7 +343,7 @@ parGrowRef = growParStich %>%
#----------Parameters to be optimized over -----------# #----------Parameters to be optimized over -----------#
#----------------------------------------------------# #----------------------------------------------------#
#Case 1 #Case A
#parameters from optimisation and XML file #parameters from optimisation and XML file
par_0 = c(tempOptGrow = 5, #from optimisation on growth curves par_0 = c(tempOptGrow = 5, #from optimisation on growth curves
kOptForFemale = 0.32, #XML - Alosa alosa kOptForFemale = 0.32, #XML - Alosa alosa
...@@ -537,6 +357,7 @@ fixedPar = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 R ...@@ -537,6 +357,7 @@ fixedPar = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 R
tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020 tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020
lengthAtHatching = 2.8,# XML - Alosa sapidissima lengthAtHatching = 2.8,# XML - Alosa sapidissima
sigmaDeltaLVonBert = .2) #XML - Alosa alosa sigmaDeltaLVonBert = .2) #XML - Alosa alosa
#optimisation #optimisation
res_Abis <- optim(par = par_0, res_Abis <- optim(par = par_0,
fn = SSEforAllMetapop, fn = SSEforAllMetapop,
...@@ -553,25 +374,105 @@ res_Abis <- optim(par = par_0, ...@@ -553,25 +374,105 @@ res_Abis <- optim(par = par_0,
Nind = 5000, Nind = 5000,
RNGseed = 1) RNGseed = 1)
#----------------------------------------------------#
#Case B
fixedPar_1 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9,
lengthAtHatching = 5.87)#mean lengthAthatching derived from Stich et al, 2020
res_B <- optim(par = par_0,
fn = SSEforAllMetapop,
#method = "BFGS",
control = list(trace = 1,
maxit = 1000),
#hessian = TRUE,
fixedPar = fixedPar_1,
regional_metapop = regional_metapop,
ogivesObs = ogivesObs,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
seasonSelected = "spring",
Nind = 5000,
RNGseed =1)
#----------------------------------------------------------#
#Case C
par_1 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 76,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45,
lengthAtHatching = 2.8,
sigmaDeltaLVonBert = .2)
fixedPar_2 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9)
#Case 2
#parameters from optimisation and XML file res_C <- optim(par = par_1,
fn = SSEforAllMetapop,
#method = "BFGS",
control = list(trace = 1,
maxit = 1000),
#hessian = TRUE,
fixedPar = fixedPar_2,
regional_metapop = regional_metapop,
ogivesObs = ogivesObs,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
seasonSelected = "spring",
Nind = 5000,
RNGseed =1)
#--------------------------------------------------
#Case D
par_2 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 47.5,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45,
lengthAtHatching = 5.87)
fixedPar_3 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9,
sigmaDeltaLVonBert = .2)
res_D <- optim(par = par_2,
fn = SSEforAllMetapop,
#method = "BFGS",
control = list(trace = 1,
maxit = 1000),
#hessian = TRUE,
fixedPar = fixedPar_3,
regional_metapop = regional_metapop,
ogivesObs = ogivesObs,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
seasonSelected = "spring",
Nind = 5000,
RNGseed =1)
#-------------------------------------------------------------------------
#case E:
par_0E = c(tempOptGrow = 5, #from optimisation on growth curves par_0E = c(tempOptGrow = 5, #from optimisation on growth curves
kOptForFemale = 0.32, #XML - Alosa alosa kOptForFemale = 0.32, #XML - Alosa alosa
kOptForMale = .21, #XML - Alosa alosa kOptForMale = .21, #XML - Alosa alosa
linfVonBert = 76, #XML - Alosa alosa linfVonBert = 76, #XML - Alosa alosa
lFirstMaturityForMale = 40, #XML - Alosa sapidissima lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45, lFirstMaturityForFemale = 45,
lengthAtHatching = 2.8)#XML - Alosa sapidissima lengthAtHatching = 2.8)#XML - Alosa sapidissima
fixedPar_E = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020 fixedPar_E = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020
tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020 tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020
# XML - Alosa sapidissima # XML - Alosa sapidissima
sigmaDeltaLVonBert = .4) sigmaDeltaLVonBert = .4)
#optimisation #optimisation
res_E <- optim(par = par_0E, res_E <- optim(par = par_0E,
fn = SSEforAllMetapop, fn = SSEforAllMetapop,
...@@ -588,10 +489,8 @@ res_E <- optim(par = par_0E, ...@@ -588,10 +489,8 @@ res_E <- optim(par = par_0E,
Nind = 5000, Nind = 5000,
RNGseed = 1) RNGseed = 1)
#---------------------------------#
#Case 3 #Case F
#parameters from optimisation and XML file
par_0F = c(tempOptGrow = 5, #from optimisation on growth curves par_0F = c(tempOptGrow = 5, #from optimisation on growth curves
kOptForFemale = 0.32, #XML - Alosa alosa kOptForFemale = 0.32, #XML - Alosa alosa
kOptForMale = .21, #XML - Alosa alosa kOptForMale = .21, #XML - Alosa alosa
...@@ -621,6 +520,8 @@ res_F <- optim(par = par_0F, ...@@ -621,6 +520,8 @@ res_F <- optim(par = par_0F,
Nind = 5000, Nind = 5000,
RNGseed = 1) RNGseed = 1)
#-------------------------------------------------------------
#Case G
par_0G = c(tempOptGrow = 5, #from optimisation on growth curves par_0G = c(tempOptGrow = 5, #from optimisation on growth curves
kOptForFemale = 0.32, #XML - Alosa alosa kOptForFemale = 0.32, #XML - Alosa alosa
...@@ -628,7 +529,6 @@ par_0G = c(tempOptGrow = 5, #from optimisation on growth curves ...@@ -628,7 +529,6 @@ par_0G = c(tempOptGrow = 5, #from optimisation on growth curves
linfVonBert = 76) #XML - Alosa alosa linfVonBert = 76) #XML - Alosa alosa
#XML - Alosa sapidissima #XML - Alosa sapidissima
fixedPar_G = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020 fixedPar_G = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020
tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020 tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020
# XML - Alosa sapidissima # XML - Alosa sapidissima
...@@ -652,99 +552,79 @@ res_G <- optim(par = par_0G, ...@@ -652,99 +552,79 @@ res_G <- optim(par = par_0G,
Nind = 5000, Nind = 5000,
RNGseed = 1) RNGseed = 1)
#--------------------------------------------------------
#case H: without ponderation of SSEmale + SSEfemale *10
par_0H = c(tempOptGrow = 5, #from optimisation on growth curves
kOptForFemale = 0.32, #XML - Alosa alosa
kOptForMale = .21, #XML - Alosa alosa
linfVonBert = 76, #XML - Alosa alosa
lFirstMaturityForMale = 40, #XML - Alosa sapidissima
lFirstMaturityForFemale = 45,
lengthAtHatching = 2.8)#XML - Alosa sapidissima
#Case 2
fixedPar_1 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9,
lengthAtHatching = 5.87)
fixedPar_H = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020
tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020
# XML - Alosa sapidissima
sigmaDeltaLVonBert = .4)
res_B <- optim(par = par_0, res_H <- optim(par = par_0H,
fn = SSEforAllMetapop, fn = SSEforAllMetapop,
#method = "BFGS", #method = "L-BFGS-B",
control = list(trace = 1, control = list(trace = 1,
maxit = 1000), maxit = 1000),
#hessian = TRUE, #hessian = TRUE,
fixedPar = fixedPar_1, fixedPar = fixedPar_H,
regional_metapop = regional_metapop, regional_metapop = regional_metapop,
ogivesObs = ogivesObs, ogivesObs = ogivesObs,
parGrowRef = parGrowRef, parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern, temperaturePattern = temperaturePattern,
seasonSelected = "spring", seasonSelected = "spring",
Nind = 5000, Nind = 5000,
RNGseed =1) RNGseed = 1)
save(res_A, res_B, res_C, res_D, res_E,res_F,res_G, file = "res_optim.Rdata") #------------------------
#Case I: without ponderation of SSEmale + SSEfemale *10
par_0I = c(tempOptGrow = 5, #from optimisation on growth curves
kOptForFemale = 0.32, #XML - Alosa alosa
kOptForMale = .21, #XML - Alosa alosa
lFirstMaturityForMale = 40, #XML - Alosa sapidissima
lFirstMaturityForFemale = 45,
lengthAtHatching = 2.8)#XML - Alosa sapidissima
#Case 3
par_1 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 76,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45,
lengthAtHatching = 2.8,
sigmaDeltaLVonBert = .2)
fixedPar_2 = c(tempMinGrow = 1.6, fixedPar_I = c(tempMinGrow = 1.6, #min temperature experienced by fish in the 11 RB from Stich et al, 2020
tempMaxGrow = 27.9) tempMaxGrow = 27.9,#max temperature experienced by fish in the 11 RB from Stich et al, 2020
# XML - Alosa sapidissima
sigmaDeltaLVonBert = .4,
linfVonBert = 50)
res_C <- optim(par = par_1, res_I <- optim(par = par_0I,
fn = SSEforAllMetapop, fn = SSEforAllMetapop,
#method = "BFGS", #method = "L-BFGS-B",
control = list(trace = 1, control = list(trace = 1,
maxit = 1000), maxit = 1000),
#hessian = TRUE, #hessian = TRUE,
fixedPar = fixedPar_2, fixedPar = fixedPar_I,
regional_metapop = regional_metapop, regional_metapop = regional_metapop,
ogivesObs = ogivesObs, ogivesObs = ogivesObs,
parGrowRef = parGrowRef, parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern, temperaturePattern = temperaturePattern,
seasonSelected = "spring", seasonSelected = "spring",
Nind = 5000, Nind = 5000,
RNGseed =1) RNGseed = 1)
par_2 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 47.5,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45,
lengthAtHatching = 5.87)
fixedPar_3 = c(tempMinGrow = 1.6, #save results from otpimisation into Rdata
tempMaxGrow = 27.9, save(res_A, res_B, res_C, res_D, res_E,res_F,res_G,res_H,res_I, file = "res_optim.Rdata")
sigmaDeltaLVonBert = .2)
res_D <- optim(par = par_2,
fn = SSEforAllMetapop,
#method = "BFGS",
control = list(trace = 1,
maxit = 1000),
#hessian = TRUE,
fixedPar = fixedPar_3,
regional_metapop = regional_metapop,
ogivesObs = ogivesObs,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
seasonSelected = "spring",
Nind = 5000,
RNGseed =1)
``` ```
```{r plot, warning = "FALSE"}
#set of optim parameters for growth ```{r local function to export simulated growth and ogives df , warning = "FALSE"}
fullParOptim = enframe(c(res_E$par,fixedPar_E)) %>%
pivot_wider()
#compute growth curves for one metapopulation #compute simulated growth curves for one metapopulation
growthCurvesForOneMetapop = function(fullPar, parGrowRef, temperaturePatternforOneMetapop, growthCurvesForOneMetapop = function(fullPar, parGrowRef, temperaturePatternforOneMetapop,
Nind = 1000, Nind = 1000,
RNGseed =1, RNGseed =1,
...@@ -794,7 +674,7 @@ growthCurvesForOneMetapop = function(fullPar, parGrowRef, temperaturePatternforO ...@@ -794,7 +674,7 @@ growthCurvesForOneMetapop = function(fullPar, parGrowRef, temperaturePatternforO
} }
#compute predicted ogives for one metapopulation
predictedOgivesForOneMetapop = function(fullPar,ogivesObs, temperaturePatternforOneMetapop, predictedOgivesForOneMetapop = function(fullPar,ogivesObs, temperaturePatternforOneMetapop,
Nind = 1000, Nind = 1000,
RNGseed =1, RNGseed =1,
...@@ -856,7 +736,7 @@ predictedOgivesForOneMetapop = function(fullPar,ogivesObs, temperaturePatternfor ...@@ -856,7 +736,7 @@ predictedOgivesForOneMetapop = function(fullPar,ogivesObs, temperaturePatternfor
} }
#compute growth and ogives for all metapopulations #compute simulated growth and ogives for all metapopulations
growthOptimForAllMetapop = list() growthOptimForAllMetapop = list()
predictedOgivesForAllMetapop = list() predictedOgivesForAllMetapop = list()
...@@ -889,10 +769,26 @@ for (metapopulation in regional_metapop){ ...@@ -889,10 +769,26 @@ for (metapopulation in regional_metapop){
growthOptim_df = rbindlist(growthOptimForAllMetapop) growthOptim_df = rbindlist(growthOptimForAllMetapop)
predictOgives_df = rbindlist(predictedOgivesForAllMetapop) predictOgives_df = rbindlist(predictedOgivesForAllMetapop)
```
```{r plot cruves from results, warning = "FALSE"}
#set of optim parameters for growth
fullParOptim = enframe(c(res_H$par,fixedPar_H)) %>%
pivot_wider()
#set of reference parameters for growth
parGrowRef = growParStich %>%
select(L0 = lengthAtHatching, K, Linf, metapop)
#plot growth curves #plot growth curves
growthOptim_df %>% growthOptim_df %>%
pivot_longer(Lfemale:L2genders,names_to = "Ltype", values_to = "L") %>% dplyr::rename("Sex_combined" = "L2genders",
group_by(age,metapop, Ltype) %>% "Male" = "Lmale",
"Female" = "Lfemale",
"Ref" = "Lref") %>%
pivot_longer(Female:Sex_combined,names_to = "Length", values_to = "L") %>%
group_by(age,metapop, Length) %>%
summarise(Lmin = min(L), summarise(Lmin = min(L),
L025 = quantile(L, .25), L025 = quantile(L, .25),
L25 = quantile(L, .25), L25 = quantile(L, .25),
...@@ -904,22 +800,22 @@ growthOptim_df %>% ...@@ -904,22 +800,22 @@ growthOptim_df %>%
#geom_ribbon(aes(ymin = Lmin, ymax = Lmax), fill = "gray30", alpha = .2) + #geom_ribbon(aes(ymin = Lmin, ymax = Lmax), fill = "gray30", alpha = .2) +
#geom_ribbon(aes(ymin = L025, ymax = L975), fill = "gray30", alpha = .2) + #geom_ribbon(aes(ymin = L025, ymax = L975), fill = "gray30", alpha = .2) +
#geom_ribbon(aes(ymin = L25, ymax = L75), fill = "gray30", alpha = .2) + #geom_ribbon(aes(ymin = L25, ymax = L75), fill = "gray30", alpha = .2) +
geom_line(aes(y = Lmed, col = Ltype)) + geom_line(aes(y = Lmed, col = Length)) +
labs(x = 'Age (y)', y = 'Length (cm)')+ labs(x = 'Age (y)', y = 'Length (cm)',color = 'Length Type')+
facet_grid(~metapop) facet_grid(~metapop)
#plot ogives #plot ogives
predictOgives_df %>% predictOgives_df %>%
pivot_longer(contains('p.'),names_to = "type", values_to = "ogives") %>% pivot_longer(contains('p.'),names_to = "type", values_to = "ogives") %>%
mutate(sex = case_when( mutate(sex = case_when(
grepl("_female", type)~ "female", TRUE ~ "male")) %>% grepl("_female", type)~ "female", TRUE ~ "male")) %>%
ggplot(aes(x=age))+ ggplot(aes(x=age))+
geom_line(aes(y = ogives, col = type)) + geom_line(aes(y = ogives, col = type)) +
labs(x = 'Age (y)', y = 'Length (cm)')+ labs(x = 'Age (y)', y = 'Length (cm)', color = "Data Type")+
facet_grid(sex~metapop)+ scale_color_manual(values =c("#FC4E07","#FC4E07","#E7B800","#E7B800"))+
guides(color =FALSE) facet_grid(.~metapop)
#guides(color =FALSE)
predictOgives_df %>% predictOgives_df %>%
...@@ -930,9 +826,9 @@ growthOptim_df %>% ...@@ -930,9 +826,9 @@ growthOptim_df %>%
ggplot(aes(x=age))+ ggplot(aes(x=age))+
geom_line(aes(y = ogives_obs), col = 'red') + geom_line(aes(y = ogives_obs), col = 'red') +
geom_line(aes(y = ogives_pred), col ="orange" )+ geom_line(aes(y = ogives_pred), col ="orange" )+
labs(x = 'Age (y)', y = 'Length (cm)')+ labs(x = 'Age (y)', y = 'Length (cm)', legend = c("Obs", "Pred"))+
facet_grid(sex~metapop)+ facet_grid(sex~metapop)
guides(color =FALSE)
``` ```
......
Supports Markdown
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