6-Growth_Calibration.Rmd 34.1 KB
Newer Older
Poulet Camille's avatar
Poulet Camille committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
---
title: "Growth Calibration"
author: "Camille POULET"
date: "23/02/2021"
output: word_document
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = FALSE, include = TRUE, 
                      warning = FALSE,
                      fig.cap = TRUE)

rm(list = ls())
```

```{r library, include=FALSE}
library(dplyr)
library(tidyr)
library(ggplot2)
21
library(ggrepel)
Poulet Camille's avatar
Poulet Camille committed
22
23
24
25
26
27
28
29
library(readr)
library(XML)
library(forcats)
library(knitr)
library(officedown)
library(flextable)
library(stringr)
library(tictoc)
patrick.lambert's avatar
patrick.lambert committed
30
library(tibble)
Poulet Camille's avatar
Poulet Camille committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

```

```{r GR3D functions and data, include=FALSE, echo = FALSE}

source("../GR3D_Rdescription/GR3Dfunction.R")
source("../GR3D_Rdescription/GR3D_NEA_env_data.R")
source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")

#source_rmd <- function(rmd_file){
#knitr::knit(rmd_file, output = tempfile())
#}

#source_rmd("NEAgrowth.Rmd")

#source(knitr::purl("NEAgrowth.Rmd", quiet=TRUE))

```

50
51
52
53
54
55
56
57
58
59
60
61
62
63
## 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.

Poulet Camille's avatar
Poulet Camille committed
64
65
66
67
```{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 %>% 
68
69
70
  distinct(wintering_offshore_basin_name,summering_offshore_basin_name) %>% 
  flextable() %>% 
  autofit()
Poulet Camille's avatar
Poulet Camille committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

#Compute an average temperature experienced by juveniles during their life at sea according to connections files i.e. wintering and summering areas.

connections_trans <- connections %>% 
  dplyr::select(-c(wintering_offshore_basin_id,
                   summering_offshore_basin_id)) %>% 
  pivot_longer(cols = c("wintering_offshore_basin_name",
                        "summering_offshore_basin_name"),
               names_to = "offshore_basin_type", 
               values_to = "offshore_basin_name") %>% 
  select(inshore_basin_name, offshore_basin_name, offshore_basin_type) %>% 
  inner_join(tempInOffshore, by = c("offshore_basin_name" = "basin_name")) %>% 
  filter(between(year,  1981, 2010)) %>% 
  mutate(offshore_basin_type = as.factor(offshore_basin_type)) %>% 
  mutate(offshore_basin_type = fct_recode(offshore_basin_type, 
                                          "summering" = "summering_offshore_basin_name",
                                          "wintering" = "wintering_offshore_basin_name"))  %>% 
  pivot_longer(cols = contains('temperature'),  names_to = "season", values_to = "temperature") %>% 
  mutate(season = str_remove(season, "_sea_surface_temperature")) 

##Select the temperature experienced by fish during the corresponding season 
#in winter
patrick.lambert's avatar
patrick.lambert committed
93
wint <- connections_trans %>% 
Poulet Camille's avatar
Poulet Camille committed
94
95
96
97
  filter(season %in% c("winter","spring"), offshore_basin_type == "wintering") %>% 
  select(inshore_basin_name,offshore_basin_name,year,season,temperature) 

#in summer
patrick.lambert's avatar
patrick.lambert committed
98
sum <- connections_trans %>% 
Poulet Camille's avatar
Poulet Camille committed
99
100
101
102
103
  filter(season %in% c("summer","fall"), offshore_basin_type == "summering") %>% 
  select(inshore_basin_name,offshore_basin_name,year,season,temperature)

#Combine the two files to build a summary dataframe
tempEffectInOffshore <- bind_rows(wint,sum) %>% 
104
105
106
  ungroup() %>% 
  arrange(inshore_basin_name, year,season) %>%
  group_by(inshore_basin_name,offshore_basin_name,season) %>% 
patrick.lambert's avatar
patrick.lambert committed
107
  summarise(temperature_O = mean(temperature), .groups = 'drop')
108
#mutate(temperatureEffect = temperatureEffect(Tmean, growPar$tempMinGrow, growPar$tempOptGrow, growPar$tempMaxGrow))
Poulet Camille's avatar
Poulet Camille committed
109
110
111
112

#delete intermediate files 
remove(sum,wint)

113
114
115
116
117
118
#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)
119
120

```
121
## Temperature in river and inshore basin
122

123
124
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.
Poulet Camille's avatar
Poulet Camille committed
125

126
```{r Build the temperature expericenced by fish during their residency in river and inshore basins, echo = FALSE, warning=FALSE}
127

128
tempEfffectInRI <- tempInriver  %>% 
129
130
131
132
133
134
135
136
137
  pivot_longer(cols = contains('temperature'),  names_to = "season", values_to = "river_temperature") %>% 
  mutate(season = str_remove(season, "_river_temperature")) %>% 
  inner_join(tempInInshore %>% 
               pivot_longer(cols = contains('temperature'),  names_to = "season", values_to = "inshore_temperature") %>% 
               mutate(season = str_remove(season, "_inshore_temperature")),
             by = c("basin_id", "basin_name", "year","season")) %>% 
  filter(between(year,  1981, 2010)) %>% 
  pivot_longer(c("river_temperature", "inshore_temperature"), names_to = "zone", values_to = "temperature") 

138
139
#TODO: code a optimiser
tempRiver <- tempEfffectInRI %>% 
patrick.lambert's avatar
patrick.lambert committed
140
  filter(season %in% c("spring","summer") & zone == "river_temperature")
141

142
tempInshore <- tempEfffectInRI %>% 
patrick.lambert's avatar
patrick.lambert committed
143
  filter(season %in% c("fall","winter") & zone == "inshore_temperature")
144
145


146
tempEffectInRI <- bind_rows(tempRiver,tempInshore) %>% 
147
148
149
  arrange(basin_name,year,season) %>% 
  select(basin_name,year, season, temperature) %>% 
  group_by(basin_name,season) %>% 
patrick.lambert's avatar
patrick.lambert committed
150
  summarize(temperature_RI = mean(temperature), .groups = 'drop') 
151
152
153
154
155
156

# tempEfffectInRiver %>% 
#   {bind_rows(
#     filter(season == "summer" & zone == "river_temperature"),
#     filter(season == "fall" & zone == "inshore_temperature"))
#     }
patrick.lambert's avatar
patrick.lambert committed
157

Poulet Camille's avatar
Poulet Camille committed
158
159
```

160
161
162
163
164
165
166
167
168
169
170
171
# 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()

```
172

Poulet Camille's avatar
Poulet Camille committed
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190

```{r, include = TRUE, fig.cap = "Growth curve without temperature effect"}

dfGrowth <- data.frame(age = seq(0,19.75,.25), season = c("spring","summer","fall","winter"), basin_name = 'reference') %>% 
  mutate(gender = 'male',
         L = vonBertalanffyGrowth(age, growPar$lengthAtHatching, growPar$linfVonBertForMale, growPar$kOptForMale)) %>% 
  bind_rows(select(.data = ., age, season, basin_name) %>%  
              mutate(gender = 'female',
                     L = vonBertalanffyGrowth(age, growPar$lengthAtHatching, growPar$linfVonBertForFemale, growPar$kOptForFemale)))

dfGrowth %>% filter(basin_name == 'reference') %>% 
  ggplot(aes(x = age)) + 
  geom_line(aes(y = L, linetype = gender, color = basin_name)) +
  labs( x = 'Age (year)', y = 'Fork length (cm)') 


```

191
192
193
194
195
196
197
198
199
200
201
## 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) %>% 
202
  pivot_wider(names_from = parameter, values_from = mean) 
203
204
205
206
207
208
209
210
211
212
213
214
215
216

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).
Poulet Camille's avatar
Poulet Camille committed
217

