Commit cb06d07d authored by Poulet Camille's avatar Poulet Camille
Browse files

Merge branch 'exploration_GR3D_process' of...

Merge branch 'exploration_GR3D_process' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D into exploration_GR3D_process
parents 73f7a2f5 3f5489ee
......@@ -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>
......
......@@ -28,10 +28,19 @@ load('SOS.rdata')
```
```{r}
source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
growParXML = data.frame(t(unlist(fishXML[[1]][["processes"]][["processesEachStep"]][["species.Grow"]]))) %>%
select(-synchronisationMode ) %>%
mutate(linfVonBertForMale = fishXML[[1]][["linfVonBertForMale"]],
linfVonBertForFemale = fishXML[[1]][["linfVonBertForFemale"]],
lengthAtHatching = fishXML[[1]][["lengthAtHatching"]],
lFirstMaturityForMale = fishXML[[1]][["lFirstMaturityForMale"]],
lFirstMaturityForFemale = fishXML[[1]][["lFirstMaturityForFemale"]]) %>%
mutate_all(~as.numeric(as.character(.)))
growPar <- Stich2020_sel %>% filter(catchment == "semelparous") %>%
growParStich <- 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,
......@@ -49,6 +58,10 @@ temperaturePattern <- growthInBasin %>%
# )) %>%
select(age, season,temperature)
```
```{r local function}
# generate a cohort of length trajectories
computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar){
ages = temperaturePattern$age
......@@ -80,19 +93,64 @@ computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, grow
return(res)
}
computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
Nind = 10,
growPar,
RNGseed =1){
set.seed(RNGseed)
ages = temperaturePattern$age
res <- expand_grid(ind = seq.int(Nind), age = ages) %>%
mutate(random = rnorm(Nind * length(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)
rnd <- res %>% filter(age == currentAge) %>% select(random) %>% unlist(use.names = FALSE)
currentL <- vonBertalanffyWithRandomVector(L = previousL,
L0 = growPar$lengthAtHatching,
Linf = growPar$Linf,
K = growPar$kOpt,
timeStepDuration = currentAge - previousAge,
randomVector = rnd,
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()
vonBertalanffyWithRandomVector = function(L, L0, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
if (sigma == 0) {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
increment = exp(mu)
}
else {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
increment = exp(randomVector * sigma + mu)
}
L = L + increment * tempEffect
return(L)
}
computeOgive = function(yearlyLengthTrajectories, lengthAtMaturity) {
ogive <- yearlyLengthTrajectories %>%
computeOgive = function(lengthTrajectories, lengthAtMaturity) {
ogive <- lengthTrajectories %>%
distinct(age) %>% arrange(age) %>%
left_join(yearlyLengthTrajectories %>% filter(L > lengthAtMaturity) %>%
left_join(lengthTrajectories %>% filter(L > lengthAtMaturity) %>%
group_by(ind) %>%
summarise(age = min(age), .groups = 'drop') %>%
group_by(age) %>%
......@@ -100,14 +158,26 @@ computeOgive = function(yearlyLengthTrajectories, lengthAtMaturity) {
by = 'age') %>%
replace_na(list(mature = 0)) %>%
mutate(immature = sum(mature) - cumsum(mature),
p = if_else(mature + immature == 0, 1, mature/(mature + immature))) %>%
p = if_else(mature + immature == 0, 0, mature/(mature + immature))) %>%
select(age, p)
return(ogive)
}
```
```{r test}
tic()
computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growParStich) %>%
filter(season == 'spring')
toc()
tic()
computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 100000, growParStich) %>%
filter(season == 'spring')
toc()
computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
Nind = 10000, growPar) %>%
......@@ -116,21 +186,18 @@ computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
```
```{r ogive varaibility}
```{r ogive variability}
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 envelop}
# envelop of growth curves
computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growPar) %>%
```{r growth curve envelope}
# envelope of growth curves
computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growParStich) %>%
group_by(age) %>%
summarise(Lmin = min(L),
L025 = quantile(L, .25),
......@@ -149,11 +216,11 @@ computeMultipleLengthTrajectories(temperaturePattern , Nind = 1000, growPar) %>%
# verify that all L < Linf
computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growPar) %>%
filter(L>growPar$Linf)
computeMultipleLengthTrajectories(temperaturePattern , Nind = 10000, growParStich) %>%
filter(L>growParStich$Linf)
```
```{r}
```{r optimisation}
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) %>%
......@@ -161,18 +228,26 @@ lengthTrajectoriesInSpring = computeMultipleLengthTrajectories(temperaturePatter
# calibrate only lengthAtMaturity
objFn_A = function(par, ogiveObs, yearlyLengthTrajectories){
SSE = ogiveObs %>%
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)) %>%
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)
}
objFn_A(par = c(lengthAtMaturity = 35),
objFn_A(par = c(lengthAtMaturity = 40),
ogiveObs,
yearlyLengthTrajectories = lengthTrajectoriesInSpring)
......@@ -191,12 +266,12 @@ computeOgive(yearlyLengthTrajectories = lengthTrajectoriesInSpring,
# ----------------------------------------------------------------------------------------------------
# calibrate lengthAtMaturity and Kopt
objFn_B = function(par, growPar, ogiveObs, temperaturePattern, Nind = 100000, season = 'spring', sigmaDeltaLVonBert = .2){
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 == season) %>%
filter(season == seasonSelect) %>%
computeOgive(lengthAtMaturity = par[1])
SSE = ogiveObs %>%
......@@ -252,4 +327,129 @@ computeMultipleLengthTrajectories(temperaturePattern,
geom_line(aes(y = Lmed)) +
labs(x = 'Age (y)', y = 'Length (cm)') +
geom_abline(slope =0, intercept = 44)
```
\ No newline at end of file
```
```{r calibration K and lengthAtMaturity by gender with ogive}
#par = res_D$par
objFn_D = function(par, fixedPar, ogivesObs, parGrowRef, temperaturePattern,
Nind = 1000, seasonSelected = 'spring',
RNGseed =1) {
fullPar <- enframe(c(par, fixedPar)) %>% pivot_wider()
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
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
femaleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern ,
Nind = 1,
growPar = parFemale %>% mutate(sigmaDeltaLVonBert=0) ,
RNGseed = RNGseed) %>%
filter(season == seasonSelected) %>%
select(age, Lfemale = L)
maleCurve <- computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern ,
Nind = 1,
growPar = parMale %>% mutate(sigmaDeltaLVonBert=0) ,
RNGseed = RNGseed) %>%
filter(season == seasonSelected) %>%
select(age, Lmale = L)
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 ) %>%
summarise(SSE = sum(SE)) %>% unlist(use.names = FALSE)
return((SSEmale + SSEfemale)*10 + SCE1growth)
}
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.01, 0.15, 0.45, 0.75, 0.92, 1, 1))
#TODO verifier si ce sont bien les parametres de Stich
parGrowRef = growParStich %>% select(L0 = lengthAtHatching, K, Linf)
par = c(tempOptGrow = 5,
kOptForFemale = 0.32,
kOptForMale = .21,
linfVonBert = 76,
lFirstMaturityForMale = 40,
lFirstMaturityForFemale = 45,
sigmaDeltaLVonBert = .2)
fixedPar = c(tempMinGrow = 1.6,
tempMaxGrow = 27.9,
lengthAtHatching = parGrowRef$L0 )
# objFn_D(par = par,
# fixedPar = fixedPar,
# ogivesObs = ogivesObs,
# parGrowRef = parGrowRef,
# temperaturePattern = temperaturePattern,
# Nind = 5000)
res_D <- optim(par = par,
fn = objFn_D,
fixedPar = fixedPar,
ogivesObs = ogivesObs,
parGrowRef = parGrowRef,
temperaturePattern = temperaturePattern,
Nind = 5000,
RNGseed =1,
control = list(trace = 1,
maxit = 1000))
```
---
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}
source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
growParXML = data.frame(t(unlist(fishXML[[1]][["processes"]][["processesEachStep"]][["species.Grow"]]))) %>%
select(-synchronisationMode ) %>%
mutate(linfVonBertForMale = fishXML[[1]][["linfVonBertForMale"]],
linfVonBertForFemale = fishXML[[1]][["linfVonBertForFemale"]],
lengthAtHatching = fishXML[[1]][["lengthAtHatching"]],
lFirstMaturityForMale = fishXML[[1]][["lFirstMaturityForMale"]],
lFirstMaturityForFemale = fishXML[[1]][["lFirstMaturityForFemale"]]) %>%
mutate_all(~as.numeric(as.character(.)))
growParStich <- 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)
```
```{r local function}
# 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)
}
computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern,
Nind = 10,
growPar,
RNGseed =1){
set.seed(RNGseed)
ages = temperaturePattern$age
res <- expand_grid(ind = seq.int(Nind), age = ages) %>%
mutate(random = rnorm(Nind * length(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)
rnd <- res %>% filter(age == currentAge) %>% select(random) %>% unlist(use.names = FALSE)
currentL <- vonBertalanffyWithRandomVector(L = previousL,
L0 = growPar$lengthAtHatching,
Linf = growPar$Linf,
K = growPar$kOpt,
timeStepDuration = currentAge - previousAge,
randomVector = rnd,
sigma = growPar$sigmaDeltaLVonBert,
tempEffect = tempEffect)
res = res %>% mutate(L = replace(L, age == ages[i], currentL) )
}
return(res)
}
vonBertalanffyWithRandomVector = function(L, L0, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
if (sigma == 0) {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
increment = exp(mu)
}
else {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
increment = exp(randomVector * sigma + mu)
}
L = L + increment * tempEffect
return(L)
}
computeOgive = function(lengthTrajectories, lengthAtMaturity) {
ogive <- lengthTrajectories %>%
distinct(age) %>% arrange(age) %>%
left_join(lengthTrajectories %>% 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, 0, mature/(mature + immature))) %>%
select(age, p)
return(ogive)
}
```
```{r test}
tic()
computeMultipleLengthTrajectories(temperaturePattern, Nind = 1000, growParStich) %>%
filter(season == 'spring')
toc()
tic()
computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 100000, growParStich) %>%
filter(season == 'spring')
toc()
computeMultipleLengthTrajectories(temperaturePattern %>% filter(age < 10) ,
Nind = 10000, growPar) %>%
filter(season == 'spring') %>%
computeOgive(lengthAtMaturity = 35)
```
```{r ogive variability}
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}
# envelope of growth curves
computeMultipleLengthTrajectories(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))