setwd("C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper")

#read in dataset at line 104 ####

#read in temperature dataset ####

temp <- read.csv("all.csv")

summary(temp$Date)
summary(temp$Time)
summary(temp$Day_F)

library(lubridate)
citation("lubridate")

# Extract date values
date <- parse_date_time(temp$Date, "m/d/y") # identify the format of the date
summary(date)

temp$date <- as.factor(date)

summary(temp$date)

#create new variable where dates are summarised into weeks ####
temp$Week <- recode(temp$date, 
                    
       "2016-12-20" = "1", "2016-12-21" = "1" , "2016-12-22" = "1" ,
       
       "2016-12-23" = "2" ,"2016-12-24" = "2", "2016-12-25" = "2", "2016-12-26" = "2", 
       "2016-12-27" = "2", "2016-12-28" = "2", "2016-12-29" = "2" ,
       
       "2016-12-30" = "3", "2016-12-31" = "3", "2017-01-01"= "3", "2017-01-02" = "3", 
       "2017-01-03" = "3", "2017-01-04" = "3", "2017-01-05" = "3", 
       
       "2017-01-06" = "4", "2017-01-07" = "4", "2017-01-08" = "4", "2017-01-09" = "4", 
       "2017-01-10" = "4", "2017-01-11" = "4", "2017-01-12" = "4", 
       
       "2017-01-13" = "5", "2017-01-14" = "5", "2017-01-15" = "5", "2017-01-16" = "5",
       "2017-01-17" = "5", "2017-01-18" = "5", "2017-01-19" = "5", 
       
       "2017-01-20" = "6", "2017-01-21" = "6", "2017-01-22" = "6", "2017-01-23" = "6",
       "2017-01-24" = "6", "2017-01-25" = "6", "2017-01-26" = "6", 
       
       "2017-01-27" = "7", "2017-01-28" = "7", "2017-01-29" = "7", "2017-01-30" = "7",
       "2017-01-31" = "7", "2017-02-01" = "7", "2017-02-02" = "7",
       
       "2017-02-03" = "8", "2017-02-04" = "8", "2017-02-05" = "8", "2017-02-06" = "8",
       "2017-02-07" = "8", "2017-02-08" = "8", "2017-02-09" = "8", 
       
       "2017-02-10" = "9", "2017-02-11" = "9", "2017-02-12" = "9", "2017-02-13" = "9",
       "2017-02-14" = "9", "2017-02-15" = "9", "2017-02-16" = "9", 
       
       "2017-02-17" = "10", "2017-02-18" = "10", "2017-02-19" = "10", "2017-02-20" = "10",
       "2017-02-21" = "10", "2017-02-22" = "10", "2017-02-23" = "10", 
       
       "2017-02-24" = "11", "2017-02-25" = "11", "2017-02-26" = "11", "2017-02-27" = "11",
       "2017-02-28" = "11", "2017-03-01" = "11", "2017-03-02" = "11", 
       
       "2017-03-03" = "12", "2017-03-04" = "12", "2017-03-05" = "12", "2017-03-06" = "12",
       "2017-03-07" = "12", "2017-03-08" = "12", "2017-03-09" = "12", 
       
       "2017-03-10" = "13", "2017-03-11" = "13", "2017-03-12" = "13", "2017-03-13" = "13",
       "2017-03-14" = "13", "2017-03-15" = "13", "2017-03-16" = "13", 
       
       "2017-03-17" = "14", "2017-03-18" = "14", "2017-03-19" = "14", "2017-03-20" = "14",
       "2017-03-21" = "14", "2017-03-22" = "14", "2017-03-23" = "14", 
       
       "2017-03-24" = "15", "2017-03-25" = "15") 
   
summary(temp$Week) 

temp$Week <- as.numeric(temp$Week)

library(dplyr)

#take the mean temperature of each week at each site ####

data_logger_temp_prev_week <- tapply(temp$Temp, list((temp$Week+1), temp$Site), mean, na.rm = T) #creates the variable of mean temperature from previous week
data_logger_temp_prev_week <- as.data.frame(data_logger_temp_prev_week)
data_logger_temp_prev_week$Week <- rownames(data_logger_temp_prev_week)
data_logger_temp_prev_week$Week <- as.numeric(data_logger_temp_prev_week$Week)