218
219
220
221
```{r Growth parameters from Stich et al, 2020, echo = FALSE, warning=FALSE}
#------------------------------------------#
#Compute L0 for Stich et al (2020)
#-----------------------------------------#
Poulet Camille's avatar
Poulet Camille committed
222
223
224
225
226
227
228

#function to compute L0 based on growth paramters from Stich et al. 
#K, Linf and t0 were provided for the three metapopulations.
#Linf is in millimeters.

vonBertalanffyLengthAtHatch = function(Linf, K, t0){
  
patrick.lambert's avatar
patrick.lambert committed
229
  L0 = Linf * (1 - exp(t0*K))
Poulet Camille's avatar
Poulet Camille committed
230
231
232
233
  
  return(L0)
}

234
235
#Results
Stich2020_sel<-Stich2020_sel %>% 
patrick.lambert's avatar
patrick.lambert committed
236
  mutate(L0_theo = vonBertalanffyLengthAtHatch(Linf, K, t0)) 
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
#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
Poulet Camille's avatar
Poulet Camille committed
258

259
growthInOffshore <- expand_grid(basin_name = nea_riverBasinFeatures$basin_name , age = seq(0,9.75,.25)) %>%
Poulet Camille's avatar
Poulet Camille committed
260
261
262
263
264
  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'
265
  )) %>% inner_join(nea_riverBasinFeatures %>% select(basin_name, lat_outlet), by ='basin_name') %>% #add latitude to identify the 3 regions 
Poulet Camille's avatar
Poulet Camille committed
266
267
268
269
270
271
272
  select(basin_name, lat_outlet, age, season)%>% 
  mutate(metapop = case_when
         (lat_outlet < 33.8 ~ "semelparous",
           lat_outlet >=33.8 & lat_outlet <= 41.28793 ~ "southern iteroparous",
           lat_outlet > 41.28793 ~ "northern iteroparous")) %>% 
  inner_join(Stich2020_sel, by = c("metapop" = "catchment")) %>% 
  group_by(basin_name,metapop) %>% 
273
  #compute the length at each age based on the growth parameters from Stich et al (2020) provided above.
patrick.lambert's avatar
patrick.lambert committed
274
  mutate(LStich = vonBertalanffyGrowth(age, L0_theo, Linf, K) / 10) %>% # in cm
Poulet Camille's avatar
Poulet Camille committed
275
  ungroup()
276

277
278
279
280
281
#------------------------------------------#
#Add temperature effect for growth curves (2020)
#------------------------------------------#

#Add the temperature time-series for offshore, river, and inshore basins
282
growthInBasin <-growthInOffshore %>% 
patrick.lambert's avatar
patrick.lambert committed
283
284
  inner_join(tempEffectInOffshore, 
             by = c("basin_name" = "inshore_basin_name", "season")) %>% 
285
  right_join(tempEffectInRI, 
patrick.lambert's avatar
patrick.lambert committed
286
287
             by = c("basin_name","season"))

288
#for each season, pick the temperature depending on the location of fish
289
290
growthInBasin <- growthInBasin %>% 
  mutate(temperature_RIO = ifelse(season == "spring" & age == 0, temperature_RI,
patrick.lambert's avatar
patrick.lambert committed
291
292
                              ifelse(season == "summer" & age == 0.25, temperature_RI,
                                     ifelse(season == "fall" & age == 0.50, temperature_RI, temperature_O)))) %>% 
293
  arrange(basin_name,age,season) %>% 
294
  select(basin_name, offshore_basin_name, metapop,age, season, temperature_O, temperature_RI, temperature_RIO, LStich)
295

Poulet Camille's avatar
Poulet Camille committed
296
297
298
#group_by(metapop,age,season) %>% 
#summarize(LStich = unique(LStich))

299
300
301
302
303
304
305
```

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}

