# Bastien Bellemin-Noel 2020

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


library(dplyr)
library(ggplot2)

rm(list = ls())
dev.off()

bud <- read.csv("budburst 16 17 18 19.csv",sep = ";", header=TRUE)
bud$year <- as.factor (bud$year)
bud$tree <- as.factor (bud$tree)

bud.2017 <- bud[bud$year=="2017",]
bud.2018 <- bud[bud$year=="2018",]
bud.2019 <- bud[bud$year=="2019",]

cat("\014")

#### 2017 ########################################################################################################################

mean.bud.2017 <- aggregate(bud.2017$mean.bud.stage, list(bud.2017$treatment, bud.2017$sp, bud.2017$day), mean)
sd.bud.2017 <- aggregate(bud.2017$mean.bud.stage, list(bud.2017$treatment, bud.2017$sp, bud.2017$day), sd)
colnames(mean.bud.2017) <- c("treatment","sp","day","mean.bud")
mean.bud.2017$sd <- sd.bud.2017$x
mean.bud.2017

ggplot(mean.bud.2017, aes(x = day, y = mean.bud))+
  geom_line(aes(linetype=sp, colour=treatment, size = sp))+
  geom_point(aes(colour="black"))+
  scale_linetype_manual(values=c("solid", "dashed"))+
  labs (title = "Mean budstage by treatment as a function of date in 2017",
        x = "Date (julian)", y = "Mean budstage (apical + north + south)" )+
  scale_color_manual(values = c("black", "orange", "dodgerblue"))+
  xlim(130, 165)+
  ylim(0.7,4)+
  scale_size_manual(values=c(1.2, 1.6))+
  theme_bw()

# model: uses the primary datasets, not the average by day e.g, 1 point/tree/day
a=glm(mean.bud.stage ~ treatment * sp + day, data= bud.2017, family= quasipoisson)
summary(a)

#### 2018 ########################################################################################################################

mean.bud.2018 <- aggregate(bud.2018$mean.bud.stage, list(bud.2018$treatment, bud.2018$sp, bud.2018$day), mean)
colnames(mean.bud.2018) <- c("treatment","sp","day","mean.bud")
mean.bud.2018

ggplot(mean.bud.2018, aes(x = day, y = mean.bud))+
  geom_line(aes(linetype=sp, colour=treatment, size = sp))+
  geom_point(aes(colour="black"))+
  scale_linetype_manual(values=c("solid", "dashed"))+
  labs (title = "Mean budstage by treatment as a function of date in 2018",
        x = "Date (julian)", y = "Mean budstage (apical + north + south)" )+
  scale_color_manual(values = c("black", "orange", "dodgerblue"))+
  xlim(130, 165)+
  ylim(0.7,4)+
  scale_size_manual(values=c(1.2, 1.6))+
  theme_bw()

# model
b=glm(mean.bud.stage ~ treatment * sp + day, data= bud.2018, family= quasipoisson)
summary(b)

#### 2019 ########################################################################################################################

# Day 140 needs to be reemoved from the dataset, mean budstage is way too high for BS t on that day. 
# see following bloc for details and graph
bud.2019.bis <- bud.2019[bud.2019$day != "140",]
mean.bud.2019 <- aggregate(bud.2019.bis$mean.bud.stage, list(bud.2019.bis$treatment, bud.2019.bis$sp, bud.2019.bis$day), mean)
colnames(mean.bud.2019) <- c("treatment","sp","day","mean.bud")
mean.bud.2019

ggplot(mean.bud.2019, aes(x = day, y = mean.bud))+
  geom_line(aes(linetype=sp, colour=treatment, size = sp))+
  geom_point(aes(colour=""))+
  scale_linetype_manual(values=c("solid", "dashed"))+
  labs (title = "Mean budstage by treatment as a function of date in 2019",
        x = "Date (julian)", y = "Mean budstage (apical + north + south)" )+
  scale_color_manual(values = c("black", "orange", "dodgerblue"))+
  xlim(130, 165)+
  ylim(0.7,4)+
  scale_size_manual(values=c(1.2, 1.6))+
  theme_bw()

# model
c=glm(mean.bud.stage ~ treatment * sp + day, data= bud.2019, family= quasipoisson)
summary(c)

#### day 140 2019 ########################################################################################################################
mean.bud.2019 <- aggregate(bud.2019$mean.bud.stage, list(bud.2019$treatment, bud.2019$sp, bud.2019$day), mean)
colnames(mean.bud.2019) <- c("treatment","sp","day","mean.bud")

library(grid)

ggplot(mean.bud.2019, aes(x = day, y = mean.bud))+
  geom_line(aes(linetype=sp, colour=treatment, size = sp))+
  geom_point(aes(colour=""))+
  scale_linetype_manual(values=c("solid", "dashed"))+
  labs (title = "error in measurements in 2019 on day 140",
        x = "Date (julian)", y = "Mean budstage (apical + north + south)" )+
  scale_color_manual(values = c("black", "orange", "dodgerblue"))+
  xlim(130, 165)+
  ylim(0.7,4)+
  scale_size_manual(values=c(1.2, 1.6))+
  geom_segment(aes(x = 136, y = 1.9, xend = 139, yend = 1.7),
                  arrow = arrow(length = unit(0.5, "cm")))+
  theme_bw()

# distribution of the bud stage by day around day 140 higlights the flaw: 
bud.2019.BS <- bud.2019[bud.2019$sp=="EPN",]
bud.2019.BS.136 <- bud.2019.BS[bud.2019.BS$day=="136",]
bud.2019.BS.140 <- bud.2019.BS[bud.2019.BS$day=="140",]
bud.2019.BS.142 <- bud.2019.BS[bud.2019.BS$day=="142",]
par(mfrow=c(3,1))
hist(bud.2019.BS.136$mean.bud.stage, col = "firebrick3", xlim = c(0,2.5),ylim = c(0,30),main ="BS mean bud stage on day 136, 2019", xlab = "Mean bud stage")
hist(bud.2019.BS.140$mean.bud.stage, col = "dodgerblue3", xlim = c(0,2.5),ylim = c(0,30),main ="BS mean bud stage on day 140, 2019", xlab = "Mean bud stage")
hist(bud.2019.BS.142$mean.bud.stage, col = "orange", xlim = c(0,2.5),ylim = c(0,30),main ="BS mean bud stage on day 142, 2019", xlab = "Mean bud stage")