--- 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(temperature_O = mean(temperature), .groups = 'drop') #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(temperature_RI = mean(temperature), .groups = 'drop') # 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)) #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) / 10) %>% # in cm ungroup() growthInBasin <-growthInOffshore %>% inner_join(tempEffectInOffshore, by = c("basin_name" = "inshore_basin_name", "season")) %>% right_join(tempEffectInFreshWat, by = c("basin_name","season")) growthInBasin <- growthInBasin %>% mutate(temperature_RIO = ifelse(season == "spring" & age == 0, temperature_RI, ifelse(season == "summer" & age == 0.25, temperature_RI, ifelse(season == "fall" & age == 0.50, temperature_RI, temperature_O)))) %>% arrange(basin_name,age,season) %>% select(basin_name, offshore_basin_name, metapop,age, season, temperature_O, temperature_RIO, LStich) #group_by(metapop,age,season) %>% #summarize(LStich = unique(LStich)) # 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'))) ``` There is 1 river ( North (MA)) in Stich et al paper that is not referenced in our database. ```{r join between Stich and GR3D database} # NEW ===== 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"), basin_name = replace(basin_name, basin_name == "Upper Chesapeake Bay", "Susquehanna")) 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() ``` ```{r range of temperature} # 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) ``` ```{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} # new ==== save(growParModifed,growthInBasin, nea_presence, Stich2020_sel, file = 'SOS.rdata' ) rm(list = ls()) load('SOS.rdata') source("../GR3D_Rdescription/GR3Dfunction.R") ``` ```{r growth optimisation function Camille, 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 #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) ``` ```{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) ``` 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 ```{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() 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') growth_optimal <- optim(c(Tmin = vecPar[1], Topt = vecPar[2], Tmax = vecPar[3], L0 = vecPar[4], Linf = vecPar[5], K = vecPar[6]), 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)) toc() # 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 #Optimized parameters growParOptim <- data.frame(value = growth_optimal$par) %>% rownames_to_column("parameter") #Rdata #save(growth_optimal, res, res_B, res_C, file = "growth_optimal.RData") #csv write_csv(growParOptim, "growParOptim.csv") ``` ```{r graph to present results} # NEW # =============================================================================== #parGrowthUsed = c(res_B$par) %>% enframe() %>% pivot_wider() parGrowthUsed = c(res_C$par, fixedPar_C) %>% enframe() %>% pivot_wider() # ================================================================================ dataVerif <- computeGrowAllBasins(data = dataCalibration, growPar = parGrowthUsed %>% mutate(sigmaDeltaLVonBert = 0) ) # 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)) save(dataVerif, growth_optimal,res, res_B, res_C, file = "growth_optimal.RData") # thermal response curve with average seasonal temperature experienced by fish dataTemperatureMetapop <- dataVerif %>% group_by(metapop, season) %>% summarise(temperature = mean(temperature), .groups = 'drop') %>% mutate(tempEffect = temperatureEffect(temperature, parGrowthUsed$tempMinGrow, parGrowthUsed$tempOptGrow, parGrowthUsed$tempMaxGrow ), 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)) %>% mutate(tempEffect = temperatureEffect(temperature, parGrowthUsed$tempMinGrow, parGrowthUsed$tempOptGrow, parGrowthUsed$tempMaxGrow )) %>% 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) ``` ```{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) ```