Commit a48b2249 authored by Lambert Patrick's avatar Lambert Patrick
Browse files

exploration air temperature water temperature

parent 42edfdf4
# # Code article V1.2 avec ajout DDL
# February 2019
library(dplyr)
library(tidyr)
library(stats)
library(lubridate)
library(purrr)
library(suncalc)
library(furrr)
library(tictoc)
library(caret)
library(readr)
library(patchwork)
library(dismo)
library(gbm)
library(vioplot)
library(fields)
library(gam)
library(voxel)
library(aplpack)
library(scales)
library(tidyverse)
rm(list=ls())
USER="Alexis"
if(USER=="Alexis"){
repertoire_travail<-getwd()
repertoire_data<-paste0(repertoire_travail,"/data/input/")
repertoire_dataGW<-paste0(repertoire_travail,"/data/GlobalWarming/")
repertoire_sorties<-paste0(repertoire_travail,"/data/output/")
}else{
repertoire_travail <- getwd()
repertoire_data <- "../data/input/"
repertoire_sorties <- "../data/output/"
}
#----------------------------------------------------------------------------
#La fonction de contour remasterisée
load(paste0(repertoire_sorties, "contourGBM.RData"))
#La fonction de calcul de dayl-length
load(paste0(repertoire_sorties, "dayLight.RData"))
#----------------------------------------------------------------------------
# Day length -------------------------------------------------------------
st <- as.Date("2019-01-01")
en <- as.Date("2019-12-31")
calendar <- seq.Date(st, en, by = "day")
dayL<-data.frame(date=calendar,
site=c(rep("Garonne", 365),rep("Dordogne",365)),
DL=c(dayLight(dates=calendar),
dayLight(dates=calendar, latitude = 44.8333 , longitude=0.35)))
dayL %>%
group_by(site) %>%
mutate(deltaDL1=ifelse(row_number() == 1, NaN, DL - lag(DL, 1)),
dataLink=format(date, "%d-%m"),
jourDeAnnee=lubridate::yday(date)) %>%
ungroup() %>%
as.data.frame()->dayL
# Day length -------------------------------------------------------------
# Data env + reproduction -------------------------------------------------------------
envCondition<-read.table(paste0(repertoire_data, "envCondition.csv"),
header=TRUE, sep=";") %>%
mutate(jourDeAnnee=lubridate::yday(date)) %>%
filter(jourDeAnnee %in% c(91:213))->data
data$date<-as.Date(data$date)
data$dataLink=base::format(data$date, "%d-%m")
repro<-read.table(paste0(repertoire_data,"reproMonitoring.csv"), sep=";", header = T)
repro$date<-as.Date(repro$date)
repro$Presence<-as.numeric(repro$Bulls>0)
data<-right_join(x=data, y=repro, by=c("lien", "date"))
# Data env + reproduction -------------------------------------------------------------
data<-left_join(x=data, y=dayL, by=c("dataLink","site"))
data<-data[,names(data) %in% c("lien", "site", 'date.x', "Temp1",
"deltaDL1.y", "Q1", "Presence")]
names(data)<-c("lien", "site", 'date', "Temp1", "Q1", "Presence", "deltaDL1")
# END du data--------------------------------------------------------------
indexes = sample(1:nrow(data),size=0.3*nrow(data))
trained_data = data[-indexes,]
test_data = data[indexes,]
# ------------------------------------------------------------------------------
# 1/
# Le modele BRT presence-absence
# ------------------------------------------------------------------------------
Variable<-c("Temp1", "Q1", "deltaDL1")
brtTempPres<-gbm.step(data=trained_data,
gbm.x = Variable, gbm.y = "Presence",
family = "bernoulli", tree.complexity = 5,
learning.rate = 0.001, bag.fraction = 0.5)
Relative_influence<-as.data.frame(summary(brtTempPres))
# ----------------------------------------------------------------------------
data$Preds<-predict.gbm(brtTempPres, data, n.trees=brtTempPres$gbm.call$best.trees, type="response")
d <- cbind(data$Presence, data$Preds)
pres <- d[d[,1]==1 , 2]
abs <- d[d[,1]==0 , 2]
e.dm<-evaluate(p=pres, a=abs)
e.dm
# Check the evaluation results
str(e.dm)
# Boxplot of presence and absence suitable values
# blue are absences, red are presences
boxplot(e.dm,col=c('blue','red'))
# Density plot of presence and absence suitable values
# blue are absences, red are presences
raster::density(e.dm)#
tr.b<-threshold(e.dm,'spec_sens')
# Get the value of the threshold=0.570914
abline(v=tr.b)
#----------------------------------------------------------------------
# Parital plot
par(mfrow=c(2,2))
a<-plot.gbm(brtTempPres,
i.var=1,
return.grid = TRUE, type="link")
b<-plot.gbm(brtTempPres,
i.var=2,
return.grid = TRUE, type="link")
c<-plot.gbm(brtTempPres,
i.var=3,
return.grid = TRUE, type="link")
a %>%
ggplot()+
geom_line(aes(x=Temp1, y=y))+
xlab("Water temperature (°C)")+
ylab("Marginal effect on logit(p)")+
theme_bw()->A
b %>%
ggplot()+
geom_line(aes(x=Q1, y=y))+
xlab("Water discharge (log scale; m3/s)")+
ylab("Marginal effect on logit(p)")+
theme_bw()->B
c %>%
ggplot()+
geom_line(aes(x=deltaDL1, y=y))+
xlab("Day-length difference (hours)")+
ylab("Marginal effect on logit(p)")+
theme_bw()->C
C+A+B+plot_layout(2,2)
#----------------------------------------------------------------------------
# 3-d dependanc plots
find.int <- gbm.interactions(brtTempPres)
find.int$interactions
find.int$rank.list
par(mfrow=c(2,2), mar=c(4,4,1,2))
# Temperature et delta day length
contourGBM(brtTempPres, 3, 1,perspective = FALSE,
smooth = "average", y.label='Water temperature (°C)',
x.label="Day-length difference (hours)")
X<-data.frame(x=data$deltaDL1,
y=data$Temp1)
aplpack::plothulls(X,fraction=1,add=TRUE, pch="", main="")
points(data$Temp1[data$Presence>0]~data$deltaDL1[data$Presence>0],
pch="+",cex=.7)
box()
#-------------------------------------------------------
# Temeprature et debit
contourGBM(brtTempPres, 1, 2,perspective = FALSE,
smooth = "average", x.label='Water temperature (°C)',
y.label="Water discharge (log scale; m3/s)")
X<-data.frame(x=data$Temp1,
y=data$Q1)
aplpack::plothulls(X,fraction=1,add=TRUE, pch="", main="")
points(data$Q1[data$Presence>0] ~ data$Temp1[data$Presence>0],
pch="+",cex=.7)
box()
#---------------------------------------------------------------
contourGBM(brtTempPres, 3, 2,perspective = FALSE,
smooth = "average", x.label="Day-length difference (hours)",
y.label="Water discharge (log scale; m3/s)")
X<-data.frame(x=data$deltaDL1,
y=data$Q1)
aplpack::plothulls(X,fraction=1,add=TRUE, pch="", main="")
points(data$deltaDL1[data$Presence>0],
data$Q1[data$Presence>0],
pch="+",cex=.7)
box()
This diff is collapsed.
This diff is collapsed.
......@@ -141,8 +141,8 @@ juvenileTolerance = c(9.8, 20, 26) # Rougier et al 2015
# juvenileTolerance = c(10.8, 20.8, 29.8) #tolerance temperature range Jatteau 2017
#
spawnerTolerance = c(10.7, 17, 25.7) - 2 # spawning temperature range Paumier et al 2019
juvenileTolerance = c(10.8, 20.8, 29.8) - 2 #tolerance temperature range Jatteau 2017
spawnerTolerance = c(10.7, 17, 25.7) # spawning temperature range Paumier et al 2019
juvenileTolerance = c(10.8, 20.8, 29.8) #tolerance temperature range Jatteau 2017
# Tmin = 8
# spawnerTolerance = c(Tmin, 20, 23)
......@@ -189,12 +189,14 @@ lines(Stot, modifiedBevertonHolt(Stot), col = 'green') # for the Garonne River
lines(c(0, ymax * exp(-sigmaZ)), c(0, ymax), col = 'green', lty = 2)
# ===================================== stock at equilibrium =============================================================
stocksAtEquilibrium = function(temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface, propFemale = 0.5, aFecundity = 135000) {
# as coded in GR3D v3
stocksAtEquilibrium = function(temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface, propFemale, aFecundity, pHoming) {
# !!!!! stock and recruit consider only females !!!!!!
objFct = function(stock, temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface, aFecundity){
survRecruit = survivingRecruit(stock, temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface, aFecundity) * propFemale
return((survRecruit - stock)^2)
objFct = function(femaleStock, temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface, aFecundity, pHoming){
survFemaleRecruit = survivingRecruit(femaleStock, temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface, aFecundity) * propFemale
survFemaleRecruit = survFemaleRecruit * pHoming
return((survFemaleRecruit - femaleStock)^2)
}
stockatequilibrium = 0
......@@ -209,22 +211,30 @@ stocksAtEquilibrium = function(temperature, juvenileTolerance, spawnerTolerance,
recruitatcrash = thermalBevertonHolt(stockatcrash * tempEffSpawner,
temperature, juvenileTolerance, surface, aFecundity) * propFemale
sigmaZatcrash = -log(stockatcrash / recruitatcrash)
Zhoming = -log(pHoming)
if (sigmaZ < sigmaZatcrash) {
# pHoming is added after the mortatly at sea
if ((sigmaZ + Zhoming) < sigmaZatcrash) {
stockattrap = optimise(objFct, c(0, stockatcrash), temperature = temperature, juvenileTolerance = juvenileTolerance, spawnerTolerance = spawnerTolerance,
sigmaZ = sigmaZ, surface = surface, aFecundity = aFecundity)$minimum
sigmaZ = sigmaZ, surface = surface, aFecundity = aFecundity, pHoming = pHoming)$minimum
stockatequilibrium = optimise(objFct, (c(1, 100) * stockatcrash), temperature = temperature, juvenileTolerance = juvenileTolerance, spawnerTolerance = spawnerTolerance,
sigmaZ = sigmaZ, surface = surface, aFecundity = aFecundity)$minimum
sigmaZ = sigmaZ, surface = surface, aFecundity = aFecundity, pHoming = pHoming)$minimum
}
}
return(c(stockatequilibrium, stockattrap, stockatcrash))
}
# ===== Parameters =================================================
temperature = seq(5,30,.1)
aFecundity = 135000 * 2
stocks = as.data.frame(t(sapply(temperature, function(temp) (stocksAtEquilibrium(temp, juvenileTolerance, spawnerTolerance, sigmaZ, surface = surfaceBasin, aFecundity = aFecundity)))))
propFemale = .5
pHoming = .75
# == simulation of spawner mortality, reprpoduction, mortality at sea, homing ==========================================
stocks = as.data.frame(t(sapply(temperature, function(temp) (stocksAtEquilibrium(temp, juvenileTolerance, spawnerTolerance, sigmaZ, surface = surfaceBasin,
propFemale = propFemale, aFecundity = aFecundity, pHoming = pHoming )))))
names(stocks) = c('atEquilibrium', 'atTrap', 'atCrash')
stocks$temperature = temperature
......@@ -241,7 +251,7 @@ stocks[stocks$temperature == tempe, ]
S = seq(0, 5e5, 1e3) # only females
recruit = thermalBevertonHolt(S * temperatureEffect(tempe, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3]),
temperature = tempe, juvenileTolerance = juvenileTolerance, surface = surfaceBasin, aFecundity = aFecundity ) * propFemale
temperature = tempe, juvenileTolerance = juvenileTolerance, surface = surfaceBasin, aFecundity = aFecundity ) * propFemale * pHoming
plot(S, recruit, type = 'l', xlab = 'Stock (# female)', ylab = 'Recruit (# of female)', main = paste(riverName, surfaceBasin))
lines(c(0, max(recruit) * exp(-sigmaZ)), c(0, max(recruit)), col = 'green', lty = 2)
......
library(openxlsx)
library(dplyr)
library(tidyr)
library(lubridate)
library(insol)
rm(list = ls())
# ===== probability to spawn from BRT Paumier et al
library(gbm)
load("BTRarticle.RData")
# ===== juvenile survival from Jatteau et al 2017 ===============================================================================
dailySurvivalEmbryoCoeff = function(temperature, modelSurvivalEmbryoCoeff, parDurationIncubation){
sel = is.na(temperature)
temperature[sel] = 0
duration = (parDurationIncubation[1] * temperature ^ parDurationIncubation[2]) # remove the round() to smooth the curve
predict_matrix = model.matrix(~ Temp + I(Temp ^ 2) + I(Temp ^ 3), data = data.frame(Temp = temperature))
surv = (1 / (1 + exp(-predict_matrix %*% modelSurvivalEmbryoCoeff))) ^ (1 / duration)
surv[sel] = NA
return(surv)
}
cumulativeSurvivalEmbryoCoeff = function(temperature, modelSurvivalEmbryoCoeff){
sel = is.na(temperature)
temperature[sel] = 0
predict_matrix = model.matrix(~ Temp + I(Temp ^ 2) + I(Temp ^ 3), data = data.frame(Temp = temperature))
surv = (1 / (1 + exp(-predict_matrix %*% modelSurvivalEmbryoCoeff)))
surv[sel] = NA
return(surv)
}
dailySurvivalLarvaeCoeff = function(temperature, modelSurvivalLarvaeCoeff) {
sel = is.na(temperature)
temperature[sel] = 0
pred_matrix = model.matrix(~Temp + I(Temp^2) + I(Temp^3), data = data.frame(Temp = temperature))
surv = exp(-1/exp(pred_matrix %*% (modelSurvivalLarvaeCoeff)))
surv[sel] = NA
return(surv)
}
computeJuvenileSurvival = function(data, parDurationIncubation, modelSurvivalEmbryoCoeff, modelSurvivalLarvaeCoeff, horizon) {
data = data %>% select(one_of(c("date", "temperature"))) %>%
mutate(temperatureEmbryo = temperature ^ -parDurationIncubation[2]) %>%
mutate(dailyEmbryoSurvival = dailySurvivalEmbryoCoeff(temperature, modelSurvivalEmbryoCoeff, parDurationIncubation) ) %>%
mutate(dailyLarvalSurvival = dailySurvivalLarvaeCoeff(temperature, modelSurvivalLarvaeCoeff)) %>%
mutate(durationIncubation = NA, juvenileSurvival = NA)
for (t in (1:(nrow(data) - (horizon + 20)))) {
# compute duration of incubation
if (!is.na(data$temperatureEmbryo[t])) {
i = 0
while (sum(data$temperatureEmbryo[t:(t + i)]) < parDurationIncubation[1]) {
i = i + 1
if (is.na(sum(data$temperatureEmbryo[t:(t + i)]))) {
i = NA
break()
}
}
}
durationIncubation = i + 1
data$durationIncubation[t] = durationIncubation
if (!is.na(durationIncubation)) {
# cumulative survival of embryo
embryoSurvival = prod(data$dailyEmbryoSurvival[t:(t + durationIncubation - 1)])
# cumulative survival of larval until horizon days of hatch
larvalSurvival = prod(data$dailyLarvalSurvival[(t + durationIncubation):(t + durationIncubation + horizon - 1)])
# cumulative juvenile survival
data$juvenileSurvival[t] = embryoSurvival * larvalSurvival
}
}
return(data)
}
parDurationIncubation = c(1124.531453, -1.834438)
modelSurvivalEmbryoCoeff = c(-1.245015E+01, 1.120331E+00 , -7.830614E-03, -6.112982E-04)
modelSurvivalLarvaeCoeff = c(-1.851093E+00, 5.132855E-01, -1.020347E-02, -7.005887E-05 )
horizon = 14
# ==== Air temperature ===================================================================================================
# données Safran (mises à disposition par Eric Sauquet)
safranDordogne = read.csv("GARDONNE_T.txt", sep = " ", skip = 1, header = FALSE, col.names = c("date", "temperature")) %>%
mutate(date = as.Date(date), site = "Dordogne") %>%
mutate(season = cut(as.numeric(format(date, '%m')), breaks = c(0, 3,6,9,12), labels = c("winter", "spring", "summer" , "fall"))) %>%
mutate(lien = paste(site, format(date, "%Y"), season, sep = "_"))
safranGaronne = read.csv("LAMAGIST_T.txt", sep = " ", skip = 1, header = FALSE, col.names = c("date", "temperature")) %>%
mutate( date = as.Date(date), site = "Garonne") %>%
mutate(season = cut(as.numeric(format(date, '%m')), breaks = c(0, 3,6,9,12), labels = c("winter", "spring", "summer" , "fall"))) %>%
mutate(lien = paste(site, format(date, "%Y"), season, sep = "_"))
# ==== water temperature + probability to spawn from Paumier et al =======================================================
# données EDF (mis à disposition par Elorri)
dataDordogne = read.xlsx( xlsxFile = 'TQ_Dordogne_Garonne.xlsx', sheet = "Dordogne_1997_2018", detectDates = TRUE) %>%
select(one_of(c("date", "temperature", "debit" ))) %>%
mutate(season = cut(as.numeric(format(date, '%m')), breaks = c(0, 3,6,9,12), labels = c("winter", "spring", "summer" , "fall"))) %>%
mutate(site = "Dordogne", lien = paste(site, format(date, "%Y"), season, sep = "_")) %>%
mutate(dayOfYear = yday(date), DL = daylength(lat = 44.8452, long = 0.4065, jd = dayOfYear, tmz = 1)[,3]) %>%
mutate(Temp1 = temperature, Q1 = log(debit)) %>%
mutate(deltaT1 = c(NA,diff(Temp1)) , deltaQ1 = c(NA,diff(Q1)), deltaDL1 = c(NA, diff(DL)))
dataDordogne$pReproduce = predict.gbm(finalBrtPres, dataDordogne, n.trees = finalBrtPres$gbm.call$best.trees, type="response")
dataDordogne$pReproduce[1] = NA
dataDordogne$juvenileSurvival = computeJuvenileSurvival(dataDordogne, parDurationIncubation, modelSurvivalEmbryoCoeff, modelSurvivalLarvaeCoeff, horizon)$juvenileSurvival
head(dataDordogne,5)
dataGaronne = read.xlsx( xlsxFile = 'TQ_Dordogne_Garonne.xlsx', sheet = "Garonne_1976_2018", detectDates = TRUE) %>%
rename(date = Date) %>% select(one_of(c("date", "temperature", "debit"))) %>%
mutate(season = cut(as.numeric(format(date, '%m')), breaks = c(0, 3,6,9,12), labels = c("winter", "spring", "summer" , "fall"))) %>%
mutate(site = "Garonne", lien = paste(site, format(date, "%Y"), season, sep = "_")) %>%
mutate(dayOfYear = yday(date), DL = daylength(lat = 44.1167, long = 0.8333, jd = dayOfYear, tmz = 1)[,3]) %>%
mutate(Temp1 = temperature, Q1 = log(debit)) %>%
mutate(deltaT1 = c(NA,diff(Temp1)) , deltaQ1 = c(NA,diff(Q1)), deltaDL1 = c(NA, diff(DL)))
dataGaronne$pReproduce = predict.gbm(finalBrtPres, dataGaronne, n.trees = finalBrtPres$gbm.call$best.trees, type="response")
dataGaronne$pReproduce[1] = NA
dataGaronne$juvenileSurvival = computeJuvenileSurvival(dataGaronne, parDurationIncubation, modelSurvivalEmbryoCoeff, modelSurvivalLarvaeCoeff, horizon)$juvenileSurvival
head(dataGaronne)
rm(finalBrtPres)
# combine data =========================================================================================================
selectedSeason = c("spring", "winter","fall", "summer")
# selectedSeason = c("spring", "winter")
# selectedSeason = c("spring")
dailyWaterTemperature = dataDordogne %>% bind_rows(dataGaronne) %>% filter(season %in% selectedSeason) %>% rename(waterTemperature = temperature)
waterData = dailyWaterTemperature %>% group_by(lien, season) %>% summarize(meanWaterTemp = mean(waterTemperature))
head(waterData)
dailyAirTemperature = safranDordogne %>% bind_rows(safranGaronne) %>% filter(season %in% selectedSeason) %>% rename(airTemperature = temperature)
airData = dailyAirTemperature %>% group_by(lien, season) %>% summarize(meanAirTemp = mean(airTemperature))
myData = waterData %>% left_join(airData, by = c("lien", "season"))
myDailyData = dailyWaterTemperature %>% left_join(dailyAirTemperature %>% select(one_of(c("date", "airTemperature", "site"))), by = c("date", "site"))
head(myDailyData)
weightedJuvSurvival = myDailyData %>% mutate(pReproduce = replace_na(pReproduce, 0)) %>% group_by(lien, season) %>% summarise(meanWaterTemp = mean(waterTemperature), meanAirTemp = mean(airTemperature),
meanSurvival = weighted.mean(x = juvenileSurvival, w = pReproduce, na.rm = TRUE))
head(weightedJuvSurvival)
temperature = seq(5,30,.1)
plot(weightedJuvSurvival$meanAirTemp, weightedJuvSurvival$meanSurvival, pch=20, xlim=range(temperature), ylim=c(0,1),
col = c("red", "blue", "yellow", "green")[as.factor(weightedJuvSurvival$season)])
legend("topleft", legend = c("winter", "spring", "summer", "fall"), col = c("red", "blue", "yellow", "green"), lty = -1, pch = 20)
embryoSurvival = cumulativeSurvivalEmbryoCoeff(temperature, modelSurvivalEmbryoCoeff)
larvalSurvival = dailySurvivalLarvaeCoeff(temperature, modelSurvivalLarvaeCoeff)^horizon
points(weightedJuvSurvivalSplash$meanWaterTemp, weightedJuvSurvivalSplash$meanSurvival, pch = 20, col = 'darkblue')
lines(temperature, embryoSurvival * larvalSurvival, col = 'red')
lines(temperature, temperatureEffect(temperature, juvenileTolerance[1], juvenileTolerance[2], juvenileTolerance[3]) * max(embryoSurvival * larvalSurvival), col = "blue")
#legend("topleft", legend = c("uniform bull", "observed splash", "juvenile tolerance", "temperature effect"),
# pch = c(20,20,-1,-1), lty = c(0,0, 1,1), col = c('black', 'green', 'red', "blue"))
# with fake day of the year ====================================
myDailyDataFake = myDailyData %>% mutate(fakeDayOfYear = case_when(season == "winter" ~ dayOfYear + 90,
season == "spring" ~ dayOfYear,
season == "summer" ~ dayOfYear - 91,
season == " winter" ~ dayOfYear - 183)) %>%
mutate(DL = case_when(site == "Garonne" ~ daylength(lat = 44.1167, long = 0.8333, jd = fakeDayOfYear, tmz = 1)[,3],
site == "Dordogne" ~ daylength(lat = 44.8452, long = 0.4065, jd = fakeDayOfYear, tmz = 1)[,3])) %>%
mutate(deltaDL1 = c(NA, diff(DL)))
myDailyDataFake$pReproduceFake = predict.gbm(finalBrtPres, myDailyDataFake, n.trees = finalBrtPres$gbm.call$best.trees, type="response")
myDailyDataFake$pReproduceFake[1] = NA
head(myDailyDataFake)
weightedJuvSurvivalFake = myDailyDataFake %>% mutate(pReproduceFake = replace_na(pReproduceFake, 0)) %>% group_by(lien, season) %>%
summarise(meanWaterTemp = mean(waterTemperature), meanAirTemp = mean(airTemperature),
meanSurvival = weighted.mean(x = juvenileSurvival, w = pReproduceFake))
plot(weightedJuvSurvivalFake$meanAirTemp, weightedJuvSurvivalFake$meanSurvival, pch = 20, xlim = range(temperature), ylim = c(0,1),
col = c("red", "blue", "yellow", "green")[as.factor(weightedJuvSurvivalFake$season)])
legend("topleft", legend = c("winter", "spring", "summer", "fall"), col = c("red", "blue", "yellow", "green"), lty = -1, pch = 20)
weightedJuvSurvivalFake %>% filter(season == "summer" & meanSurvival <.45)
myDailyDataFake %>% filter(between(date, as.Date("2003-06-25"), as.Date("2003-07-10")), site == 'Dordogne' )
# calibration ============================================================================================
objFun = function(par, observationData) {
SE = (observationData$meanSurvival - par[4] * Rosso2(observationData$meanAirTemp, par))^2
return(sum(SE, na.rm = TRUE))
}
objFunLehman = function(par, observationData) {
SE = (observationData$meanSurvival - par[4] * modifiedLehman(observationData$meanAirTemp, par))^2
return(sum(SE, na.rm = TRUE))
}
data = weightedJuvSurvivalFake %>% drop_na()
resRosso2 = nlm(objFun, p = c(10,17,25,.6), observationData = data)
resLehman = nlm(objFunLehman, p = c(10,17,25,.6), observationData = data)
print(resRosso2)
plot(data$meanAirTemp, data$meanSurvival, pch=20, xlim = range(temperature), ylim=c(0,1),
col = c("red", "blue", "yellow", "green")[as.factor(data$season)])
legend("topleft", legend = c("winter", "spring", "summer", "fall"), col = c("red", "blue", "yellow", "green"), lty = -1, pch = 20)
# lines(temperature, Rosso2(temperature, juvenileTolerance) * max(embryoSurvival * larvalSurvival), col = "blue", lty =2)
lines(temperature, Rosso2(temperature, resRosso2$estimate[1:3]) * resRosso2$estimate[4], col = "blue")
lines(temperature, modifiedLehman(temperature, resLehman$estimate[1:3]) * resLehman$estimate[4], col ='cyan')
library(dplyr)
library(tidyr)
library(survival)
library(lme4)
#library(pracma)
......@@ -28,7 +29,7 @@ dailySurvivalLarvae = function(temperature, modelSurvivalLarvae) {
return(surv)
}
# with lme coeff
# directly with coeff
dailySurvivalEmbryoCoeff = function(temperature, modelSurvivalEmbryoCoeff, parDurationIncubation){
sel = is.na(temperature)
temperature[sel] = 0
......@@ -39,6 +40,16 @@ dailySurvivalEmbryoCoeff = function(temperature, modelSurvivalEmbryoCoeff, parDu
return(surv)
}
cumulativeSurvivalEmbryoCoeff = function(temperature, modelSurvivalEmbryoCoeff){
sel = is.na(temperature)