patrick.lambert's avatar
patrick.lambert committed
306
307
308
309
310
311
312
# 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')))

313
#rename the basin_name to fit with our list of river_basin name
patrick.lambert's avatar
patrick.lambert committed
314
315
316
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"),
317
         basin_name = replace(basin_name, basin_name == "Upper Chesapeake Bay", "Susquehanna"))
patrick.lambert's avatar
patrick.lambert committed
318
319
320
321
322
323
324
325
326

Stich2020_sel %>% filter(!catchment %in% c('coastwide', 'northern iteroparous', 'semelparous', 'southern iteroparous' )) %>%
  select(catchment, basin_name) %>% 
  left_join(growthInBasin %>% distinct(basin_name, metapop), by = c("basin_name") ) %>% 
  arrange(metapop) %>% 
  flextable() %>% 
  autofit()

```
327
328
329
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.
patrick.lambert's avatar
patrick.lambert committed
330

331
```{r range of temperature, include = TRUE}
patrick.lambert's avatar
patrick.lambert committed
332
333
334
335
336
337
338
339
340
341
342
343
344
# NEW =====
growthInBasin %>% 
  group_by(obs_1900_1950) %>% 
  summarise(Tmin = min(temperature_RIO), Tmax = max(temperature_RIO), .groups = 'drop')

growthInBasin %>% 
  inner_join(Stich2020_sel %>% select(basin_name), by ='basin_name') %>% 
  group_by(obs_1900_1950) %>% 
  summarise(Tmin = min(temperature_RIO), Tmax = max(temperature_RIO), .groups = 'drop')

growthInBasin %>% 
  filter(temperature_RIO < 1.6) %>% 
  distinct(basin_name, obs_1900_1950)
Poulet Camille's avatar
Poulet Camille committed
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376

#----------------------------- NEW GRAPH WITH TEMP EFFECT IN THE 11 RIVERS FROM STICH
growthInBasin_2 = growthInBasin %>% 
  mutate(temp_ref = case_when(
    age == 0.00|age == 0.25 ~ "river", 
    age == 0.50 ~ "inshore",
    TRUE~ as.character("offshore"))) %>% 
  arrange(age,season) 

