6-Growth_Calibration.Rmd 32.2 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
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
## 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) %>% 
  pivot_wider(names_from = parameter, values_from = mean) %>% 

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
```

```{r Growth in GR3D with temperature effect, echo = FALSE, warning=FALSE}

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

351
growParUnisex <- growPar %>% 
Poulet Camille's avatar
Poulet Camille committed
352
353
354
355
356
357
  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 %>% 
358
359
#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
360
361
362
363
#select(basin_name,lat_outlet, metapop,age,season,LStich, temperatureEffect)

```

patrick.lambert's avatar
patrick.lambert committed
364
365
366

```{r}
# new ====
367
save(growParUnisex,growthInBasin, nea_presence, Stich2020_sel,  file  = 'SOS.rdata' )
patrick.lambert's avatar
patrick.lambert committed
368
369
370
371
372
373
374
375
376

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



```{r growth optimisation function Camille, echo = FALSE, warning=FALSE}
Poulet Camille's avatar
Poulet Camille committed
377
378
379
380
381
382
383
384
#vector on paramters to get optimized 
#par[1] = Tmin
#par[2] = Topt
#par [3] = Tmax
#par[4] = L0
#par[5] = Linf
#par[6] = k

385
386
#Return data
computeGrowOpt = function(par,data){
Poulet Camille's avatar
Poulet Camille committed
387
388
389
390
391
392
  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 %>% 
393
394
395
396
397
398
399
    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
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
  
  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)) 
    
  }
  
418
419
420
421
422
423
424
425
426
  
  return(data)
}

#computeGrowOpt(vecPar,growthInOffshore)

#Function to be optimized 
computeSSE = function(par,data){
  
patrick.lambert's avatar
patrick.lambert committed
427
428
  data = computeGrowOpt(par,data)
  
Poulet Camille's avatar
Poulet Camille committed
429
  SSE = data %>% 
430
431
432
433
434
435
    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
436
437
438
439
440
  #ungroup() %>% 
  #select(sumSquare) 
  #summarize(sum = sum(sumSquare))
  
  #return(SSE = as.vector(SSE$sum))
441
  return(SSE)
Poulet Camille's avatar
Poulet Camille committed
442
443
}

444
445
#computeSSE (vecPar, growthInOffshore)

patrick.lambert's avatar
patrick.lambert committed
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
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
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
```

```{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)
}

#  avec Topt, Tmin = Topt - espsilonMinus, Tmax = Topt + epsilonPlus 
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)
}

# équivalent de la fonction computeSSE (un peu plus rapide)
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)
}

# équivalent de la fonction objFn_B mais avec la possibilité de fixer des parametres
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
vecPar = growParModifed %>% 
  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
595
596
597

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
598
599
600
601
602
603
604
605
606
607
608
609
610
611
```{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
612

patrick.lambert's avatar
patrick.lambert committed
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640

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
641
growth_optimal <- optim(c(Tmin = vecPar[1],
642
643
644
645
646
                          Topt = vecPar[2],
                          Tmax = vecPar[3],
                          L0 = vecPar[4],
                          Linf = vecPar[5],
                          K = vecPar[6]),
patrick.lambert's avatar
patrick.lambert committed
647
648
649
650
651
652
653
                        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
654
655
656

toc()

patrick.lambert's avatar
patrick.lambert committed
657
658
659
660
661
662
663
664
665
666
667
668
# 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


669
#Optimized parameters 
Poulet Camille's avatar
Poulet Camille committed
670
671
672
673
growParOptim <- data.frame(value = growth_optimal$par) %>% 
  rownames_to_column("parameter")

#Rdata
Poulet Camille's avatar
Poulet Camille committed
674
#save(growth_optimal, res, res_B, res_C, file = "growth_optimal.RData")
Poulet Camille's avatar
Poulet Camille committed
675
676
677
678
679
680

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

```

patrick.lambert's avatar
patrick.lambert committed
681
682
683
684
685

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

dataVerif <- computeGrowAllBasins(data = dataCalibration, 
690
691
                                  growPar = growParOptim %>% 
                                  mutate(sigmaDeltaLVonBert = 0))
patrick.lambert's avatar
patrick.lambert committed
692
693
694
695
696
697
698
699
700

# 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
701
save(dataVerif, growth_optimal,res, res_B, res_C, file = "growth_optimal.RData")
patrick.lambert's avatar
patrick.lambert committed
702
703
704
705
706

# thermal response curve with average seasonal temperature experienced by fish 
dataTemperatureMetapop <- dataVerif %>% 
  group_by(metapop, season) %>% 
  summarise(temperature = mean(temperature), .groups = 'drop') %>% 
707
  mutate(tempEffect = temperatureEffect(temperature, growParOptim$tempMinGrow, growParOptim$tempOptGrow, growParOptim$tempMaxGrow ),
patrick.lambert's avatar
patrick.lambert committed
708
709
710
711
712
         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)) %>% 
713
  mutate(tempEffect = temperatureEffect(temperature, growParOptim$tempMinGrow, growParOptim$tempOptGrow, growParOptim$tempMaxGrow )) %>% 
patrick.lambert's avatar
patrick.lambert committed
714
715
716
717
718
719
  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
720
721
```{r draw growth curves}

Poulet Camille's avatar
Poulet Camille committed
722
#load("growth_optimal.RData")
Poulet Camille's avatar
Poulet Camille committed
723

724
growthCurves = computeGrowOpt(growth_optimal$par, growthInOffshore)
725
save(growthCurves, file = "growthCurves.RData")
Poulet Camille's avatar
Poulet Camille committed
726
727
728
729
730
731
732

#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)") +
733
734
735
  #scale_x_continuous(limits = c(0, 10)) +
  geom_hline(yintercept = 45, color = "pink", linetype =1)+ #female
  geom_hline(yintercept = 40, color = "lightblue", linetype = 1)+
736
737
  geom_vline(aes(age), xintercept = 8, linetype = 1)+
  facet_wrap(.~ metapop, ncol =1)
Poulet Camille's avatar
Poulet Camille committed
738
739
740
741
742

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

growthCurves %>% 
743
  ggplot(aes(x = temperature)) + 
Poulet Camille's avatar
Poulet Camille committed
744
745
746
747
748
749
750
751
752
753
  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,
754
755
756
                                                       growParOptim$Tmin,
                                                       growParOptim$Topt,
                                                       growParOptim$Tmax))
Poulet Camille's avatar
Poulet Camille committed
757
758
759
760
761
762
763
764

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

```

765
#Growth Per Gender 
Poulet Camille's avatar
Poulet Camille committed
766

767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
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
785
786
787



788
789
790
791
792
parGrow = data.frame(Linf = 45.4,
                     K = 0.49,
                     L0 = 4.4,
                     LatMaturity = 40 ,
                     diffMeanAgeAtMaturity = 1.5) 
Poulet Camille's avatar
Poulet Camille committed
793

794
795
796
797
798
799
800
801
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
802

803
804
805
objFUN(par = c(Kfemale = 0.45, Kmale = .52),
       ages = ages, 
       parGrow = parGrow)
Poulet Camille's avatar
Poulet Camille committed
806

807
808
809
810
res = optim(par = c(Kfemale = 0.45, Kmale = .52), 
            fn = objFUN,
            ages = ages, 
            parGrow = parGrow)
Poulet Camille's avatar
Poulet Camille committed
811
812


813
814
815
816
817
818
819
820
821
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
822
823


824
825
826
827
828
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
829

830
```
Poulet Camille's avatar
Poulet Camille committed
831
832