--- 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) library(readr) library(XML) library(forcats) library(knitr) library(officedown) library(flextable) library(stringr) library(tictoc) library(tibble) ``` ```{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)) ``` ```{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 %>% distinct(wintering_offshore_basin_name,summering_offshore_basin_name) #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 wint<-connections_trans %>% filter(season %in% c("winter","spring"), offshore_basin_type == "wintering") %>% select(inshore_basin_name,offshore_basin_name,year,season,temperature) #in summer sum<-connections_trans %>% 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) %>% ungroup() %>% arrange(inshore_basin_name, year,season) %>% group_by(inshore_basin_name,offshore_basin_name,season) %>% summarise(Tmean = mean(temperature)) #mutate(temperatureEffect = temperatureEffect(Tmean, growPar$tempMinGrow, growPar$tempOptGrow, growPar$tempMaxGrow)) #delete intermediate files remove(sum,wint) ``` ```{r Build the temperature expericenced by fish during their residency in river and inshore basins, echo = FALSE, warning=FALSE} #TODO #Contruire un tableau, du même style que celui-ci avec les temperatures du natal basins #season = c("spring","summer","fall","winter") #zone = c("river_tempertaure","inshore_temperature") tempEfffectInRiver <- tempInriver %>% 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") #Optimzation tempRiver <- tempEfffectInRiver %>% filter(season %in% c("spring","summer") & zone == "river_temperature") tempInshore <- tempEfffectInRiver %>% filter(season %in% c("fall","winter") & zone == "inshore_temperature") tempEffectInFreshWat <- bind_rows(tempRiver,tempInshore) %>% arrange(basin_name,year,season) %>% select(basin_name,year, season, temperature) %>% group_by(basin_name,season) %>% summarize(Tmean_RIB = mean(temperature)) # tempEfffectInRiver %>% # {bind_rows( # filter(season == "summer" & zone == "river_temperature"), # filter(season == "fall" & zone == "inshore_temperature")) # } ``` The growth parameters used in GR3D give the following curves when random and temperature effects are not considered. ```{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)') ``` To optimize the growth parameters included in GR3D for the Amercian shad, we used values gave by Stich et al, 2020. ```{r Growth from Stich et al, 2020, echo = FALSE, warning=FALSE} #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){ L0 = Linf * (1-exp(t0*K)) return(L0) } #vonBertalanffyTry(454,0.49,-0.21) #parameters from Stich 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) %>% mutate(L0_theo = vonBertalanffyLengthAtHatch(Linf,K,t0)) %>% mutate_at(vars(c(Linf,L0_theo)), funs(./10)) #from mm to cm #create a data frame and compute L at age with values from Stich growthInOffshore <- expand_grid(basin_name = nea_riverBasinFeatures$basin_name , age = seq(0,9.75,.25)) %>% 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' )) %>% inner_join(nea_riverBasinFeatures %>% select(basin_name, lat_outlet), by ='basin_name') %>% 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) %>% mutate(LStich = vonBertalanffyGrowth(age, L0_theo, Linf, K)) %>% ungroup() growthInOffshore <-growthInOffshore %>% inner_join(tempEffectInOffshore, by = c("basin_name" ="inshore_basin_name", "season")) %>% right_join(tempEffectInFreshWat, by = c("basin_name","season")) growthInOffshore<-growthInOffshore %>% mutate(temperature = ifelse (season == "spring" & age == 0, Tmean_RIB, ifelse(season == "summer" & age == 0.25, Tmean_RIB, ifelse(season == "fall" & age == 0.50, Tmean_RIB, Tmean)))) %>% arrange(basin_name,age,season) %>% select(basin_name,offshore_basin_name,metapop,age, season,temperature,LStich) #group_by(metapop,age,season) %>% #summarize(LStich = unique(LStich)) ``` ```{r Growth in GR3D with temperature effect, echo = FALSE, warning=FALSE} #define parameters for both sex based on value from the XML file growParModifed <- growPar %>% 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 %>% #inner_join(tempEffectInOffshore %>% #select(inshore_basin_name,season,temperatureEffect), by =c("basin_name" ="inshore_basin_name", "season" ="season")) %>% #select(basin_name,lat_outlet, metapop,age,season,LStich, temperatureEffect) ``` ```{r growth optimisation, echo = FALSE, warning=FALSE} #vector on paramters to get optimized #par[1] = Tmin #par[2] = Topt #par [3] = Tmax #par[4] = L0 #par[5] = Linf #par[6] = k vecPar = growParModifed %>% select(tempMinGrow,tempOptGrow,tempMaxGrow,lengthAtHatching,Linf, kOpt) %>% pivot_longer(tempMinGrow:kOpt, names_to = "parameter",values_to = "value") vecPar = as.vector(vecPar$value) #Return data computeGrowOpt = function(par,data){ 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 %>% 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() 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)) } return(data) } #computeGrowOpt(vecPar,growthInOffshore) #Function to be optimized computeSSE = function(par,data){ data = computeGrowOpt(par,data) SSE = data %>% group_by(metapop) %>% mutate(W = 1/n()) %>% ungroup() %>% transmute(squareError = W * ((L - LStich)^2)) %>% summarise(sumSquare = sum(squareError)) %>% unlist() #ungroup() %>% #select(sumSquare) #summarize(sum = sum(sumSquare)) #return(SSE = as.vector(SSE$sum)) return(SSE) } #computeSSE (vecPar, growthInOffshore) tic() growth_optimal <- optim(c(Tmin = vecPar[1], Topt = vecPar[2], Tmax = vecPar[3], L0 = vecPar[4], Linf = vecPar[5], K = vecPar[6]), computeSSE, data = growthInOffshore, lower = c(0,0,0,0,0,0), method = "L-BFGS-B", control = list(trace = TRUE, parscale = c(1,1,1,0.1,1,0.01))) toc() #Optimized parameters growParOptim <- data.frame(value = growth_optimal$par) %>% rownames_to_column("parameter") #Rdata save(growth_optimal, file = "growth_optimal.RData") #csv write_csv(growParOptim, "growParOptim.csv") ``` ```{r draw growth curves} load("growth_optimal.RData") growthCurves = computeGrowOpt(growth_optimal$par, growthInOffshore) save(growthCurves, file = "growthCurves.RData") #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)") + #scale_x_continuous(limits = c(0, 10)) + geom_hline(yintercept = 45, color = "pink", linetype =1)+ #female geom_hline(yintercept = 40, color = "lightblue", linetype = 1)+ geom_vline(aes(age), xintercept = 8, linetype = 1)+ facet_wrap(.~ metapop, ncol =1) #temperature effect temp <- seq(0,45,0.1) growthCurves %>% ggplot(aes(x = temperature)) + 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, growParOptim$Tmin, growParOptim$Topt, growParOptim$Tmax)) df %>% ggplot(aes(x = temperature)) + geom_line(aes(y = temperatureEffect)) + geom_point(aes(y = temperatureEffect), col = 'red', alpha = 0.4) ``` ```{r growth try, echo = FALSE, warning=FALSE} #k depending on water temperature in ofsshore basins kTempEffect = function(kopt, temp, Tmin, Topt, Tmax){ kTemp <- kopt * temperatureEffect(temp, Tmin, Topt, Tmax) return(kTemp)} #Generic function vonBertalanffyGrowthWithTempEffect = function(temp, age, L0, Linf, Kopt, Tmin, Topt, Tmax){ koptTemp = kTempEffect(Kopt, temp, Tmin, Topt, Tmax) t0 = log(1 - L0 / Linf) / koptTemp return(Linf * (1 - exp(-koptTemp * (age - t0)))) } ```