GR3Dfunction.R 13.5 KB
Newer Older
1
2
3
#----------------------------------------------
# Temperature effect
#----------------------------------------------
4
5
# temperature effect in GR3D
# from Rosso
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

temperatureEffect = function(tempWater, Tmin, Topt, Tmax, effectFunction = Rosso){
  return(effectFunction(tempWater, Tmin, Topt, Tmax))
}



modifiedLehman = function(temperature, Tmin, Topt, Tmax){
  # Svirezhev, Yu. M., Krysanova, V. P. & Voinov, A. A. (1984) Mathematical modelling of a fish pond ecosystem. 
  #   Ecological Modelling 21, 315–337.
  
  # rate(Tmin or Tmax) = exp(-4.6) = 1 %
  
  rate = rep(0, length(temperature))
  rate[temperature < Topt] = exp(-4.6 * ((Topt - temperature[temperature < Topt]) / (Topt - Tmin))^4)
  rate[temperature >= Topt] = exp(-4.6 * ((temperature[temperature >= Topt] - Topt) / (Tmax - Topt))^4)
  
  return(rate)
}

Rosso = function(tempWater, Tmin, Topt, Tmax){
27
28
29
30
31
32
  response = (tempWater - Tmin) * (tempWater - Tmax) / ((tempWater - Tmin) * (tempWater - Tmax) - (tempWater - Topt)^2)
  
  response[tempWater <= Tmin | tempWater >= Tmax] = 0
  return(response)
}

33
34
thermalRange = function(pct = 0.8, Tmin, Topt, Tmax,...){
  lower = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax,...) - pct,
patrick.lambert's avatar
patrick.lambert committed
35
           interval = c(Tmin, Topt))$root
36
  upper = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax,...) - pct,
patrick.lambert's avatar
patrick.lambert committed
37
38
39
40
                  interval = c(Topt, Tmax))$root
  return(c(lower = lower,upper = upper))
}

41
# thermalRange(Tmin = 3, Topt = 17, Tmax = 27, effectFunction = modifiedLehman )
42

43
44
45
#----------------------------------------------
# Growth simulation
#---------------------------------------------
46
47
48
49
50
51
52
53
54
55
56
57
58
# von Bertalaffy : L = f(age)
vonBertalanffyGrowth = function(age, L0, Linf, K){
  t0 = log(1 - L0 / Linf) / K
  return(Linf * (1 - exp(-K * (age - t0))))
}

# von Bertalanffy inverse  : age = f(L)
# with L0 rather the usual t0 parameter
vonBertalanffyInverse = function(L, L0, Linf, K){
  t0 = log(1 - L0/Linf)/K
  return(t0 - log(1 - L/Linf)/K)
}

Poulet Camille's avatar
Poulet Camille committed
59

60
61
# von Bertalanffy increment
# pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche
patrick.lambert's avatar
patrick.lambert committed
62

patrick.lambert's avatar
patrick.lambert committed
63
vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){
64
65
66
67
  tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26)
  L = matrix(nrow = nStep + 1)
  L[1] = L0
  for (i in 1:nStep) {
patrick.lambert's avatar
patrick.lambert committed
68
    mu = log((Linf - L[i]) * (1 - exp(-K * timeStepDuration))) - sigma * sigma / 2
69
70
71
72
73
74
75
76
77
    increment = exp(rnorm(1, mu, sigma))
    if (withTempEffect) {
      increment = increment * tempEffect[((i - 1) %% 4) + 1]
    }
    L[i + 1] = L[i] + increment
  }
  return(L)
}

patrick.lambert's avatar
patrick.lambert committed
78
vonBertalanffyWithNextIncrement = function(L,  Linf, K, timeStepDuration, sigma, tempEffect ){
79
  if (sigma == 0) {
patrick.lambert's avatar
patrick.lambert committed
80
    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
81
82
83
    increment = exp(mu)
  }
  else {
patrick.lambert's avatar
patrick.lambert committed
84
    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
85
    increment = exp(rnorm(length(mu), mu, sigma))
86
87
88
89
90
  }
  L = L + increment * tempEffect
  return(L)
}

patrick.lambert's avatar
patrick.lambert committed
91
# generate a cohort of length trajectories
92
computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar, ...){
patrick.lambert's avatar
patrick.lambert committed
93
94
95
96
97
98
99
100
  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, 
101
                                                 growPar$tempMaxGrow, ...), 
patrick.lambert's avatar
patrick.lambert committed
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
           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, 
                                                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, 
124
                                                           RNGseed =1, ...){
patrick.lambert's avatar
patrick.lambert committed
125
126
127
128
129
130
131
132
133
134
135
136
  set.seed(RNGseed)
  ages = temperaturePattern %>% 
    distinct(age) %>% 
    unlist(use.names = FALSE)
  
  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, 
137
                                                 growPar$tempMaxGrow, ...), 