growthInBasin_2 = growthInBasin_2 %>% 
  mutate(temp_ref_basin = case_when(
   temp_ref =="offshore" & season =="summer" ~ "summering_offshore",
   temp_ref =="offshore" & season =="fall" ~ "summering_offshore",
   temp_ref =="offshore" & season =="winter" ~ "wintering_offshore",
   temp_ref == "offshore" & season =="spring" ~ "wintering_offshore",
   temp_ref == "inshore" ~ "inshore",
   temp_ref =="river"~ "river")) 

growthInBasin_2 %>% 
  inner_join(Stich2020_sel %>% select(basin_name), by ='basin_name') %>% 
  filter(obs_1900_1950 == 1) %>% 
  inner_join(nea_riverBasinFeatures %>% select(basin_name,basin_id,lat_outlet), by = c('basin_name' = 'basin_name')) %>% 
  mutate(basin_name = fct_reorder(basin_name,lat_outlet)) %>% 
  ggplot(aes(x = temperature_RIO, y = basin_name), show.legend = FALSE) +
  geom_path(show.legend = FALSE) + 
  geom_label((aes(label = str_to_title(str_sub(season,1,2)), fill = temp_ref_basin)), colour = "white",size = 2, show.legend = TRUE)+
  scale_fill_manual(values= c("lightblue","blue","orange","purple"), levels(growthInBasin_2$temp_ref_basin))+
  labs(title = '', x = "Average seasonal temperature (°C) ", y = 'River Basins')
  



Poulet Camille's avatar
Poulet Camille committed
377
378
```

379
380
As Stich et al 2020 provided a growth curve for both sex combined, we draw a parameters set for unisex based on the growth parameters for males and females from the XML file.

Poulet Camille's avatar
Poulet Camille committed
381
382
383
384
```{r Growth in GR3D with temperature effect, echo = FALSE, warning=FALSE}

#define parameters for both sex based on value from the XML file 

385
growParUnisex <- growPar %>% 
Poulet Camille's avatar
Poulet Camille committed
386
387
388
389
390
391
  transmute(tempMinGrow,tempMaxGrow,tempOptGrow, sigmaDeltaLVonBert,lengthAtHatching,
            kOpt = mean(c(kOptForFemale,kOptForMale)),
            Linf = mean(c(linfVonBertForFemale,linfVonBertForMale))) 

#merge the temperature effect to compute the relative VBGF 
#growthInOffshore<-growthInOffshore %>% 
392
393
#inner_join(tempEffectInOffshore %>% 
#select(inshore_basin_name,season,temperatureEffect), by =c("basin_name" ="inshore_basin_name", "season" ="season")) %>% 
Poulet Camille's avatar
Poulet Camille committed
394
395
396
397
#select(basin_name,lat_outlet, metapop,age,season,LStich, temperatureEffect)

```

patrick.lambert's avatar
patrick.lambert committed
398
399
400

```{r}
# new ====
401
#save all file  that will be needed for optimization process as SOS data
402
save(growParUnisex,growthInBasin, nea_presence, Stich2020_sel,  file  = 'SOS.rdata' )
patrick.lambert's avatar
patrick.lambert committed
403
404
405
406
407
408

rm(list = ls())
load('SOS.rdata')
source("../GR3D_Rdescription/GR3Dfunction.R")
```

409
Based on growth curves from Stich et al, 2020, we looked for growth parameters (L0, Linf, Kopt) and the thermal range regulating growth, that minimize the square distance between the mean growth curve derived from Stich et al, 2020 and the growth curve to implemented into the GR3D model for each meta-population. 
patrick.lambert's avatar
patrick.lambert committed
410
411

```{r growth optimisation function Camille, echo = FALSE, warning=FALSE}
Poulet Camille's avatar
Poulet Camille committed
412
413
414
415
416
417
418
419
#vector on paramters to get optimized 
#par[1] = Tmin
#par[2] = Topt
#par [3] = Tmax
#par[4] = L0
#par[5] = Linf
#par[6] = k

