# Script for Calculating Respiration Rate
# Modified from N. Miller code (lab protocols Dropbox folder)
# Modified from E. Armstrong 
# J. VanWert 13 Oct 2016

rm(list=ls()) #clear workspace
setwd("~/R/Elysia_stylifera")


library(xlsx)
library(ggplot2)


##### Define Functions ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### 


## Summarizes data.
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}


##### Load in Data ##### ##### ##### ##### #####


A1 = read.xlsx("Amb_Trial1_Day1.xlsx", sheetName="Data")
A1$Date = as.POSIXct(A1$Date, format ="%d.%m.%y %H:%M:%S")
A2 = read.xlsx("Amb_Trial2_Day1.xlsx", sheetName="Data")
A2$Date = as.POSIXct(A2$Date, format ="%d.%m.%y %H:%M:%S")
A3 = read.xlsx("Amb_Trial3_Day1.xlsx", sheetName="Data")
A3$Date = as.POSIXct(A3$Date, format ="%d.%m.%y %H:%M:%S")
H1 = read.xlsx("Heat_Trial1_Day1.xlsx", sheetName="Data")
H1$Date = as.POSIXct(H1$Date, format ="%d.%m.%y %H:%M:%S")
H1S = read.xlsx("Heat_Trial1_Day5.xlsx", sheetName="Data")
H1S$Date = as.POSIXct(H1S$Date, format ="%d.%m.%y %H:%M:%S")
H2 = read.xlsx("Heat_Trial2_Day1.xlsx", sheetName="Data")
H2$Date = as.POSIXct(H2$Date, format ="%d.%m.%y %H:%M:%S")
H2S = read.xlsx("Heat_Trial2_Day5.xlsx", sheetName="Data")
H2S$Date = as.POSIXct(H2S$Date, format ="%d.%m.%y %H:%M:%S")
H3 = read.xlsx("11_15_16_HTrial_R3_Oxygen.xlsx", sheetName="Data")
H3$Date = as.POSIXct(H3$Date, format ="%d.%m.%y %H:%M:%S")
H3S = read.xlsx("11_18_16_HTrialR3_day4_Oxygen.xlsx", sheetName="Data")
H3S$Date = as.POSIXct(H3S$Date, format ="%d.%m.%y %H:%M:%S")
A1S = read.xlsx("Amb_Trial1_Day5.xlsx", sheetName="Data")
A1S$Date = as.POSIXct(A1S$Date, format ="%d.%m.%y %H:%M:%S")
A2S = read.xlsx("Amb_Trial2_Day5.xlsx", sheetName="Data")
A2S$Date = as.POSIXct(A2S$Date, format ="%d.%m.%y %H:%M:%S")
A3S = read.xlsx("Amb_Trial3_Day5.xlsx", sheetName="Data")
A3S$Date = as.POSIXct(A3S$Date, format ="%d.%m.%y %H:%M:%S")



#subset of slope
A1Slope_subset = A1[!(A1$Sec<2*60) & !(A1$Sec>5*60),]
A2Slope_subset = A2[!(A2$Sec<2*60) & !(A2$Sec>5*60),]
A3Slope_subset = A3[!(A3$Sec<2*60) & !(A3$Sec>5*60),]
H1Slope_subset = H1[!(H1$Sec<2*60) & !(H1$Sec>5*60),]
H1SSlope_subset = H1S[!(H1S$Sec<2*60) & !(H1S$Sec>5*60),]
H2Slope_subset = H2[!(H2$Sec<0*60) & !(H2$Sec>5*60),] #started 4 minutes late
H2SSlope_subset = H2S[!(H2S$Sec<2*60) & !(H2S$Sec>5*60),] 
H3Slope_subset = H3[!(H3$Sec<0*60) & !(H3$Sec>215),] #light on during reading at 4 min.. move over 
H3SSlope_subset = H3S[!(H3S$Sec<8*60) & !(H3S$Sec>10*60),] #light temporarily broke, shift time over
A1SSlope_subset = A1S[!(A1S$Sec<2*60) & !(A1S$Sec>5*60),]
A2SSlope_subset = A2S[!(A2S$Sec<2*60) & !(A2S$Sec>5*60),]
A3SSlope_subset = A3S[!(A3S$Sec<2*60) & !(A3S$Sec>5*60),]