patrick.lambert's avatar
patrick.lambert committed
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
           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, 
                                               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)
}

#computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 10, 
#growPar = growParUnisex,
# RNGseed = 1)

vonBertalanffyWithRandomVector = function(L, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
  if (sigma == 0) {
    mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))),-Inf)
    #mu =  log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration)))
    increment = exp(mu)
  }
  else {
    mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2, -Inf)
    #      mu = log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2
    increment = exp(randomVector * sigma + mu)
  }
  L = pmin(Linf, L + increment) 
  return(L)
}

patrick.lambert's avatar
patrick.lambert committed
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
computeOgive = function(lengthTrajectories, lengthAtMaturity){
  ogive <- lengthTrajectories %>%
    group_by(age) %>%
    summarise(nTotal = n()) %>%
    left_join(lengthTrajectories %>%  group_by(age) %>%
                filter(L >= lengthAtMaturity) %>%
                summarise(taller = n(), .groups = 'drop'),
              by = 'age') %>%
    replace_na(list(taller = 0)) %>%
    mutate(mature = c(0,diff(taller)), 
           immature = nTotal - cumsum(mature),
           p = if_else(mature + immature > 0 , mature / (mature + immature), 1)) %>% 
    select(age, p)
  return(ogive)
}
patrick.lambert's avatar
patrick.lambert committed
195

patrick.lambert's avatar
patrick.lambert committed
196
197
# deprecated
computeOgive3 = function(lengthTrajectories, lengthAtMaturity) {
patrick.lambert's avatar
patrick.lambert committed
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
  ogive <-
    lengthTrajectories %>% 
    group_by(age) %>% 
    summarise(nTotal=n()) %>% 
    left_join(lengthTrajectories %>%  group_by(age) %>% 
                filter(L > lengthAtMaturity) %>% 
                summarise(mature = n(), .groups = 'drop'), 
              by = 'age') %>% 
    replace_na(list(mature = 0)) %>% 
    mutate(p =  mature/nTotal) %>% 
    select(age, p)
  
  return(ogive)
} 

213
214
215
216
217
218
219
220
221
222
223
# -------------------------------------------------------
# Dispersal      
# see see (Chapman et al., 2007; Holloway et al., 2016)
# -------------------------------------------------------
logitKernel =  function(distance, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance,
                        basinSurface = 0, alpha2 = 0, meanBasinSurface = 0, standardDeviationBasinSurface = 1,
                        fishLength = 0, alpha3 = 0, meanLengthAtRepro = 0, standardDeviationLengthAtRepro = 1) {
  logitW =  alpha0 - alpha1 * (distance - meanInterDistance) / standardDeviationInterDistance +
    alpha2 * (basinSurface - meanBasinSurface) / standardDeviationBasinSurface +
    alpha3 * (fishLength - meanLengthAtRepro) / standardDeviationInterDistance;
  
224
225
226
227
228
229
230
231
232
233
234
235
236
  W = 1 / (1+ exp(-logitW))
  return(W)
}
# extended negative exponential kernel
eneKernel =  function(distance, alpha_D, beta_D){
  W = exp(-alpha_D * distance ^ beta_D)
  return(W)
}

strayerSurvival = function(distance, m ){
  # m: extra mortality rate for strayers (L-1)
  return(exp(-m * distance))
}
237
238
239
240
241

#-------------------------------------------------------
# Spwaner survival before reproduction 
# Dome-shape curve with temperature effect
#-------------------------------------------------------
242
spawnerSurvivalPreReproduction <- function(Triver,Tmin, Topt, Tmax, parSurv){
243
244
  
  #River survival
245
  survProbOptRiver =  parSurv$survProbOptRiver
246
  
Poulet Camille's avatar
Poulet Camille committed
247
  spRiver = survProbOptRiver * temperatureEffect(Triver, Tmin, Topt, Tmax)
248
249
250
251
252
253
254
255
256
257
258
259
  
  #StockAfterSurv = S * SpRiver
  
  return(spRiver)
  
}