420
421
#Return data
computeGrowOpt = function(par,data){
Poulet Camille's avatar
Poulet Camille committed
422
423
424
425
426
427
  data <- data %>% arrange(age, basin_name) %>% 
    mutate(L = if_else(age == 0, par[4], 0))
  ages <- data %>% select(age) %>% distinct() %>% arrange(age) %>% unlist(use.names = FALSE)   
  
  
  data <- data %>% 
428
429
430
431
432
433
434
    arrange(basin_name,age, season) %>% 
    group_by(basin_name,age,season) %>% 
    mutate(temperatureEffect = temperatureEffect(temperature, 
                                                 Tmin = par[1],
                                                 Topt = par[2],
                                                 Tmax = par[3])) %>% 
    ungroup()
Poulet Camille's avatar
Poulet Camille committed
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
  
  for (i in 2:length(ages)) {
    previousAge  = ages[i - 1]
    currentAge = ages[i]
    tempEffect = data %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE)
    previousL <- data  %>% filter(age == previousAge) %>% select(L) %>% unlist(use.names = FALSE)
    
    currentL <- vonBertalanffyWithNextIncrement(L = previousL, 
                                                L0 = par[4], 
                                                Linf = par[5], 
                                                K = par[6],
                                                timeStepDuration = currentAge - previousAge, 
                                                sigma = 0, 
                                                tempEffect = tempEffect)
    data = data %>% mutate(L =replace(L, age == ages[i], currentL)) 
    
  }
  
453
454
455
456
457
458
459
460
461
  
  return(data)
}

#computeGrowOpt(vecPar,growthInOffshore)

#Function to be optimized 
computeSSE = function(par,data){
  
patrick.lambert's avatar
patrick.lambert committed
462
463
  data = computeGrowOpt(par,data)
  
Poulet Camille's avatar
Poulet Camille committed
464
  SSE = data %>% 
465
466
467
468
469
470
    group_by(metapop) %>% 
    mutate(W = 1/n()) %>% 
    ungroup() %>% 
    transmute(squareError = W * ((L - LStich)^2)) %>% 
    summarise(sumSquare = sum(squareError)) %>% 
    unlist()
Poulet Camille's avatar
Poulet Camille committed
471
472
473
474
475
  #ungroup() %>% 
  #select(sumSquare) 
  #summarize(sum = sum(sumSquare))
  
  #return(SSE = as.vector(SSE$sum))
476
  return(SSE)
Poulet Camille's avatar
Poulet Camille committed
477
478
}

479
480
#computeSSE (vecPar, growthInOffshore)

patrick.lambert's avatar
patrick.lambert committed
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
```

```{r growth optimisation function Patrick, echo = FALSE, warning=FALSE}
computeGrowAllBasins = function(data, growPar){
  data <- data %>% arrange(age, basin_name) %>% 
    mutate(L = if_else(age == 0, growPar$lengthAtHatching, 0),
           temperatureEffect = temperatureEffect(temperature, growPar$tempMinGrow, growPar$tempOptGrow, growPar$tempMaxGrow))
  ages <- data %>% distinct(age) %>% arrange(age) %>% unlist(use.names = FALSE)   
  
  for (i in 2:length(ages)) {
    previousAge  = ages[i - 1]
    currentAge = ages[i]
    tempEffect = data %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE)
    previousL <- data %>% 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)
    data = data %>% mutate(L =replace(L, age == ages[i], currentL) ) 
  }
  return(data)
}


par2parGrowth = function(par){
  parGrowth = data.frame(tempMinGrow = unname(par['tempOptGrow'] - par['epsilonMinus']),
                         tempOptGrow = unname(par['tempOptGrow']),
                         tempMaxGrow = unname(par['tempOptGrow'] + par['epsilonPlus']),
                         lengthAtHatching = unname(par['lengthAtHatching']),
                         Linf = unname(par['Linf']),
                         kOpt = unname(par['kOpt']), 
                         sigmaDeltaLVonBert = 0 )
  return(parGrowth)
}

520
#  with Topt, Tmin = Topt - espsilonMinus, Tmax = Topt + epsilonPlus 
patrick.lambert's avatar
patrick.lambert committed
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
objFn = function(par, data) {
  growPar = par2parGrowth(par)
  data = computeGrowAllBasins(data, growPar) 
  
  SSE <-  data %>% 
    # fit only on one-season L to avoid smoothing of seasonal growth
    # filter(season == 'spring') %>%
    # mutate(age_year = floor(age)) %>% 
    # group_by(basin_name, metapop, age_year) %>% 
    # summarise( L = mean(L), LStich = mean(LStich), .groups = 'drop') %>% 
    # on pondère chaque basin versant pour garder des points équivalents entre les metapop
    group_by(metapop) %>% mutate(W = 1/n()) %>% ungroup() %>% 
    transmute(squareError = W * ((L - LStich)^2)) %>% summarise(sum(squareError)) %>% unlist(use.names = FALSE)
  return(SSE)
}