## Find Slope of Best Fit Lines For Respiration Data ## ## ## ##
##Analysis of Slopes (O2 Consumption Rates)


#Calculates the slope of Oxy by Sec for Each WellRD in VHAData and returns Dataframe with WellRd, y-Intercept, and Slope (Sec)
library(plyr)

ASlopeData1 <- ddply(A1Slope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})
ASlopeData1= ASlopeData1[-c(3),] #slug too big; started in low O2 state


ASlopeData2 <- ddply(A2Slope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})
ASlopeData2 = ASlopeData2[-c(9),] #B5 slug too large and at anaerobic state (<100 mmol conc O2 entire time)
ASlopeData2 = ASlopeData2[-c(5),] #A6 not included in repeated measures ANOVA bc missing from starved ASSlopeData2 (slug escape)


ASlopeData3 <- ddply(A3Slope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})
ASlopeData3= ASlopeData3[-c(6),] #A3 missing from starved data
ASlopeData3= ASlopeData3[-c(8),] #B3 missing from starved data


ASSlopeData1 <- ddply(A1SSlope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})
ASSlopeData1= ASSlopeData1[-c(3),] #slug too big; started in low O2 state

ASSlopeData2 <- ddply(A2SSlope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})

ASSlopeData3 <- ddply(A3SSlope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})

HSlopeData1 <- ddply(H1Slope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})

HSlopeData1= HSlopeData1[-c(9),] #B3 outlier

HSSlopeData1 <- ddply(H1SSlope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})


HSlopeData2 <- ddply(H2Slope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})

HSSlopeData2 <- ddply(H2SSlope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})

HSlopeData3 <- ddply(H3Slope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})
HSlopeData3= HSlopeData3[-c(5),] #A5 missing from starved trial, escaped


HSSlopeData3 <- ddply(H3SSlope_subset, c("WellID", "WellRd", "Mass","Treatment", "Overall_Treatment", "Temp", "Starved"), function(x) {
  model <- lm(Oxy ~ Sec, data = x) #lm(LogOxy ~ Sec, data = x)
  coef(model)
})



############################################################################################################################




# Subtract blank (A1) from each slope and divide by mass 

ASlopeData1$TrueSlope = (ASlopeData1$Sec - ASlopeData1[1,9])/ASlopeData1$Mass
ASlopeData2$TrueSlope = (ASlopeData2$Sec - ASlopeData2[1,9])/ASlopeData2$Mass
ASlopeData3$TrueSlope = (ASlopeData3$Sec - ASlopeData3[1,9])/ASlopeData3$Mass
ASSlopeData1$TrueSlope = (ASSlopeData1$Sec - ASSlopeData1[1,9])/ASSlopeData1$Mass
ASSlopeData2$TrueSlope = (ASSlopeData2$Sec - ASSlopeData2[1,9])/ASSlopeData2$Mass
ASSlopeData3$TrueSlope = (ASSlopeData3$Sec - ASSlopeData3[1,9])/ASSlopeData3$Mass
HSlopeData1$TrueSlope = (HSlopeData1$Sec - HSlopeData1[1,9])/HSlopeData1$Mass
HSSlopeData1$TrueSlope = (HSSlopeData1$Sec - HSSlopeData1[1,9])/HSSlopeData1$Mass
HSlopeData2$TrueSlope = (HSlopeData2$Sec - HSlopeData2[1,9])/HSlopeData2$Mass
HSSlopeData2$TrueSlope = (HSSlopeData2$Sec - HSSlopeData2[1,9])/HSSlopeData2$Mass
HSlopeData3$TrueSlope = (HSlopeData3$Sec - HSlopeData3[1,9])/HSlopeData3$Mass
HSSlopeData3$TrueSlope = (HSSlopeData3$Sec - HSSlopeData3[1,9])/HSSlopeData3$Mass



