# Bastien Bellemin-Noel 2020

# This scripts corresponds to the "development time" and "pupal mass" analysis in the article 
# "Improved performance of the eastern spruce budworm on black spruce as warming temperatures disrupt phenological defenses" 

rm(list = ls())
dev.off()
cat("\014")

library(dplyr)
library(ggplot2)
library(car)
library(lme4)
library(nlme)
library(lmerTest)
library(AICcmodavg)
source("LME_diagnostics.R")

pup <- read.csv("pup.csv")
pup$year <- as.factor(pup$year)
pup$tree <- as.factor(pup$tree)
pup$tface <- as.factor(pup$tface)

#################################################### figures  #################################################
pup_m <- pup[pup$sex=="M",]
pup_f <- pup[pup$sex=="F",]

# Boxplot males on 3 years
ggplot(pup_m, aes(x = sp, y = mass, fill = treatment))+
  geom_boxplot()+
  ylim(0.03,0.140)+
  facet_grid(.~year,scales="free")+
  labs (title = "Males mass", x = "", y = "mass (g)" )+
  scale_fill_manual(values = c("firebrick","moccasin"))+
  theme_bw()

# Boxplot females on 3 years
ggplot(pup_f, aes(x = sp, y = mass, fill = treatment))+
  geom_boxplot()+
  ylim(0.03,0.140)+
  facet_grid(.~year,scales="free")+
  labs (title = "Female mass", x = "", y = "mass (g)" )+
  scale_fill_manual(values = c("pink2","palegoldenrod"))+
  theme_bw()

# Boxplots for individual variables 
par(mfrow = c(2,2))
plot(pup$sex, pup$mass)
plot(pup$treatment, pup$mass)
plot(pup$sp, pup$mass)
plot(pup$year, pup$mass)

#### global model ####

fit <- lme(mass ~ treatment * sp * sex*  year, random= ~1|tface/tree, data=pup)

diag_plots(fit, level = 2)
res <-residuals(fit, level = 2, type = "pearson")

plot(pup$sex, res)
plot(pup$treatment, res)
plot(pup$sp, res)
plot(pup$year, res)
dev.off()

# weighted variance by sex
fit_varsex <- lme(mass ~ treatment * sp * sex*  year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|sex))
diag_plots(fit_varsex, level = 2)

# weighted variance by year
fit_varyr <- lme(mass ~ treatment * sp * sex*  year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year))
# by both year and sex
fit_varyrsex <- lme(mass ~ treatment * sp * sex*  year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex))

# AICc score comparison
AICc(fit_varsex)
AICc(fit_varyr)
AICc(fit_varyrsex)
# varaince structure according to year and sex seems best

diag_plots(fit_varyrsex, level = 2)
res <-residuals(fit_varyrsex, level = 2, type = "pearson")

range(pup$mass) # rescale to get bigger values
pup$mass2 <- pup$mass*1000

# Selection of fixed effects 
fit_var1 <- lme(mass2 ~ treatment * sp * sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")
fit_var2 <- lme(mass2 ~ treatment * sp + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")
fit_var3 <- lme(mass2 ~ treatment + sp + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")
fit_var4 <- lme(mass2 ~ treatment + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")
fit_var5 <- lme(mass2 ~ treatment * year + sex + sp, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")
fit_var6 <- lme(mass2 ~ year + sex + sp, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")
fit_var7 <- lme(mass2 ~ year + sex, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")
fit_var8 <- lme(mass2 ~ treatment + sp + sex + year + treatment * sp + treatment * sex + treatment * year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "ML")

aictab(list(fit_var1,fit_var2,fit_var3,fit_var4,fit_var5,fit_var6,fit_var7,fit_var8))
# model that provides the best fit is without treatment but with sp (6)
# 2nd best is without sp and treatment (7)
# further proof that they have no influence ? 
# 3rd is without sp and no interactions (4)

#### final model mass ####
pup$mass2 <- pup$mass*1000
fit_var3REML <- lme(mass2 ~ treatment + sp + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex), method = "REML")
diag_plots(fit_var3REML, level = 2)
summary(fit_var3REML)


#################### Development time analysis ####################################
fit <- lme(date.picked ~ treatment + sp + sex + year, random= ~1|tface/tree, data=pup, method = "ML")

diag_plots(fit, level = 2)
res <-residuals(fit, level = 2, type = "pearson")

par(mfrow = c(2,2))
plot(pup$sex, res)
plot(pup$treatment, res)
plot(pup$sp, res)
plot(pup$year, res)
dev.off()

# weighted variance by sex
fit_varsex <- lme(date.picked ~ treatment + sp + sex+  year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|sex))
diag_plots(fit_varyr, level = 2)

# weighted variance by year
fit_varyr <- lme(date.picked ~ treatment + sp + sex+  year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year))
# by both year and sex
fit_varyrsex <- lme(date.picked ~ treatment + sp + sex+  year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year*sex))

# AICc score comparison
AICc(fit_varsex)
AICc(fit_varyr)
AICc(fit_varyrsex)

fit_var1 <- lme(date.picked ~ treatment * sp * sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "ML")
fit_var2 <- lme(date.picked ~ treatment * sp + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "ML")
fit_var3 <- lme(date.picked ~ treatment + sp + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "ML")
fit_var4 <- lme(date.picked ~ treatment + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "ML")
fit_var5 <- lme(date.picked ~ treatment * year + sex + sp, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "ML")
fit_var6 <- lme(date.picked ~ year + sex + sp, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "ML")
fit_var7 <- lme(date.picked ~ year + sex, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "ML")

aictab(list(fit_var1,fit_var2,fit_var3,fit_var4,fit_var5,fit_var6,fit_var7))

#### final model dvp ####
fit.REML <- lme(date.picked ~ treatment + sp + sex + year, random= ~1|tface/tree, data=pup, weights = varIdent(form = ~ 1|year), method = "REML")
diag_plots(fit.REML, level = 2)
summary(fit.REML)