537
# same as computeSSE but faster
patrick.lambert's avatar
patrick.lambert committed
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
objFn_B = function(par, data) {
  growPar = enframe(par) %>% 
    pivot_wider() %>% mutate(sigmaDeltaLVonBert = 0)
  data = computeGrowAllBasins(data, growPar) 
  
  SSE <-  data %>% 
    # fit only on one-season L to avoid smoothing of seasonal growth
    # filter(season == 'spring') %>%
    # mutate(age_year = floor(age)) %>% 
    # group_by(basin_name, metapop, age_year) %>% 
    # summarise( L = mean(L), LStich = mean(LStich), .groups = 'drop') %>% 
    # on pondère chaque basin versant pour garder des points équivalents entre les metapop
    group_by(metapop) %>% mutate(W = 1/n()) %>% ungroup() %>% 
    transmute(squareError = W * ((L - LStich)^2)) %>% summarise(sum(squareError)) %>% unlist(use.names = FALSE)
  return(SSE)
}

555
# same as objFn_B with a vector for starting initaila values
patrick.lambert's avatar
patrick.lambert committed
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
objFn_C = function(par, data, fixedPar) {
  growPar = enframe(c(par, fixedPar)) %>% pivot_wider() %>% 
    mutate(sigmaDeltaLVonBert = 0)
  data = computeGrowAllBasins(data, growPar) 
  
  SSE <-  data %>% 
    # fit only on one-season L to avoid smoothing of seasonal growth
    # filter(season == 'spring') %>%
    # mutate(age_year = floor(age)) %>% 
    # group_by(basin_name, metapop, age_year) %>% 
    # summarise( L = mean(L), LStich = mean(LStich), .groups = 'drop') %>% 
    # on pondère chaque basin versant pour garder des points équivalents entre les metapop
    group_by(metapop) %>% mutate(W = 1/n()) %>% ungroup() %>% 
    transmute(squareError = W * ((L - LStich)^2)) %>% summarise(sum(squareError)) %>% unlist(use.names = FALSE)
  return(SSE)
}

```

```{r starting parameters for optimisation}
# starting point from XML
577
vecPar = growParUnisex %>% 
patrick.lambert's avatar
patrick.lambert committed
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
  select(tempMinGrow, tempOptGrow, tempMaxGrow, lengthAtHatching, Linf, kOpt) %>% 
  pivot_longer(tempMinGrow:kOpt, names_to = "parameter", values_to = "value") 
vecPar = as.vector(vecPar$value)

par0 =  c(tempOptGrow = vecPar[2], 
          epsilonMinus = vecPar[2] - vecPar[1], 
          epsilonPlus =  vecPar[3] - vecPar[2], 
          lengthAtHatching = vecPar[4], 
          Linf = vecPar[5], 
          kOpt = vecPar[6])

par0_B = c(tempMinGrow = vecPar[1], 
           tempOptGrow = vecPar[2], 
           tempMaxGrow = vecPar[3], 
           lengthAtHatching = vecPar[4], 
           Linf = vecPar[5], 
           kOpt = vecPar[6])

par0_C = par0_B[c('tempOptGrow', 'lengthAtHatching', 'Linf', 'kOpt')]
#fixedPar_C =  par0_B[c('tempMinGrow', 'tempMaxGrow')]
fixedPar_C =  c(tempMinGrow = 1.6,  tempMaxGrow = 27.9)

# # another starting point
# par0 =  c(tempOptGrow = 7, 
#           epsilonMinus = 5, 
#           epsilonPlus = 10,
#           lengthAtHatching = 5, 
#           Linf = 45, 
#           kOpt = .5)
# 
# par0_B = par2parGrowth(par0) %>% 
#   select(-sigmaDeltaLVonBert) %>% 
#   unlist()
# 
# par0_C = par0_B[c('tempOptGrow', 'lengthAtHatching', 'Linf', 'kOpt')]
# #fixedPar_C =  par0_B[c('tempMinGrow', 'tempMaxGrow')]
# fixedPar_C =  c(tempMinGrow = 1.6,  tempMaxGrow = 27.9)
# 
# vecPar = par2parGrowth(par0) %>% 
#   select(-sigmaDeltaLVonBert) %>% 
#   pivot_longer(tempMinGrow:kOpt, names_to = "parameter",values_to = "value") 
# vecPar = as.vector(vecPar$value)

```


```{r data for calibration}
dataCalibration = growthInBasin %>% 
  inner_join(Stich2020_sel %>% select(basin_name), by ='basin_name') %>% 
  mutate(temperature = temperature_RIO)
```

patrick.lambert's avatar
patrick.lambert committed
630

631
<!-- Avec uniqument le jeu de Stich cela part en vrille avec objFn_B et computeSSE ( sans donner les mémes resultats argh) avec Topt < Tmin. On a des résulats corrects avec objFn ( avec Topt +- epsilon). Mais comme on va utiliser l'option en fixant Tmin et Tmax -->
patrick.lambert's avatar
patrick.lambert committed
632

patrick.lambert's avatar
patrick.lambert committed
633
634
635
636
637
638
639
640
641
642
643
644
645
646
```{r optimisation run}

