##### Calculating the Tprefs data for analyses and for calculating Db and De

# Packages ----------------------------------------------------------------
library(nlme)
library(ggplot2)
library(MASS)
library(MuMIn)
library(Rmisc)
library(stringr)

# LOADING THE ALL_TPREFS.csv dataset : tpref

# Calibration Tpref --------------------------------------------------------
tpref$TEMP <- -4.50 + 1.17*tpref$TEMP


# For analyses of the Tpref variation in ROB ------------------------------------------
tprefROB <- tpref[tpref$POP=="ROB",]

# For comparison of Tpref between pops ----------------
tprefPOP <- tpref[tpref$AGE=="A",]

### THESE TWO DATASETS HAVE TO SAVED



# TPREFS AT THE POP AND INDIV CATEGORY SCALE ------------------------------

# Deleting outliers data (see methods) -------------------------------------
tpref <- tpref[tpref$TEMP>25 & tpref$TEMP<41,]


# Loop for the calculation --------
POP <- unique(tpref$POP)
tpref$DATE <- as.POSIXlt(tpref$DATE,format=c("%Y-%m-%d"))
tpref$YEAR <- tpref$DATE$year+1900
year <- unique(tpref$YEAR)
count <- data.frame(xtabs(~POP+YEAR,tpref))
tpref2 <- count[count$Freq>0,][,1:2]
tpref2$SEX <- "F"
tpref3 <- tpref2
tpref3$SEX <- "M"
tpref2 <- rbind(tpref2,tpref3)
rm(tpref3)
tpref2$AGE <- "A"
tpref2 <- rbind(tpref2,c("ROB",2017,"F","SA"),c("ROB",2017,"M","SA"),c("ROB",2018,"F","SA"),c("ROB",2018,"M","SA"))
tpref2$QUANTl <- NA
tpref2$QUANTm <- NA
tpref2$QUANTh <- NA
tpref2$MEAN <- NA
for (l in 1:dim(tpref2)[1]) {
  p <- as.character(tpref2$POP[l])
  y <- as.character(tpref2$YEAR[l])
  s <- as.character(tpref2$SEX[l])
  a <- as.character(tpref2$AGE[l])
  tpref_subset <- tpref[tpref$POP==p & tpref$YEAR==y & tpref$SEX==s & tpref$AGE==a,]
  tpref2$QUANTl[l] <- as.numeric(quantile(tpref_subset$TEMP,na.rm=T)[2])
  tpref2$QUANTm[l] <- as.numeric(quantile(tpref_subset$TEMP,na.rm=T)[3])
  tpref2$QUANTh[l] <- as.numeric(quantile(tpref_subset$TEMP,na.rm=T)[4])
  tpref2$MEAN[l] <- as.numeric(mean(tpref_subset$TEMP,na.rm=T))
}