#combine slugs of each treatment
ASlopeDataAll_Real = rbind(ASlopeData1, ASlopeData2, ASlopeData3)
ASSlopeDataAll_Real = rbind(ASSlopeData1, ASSlopeData2, ASSlopeData3)
HSlopeDataAll_Real = rbind(HSlopeData1, HSlopeData2, HSlopeData3)
HSSlopeDataAll_Real = rbind(HSSlopeData1, HSSlopeData2, HSSlopeData3)
AmbHighStarvedData = rbind (ASlopeDataAll_Real, HSlopeDataAll_Real, HSSlopeDataAll_Real, ASSlopeDataAll_Real)




#amb day 1, amb day 5, high temp day 1, high temp day 5
#multiply all slopes by -1, divide by 2 to convert to minutes
AmbHighStarvedData$NewTrueSlope = AmbHighStarvedData$TrueSlope*-1
AmbHighStarvedData$RealNewTrueSlope = AmbHighStarvedData$NewTrueSlope/2



#ANOVA and TUKEY 

justdark = subset(AmbHighStarvedData, Overall_Treatment %in% c("Dark"))
require(nlme)
require(multcomp)
justdark$Starved = as.factor(justdark$Starved)
justdark$Temp = as.factor(justdark$Temp)
Lme.dark <- lme(RealNewTrueSlope~Temp*Starved, random=~1|WellID, data=justdark)
summary(Lme.dark)
summary(glht(Lme.dark, linfct=mcp(Temp="Tukey")))

##USED THIS TEST##
#DARK
justdark$together <- interaction(justdark$Temp, justdark$Starved)
model <- lme(RealNewTrueSlope~together, random =~1|WellID/Starved, data=justdark)
summary(model)
summary(glht(model, linfct=mcp(together="Tukey")))


AmbHighStarvedData$Temp = as.factor(AmbHighStarvedData$Temp) #make temperature a factor
summdata = summarySE(justdark, measurevar="RealNewTrueSlope", groupvars=c("Temp","Starved"))


#LIGHT
justlight = subset(AmbHighStarvedData, Overall_Treatment %in% c("Light"))
require(nlme)
require(multcomp)
justlight$Starved = as.factor(justlight$Starved)
justlight$Temp = as.factor(justlight$Temp)
justlight$together <- interaction(justlight$Temp, justlight$Starved)
modellight <-lme(RealNewTrueSlope~together, random =~1|WellID/Starved, data=justlight)
summary(modellight)
summary(glht(modellight, linfct=mcp(together="Tukey")))

#PLOT STARVATION EFFECT
#DARK
dodge= position_dodge(0.2)
plotint = ggplot(summdata, aes(x=Starved, y=RealNewTrueSlope, ymin=RealNewTrueSlope-ci, ymax=RealNewTrueSlope+ci, group=Temp, shape=Temp))+  
  geom_pointrange(size=1.1, position=dodge)+
  geom_line(position=dodge)+
  scale_x_discrete(labels=c("1", "5")) + 
  xlab("Starvation Day") +
  ylab(expression(paste("Oxygen consumption rate (",mu,mol," ", O[2]," ", g^-1, min^-1,")", sep="")))+
  theme_bw() + theme(text = element_text(size=15, face="bold")) +
  theme(panel.grid.major.x = element_blank(), panel.grid.major.y = element_line())+ 
  labs(shape = "Temperature (°C)") +
  annotate("text", x = 1.44, y = 4.4, label = "** p <0.01", size = 6) 
plotint

#LIGHT

justlight = subset(AmbHighStarvedData, Overall_Treatment %in% c("Light"))
summlightdata = summarySE(justlight, measurevar="RealNewTrueSlope", groupvars=c("Temp","Starved"))