#reshape the data so Site is one column ####

library(reshape2)
logger_final <- melt(data_logger_temp_prev_week, id.vars=6)
logger_final

library(dplyr)
logger_final_2 <- logger_final %>% 
  rename(
    Site = variable,
    Temperature = value,
  )

library(plyr)
logger_final_2$Site <-  recode_factor(logger_final_2$Site, "CorA"="Cor A", "CorB"="Cor B", "CorC" = "Cor C", "CorD" = "Cor D", "Hum" = "Humble")
summary(logger_final_2)
logger_final_2$Week <- as.factor(logger_final_2$Week)
logger_final_2$Site <- as.factor(logger_final_2$Site)
summary(logger_final_2)

#write dataset to csv ####
write.csv(logger_final_2, "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/logger_final_2.csv")

#read in csv ####

logger_final_2 <- read.csv("logger_final_2.csv")
library(plyr)
library(dplyr)
logger_final_2 <- logger_final_2 %>% select(-X)

#install.packages("performance") for model comparisons
library(performance)


#install.packages("ggpmisc")


library(ggpmisc) # for stats on line graphs


# protein dataset ####

dframe1 <-read.csv("protein.csv") #choose file - Protein.csv
dframe1 #have a look at the data, if you want to
summary(dframe1) # summarise the data. Check your categorical and cont vars
dframe1$Site <- as.factor(dframe1$Site)

protein_temp <- right_join(dframe1, logger_final_2, by =c("Week", "Site"), na.rm = TRUE)
protein_temp <- na.omit(protein_temp)

protein_temp$Temperature <- as.numeric(protein_temp$Temperature)
na.omit(protein_temp)
summary(protein_temp$Site)

protein_temp<- subset(protein_temp, Site == "Cor A" | Site == "Cor B" | Site == "Cor C" | Site == "Humble")
summary(protein_temp$Site)

summary(protein_temp$Date ~ protein_temp$Week)
protein_temp$NumWeek <- as.numeric(protein_temp$Week)#convert variable Week into a continuous variable
protein_temp$NumBlock <- as.numeric(protein_temp$Block)
protein_temp$Block <- as.factor(protein_temp$Block)

hist(protein_temp$protein) #investigate Protein.content distribution
plot(protein_temp$NumBlock) #check Block has been converted to categorical
plot(protein_temp$Block)

library(ggplot2)
my.formula = y ~ x

ggplot(protein_temp, aes(Site,protein)) + geom_point() + geom_boxplot(notch = TRUE)#plot protein content variation across weeks
plot_protein_temp <- ggplot(protein_temp, aes(Temperature,protein, colour = Site)) + 
  geom_point() + geom_smooth(method = "lm", se = TRUE, formula = my.formula) + 
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +         
  labs(x = "Temperature (°C)", y = expression("Protein content" ~(ug ~ mg^2))) #plot protein contect variation against temperature
plot_protein_temp


my.formula <- y ~ x

library(mgcv)

# the k value in a gam analysis sets how "wiggly" the line is, 
# so you can either choose this manually and compare between models using both the AIC value, gam.check and residuals
# with different k-values (8 is a good place to start) or make k = -1, 
# where the model will automatically select the best K. 

 #remember to select appropriate family for the data you have 
library(mgcv)

model1_null <- gam(protein ~ 1, method = "GCV.Cp",
                na.action = na.exclude, # excludes all NAs
                family=Gamma (link = "identity"), data = protein_temp)

model1 <- gam(protein ~ s(Week, by = Site, k=-1) + Temperature +
              + s(Block, bs = "re") + Total.Wieght + Site, method = "GCV.Cp",
              na.action = na.exclude, # excludes all NAs
              family=Gamma (link = "identity"), data = protein_temp)

model2 <- gam(protein ~ s(NumWeek, k=-1) + Temperature:Site + 
              + s(Block, bs = "re") + Total.Wieght + Site,
              na.action = na.exclude, # excludes all NAs
              family=Gamma (link = "identity"), data = protein_temp) #remember to select appropriate family for the data you have 

plot(protein ~ Total.Wieght, data = protein_temp)