#----------------------------------------------
# Stock Recruitement relationship
#See Rougier et al, 2014; 2015 and Jatteau et al., 2017
#----------------------------------------------
#survival of larvae 14dph from Jatteau et al, 2017
260
261
262
263
264
265
266
267
268
269
stockRecruitementRelationship <- function(Triver,Tmin, Topt, Tmax, survivalStock, parRep, parSurv, surfaceWatershed = 84810) {
  
  #surfaceWatershed = 84810 Surface de la Garonne 
  lambda = parRep$lambda
  deltaT = parRep$deltaT
  survOptRep =  parRep$survOptRep
  eta = parRep$eta #simule l'effet allee
  ratioS95_S50 = parRep$ratioS95_S50
  alpha =  parRep$a #species fecundity 
  survProbOptRiver =  parSurv$survProbOptRiver
270
271
272
273
  periodAtSea = 5 - deltaT
  
  ###--------------------- SR relationship -----------------
  
Poulet Camille's avatar
Poulet Camille committed
274
  #parametre c de la RS de BH int?grant un effet du BV consid?r? 
275
276
  cj = lambda/surfaceWatershed
  
Poulet Camille's avatar
Poulet Camille committed
277
  #parametre b repr?sentant la mortalit? densit? d?pendante de la RS de BH int?grant un effet de la temperature
278
279
280
281
282
  # bj = (-(1/deltaTrecruitement))*
  #   log(survOptRep * temperatureEffect(temp, 9.8, 20.0, 26.0))
  
  bj = -log(survOptRep * temperatureEffect(Triver, Tmin , Topt , Tmax)) / deltaT
  
Poulet Camille's avatar
Poulet Camille committed
283
  #parametre a (f?condit? de l'esp?ce) de la RS de BH int?grant un effet de la temperature
284
285
286
  alphaj = (bj * exp(-bj * deltaT)) / (cj * (1 - exp(-bj * deltaT)))
  alphaj[is.na(alphaj)] <- 0
  
Poulet Camille's avatar
Poulet Camille committed
287
  #Bj = param?tre de la relation SR int?grant l'effet de la temp?rature 
288
289
  betaj = bj/(alpha*cj*(1 - exp(-bj*deltaT)))
  
Poulet Camille's avatar
Poulet Camille committed
290
  #p = proportion de g?niteurs participant ? la reproduction en focntion de la quantit? de g?niteur total
291
292
  #p = 1/(1+exp(-log(19)*(S-n)/(Ratio*surfaceWatershed)))
  
293
  S95 = eta * surfaceWatershed
294
295
296
297
  S50 = S95/ratioS95_S50
  
  p = 1/(1 + exp(-log(19)*(survivalStock - S50)/(S95 - S50)))
  
Poulet Camille's avatar
Poulet Camille committed
298
  #relation Stock Recrutement ie calcul le nombre de recrues en fonction du nombre de g?niteurs et de la T en int?grant l'effet Allee 
299
300
301
  
  #R0 = aj * S * p 
  
302
  alleeEffect = 1/(1 + exp(-log(19)*(survivalStock - eta/ratioS95_S50*surfaceWatershed)/(eta*surfaceWatershed - eta/ratioS95_S50*surfaceWatershed)))
303
304
  
  
Poulet Camille's avatar
Poulet Camille committed
305
  Rj = (alphaj * survivalStock * alleeEffect)/(betaj + survivalStock * alleeEffect)
306
307
308
309
310
311
312
313
314
315
316
317
  
  #Rj = ((aj * S) * p)/(Bj +S * p)
  
  stockRecruitement = as.vector(Rj)
  
  return(stockRecruitement)
}


# ----------------------------------------------
# Spawner survival after reproduction 
#----------------------------------------------
318
319
320
321
322
323
324
325
326
327
328
#Logit with temperature effect depending on a Triver 
spawnerSurvivalPostReproductionTempRiver <- function(Triver, coeffa, coeffb){
  
  spRiverPostSpawn = (coeffb/ (1 + exp(coeffa*(Triver))))
  
  return(spRiverPostSpawn)
  
}

#Logit with temperature effect depending on the difference between the river temperature and Tref
spawnerSurvivalPostReproductionTempRef <- function(Triver,Tref, coeffa, coeffb){
329
330
331
332
333
334
335
  
  spRiverPostSpawn = (coeffb/ (1 + exp(coeffa*(Triver-Tref))))
  
  return(spRiverPostSpawn)
  
}

336
337
338
339
340
341
342
#Logit for survival after reproduction with log19
#similar to spawnerSurvivalPostReproductionWithTempRef which include a

logit2 = function(Triver, Tref,minTempForIteroparity){
  return( 1/ (1+exp((log(19)/(minTempForIteroparity-Tref))*(Triver-Tref))))
}

343
#Dome-shape curve with temperature effect 
344
spawnerSurvivalPostReproductionWithBellCurve <- function(Triver, Tmin, Topt, Tmax, coeffb){
345
346
347
348
349
350
351
352
353
354
355
  
  #P1: 
  #SpRiverPostSpawn = probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit*(TemperatureEffect(Triver, Tmin, Topt, Tmax))))
  
  #P2: 
  #SpRiverPostSpawn = probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit))
  #SpRiverPostSpawnTempEffect = (probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit))) * TemperatureEffect(Triver, Tmin, Topt, Tmax)
  
  #P3: 
  spRiverPostSpawn = coeffb * temperatureEffect(Triver, Tmin, Topt, Tmax)
  
356
  #stockAfterSpawn = S * spRiverPostSpawn 
357
358
359
360
361
362
  
  return(spRiverPostSpawn)
}



363