diff --git a/exploration/NEA_calibration_offline/1-Calibration_Thermal_Curves.Rmd b/exploration/NEA_calibration_offline/1-Calibration_Thermal_Curves.Rmd index 9ca6ee700e0d0bbf939a322b86ccea45a11fba70..0a4f11a91190334762fa8118d43f11a068c17414 100644 --- a/exploration/NEA_calibration_offline/1-Calibration_Thermal_Curves.Rmd +++ b/exploration/NEA_calibration_offline/1-Calibration_Thermal_Curves.Rmd @@ -31,7 +31,7 @@ library(tidymv) ``` -```{r upload data, warning=FALSE} +```{r upload data, warning= FALSE} source("../GR3D_Rdescription/GR3Dfunction.R") source("../GR3D_Rdescription/GR3D_NEA_env_data.R") @@ -86,7 +86,7 @@ rectTemperature = function(Tmin, Topt, Tmax, pct = 0.8, col ='red'){ ``` -```{r Observations of presence/absence according to temperature, include = TRUE, fig.cap= 'Temperature experienced by spawners during their river stage in spring'' } +```{r Observations of presence/absence according to temperature, include = TRUE, fig.cap= 'Temperature experienced by spawners during their river stage in spring' } #plot of temperature ranges obs_spg <- nea_presence_temp %>% @@ -149,7 +149,12 @@ rectTempRange = rectTemperature(Tmin = survivePar$TminSurv, #------- Compute spawners survival in river before reproduction --------# -SpawnerSurvivalPreRep <- spawnerSurvivalPreReproduction(nea_presence_temp$spring_river_temperature, survivePar$TminSurv,survivePar$ToptSurv, survivePar$TmaxSurv) +SpawnerSurvivalPreRep <- spawnerSurvivalPreReproduction( + Triver = nea_presence_temp$spring_river_temperature, + Tmin = survivePar$TminSurv, + Topt = survivePar$ToptSurv, + Tmax = survivePar$TmaxSurv, + parSurv = survivePar) #------- Plot results from GAM with ggplot------------------------------# #with ggplot @@ -162,7 +167,7 @@ gam_gp_spg <- plotGAM(gam_spg, smooth.cov = "spring_river_temperature")+ natural_pred <-pred_spg %>% ggplot(aes(spring_river_temperature))+ geom_line(aes(y = natural_fit))+ - geom_line(aes(y = natural_fit_80_pct), col = 'red', alpha = 0.3)+ + #geom_line(aes(y = natural_fit_80_pct), col = 'red', alpha = 0.3)+ theme_bw()+ annotate("rect", xmin = rectTempRange$xmin, xmax = rectTempRange$xmax, ymin = rectTempRange$ymin, ymax = rectTempRange$ymax, alpha = rectTempRange$alpha, fill = rectTempRange$col) @@ -189,7 +194,7 @@ plot_grid(obs_spg, gam_gp_spg, natural_pred, labels=c("A","B","C"), ncol = 3, nr ``` -```{r optimization, include = TRUE } +```{r optimization - OBSOLETE, include = TRUE } #Found Tmin and Tmax to have a survival of 80% accross these ranges @@ -223,87 +228,287 @@ fixedPar_4 = c(tempMinSurv = 8, tempMaxSurv = 26) spSurvival_df = nea_presence_temp -spSurvival_df %>% +survivePar_0 = spSurvival_df %>% group_by(obs_1900_1950) %>% summarise(Tmin = min(spring_river_temperature), Tmax = max(spring_river_temperature), .groups='drop') %>% - #filter(obs_1900_1950 == 1) %>% + filter(obs_1900_1950 == 1) + +objectiveFunction = function(par, data, fixedPar){ - objectiveFunction = function(par, data, fixedPar){ - - df <- data %>% - group_by(basin_id, basin_name,year) %>% - mutate(spSurvival = spawnerSurvivalPreReproduction(spring_river_temperature, - Tmin = fixedPar['tempMinSurv'], - Topt = par ['tempOptSurv'], - Tmax =fixedPar['tempMaxSurv'])) - tempRange <- df %>% - group_by(obs_1900_1950) %>% - summarise(Tmin = min(spring_river_temperature), - Tmax = max(spring_river_temperature), .groups='drop') %>% - filter(obs_1900_1950 == 1) - - - df<-df %>% - filter(obs_1900_1950 == 1) %>% - group_by(basin_id,basin_name) %>% - mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.80, fifelse (spring_river_temperature == tempRange$Tmax, 0.80, - spSurvival))) %>% - select(basin_id,basin_name,spring_river_temperature,spSurvival,spSurvival_opt,obs_1900_1950) %>% - arrange(spring_river_temperature,spSurvival) - #mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.8 * spSurvival, fifelse (spring_river_temperature == tempRange$Tmax, 0.8 * spSurvival, - #spSurvival))) - - - SSE = df %>% - group_by(basin_id,basin_name) %>% - filter(spSurvival != spSurvival_opt) %>% - mutate(SSE = (spSurvival - spSurvival_opt)^2) %>% - ungroup() %>% - summarise(sum(SSE), .groups = 'drop') %>% - unlist(use.names = FALSE) - - return(SSE) - } + df <- data %>% + group_by(basin_id, basin_name,year) %>% + mutate(spSurvival = spawnerSurvivalPreReproduction(spring_river_temperature, + Tmin = fixedPar['tempMinSurv'], + Topt = par ['tempOptSurv'], + Tmax =fixedPar['tempMaxSurv'], + parSurv = sruvivePar)) + tempRange <- df %>% + group_by(obs_1900_1950) %>% + summarise(Tmin = min(spring_river_temperature), + Tmax = max(spring_river_temperature), .groups='drop') %>% + filter(obs_1900_1950 == 1) + + + df<-df %>% + filter(obs_1900_1950 == 1) %>% + group_by(basin_id,basin_name) %>% + mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.80, fifelse (spring_river_temperature == tempRange$Tmax, 0.80, + spSurvival))) %>% + select(basin_id,basin_name,spring_river_temperature,spSurvival,spSurvival_opt,obs_1900_1950) %>% + arrange(spring_river_temperature,spSurvival) + #mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.8 * spSurvival, fifelse (spring_river_temperature == tempRange$Tmax, 0.8 * spSurvival, + #spSurvival))) + + + SSE = df %>% + group_by(basin_id,basin_name) %>% + filter(spSurvival != spSurvival_opt) %>% + mutate(SSE = (spSurvival - spSurvival_opt)^2) %>% + ungroup() %>% + summarise(sum(SSE), .groups = 'drop') %>% + unlist(use.names = FALSE) + + return(SSE) +} -#objectiveFunction(par_0,spSurvival_df) -res <- optim(par = par_1, - fn = objectiveFunction, - data = spSurvival_df) -#method = "L-BFGS-B", -#lower = c(-10,14,26), -#upper = c(10,21,30)) +``` -res_3 <- optim(par = par_ini_3, - fn = objectiveFunction, - data = spSurvival_df, - fixedPar =fixedPar_3, - #method = "L-BFGS-B", - #lower = c(8), - upper = c(21)) -res_4 <- optim(par = par_ini_3, - fn = objectiveFunction, - data = spSurvival_df, - fixedPar =fixedPar_4) -#method = "L-BFGS-B", -#lower = c(8), +```{r optimization, include = TRUE } -#--------------------------------------# -#----------- Plots curves --------------# -#---------------------------------------# +#-----------------------------------# +#Initial values to be optimised over# +#-----------------------------------# -spSurvival_df <- spSurvival_df %>% +#Case 1: optimisation on 3 parameters (TminSurv, ToptSurv, TmaxSurv) +#with starting values for Tmin and Tmax from A.alosa +par_0 = c(TminSurv = survivePar$TminSurv, + #ToptSurv = survivePar$ToptSurv, + TmaxSurv = survivePar$TmaxSurv) + + +c(TminSurv = survivePar$TminSurv, + #ToptSurv = survivePar$ToptSurv, + TmaxSurv = survivePar$TmaxSurv) + +#Case 1B: optimisation on 3 parameters (TminSurv, ToptSurv, TmaxSurv) +#with starting values for Tmin and Tmax corresponding to temeperatures ranges experienced by fish in river in spring over 1900-1950 + +tempRange <- nea_presence_temp %>% + group_by(obs_1900_1950) %>% + summarise(Tmin = min(spring_river_temperature), + Tmax = max(spring_river_temperature), .groups='drop') %>% + filter(obs_1900_1950 == 1) + + +par_1 = c(TminSurv = tempRange$Tmin, + ToptSurv = survivePar$ToptSurv, + TmaxSurv = tempRange$Tmax) + + + +#Case 2: optimisation on ToptSurv +#with starting values for Tmin and Tmax corresponding to temeperatures ranges experienced by fish in river in spring over 1900-1950 + +par_2 = par_0 ['ToptSurv'] +fixedPar_2 = par_1[c('TminSurv', 'TmaxSurv')] + + +#Case 3: optimisation on TminSurv and TmaxSurv +#with starting values for ToptSurv corresponding to A.alosa +par_3 = par_1[c('TminSurv', 'TmaxSurv')] +fixedPar_3 = par_0[c('ToptSurv')] + +#-----------------------------------# +#Ojectives function# +#-----------------------------------# + +survivalObs = tibble(temperature = c(4.15, 25.3), + survival = c(0.5, 0.5)) + + +objectiveFunction_A = function(par, survObs, lambda = 2/3){ + + Tmin = par['TminSurv'] + Tmax = par['TmaxSurv'] + Topt = Tmin + lambda *(Tmax - Tmin) + + SSE = sum((temperatureEffect(tempWater = survivalObs$temperature, + Tmin = Tmin, + Topt = Topt, + Tmax = Tmax) - survObs$survival)^2) + + # SSE_inf = (temperatureEffect(tempWater = survivalObs$temperature[1], + # Tmin = Tmin, + # Topt = Topt, + # Tmax = Tmax) - survObs$survival[1])^2 + # + # SSE_sup <- (temperatureEffect(tempWater = survivalObs$temperature[2], + # Tmin = Tmin, + # Topt = Topt, + # Tmax = Tmax) - survObs$survival[2])^2 + # + # SSE_tot <- unname(SSE_inf + SSE_sup) + + return(SSE) +} + + +objectiveFunction_A(par = par_0, survObs = survivalObs, lambda = 2/3) + +lambda = 0.75 + +res_0A <- optim(par = c(TminSurv = 0, + TmaxSurv = 40), + fn = objectiveFunction_A, + survObs = survivalObs, + lambda = lambda, + control = list(trace = 1)) + + +temperature = seq(-10,40,0.1) +TminSurv = res_0A$par['TminSurv'] +TmaxSurv = res_0A$par['TmaxSurv'] +ToptSurv = TminSurv + lambda *(TmaxSurv - TminSurv) + + +ggplot() + + geom_line(aes(temperature, temperatureEffect (temperature, res_0A$par['TminSurv'], + Topt = ToptSurv, + res_0A$par['TmaxSurv']))) + + #coord_cartesian(ylim(0,1))+ + #scale_y_continuous(limits= c(0,1))+ + geom_vline(xintercept = 4.15)+ + geom_vline(xintercept = 25.3)+ + geom_hline(yintercept = 0.8) + + + + + +res_1A <- optim(par = par_1, + fn = objectiveFunction_A, + tempMinGAM = 4.15, + tempMaxGAM = 22) + + +objectiveFunction_Test = function(par, Tmin = 0, survObs){ + + Topt = par['ToptSurv'] + Tmax = par['TmaxSurv'] + Tmin = Tmin + + SSE = sum((temperatureEffect(tempWater = survivalObs$temperature, + Tmin = Tmin, + Topt = Topt, + Tmax = Tmax) - survObs$survival)^2) + + # SSE_inf = (temperatureEffect(tempWater = survivalObs$temperature[1], + # Tmin = Tmin, + # Topt = Topt, + # Tmax = Tmax) - survObs$survival[1])^2 + # + # SSE_sup <- (temperatureEffect(tempWater = survivalObs$temperature[2], + # Tmin = Tmin, + # Topt = Topt, + # Tmax = Tmax) - survObs$survival[2])^2 + # + # SSE_tot <- unname(SSE_inf + SSE_sup) + + return(SSE) +} + +Tmin = 0 + +res_test <- optim(par = c(ToptSurv = 17, + TmaxSurv = 26), + fn = objectiveFunction_Test, + Tmin = Tmin, + survObs = survivalObs, + control = list(trace = 1)) + + +ggplot() + + geom_line(aes(temperature, temperatureEffect(temperature, + Topt = res_test$par['ToptSurv'], + Tmax = res_test$par['TmaxSurv'], + Tmin = Tmin))) + + #coord_cartesian(ylim(0,1))+ + #scale_y_continuous(limits= c(0,1))+ + geom_vline(xintercept = 4.15)+ + geom_vline(xintercept = 25.3)+ + geom_hline(yintercept = 0.8) + + + + + + + +objectiveFunction_B = function(par, fixedPar){ + + SSE_inf = (temperatureEffect(tempWater = fixedPar['TminSurv'], + Tmin = fixedPar['TminSurv'], + Topt = par['ToptSurv'], + Tmax = fixedPar['TmaxSurv'])-0.8)^2 + + SSE_sup <- (temperatureEffect(tempWater = fixedPar['TmaxSurv'], + Tmin = fixedPar['TminSurv'], + Topt = par['ToptSurv'], + Tmax = fixedPar['TmaxSurv']) - 0.8)^2 + + SSE_tot <- unname(SSE_inf + SSE_sup) + + return(as.vector(SSE_tot)) +} + +res_2B <- optim(par = par_2, + fn = objectiveFunction_B, + fixedPar = fixedPar_2) + + + +objectiveFunction_C = function(par, tempMinGAM,tempMaxGAM){ + + SSE_inf = (temperatureEffect(tempWater = tempMinGAM, + Tmin = par['TminSurv'], + Topt = par['ToptSurv'], + Tmax = fixedPar['TmaxSurv'])-0.8)^2 + + SSE_sup <- (temperatureEffect(tempWater = tempMaxGAM, + Tmin = par['TminSurv'], + Topt = par['ToptSurv'], + Tmax = par['TmaxSurv']) - 0.8)^2 + + ToptObs = par['ToptSurv'] + ToptTar = 2/3* (par['TmaxSurv'] - par['TminSurv']) + SSE_opt <- (ToptObs - ToptTar)^2 + + + SSE_tot <- unname(SSE_inf + SSE_sup +SSE_opt) + + return(as.vector(SSE_tot)) +} + +res_2C <- optim(par = par_2, + fn = objectiveFunction_C, + tempMinGAM = 4.15, + tempMaxGAM = 25.3) + +spSurvival_df <- nea_presence_temp %>% group_by(basin_id, basin_name,year) %>% mutate(spSurvivalOpt = spawnerSurvivalPreReproduction(spring_river_temperature, - Tmin = res$par['tempMinSurv'], - Topt = res$par['tempOptSurv'], - Tmax = res$par['tempMaxSurv']), + Tmin = res_1A$par['TminSurv'], + Topt = res_1A$par['ToptSurv'], + Tmax = res_1A$par['TmaxSurv'], + parSurv = survivePar), spSurvival = spawnerSurvivalPreReproduction(spring_river_temperature, Tmin = survivePar$TminSurv, Topt = survivePar$ToptSurv, - Tmax = survivePar$TmaxSurv)) %>% + Tmax = survivePar$TmaxSurv, + parSurv = survivePar)) %>% select(basin_id,basin_name,year,spring_river_temperature,obs_1900_1950,spSurvival, spSurvivalOpt) @@ -312,10 +517,6 @@ spSurvival_df %>% geom_line(aes(y= spSurvival), col = "orange", alpha = 0.3) + geom_line(aes(y= spSurvivalOpt), col = "red", alpha = 0.3)+ theme_bw() - -#------------------------ - - ``` @@ -324,9 +525,6 @@ spSurvival_df %>% After hatching, juveniles spend their first summer in freshwater before emigrating seaward, in the inshore basin to spend fall. Considering that spawner reproduce in spring and a larvae survival 14dph, we used the spring riverine temperature for thermal range optimisation. - -TO COMPLETE BY THERMAL RANGE - ```{r optimisation with gam } S = 100 #number of spawners arriving at spawning grounds @@ -400,7 +598,7 @@ spSurvival_df %>% Tmax = max(spring_river_temperature), .groups='drop') %>% #filter(obs_1900_1950 == 1) %>% -objectiveFunctionLarv = function(par, data, fixedPar){ + objectiveFunctionLarv = function(par, data, fixedPar){ df <- data %>% group_by(basin_id, basin_name,year) %>% @@ -409,35 +607,35 @@ objectiveFunctionLarv = function(par, data, fixedPar){ Topt = par ['tempOptRep'], Tmax = par['tempMaxRep'], survivalStock = spawnerSurvivalPreReproduction(spring_river_temperature, survivePar$TminSurv,survivePar$ToptSurv, survivePar$TmaxSurv))) - - -tempRange <- df %>% - group_by(obs_1900_1950) %>% - summarise(Tmin = min(spring_river_temperature), - Tmax = max(spring_river_temperature), .groups='drop') %>% - filter(obs_1900_1950 == 1) - - -df<-df %>% - filter(obs_1900_1950 == 1) %>% - group_by(basin_id,basin_name) %>% - mutate(spLarvSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.80, fifelse (spring_river_temperature == tempRange$Tmax, 0.80, - spLarvSurvival))) %>% - select(basin_id,basin_name,spring_river_temperature,spLarvSurvival,spLarvSurvival_opt,obs_1900_1950) %>% - arrange(spring_river_temperature,spLarvSurvival) -#mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.8 * spSurvival, fifelse (spring_river_temperature == tempRange$Tmax, 0.8 * spSurvival, - #spSurvival))) - - -SSE = df %>% - group_by(basin_id,basin_name) %>% - filter(spLarvSurvival != spLarvSurvival_opt) %>% - mutate(SSE = (spLarvSurvival - spLarvSurvival_opt)^2) %>% - ungroup() %>% - summarise(sum(SSE), .groups = 'drop') %>% - unlist(use.names = FALSE) - -return(SSE) + + + tempRange <- df %>% + group_by(obs_1900_1950) %>% + summarise(Tmin = min(spring_river_temperature), + Tmax = max(spring_river_temperature), .groups='drop') %>% + filter(obs_1900_1950 == 1) + + + df<-df %>% + filter(obs_1900_1950 == 1) %>% + group_by(basin_id,basin_name) %>% + mutate(spLarvSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.80, fifelse (spring_river_temperature == tempRange$Tmax, 0.80, + spLarvSurvival))) %>% + select(basin_id,basin_name,spring_river_temperature,spLarvSurvival,spLarvSurvival_opt,obs_1900_1950) %>% + arrange(spring_river_temperature,spLarvSurvival) + #mutate(spSurvival_opt = fifelse(spring_river_temperature == tempRange$Tmin, 0.8 * spSurvival, fifelse (spring_river_temperature == tempRange$Tmax, 0.8 * spSurvival, + #spSurvival))) + + + SSE = df %>% + group_by(basin_id,basin_name) %>% + filter(spLarvSurvival != spLarvSurvival_opt) %>% + mutate(SSE = (spLarvSurvival - spLarvSurvival_opt)^2) %>% + ungroup() %>% + summarise(sum(SSE), .groups = 'drop') %>% + unlist(use.names = FALSE) + + return(SSE) } #objectiveFunction(par_0,spSurvival_df) @@ -491,11 +689,44 @@ spSurvival_df %>% ``` + + +#Spawner survival after the reproduction + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +#OLD FILES ```{r Observations of presence/absence according to temperature }