plot(Temperature ~ Week, data = protein_temp)
hist(protein_temp$Temperature)

anova.gam(model1) #adding this in due to reviewer comments

compare_performance(model1_null, model1, model2, rank = TRUE) #comparing the impact of smoothing Week seperated by site
# or not seperated - model 1 is better!

gam.check(model1)
compare_performance(model1) #vary family and link function to determine best model fit
anova.gam(model1) #shows us the impact of Site and Weight on Protein as a whole,
# rather than seperated 

summary(model1) #temperature is not sig.

summary_model <- summary(model1)
write.csv(summary_model$p.table, 
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/protein_gam.csv")
write.csv(summary_model$s.table,
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/protein_gam_smooth.csv")
plot(model1)

gam.fit(model1)
gam.check(model1)
summary.gam(model1)

residuals <- resid(model1, type = "pearson", na.rm = T)
sresid <- (residuals - mean(residuals, na.rm = T))/sd(residuals, na.rm = T)
plot(sresid ~ model1$fitted.values) # you are looking for no obvious pattern/trend
hist(sresid) # looking for a normal distribution of residuals

require("nlme")

#### generation of prediction plots for model ####
plot(model1)

pdat <- expand.grid(Week = seq(min(protein_temp$Week, na.rm = T),max(protein_temp$Week, na.rm = T), length = 201), #variable of interest with appropriate spacing between points, 1 in this case
                    Total.Wieght = mean(protein_temp$Total.Wieght, na.rm = T),
                    Site = c("Cor A", "Cor B", "Cor C", "Humble"),
                    Temperature = mean(protein_temp$Temperature, na.rm = T),
                    Block = sample("Block", 1)) #generate prediction grid using
#across all weeks, with mean weight, across all sites and randomised Block#

preds <- predict(model1, pdat, type="response", se.fit=T) 

pred <- predict(model1, newdata = pdat, type = 'link', se.fit = TRUE,
                exclude = "s(Block)")

pred_df <- cbind(pdat, as.data.frame(pred))

## create confidence interval and back transform
pred_df <- transform(pred_df,
                     fitted_response = fit,
                     fitted_upper    = pred_df$fit + pred_df$se.fit,
                     fitted_lower    = pred_df$fit - pred_df$se.fit)


library('ggplot2')
require(RColorBrewer)
brewer.pal(5, "Set1")

p1 <- ggplot(pred_df, aes(x = Week, y = fitted_response)) +
  geom_ribbon(aes(fill = Site, ymin = fitted_lower,
                  ymax = fitted_upper), alpha = 0.25) +
  geom_line(aes(colour = Site), size = 1) + 
  theme_classic() + scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  labs(x = "Week", y = expression("Protein content" ~(ug ~ mg^2))) +
  geom_point(data = protein_temp, aes(x = NumWeek, y = protein, col = Site), alpha = 0.2, size = 1) + theme_classic(base_size = 20)
p1

# Carbohydrate dataset ####

dframe2 <-read.csv("carbs.csv") #choose file - carbs.csv
dframe2 #have a look at the data, if you want to
summary(dframe2) # summarise the data. Check your categorical and cont vars

carb_temp <- right_join(dframe2, logger_final_2, by =c("Week", "Site"), na.rm = TRUE)
carb_temp <- na.omit(carb_temp)

carb_temp$Temperature <- as.numeric(carb_temp$Temperature)
na.omit(carb_temp$Sample..)
summary(carb_temp)

summary(carb_temp$Site)

carb_temp<- subset(carb_temp, Site == "Cor A" | Site == "Cor B" | Site == "Cor C" | Site == "Humble")
summary(carb_temp$Site)

summary(carb_temp$Date ~ carb_temp$Week)
carb_temp$NumWeek <- as.numeric(carb_temp$Week)#convert variable Week into a continuous variable
carb_temp$NumBlock <- as.numeric(carb_temp$Block)
carb_temp$Block <- as.factor(carb_temp$Block)

hist(carb_temp$carbs) #investigate carb.content distribution
plot(carb_temp$NumBlock) #check Block has been converted to categorical
plot(carb_temp$Block)

library(ggplot2)
ggplot(carb_temp, aes(Week,carbs, colour = Site)) + geom_point() + geom_smooth(method = "lm")#plot carb content variation across weeks
plot_carb_temp <- ggplot(carb_temp, aes(Temperature,carbs, colour = Site)) +
  geom_point() + geom_smooth(method = "lm", se = TRUE, formula = my.formula) + 
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +         
  geom_point() + labs (x = "Temperature (°C)", y = "Carbohydrate content" ~(ug ~ mg^2)) #plot carb contect variation against temperature

plot_carb_temp

library(mgcv)

model3_null <- gam(carbs ~ 1, method = "GCV.Cp",
                  na.action = na.exclude, # excludes all NAs
                  family=gaussian (link = "identity"), data = carb_temp) #remember to select appropriate family for the data you have 

model3 <- gam(carbs ~ s(Week, by = Site, k=-1) + Temperature +
              + s(Block, bs = "re") + Weight + Site, method = "GCV.Cp",
              na.action = na.exclude, # excludes all NAs
              family=gaussian (link = "identity"), data = carb_temp) #remember to select appropriate family for the data you have 

model4 <- gam(carbs ~ s(Week, k=-1) + Temperature +
              + s(Block, bs = "re") + Weight + Site,
              na.action = na.exclude, # excludes all NAs
              family=gaussian (link = "identity"), data = carb_temp) #remember to select appropriate family for the data you have 



compare_performance(model3_null, model3, model4, rank = TRUE) # model3 is better


anova.gam(model3, model4)
compare_performance(model3) #vary family and link function to determine best model fit


summary(model3) #temp not sig.
anova.gam(model3)

summary_model <- summary(model3)
write.csv(summary_model$p.table, 
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/carb_gam.csv")
write.csv(summary_model$s.table,
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/carb_gam_smooth.csv")


plot(model3)
gam.check(model3)


residuals <- resid(model3, type = "pearson", na.rm = T)
sresid <- (residuals - mean(residuals, na.rm = T))/sd(residuals, na.rm = T)
plot(sresid ~ model3$fitted.values) # you are looking for no obvious pattern/trend
hist(sresid) # looking for a normal distribution of residuals

require("nlme")

#### generation of prediction plots for model ####
plot(model3)

pdat2 <- expand.grid(Week = seq(min(carb_temp$Week, na.rm = T),max(carb_temp$Week, na.rm = T), length = 201),
                     Temperature = mean(carb_temp$Temperature, na.rm = T),#variable of interest with appropriate spacing between points, 1 in this case
                    Weight = mean(carb_temp$Weight, na.rm = T),
                    Site = c("Cor A", "Cor B", "Cor C", "Humble"),
                    Block = sample("Block", 1)) #generate prediction grid using
#across all weeks, with mean weight, across all sites and randomised Block#

preds2 <- predict(model3, pdat2, type="response", se.fit=T) 

pred2 <- predict(model3, newdata = pdat2, type = 'link', se.fit = TRUE,
                exclude = "s(Block)")

pred_df2 <- cbind(pdat2, as.data.frame(pred2))

## create confidence interval and back transform
pred_df2 <- transform(pred_df2,
                     fitted_response = fit,
                     fitted_upper    = pred_df2$fit + pred_df2$se.fit,
                     fitted_lower    = pred_df2$fit - pred_df2$se.fit)


library('ggplot2')
require(RColorBrewer)
brewer.pal(5, "Set1")


p2 <- ggplot(pred_df2, aes(x = Week, y = fitted_response)) +
  geom_ribbon(aes(fill = Site, ymin = fitted_lower,
                  ymax = fitted_upper), alpha = 0.25) +
  geom_line(aes(colour = Site), size = 1) + 
  theme_classic() + scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  labs (x = "Week", y = "Carbohydrate content" ~(ug ~ mg^2))+
  geom_point(data = carb_temp, aes(x = NumWeek, y = carbs, col = Site), alpha = 0.2, size = 1) + theme_classic(base_size = 20)
p2

# Lipid dataset ####

dframe3 <-read.csv("lipid.csv") #choose file - lipid.csv
dframe3 #have a look at the data, if you want to
summary(dframe3) # summarise the data. Check your categorical and cont vars

lipid_temp <- right_join(dframe3, logger_final_2, by =c("Week", "Site"), na.rm = TRUE)
lipid_temp <- na.omit(lipid_temp)

lipid_temp$Temperature <- as.numeric(lipid_temp$Temperature)
na.omit(lipid_temp$Sample..)
summary(lipid_temp)

summary(lipid_temp$Site)

lipid_temp<- subset(lipid_temp, Site == "Cor A" | Site == "Cor B" | Site == "Cor C" | Site == "Humble")
summary(lipid_temp$Site)

summary(lipid_temp$Date ~ lipid_temp$Week)
lipid_temp$NumWeek <- as.numeric(lipid_temp$Week)#convert variable Week into a continuous variable
lipid_temp$NumBlock <- as.numeric(lipid_temp$Block)
lipid_temp$Block <- as.factor(lipid_temp$Block)

hist(lipid_temp$lipid) #investigate lipid.content distribution
plot(lipid_temp$NumBlock) #check Block has been converted to categorical
plot(lipid_temp$Block)

library(ggplot2)
ggplot(lipid_temp, aes(Week,lipid, colour = Site)) + geom_point() + geom_smooth(method = "lm")#plot lipid content variation across weeks
plot_lipid_temp <- ggplot(lipid_temp, aes(Temperature,lipid, colour = Site)) + 
  geom_point() + geom_point() + geom_smooth(method = "lm", se = TRUE, formula = my.formula) + 
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) + labs (x = "Temperature (°C)", y = "Lipid content" ~(ug ~ mg^2)) #plot lipid contect variation against temperature

plot_lipid_temp
library(mgcv)

model5_null <- gam(lipid ~ 1,
                   na.action = na.exclude, # excludes all NAs
                   family=Gamma (link = "inverse"), data = lipid_temp)

model5 <- gam(lipid ~ s(NumWeek, by = Site, k=-1) + Temperature +
              + s(Block, bs = "re") + Weight + Site,
              na.action = na.exclude, # excludes all NAs
              family=Gamma (link = "inverse"), data = lipid_temp) #remember to select appropriate family for the data you have 

model6 <- gam(lipid ~ s(Week, k=-1) + Temperature  +
              + s(Block, bs = "re") + Weight + Site,
              na.action = na.exclude, # excludes all NAs
              family=Gamma (link = "inverse"), data = lipid_temp) #remember to select appropriate family for the data you have 


compare_performance(model5_null,model5, model6, rank = TRUE) # model 6 is better - so we will NOT smooth - but it's a close one!
#week by site - keeping the overall model more simple
anova.gam(model6, model7)

compare_performance(model6) #vary family and link function to determine best model fit
summary(model6) #temp not sig. 
anova.gam(model6)

summary_model <- summary(model6)
write.csv(summary_model$p.table, 
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/lipid_gam.csv")
write.csv(summary_model$s.table,
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/lipid_gam_smooth.csv")
plot(model2)
gam.check(model2)

residuals <- resid(model6, type = "pearson", na.rm = T)
sresid <- (residuals - mean(residuals, na.rm = T))/sd(residuals, na.rm = T)
plot(sresid ~ model6$fitted.values) # you are looking for no obvious pattern/trend
hist(sresid) # looking for a normal distribution of residuals

#### plots for model ####

pdat3 <- expand.grid(Week = seq(min(lipid_temp$Week, na.rm = T),max(lipid_temp$Week, na.rm = T), length = 201), #variable of interest with appropriate spacing between points, 1 in this case
                    Temperature = mean(lipid_temp$Temperature, na.rm = T),
                     Weight = mean(lipid_temp$Weight, na.rm = T),
                    Site = c("Cor A", "Cor B", "Cor C", "Humble"),
                    Block = sample("Block", 1)) #
 #generate prediction grid using
#across all weeks, with mean weight, across all sites

preds3 <- predict(model6, pdat3, type="response", se.fit=T) 

pred3 <- predict(model6, newdata = pdat3, type = 'link', se.fit = TRUE,
                exclude = "s(Block)")

pred_df3 <- cbind(pdat3, as.data.frame(preds3))

## create confidence interval and back transform
pred_df3 <- transform(pred_df3,
                     fitted_response = fit,
                     fitted_upper    = pred_df3$fit + pred_df3$se.fit,
                     fitted_lower    = pred_df3$fit - pred_df3$se.fit)


lipid <- ggplot(pred_df3, aes(factor(Site), fit, fill = Site)) + 
  geom_boxplot(alpha = 0.5) + 
  theme_classic() + scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  labs (x = "Site", y = "Lipid content" ~(ug ~ mg^2))+
  geom_point(data = lipid_temp, aes(x = Site, y = lipid, col = Site), alpha = 0.2, size = 1)  + theme_classic(base_size = 20)
lipid


pdat4 <- expand.grid(Weight = seq(min(lipid_temp$Weight, na.rm = T),max(lipid_temp$Weight, na.rm = T), length = 201), #variable of interest with appropriate spacing between points, 1 in this case
                    Week = 16,
                    Temperature = mean(lipid_temp$Temperature, na.rm = T),
                    Site = "Cor A",
                    Block = sample("Block", 1)) #
#generate prediction grid using
#across all weeks, with mean weight, across all sites

preds4 <- predict(model6, pdat4, type="response", se.fit=T) 

pred4 <- predict(model6, newdata = pdat4, type = 'link', se.fit = TRUE,
                exclude = "s(Block)")

pred_df4 <- cbind(pdat4, as.data.frame(preds4))
pred_df4 <- transform(pred_df4,
                      fitted_response = fit,
                      fitted_upper    = pred_df4$fit + pred_df4$se.fit,
                      fitted_lower    = pred_df4$fit - pred_df4$se.fit)


weight <- ggplot(pred_df4, aes(x = Weight, y = fitted_response)) +
  geom_ribbon(aes(ymin = fitted_upper,
                  ymax = fitted_lower), alpha = 0.25, fill = "purple") +
  geom_line(aes(), size = 1, col = "purple") + 
  theme_classic() + scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  labs (x = "Mass (mg)", y = "Lipid content" ~(ug ~ mg^2))+
  geom_point(data = lipid_temp, aes(x = Weight, y = lipid), alpha = 0.2, size = 1) + theme_classic(base_size = 20)
weight

library(ggpubr)
figure <- ggarrange(lipid, weight,
                    labels = c("A", "B"),
                    ncol = 1, nrow = 2)
figure

# MR dataset ####

dframe4 <-read.csv("MR.csv") #choose file - MR.csv
dframe4 #have a look at the data, if you want to
summary(dframe4) # summarise the data. Check your categorical and cont vars

MR_temp <- right_join(dframe4, logger_final_2, by =c("Week", "Site"), na.rm = TRUE)
MR_temp <- na.omit(MR_temp)

MR_temp$Temperature <- as.numeric(MR_temp$Temperature)
na.omit(MR_temp$Sample..)
summary(MR_temp)

summary(MR_temp$Site)

MR_temp<- subset(MR_temp, Site == "Cor A" | Site == "Cor B" | Site == "Cor C" | Site == "Humble")
summary(MR_temp$Site)


summary(MR_temp$Date ~ MR_temp$Week)
MR_temp$NumWeek <- as.numeric(MR_temp$Week)#convert variable Week into a continuous variable


hist(MR_temp$MR) #investigate MR.content distribution

library(ggplot2)
ggplot(MR_temp, aes(Week,MR, colour = Site)) + geom_point() + geom_smooth(method = "lm")#plot MR content variation across weeks
plot_MR_temp <- ggplot(MR_temp, aes(Temperature,MR, colour = Site)) + ylim(0,5) +
  geom_point() + geom_point() + geom_smooth(method = "lm", se = TRUE, formula = my.formula) + 
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +labs (x = "Temperature (°C)", y = "Metabolic Rate" ~(ul ~O[2] ~ h^-1 ~mg^-1)) #plot MR contect variation against temperature

plot_MR_temp

library(mgcv)
model7_null <- gam(MR ~ 1,
                   na.action = na.exclude, # excludes all NAs
                   family=Gamma (link = "inverse"), data = MR_temp)


model7 <- gam(MR ~ s(Week, by = Site, k=-1) + Temperature +
              + Weight + Site,
              na.action = na.exclude, # excludes all NAs
              family=Gamma (link = "inverse"), data = MR_temp) #remember to select appropriate family for the data you have 

model8 <- gam(MR ~ s(Week, k=-1) +
                + Weight + Site,
              na.action = na.exclude, # excludes all NAs
              family=Gamma (link = "inverse"), data = MR_temp) #remember to select appropriate family for the data you have 

compare_performance(model7_null,model7, model8, rank = TRUE) # model 7 is much better - so we WILL smooth 
anova.gam(model7, model8)

#week by site 
compare_performance(model7, model7.1, rank = TRUE) #vary family and link function to determine best model fit
summary(model7) #temp not sig. 
anova.gam(model7)

summary_model <- summary(model7)
write.csv(summary_model$p.table, 
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/MR_gam.csv")
write.csv(summary_model$s.table,
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/MR_gam_smooth.csv")


plot(model7)
gam.check(model7)

residuals <- resid(model7, type = "pearson", na.rm = T)
sresid <- (residuals - mean(residuals, na.rm = T))/sd(residuals, na.rm = T)
plot(sresid ~ model7$fitted.values) # you are looking for no obvious pattern/trend
hist(sresid) # looking for a normal distribution of residuals

#### generation of prediction plots for model ####
plot(model1)

pdat5 <- expand.grid(Week = seq(min(MR_temp$Week, na.rm = T),max(MR_temp$Week, na.rm = T), length = 201), #variable of interest with appropriate spacing between points, 1 in this case
                    Weight = mean(MR_temp$Weight, na.rm = T),
                    Temperature = mean(MR_temp$Temperature, na.rm = T),
                    Site = c("Cor A", "Cor B", "Cor C", "Humble")
                    ) #generate prediction grid using
#across all weeks, with mean weight, across all sites

preds5 <- predict(model7, pdat5, type="response", se.fit=T) 

pred5 <- predict(model7, newdata = pdat5, type = 'link', se.fit = TRUE,
                exclude = "s(Block)")

pred_df5 <- cbind(pdat5, as.data.frame(preds5))

## create confidence interval and back transform
pred_df5 <- transform(pred_df5,
                     fitted_response = fit,
                     fitted_upper    = pred_df5$fit + pred_df5$se.fit,
                     fitted_lower    = pred_df5$fit - pred_df5$se.fit)



library('ggplot2')
require(RColorBrewer)
brewer.pal(5, "Set1")

p4 <- ggplot(pred_df5, aes(x = Week, y = fitted_response)) +
  geom_ribbon(aes(fill = Site, ymin = fitted_upper,
                  ymax = fitted_lower), alpha = 0.25) +
  geom_line(aes(colour = Site), size = 1) + 
  theme_classic() + scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  labs (x = "Week", y = "Metabolic Rate" ~(ul ~O[2] ~ h^-1 ~mg^-1))+
  geom_point(data = dframe4, aes(x = Week, y = MR, col = Site), alpha = 0.2, size = 1) + theme_classic(base_size = 20)
p4




# Weight dataset uses MR dataset####

plot(Weight ~ Site, data=MR_temp, notch = TRUE) #plot MR content variation across weeks
plot(Weight ~ Week, data=MR_temp)
hist(MR_temp$Weight)
plot(MR_temp$Weight)

MR_temp$Weight <-as.numeric(MR_temp$Weight)

library(ggplot2)
plot_weight_temp <- ggplot(MR_temp, aes(Temperature, Weight, colour = Site)) + ylim(0,0.65) +
  
  geom_point() + geom_smooth(method = "lm", se = TRUE, formula = my.formula) + 
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) + geom_point() + labs (x = "Temperature (°C)", y = "Mass (mg)")
plot_weight_temp

library(mgcv)

model9_null <- gam(Weight ~ 1,
                   na.action = na.exclude, # excludes all NAs
                   family=Gamma(link = "identity"), data = MR_temp) 

model9 <- gam(Weight ~ s(Week, by = Site, k=-1) + Temperature +
                Site,
              na.action = na.exclude, # excludes all NAs
              family=Gamma(link = "identity"), data = MR_temp) #remember to select appropriate family for the data you have 

model10 <- gam(Weight ~ s(Week, k=-1) + Temperature +
              + Site,
              na.action = na.exclude, # excludes all NAs
              family=Gamma (link = "identity"), data = MR_temp) #remember to select appropriate family for the data you have 

compare_performance(model9_null, model9, model10, rank = TRUE) # model 9 is much better - so we will  smooth 

anova.gam(model9, model10)

#week by site 
compare_performance(model9) #vary family and link function to determine best model fit

summary(model9) #temperature has a significantly negative correlation with overall weight i.e. as temperature goes up, weight goes down

summary_model <- summary(model9)
write.csv(summary_model$p.table, 
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/mass_gam.csv")
write.csv(summary_model$s.table,
          file = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/mass_gam_smooth.csv")

anova.gam(model9)
plot(model9)
gam.check(model9)

residuals <- resid(model9, type = "pearson", na.rm = T)
sresid <- (residuals - mean(residuals, na.rm = T))/sd(residuals, na.rm = T)
plot(sresid ~ model9$fitted.values) # you are looking for no obvious pattern/trend
hist(sresid) # looking for a normal distribution of residuals

#### generation of prediction plots for model ####
plot(model9)

pdat6 <- expand.grid(Week = seq(min(MR_temp$Week, na.rm = T),max(MR_temp$Week, na.rm = T)), #variable of interest with appropriate spacing between points, 1 in this case
                    Site = c("Cor A", "Cor B", "Cor C", "Humble"),
                    Temperature = mean(MR_temp$Temperature, na.rm = T)) #generate prediction grid using

#across all weeks, with mean weight, across all sites

preds6 <- predict(model9, pdat6, type="response", se.fit=T) 

pred6 <- predict(model9, newdata = pdat6, type = 'link', se.fit = TRUE)

pred_df6 <- cbind(pdat6, as.data.frame(preds6))

## create confidence interval and back transform
pred_df6 <- transform(pred_df6,
                     fitted_response = fit,
                     fitted_upper    = pred_df6$fit + pred_df6$se.fit,
                     fitted_lower    = pred_df6$fit - pred_df6$se.fit)

library('ggplot2')
require(RColorBrewer)
brewer.pal(5, "Set1")

p5 <- ggplot(pred_df6, aes(x = Week, y = fitted_response)) +
  geom_ribbon(aes(fill = Site, ymin = fitted_lower,
                  ymax = fitted_upper), alpha = 0.25) +
  geom_line(aes(colour = Site), size = 1) + 
  theme_classic() + scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  labs (x = "Week", y = "Mass (mg)")+
  geom_point(data = MR_temp, aes(x = Week, y = Weight, col = Site), alpha = 0.2, size = 1) + theme_classic(base_size = 20)
p5


#install.packages("ggpubr")
library(ggpubr)

figure_minus_corD <- ggarrange(p5, p4, p2, lipid, weight, p1,
                   ncol = 2, nrow = 3, labels = "AUTO", common.legend = TRUE,font.label = list(size = 30, face = "bold", color ="black") )
figure_minus_corD

ggsave(
  "figure_s2_minus_CorD.tiff",
  plot = figure_minus_corD,
  device = "tiff",
  path = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/",
  scale = 1,
  width = 450,
  height = 450,
  units = c("mm"),
  dpi = 600,
  limitsize = TRUE)

temp_week_plot <- ggplot(protein_temp, aes(Week, Temperature, colour = Site)) + geom_point() + geom_line()
temp_week_plot

temp_variables <- ggarrange(temp_week_plot, plot_protein_temp, plot_carb_temp, plot_lipid_temp, plot_MR_temp, plot_weight_temp,
                            labels = "AUTO", common.legend = TRUE,font.label = list(size = 16, face = "bold", color ="black"))
temp_variables

ggsave(
  "temp_variables.tiff",
  plot = temp_variables,
  device = "tiff",
  path = "C:/Users/jackj/OneDrive - University of Kentucky/PhD_research/gam_analysis_for_drew_paper/new_stats_results/",
  scale = 1,
  width = 450,
  height = 350,
  units = c("mm"),
  dpi = 600,
  limitsize = TRUE)