# to clear the workspace
rm(list =ls())

#load libraries
library(multcomp)
library(xml2)
library(grid)
library(gridExtra)
library(RColorBrewer)
library(tidyverse)
library(dplyr)
library(plyr)
library(stats)
library(MASS)

# load data
setwd("C:/Users/awvater/Desktop")
ConditionIndexData <- read.csv("C:/Users/awvater/Desktop/Figure_3data.csv")
View(ConditionIndexData)

#-------Day0--------b
Day0 <- subset(ConditionIndexData, Day == 0)
summary(Day0)

LMDay0 <- lm(CI.g.cm.3. ~ C_T_S, data = Day0)

ComparisonsDay0 <- glht (LMDay0, mcp(C_T_S = "Tukey"))
summary(ComparisonsDay0)

# show q-q plot before windsorization
qqnorm(resid(LMDay0))
shapiro.test (resid(LMDay0))

# --------------------- Code to Windsorize Day 0 ----------------------
res = resid(LMDay0)

# get 1% and 99% for residual
q <- quantile(res, probs=c(.01, .99), na.rm = TRUE)

# perform windsorization 
for (i in 1:length(res)) {
  if (res[i] < q[1]) {
    Day0$CI.g.cm.3.[i] <- Day0$CI.g.cm.3.[i] - res[i] + q[1]
  } else if (res[i] > q[2]) {
    Day0$CI.g.cm.3.[i] <- Day0$CI.g.cm.3.[i] - res[i] + q[2]
  }
}

# --------------------------------------------------------------------

# create new linear model with windsorized data
LMDay0 <- lm(CI.g.cm.3. ~ C_T_S, data = Day0)

# look at q-q norm and Shapiro test windsorized data
qqnorm(resid(LMDay0))
shapiro.test (resid(LMDay0))

ComparisonsDay0 <- glht (LMDay0, mcp(C_T_S = "Tukey"))
summary(ComparisonsDay0)

#------------Day161------------------------------------------

Day161 <- subset(ConditionIndexData, Day == 161)
summary(Day161)

LMDay161 <- lm(CI.g.cm.3. ~ C_T_S, data = Day161)

ComparisonsDay161 <- glht (LMDay161, mcp(C_T_S = "Tukey"))
summary(ComparisonsDay161)

# show q-q plot before windsorization
qqnorm(resid(LMDay161))
shapiro.test (resid(LMDay161))

# --------------------- Code to Windsorize Day 161 ----------------------
res = resid(LMDay161)

# get 1% and 99% for residual
q <- quantile(res, probs=c(.01, .99), na.rm = TRUE)

# perform windsorization 
for (i in 1:length(res)) {
  if (res[i] < q[1]) {
    Day161$CI.g.cm.3.[i] <- Day161$CI.g.cm.3.[i] - res[i] + q[1]
  } else if (res[i] > q[2]) {
    Day161$CI.g.cm.3.[i] <- Day161$CI.g.cm.3.[i] - res[i] + q[2]
  }
}

# --------------------------------------------------------------------

# create new linear model with windsorized data
LMDay161 <- lm(CI.g.cm.3. ~ C_T_S, data = Day161)

# look at q-q norm and Shapiro test windsorized data
qqnorm(resid(LMDay161))
shapiro.test (resid(LMDay161))

#----ANAOVA-----
ComparisonsDay161 <- glht (LMDay161, mcp(C_T_S = "Tukey"))
summary(ComparisonsDay161)
#-----------------------Day343----------------------


Day343 <- subset(ConditionIndexData, Day == 343)
summary(Day343)

LMDay343 <- lm(CI.g.cm.3. ~ C_T_S, data = Day343)

ComparisonsDay343 <- glht (LMDay343, mcp(C_T_S = "Tukey"))
summary(ComparisonsDay343)

# show q-q plot before windsorization
qqnorm(resid(LMDay343))
shapiro.test (resid(LMDay343))

# --------------------- Code to Windsorize Day 343 ----------------------
res = resid(LMDay343)

# get 1% and 99% for residual
q <- quantile(res, probs=c(.01, .99), na.rm = TRUE)

# perform windsorization 
for (i in 1:length(res)) {
  if (res[i] < q[1]) {
    Day343$CI.g.cm.3.[i] <- Day343$CI.g.cm.3.[i] - res[i] + q[1]
  } else if (res[i] > q[2]) {
    Day343$CI.g.cm.3.[i] <- Day343$CI.g.cm.3.[i] - res[i] + q[2]
  }
}

# --------------------------------------------------------------------

# create new linear model with windsorized data
LMDay343 <- lm(CI.g.cm.3. ~ C_T_S, data = Day343)

# look at q-q norm and Shapiro test windsorized data
qqnorm(resid(LMDay343))
shapiro.test (resid(LMDay343))

#----ANAOVA-----
ComparisonsDay343 <- glht (LMDay343, mcp(C_T_S = "Tukey"))
summary(ComparisonsDay343)


# ------- construct data frame of the conditioned data -----
dtemp <- rbind(Day0, Day161)
data <- rbind(dtemp,Day343)

datasum <- ddply(data, c("Day", "C_T_S"), summarise,
                 N = length(CI.g.cm.3.), 
                 mean=mean(CI.g.cm.3.),
                 sd=sd(CI.g.cm.3.),
                 se=sd/sqrt(N) )

pd <- position_dodge(2) # move them .05 to the left and right

ggplot(datasum, aes(x=Day, y=mean, colour=C_T_S, group=C_T_S)) + 
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), colour="black", width=5) +
  geom_line(position=pd) +
  geom_point(position=pd, size=3, shape=21, fill="white") +
  xlab("Days") + 
  ylab("Condition Index (total weight/length^3") +
  scale_colour_hue(name=" Experimental Treatment Groups",
                   labels = c("Ca. Xc - exposed, Ambient, Red",
                              "Ca. Xc - exposed, Ambient, White",
                              "Ca. Xc - exposed, Heated, Red",
                              "Ca. Xc - exposed, Heated, White",
                              "Naive, Heated, Red",
                              "Naive, Heated, White"))


Mean <- aggregate( CI.g.cm.3.~C_T_S, Day0, mean, na.rm=TRUE )


ComparisonsDay343