dodge= position_dodge(0.2)
plotint = ggplot(summlightdata, aes(x=Starved, y=RealNewTrueSlope, ymin=RealNewTrueSlope-ci, ymax=RealNewTrueSlope+ci, group=Temp, shape=Temp))+  
  geom_pointrange(size=1.1, position=dodge)+
  geom_line(position=dodge)+
  scale_x_discrete(labels=c("1", "5")) + 
  xlab("Starvation Day") +
  ylab(expression(paste("Oxygen consumption rate (",mu,mol," ", O[2]," ", g^-1, min^-1,")", sep="")))+
  theme_bw() + theme(text = element_text(size=15, face="bold")) +
  theme(panel.grid.major.x = element_blank(), panel.grid.major.y = element_line())+ 
  labs(shape = "Temperature (°C)")+
  annotate("text", x = 2.1, y = 5.5, label = "** p <0.01", size = 6)
plotint


##JUST DARK UNDER 4 CONDITIONS
fig2 = ggplot(subset(AmbHighStarvedData, Treatment %in% c("Amb_Dark", "S_Amb_Dark", "HDark", "SHDark"))) + 
  aes(x=Treatment, y=RealNewTrueSlope, fill=factor(Temp))+
  xlab("Starvation Period") +
  ylab(expression(paste("Oxygen consumption rate (",mu,mol," ", O[2]," ", g^-1, min^-1,")", sep="")))+
  scale_x_discrete(limits =c("Amb_Dark", "S_Amb_Dark", "HDark", "SHDark"), labels=c("Day 1", "Day 5","Day 1", "Day 5")) +
  theme_bw() +
  scale_fill_manual(values=c("lightgrey", "grey36"), guide = guide_legend(title = "Temperature  (°C)"))+
  theme(text = element_text(size=15, face="bold")) +
  theme(panel.grid.major.x = element_blank(), panel.grid.major.y = element_line())+
  geom_boxplot()
fig2
#ggtitle("Oxygen consumption under starvation conditions")


#DARK MINUS LIGHT 
justlight = subset(AmbHighStarvedData, Overall_Treatment %in% c("Light"))
justdark = subset(AmbHighStarvedData, Overall_Treatment %in% c("Dark"))

darklight = rbind(justlight, justdark)
summdata = summarySE(darklight, measurevar="RealNewTrueSlope", groupvars=c("Temp","Starved", "Treatment")) #will subtract means (dark - light)

#Create data frame with calculated diff in mean and new CI (google Error Propagation- msu website)
temptreatment = c("28", "28", "30", "30")  
starvationday = c("1", "5", "1", "5")
meandiff = c(4.083936,	1.3688681,	1.7643514,	0.8618287)
cia = sqrt(2.6060736^2 + 0.8725253^2)
cias = sqrt(1.7356738^2 + 1.0163853^2)
ciah = sqrt(1.8022785^2 + 1.1734258^2)
ciahs = sqrt(1.5199312^2 + 1.7833310^2)
newci = c(2.748, 2.011, 2.151, 2.343)
diffdata = data.frame (temptreatment, starvationday, meandiff, newci)


#PLOT RELATIVE NET PHOTOSYNTHESIS
fig3 = ggplot(diffdata, aes(x=starvationday, y=meandiff, ymin = meandiff - newci, ymax = meandiff + newci, group=temptreatment, shape=temptreatment)) +
  geom_pointrange(size=1, position=dodge)+
  geom_line(position=dodge)+
  geom_hline(aes(yintercept=0), color="red", linetype="dashed")+
  # scale_x_discrete(labels=c("Day 1", "Day 5")) + 
  xlab("Starvation Day") +
  ylab(expression(paste("Relative net photosynthesis (",mu,mol," ", O[2]," ", g^-1, min^-1,")", sep="")))+
  theme_bw() + theme(text = element_text(size=15, face="bold")) +
  #theme(panel.grid.major.x = element_blank(), panel.grid.major.y = element_line())
  labs(shape = "Temperature (°C)")+
  scale_color_manual(labels = c("28", "30"), values = c("black", "black"))
fig3