tic('objFn')
#  avec Topt, Tmin = Topt - espsilonMinus, Tmax = Topt + epsilonPlus 
res <- optim(par = par0, 
             fn = objFn, 
             data = dataCalibration,
             # method = "L-BFGS-B",
             # lower = c(2, 2, 2, 2, 30, 0.2),
             # upper = c(20, 20, 20, 20, 100, 2),
             control = list(trace = 0, 
                            #parscale = c(1,1,1,1,1,0.1),
                            maxit = 1000))
toc()
Poulet Camille's avatar
Poulet Camille committed
647

patrick.lambert's avatar
patrick.lambert committed
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675

tic('objFn_B')
# équivalent de la fonction computeSSE (un peu plus rapide)
res_B <- optim(par = par0_B, 
               fn = objFn_B, 
               data = dataCalibration,
               # method = "L-BFGS-B",
               # lower = c(0,2,10,2,30,0.2),
               # upper = c(5,10,30,20,100,2),
               control = list(trace = 0, 
                              maxit = 1000))
toc()


tic('objFn_C')
# équivalent de la fonction objFn_B mais avec la possibilité de fixer des parametres
res_C <- optim(par = par0_C, 
               fn = objFn_C, 
               data = dataCalibration,
               fixedPar = fixedPar_C,
               # method = "L-BFGS-B",
               # lower = c(0,2,10,2,30,0.2),
               # upper = c(5,10,30,20,100,2),
               control = list(trace = 0, 
                              maxit = 1000))
toc()

tic('computeSSE')
Poulet Camille's avatar
Poulet Camille committed
676
growth_optimal <- optim(c(Tmin = vecPar[1],
677
678
679
680
681
                          Topt = vecPar[2],
                          Tmax = vecPar[3],
                          L0 = vecPar[4],
                          Linf = vecPar[5],
                          K = vecPar[6]),
patrick.lambert's avatar
patrick.lambert committed
682
683
684
685
686
687
688
                        fn = computeSSE, 
                        data = dataCalibration,
                        # lower = c(0,0,0,0,0,0),
                        # method = "L-BFGS-B",
                        control = list(trace = 0, 
                                       #parscale = c(1,1,1,0.1,1,0.01), 
                                       maxit = 10000))
Poulet Camille's avatar
Poulet Camille committed
689
690
691

toc()

patrick.lambert's avatar
patrick.lambert committed
692
693
694
695
696
697
698
699
700
701
702
703
# results
par2parGrowth(res$par)
res_B$par
c(res_C$par, fixedPar_C)
growth_optimal$par

res$value
res_B$value
res_C$value
growth_optimal$value


704
#Optimized parameters 
Poulet Camille's avatar
Poulet Camille committed
705
706
707
708
growParOptim <- data.frame(value = growth_optimal$par) %>% 
  rownames_to_column("parameter")

#Rdata
Poulet Camille's avatar
Poulet Camille committed
709
#save(growth_optimal, res, res_B, res_C, file = "growth_optimal.RData")
Poulet Camille's avatar
Poulet Camille committed
710
711
712
713
714
715

#csv
write_csv(growParOptim, "growParOptim.csv")

```

patrick.lambert's avatar
patrick.lambert committed
716
717
718
719
720

```{r graph to present results}
# NEW 
# ===============================================================================
#parGrowthUsed = c(res_B$par) %>% enframe() %>% pivot_wider()
721
growParOptim = c(res_C$par, fixedPar_C) %>% enframe() %>% pivot_wider()
patrick.lambert's avatar
patrick.lambert committed
722
723
724
# ================================================================================

dataVerif <- computeGrowAllBasins(data = dataCalibration, 
725
726
                                  growPar = growParOptim %>% 
                                  mutate(sigmaDeltaLVonBert = 0))
patrick.lambert's avatar
patrick.lambert committed
727
728
729
730
731
732
733
734
735

# graph by meta pop
dataVerif %>% group_by(metapop, age) %>% 
  summarise(L = mean(L), LStich = mean(LStich), .groups = 'drop') %>% 
  mutate(metapop = factor(metapop, levels = c('northern iteroparous',  'southern iteroparous', 'semelparous'))) %>% 
  ggplot(aes(x = age, color = metapop)) + 
  geom_line(aes(y = LStich)) +
  geom_line(aes(y = L)) 

