Commit 002ef709 authored by patrick.lambert's avatar patrick.lambert
Browse files

first attempt to calibrate maturation

parent 59f49844
---
date: "`r Sys.Date()`"
author: "Your Name"
title: "officedown template"
output:
officedown::rdocx_document:
mapstyles:
Normal: ['First Paragraph']
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.cap = TRUE)
library(officedown)
library(officer)
library(tidyverse)
library(tictoc) # see timeR
```
```{r}
rm(list = ls())
source("../GR3D_Rdescription/GR3Dfunction.R")
load('SOS.rdata')
```
```{r}
growPar <- Stich2020_sel %>% filter(catchment == "semelparous") %>%
select(Linf, K, L0_theo) %>%
mutate(tempMinGrow = 1.6, tempOptGrow = 5, tempMaxGrow = 27.9,
Linf =Linf/10, lengthAtHatching= L0_theo/10, kOpt = K,
sigmaDeltaLVonBert = 0.4)
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(age, season,temperature)
# generate a cohort of length trajectories
computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar){
ages = temperaturePattern$age
res <- expand_grid(ind = seq.int(Nind), age = ages) %>%
inner_join(temperaturePattern, by = 'age') %>%
arrange(age, ind) %>%
mutate(temperatureEffect = temperatureEffect(temperature,
growPar$tempMinGrow,
growPar$tempOptGrow,
growPar$tempMaxGrow),
L = if_else(age == 0, growPar$lengthAtHatching, 0))
for (i in 2:length(ages)) {
previousAge = ages[i - 1]
currentAge = ages[i]
tempEffect = res %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE)
previousL <- res %>% filter(age == previousAge) %>% select(L) %>% unlist(use.names = FALSE)
currentL <- vonBertalanffyWithNextIncrement(L = previousL,
L0 = growPar$lengthAtHatching,
Linf = growPar$Linf,
K = growPar$kOpt,
timeStepDuration = currentAge - previousAge,
sigma = growPar$sigmaDeltaLVonBert,
tempEffect = tempEffect)
res = res %>% mutate(L =replace(L, age == ages[i], currentL) )
}
return(res)
}
tic()
computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growPar) %>%
filter(season == 'spring')
toc()
computeOgive = function(yearlyLengthTrajectories, lengthAtMaturity) {
ogive <- yearlyLengthTrajectories %>%
distinct(age) %>% arrange(age) %>%
left_join(yearlyLengthTrajectories %>% filter(L > lengthAtMaturity) %>%
group_by(ind) %>%
summarise(age = min(age), .groups = 'drop') %>%
group_by(age) %>%
summarise(mature = n(), .groups = 'drop'),
by = 'age') %>%
replace_na(list(mature = 0)) %>%
mutate(immature = sum(mature) - cumsum(mature),
p = if_else(mature + immature == 0, 1, mature/(mature + immature))) %>%
select(age, p)
return(ogive)
}
computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
Nind = 10000, growPar) %>%
filter(season == 'spring') %>%
computeOgive(lengthAtMaturity = 35)
```
```{r ogive varaibility}
replicate(10, computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
Nind = 100000, growPar) %>%
filter(season == 'spring') %>%
computeOgive(lengthAtMaturity = 35) %>% select(p) %>% unlist(use.names = FALSE))
```
```{r growth curve envelop}
# envelop of growth curves
computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growPar) %>%
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))
# verify that all L < Linf
computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growPar) %>%
filter(L>growPar$Linf)
```
```{r}
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)
}
objFn_A(par = c(lengthAtMaturity = 35),
ogiveObs,
yearlyLengthTrajectories = lengthTrajectoriesInSpring)
res_A <- optim(par = c(lengthAtMaturity = 35),
fn = objFn_A,
ogiveObs = ogiveObs,
yearlyLengthTrajectories = lengthTrajectoriesInSpring,
method = 'Brent',
lower = 0,
upper = growPar$Linf)
computeOgive(yearlyLengthTrajectories = lengthTrajectoriesInSpring,
lengthAtMaturity = res_A$par) %>%
inner_join(ogiveObs, by = 'age', suffix = c(".pred", ".obs"))
# ----------------------------------------------------------------------------------------------------
# calibrate lengthAtMaturity and Kopt
objFn_B = function(par, growPar, ogiveObs, temperaturePattern, Nind = 100000, season = 'spring', sigmaDeltaLVonBert = .2){
growPar$K = par[2]
growPar$sigmaDeltaLVonBert = sigmaDeltaLVonBert
ogivePredict <- computeMultipleLengthTrajectories(temperaturePattern , Nind = Nind, growPar) %>%
filter(season == season) %>%
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')
computeMultipleLengthTrajectories(temperaturePattern,
Nind = 100000,
growPar %>% mutate(K = res_B$par[2])) %>%
filter(season == 'spring') %>%
computeOgive(lengthAtMaturity = res_B$par[1]) %>%
inner_join(ogiveObs, by = 'age', suffix = c('.pred', '.obs'))
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)
```
\ No newline at end of file
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