Commit 58dc246c authored by Poulet Camille's avatar Poulet Camille
Browse files

Correction in Ogives

parent d6fc7b15
......@@ -15,6 +15,7 @@ library(officedown)
library(officer)
library(tidyverse)
library(tictoc) # see timeR
library(data.table)
......@@ -27,6 +28,7 @@ source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
#load("growth_optimal.RData")
load('SOS.rdata')
load('res_optim.rdata')
```
```{r load growth parameters and tempertaure pattern }
......@@ -42,7 +44,7 @@ growParStich <- Stich2020_sel %>%
select(catchment, Linf, K, L0_theo) %>%
rename ("metapop" = "catchment") %>%
mutate(tempMinGrow = 1.6, tempOptGrow = 5, tempMaxGrow = 27.9,
Linf =Linf/10, lengthAtHatching= L0_theo/10, kOpt = K,
Linf = Linf/10, lengthAtHatching= L0_theo/10, kOpt = K,
sigmaDeltaLVonBert = 0.4)
temperaturePattern <- growthInBasin %>%
......@@ -133,6 +135,9 @@ computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
return(res)
}
#computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 10,
#growPar = growParUnisex,
# RNGseed = 1)
......@@ -168,7 +173,6 @@ computeOgive = function(lengthTrajectories, lengthAtMaturity) {
return(ogive)
}
```
```{r test }
......@@ -199,6 +203,24 @@ replicate(10, computeMultipleLengthTrajectories(temperaturePattern %>% filter(a
filter(season == 'spring') %>%
computeOgive(lengthAtMaturity = 35) %>% select(p) %>% unlist(use.names = FALSE))
toc()
#----------------------------------------
replicate(10, computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
Nind = 100000, growPar) %>%
filter(season == 'spring') %>%
computeOgive(lengthAtMaturity = 35) %>% select(p) %>% unlist(use.names = FALSE))
#plot results
growParFemale = growPar %>% pivot_longer(1:11,names_to = "name",values_to = 'par') %>%
tempMinGrow
par_maturity_optim= enframe(c(res_A$par,fixedPar)) %>%
pivot_wider() %>%
select(contains("Maturity"))
```
```{r growth curve envelope}
......@@ -475,6 +497,7 @@ SSEforOneMetapop = function(fullPar, ogivesObs, parGrowRef, temperaturePatternf
select(age, p = p_female) %>%
inner_join(ogivePredictFemale,
by = 'age', suffix = c('.obs','.pred')) %>%
#add penality
mutate(squareError = (p.obs - p.pred)^2) %>%
summarise(SSE = sum(squareError),) %>%
unlist(use.names = FALSE)
......@@ -533,7 +556,7 @@ SSEforOneMetapop = function(fullPar, ogivesObs, parGrowRef, temperaturePatternf
#Fonction to be optimzed over
SSEforAllMetapop = function(par, fixedPar, regional_metapop,ogiveObs, parGrowRef, temperaturePattern, Nind, seasonSelected, RNGseed = 1){
SSEforAllMetapop = function(par, fixedPar, regional_metapop,ogivesObs, parGrowRef, temperaturePattern, Nind, seasonSelected, RNGseed = 1){
fullPar = enframe(c(par, fixedPar)) %>% pivot_wider()
......@@ -563,16 +586,17 @@ SSEforAllMetapop = function(par, fixedPar, regional_metapop,ogiveObs, parGrowRef
#Verification
SSEforAllMetapop(regional_metapop = regional_metapop,
par = par_0,
tic()
SSEforAllMetapop(par = par_0,
fixedPar = fixedPar,
ogiveObs = ogiveObserved,
regional_metapop = regional_metapop,
ogivesObs = ogivesObs,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
Nind = 1000,
seasonSelected = "spring",
RNGseed = 1)
toc()
# computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern = temperaturePattern %>% filter(metapop == "semelparous"),
# Nind = 10,
......@@ -594,7 +618,7 @@ SSEforAllMetapop(regional_metapop = regional_metapop,
# p_male = c(0, 0, 0.07, 0.29, 0.60, 0.86, 1, 1, 1))
# for 3 metapopulations
ogivesObserved = expand_grid(age = seq.int(9), metapop = regional_metapop) %>%
ogivesObs = expand_grid(age = seq.int(9), metapop = regional_metapop) %>%
arrange(metapop, age) %>%
group_by(metapop) %>%
mutate(p_female = fifelse(
......@@ -625,55 +649,182 @@ par_0 = c(tempOptGrow = 5, #from optimisation on growth curves
kOptForMale = .21, #XML - Alosa alosa
linfVonBert = 76, #XML - Alosa alosa
lFirstMaturityForMale = 40, #XML - Alosa sapidissima
lFirstMaturityForFemale = 45,#XML - Alosa sapidissima
sigmaDeltaLVonBert = .2) #XML - Alosa alosa
lFirstMaturityForFemale = 45)#XML - Alosa sapidissima
fixedPar = 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
lengthAtHatching = 2.8)# XML - Alosa sapidissima
lengthAtHatching = 2.8,# XML - Alosa sapidissima
sigmaDeltaLVonBert = .2) #XML - Alosa alosa
#optimisation
res_A <- optim(par = par_0,
res_Abis <- optim(par = par_0,
fn = SSEforAllMetapop,
method = "BFGS",
method = "L-BFGS-B",
control = list(trace = 1,
maxit = 1000),
hessian = TRUE,
#hessian = TRUE,
fixedPar = fixedPar,
regional_metapop = regional_metapop,
ogivesObs = ogivesObserved,
ogivesObs = ogivesObs,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
seasonSelected = "spring",
Nind = 5000,
RNGseed =1)
RNGseed = 1)
#plot results
par_maturity_optim= enframe(c(res_A$par,fixedPar)) %>%
pivot_wider() %>%
select(contains("Maturity"))
par_grow_optim_female = enframe(c(res_A$par,fixedPar)) %>%
pivot_wider() %>%
select(-c('kOptForMale','lFirstMaturityForMale'))
#growth curves
growthCurves_female = computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , Nind = 10000,
par_grow_optim_female %>%
pivot_wider(), RNGseed = 1)
computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern , Nind = 1000, growParStich) %>%
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 = c(40.3, 42))
#growth curve according to metapopulations
growthCurves %>%
ggplot(aes(x = age)) +
geom_line(aes(y = L, col = metapop), linetype = "dashed") +
geom_line(aes(y = LStich, col = metapop))+
labs( x = 'Age (year)', y = 'Fork length (cm)', color = "Metapopulation (Stich)") +
#scale_x_continuous(limits = c(0, 10)) +
geom_hline(yintercept = 45, color = "pink", linetype =1)+ #female
geom_hline(yintercept = 40, color = "lightblue", linetype = 1)+
geom_vline(aes(age), xintercept = 8, linetype = 1)+
facet_wrap(.~ metapop, ncol =1)
#temperature effect
temp <- seq(0,45,0.1)
growthCurves %>%
ggplot(aes(x = temperature)) +
geom_line(aes(y = temperatureEffect)) +
geom_jitter(aes(y = temperatureEffect, col = metapop))
#Case 2
fixedPar_1 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9,
lengthAtHatching = 5.87)
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)
save(res_A, res_B, res_C, res_D,file = "res_optim.Rdata")
#Case 3
par_1 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 76,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45,
sigmaDeltaLVonBert = .2,
lengthAtHatching = 2.8)
lengthAtHatching = 2.8,
sigmaDeltaLVonBert = .2)
fixedPar_1 = c(tempMinGrow = 1.6,
fixedPar_2 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9)
#Case 3
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)
par_2 = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 76,
linfVonBert = 47.5,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45)
lFirstMaturityForFemale = 45,
lengthAtHatching = 5.87)
fixedPar_2 = c(tempMinGrow = 1.6,
fixedPar_3 = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9,
sigmaDeltaLVonBert = .2,
lengthAtHatching = 2.8)
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)
```
......
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