Poulet Camille's avatar
Poulet Camille committed
736
save(dataVerif, growth_optimal,res, res_B, res_C, file = "growth_optimal.RData")
patrick.lambert's avatar
patrick.lambert committed
737
738
739
740
741

# thermal response curve with average seasonal temperature experienced by fish 
dataTemperatureMetapop <- dataVerif %>% 
  group_by(metapop, season) %>% 
  summarise(temperature = mean(temperature), .groups = 'drop') %>% 
742
  mutate(tempEffect = temperatureEffect(temperature, growParOptim$tempMinGrow, growParOptim$tempOptGrow, growParOptim$tempMaxGrow ),
patrick.lambert's avatar
patrick.lambert committed
743
744
745
746
747
         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)) %>% 
748
  mutate(tempEffect = temperatureEffect(temperature, growParOptim$tempMinGrow, growParOptim$tempOptGrow, growParOptim$tempMaxGrow )) %>% 
patrick.lambert's avatar
patrick.lambert committed
749
750
751
752
753
754
  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) 

```

Poulet Camille's avatar
Poulet Camille committed
755
756
```{r draw growth curves}

Poulet Camille's avatar
Poulet Camille committed
757
#load("growth_optimal.RData")
Poulet Camille's avatar
Poulet Camille committed
758

759
growthCurves = computeGrowOpt(growth_optimal$par, growthInOffshore)
760
save(growthCurves, file = "growthCurves.RData")
Poulet Camille's avatar
Poulet Camille committed
761
762
763
764
765
766
767

#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)") +
768
769
770
  #scale_x_continuous(limits = c(0, 10)) +
  geom_hline(yintercept = 45, color = "pink", linetype =1)+ #female
  geom_hline(yintercept = 40, color = "lightblue", linetype = 1)+
771
772
  geom_vline(aes(age), xintercept = 8, linetype = 1)+
  facet_wrap(.~ metapop, ncol =1)
Poulet Camille's avatar
Poulet Camille committed
773
774
775
776
777

#temperature effect 
temp <- seq(0,45,0.1)

growthCurves %>% 
778
  ggplot(aes(x = temperature)) + 
Poulet Camille's avatar
Poulet Camille committed
779
780
781
782
783
784
785
786
787
788
  geom_line(aes(y = temperatureEffect)) +
  geom_jitter(aes(y = temperatureEffect, col = metapop))


#long fromat for growParOptim
growParOptim <-growParOptim %>% 
  pivot_wider(names_from = parameter, values_from = value)

df <- data.frame(temperature = temp,
                 temperatureEffect = temperatureEffect(temp,
789
790
791
                                                       growParOptim$Tmin,
                                                       growParOptim$Topt,
                                                       growParOptim$Tmax))
Poulet Camille's avatar
Poulet Camille committed
792
793
794
795
796
797
798
799

df %>% 
  ggplot(aes(x = temperature)) + 
  geom_line(aes(y = temperatureEffect)) +
  geom_point(aes(y = temperatureEffect), col = 'red', alpha = 0.4)

```

800
#Growth Per Gender 
Poulet Camille's avatar
Poulet Camille committed
801

802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
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)

```{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)
  
  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
  
  return(unname(SCE1 + SCE2))
}
Poulet Camille's avatar
Poulet Camille committed
820
821
822



823
824
825
826
827
parGrow = data.frame(Linf = 45.4,
                     K = 0.49,
                     L0 = 4.4,
                     LatMaturity = 40 ,
                     diffMeanAgeAtMaturity = 1.5) 
Poulet Camille's avatar
Poulet Camille committed
828

829
830
831
832
833
834
835
836
parGrow = parGrowthUsed %>% 
  select(lengthAtHatching,Linf,kOpt) %>% 
  rename(c("L0" = "lengthAtHatching",
           "K" ="kOpt")) %>% 
  mutate(Lmaturity = 40,
         diffMeanAgeAtMaturity = 0.5)
         
ages = seq_len(12)
Poulet Camille's avatar
Poulet Camille committed
837

838
839
840
objFUN(par = c(Kfemale = 0.45, Kmale = .52),
       ages = ages, 
       parGrow = parGrow)
Poulet Camille's avatar
Poulet Camille committed
841

842
843
844
845
res = optim(par = c(Kfemale = 0.45, Kmale = .52), 
            fn = objFUN,
            ages = ages, 
            parGrow = parGrow)
Poulet Camille's avatar
Poulet Camille committed
846
847


848
849
850
851
852
853
854
855
856
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) 
Poulet Camille's avatar
Poulet Camille committed
857
858


859
860
861
862
863
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))
         
Poulet Camille's avatar
Poulet Camille committed
864

865
```
Poulet Camille's avatar
Poulet Camille committed
866
867