# Peer J
# Kinsey M Brock et al. 
# Code for I. Color, II. Bite Force, and III. Chemical Analyses

# I. COLOR ##############################################################################################

# Female Color Data Prep - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  #
library(pavo)
set.seed(1612217)
# Read in the data:
specs_OGfem <- getspec("~/path_name/2017_morph_fem_spec", ext = "txt", decimal = ",", subdir = FALSE, subdir.names = FALSE)

specs.orangefem <- subset(specs_OGfem, "_OO")
head(specs.orangefem)[1:8]

specs.whitefem <- subset(specs_OGfem, "_WW")
head(specs.whitefem)[1:8]

specs.yellowfem <- subset(specs_OGfem, "_YY")
head(specs.yellowfem)[1:8]

mu.femspecs <- aggspec(specs_OGfem, by = 6, FUN = mean)
mu.femspecs[1:5, 1:5]
dim(mu.femspecs) 

femmorph <- gsub("\\d", '', names(mu.femspecs))[-1]

mu.femmorphspec <- aggspec(mu.femspecs, by = femmorph, FUN = mean)
round(mu.femmorphspec[1:5, ], 2)

femmorphspec.sm <- procspec(mu.femmorphspec, opt = "smooth", span = 0.2)

femspec.sm <- procspec(mu.femspecs, opt = "smooth", span = 0.2)

femspec.sm.min <- procspec(femspec.sm, opt = "min")

femspec.bin <- procspec(mu.femspecs, opt = c("bin", "center"))

colr <- spec2rgb(mu.femspecs)
wls <- as.numeric(colnames(femspec.bin))

aggplot(femspec.sm, femmorph, 
        FUN.center = median,
        ylim = c(0, 100),
        xlim = c(300, 700),
        alpha = 0.4, legend = TRUE, lcol = "black", lty = (values=c("solid", "dotted", "longdash")), shadecol = (values=c("yellow", "gray", "darkorange")))

write.csv(t(as.data.frame(summary(femspec.sm, subset = FALSE, wlmin = NULL,
                                  wlmax = NULL))), file="dfacol_femmale2017.csv") 

#  Discriminant Function Analysis of Morph Types (Females 2017) - - - - - - - - - - - - - - - - - - - - #
library(tidyverse)
library(caret)

data.col2017f <- read.csv("dfacol_female2017.csv")

#Identify and remove colinear variables:
col.model <- lm(morph ~ B1 + B2 + B3 + S1U + S1V + S1B + S1G + S1Y + S1R + S2 + S3 + S4 + 
                  S5 + S6 + S7 + S8 + S9 + S10 + H1 + H2 + H3 + H4 + H5, data = data.col2017f)
summary(col.model)

X<-data.col2017f[,2:23]
library(GGally)
ggpairs(X)

data <- read.csv("dfa_female2017.csv")

library(pander)
library(magrittr)
library(dplyr)

data %>% head %>% pander

#data w descriptives
library(tables)
tabular( B2 + B3 + S1U + S1Y + S1B + S1G + S1R + S6 + S7 + H1 
         + H2 ~ Heading() * morph * (Heading(Mean) * mean + Heading(SD) *sd) * Format(sprintf("%.2f")), 
         data = data) %>% pander(style = "multiline")

#lda
library(MASS)
m <- lda(morph ~ B2 + B3 + S1U + S1Y + S1B + S1G + S1R + S6 + S7 + H1, data)
m

#observed vs. predicted morph categories
p <- predict(m)
freqtable <- table(p$class, data$morph) 
rownames(freqtable) <- paste0("Predicted ", data$morph %>% levels)
freqtable %>% addmargins %>% pander("Observed vs. Predicted Frequencies")

#proportions
prop.table(freqtable) %>% addmargins %>% pander("Proportions")

#discriminant functions
data.frame(p$posterior, 
           ObservedGroup = data$morph, 
           PredictedGroup = p$class, p$x) %>% 
  head(20) %>% pander

data$LDA1 <- p$x[,1]
data$LDA2 <- p$x[,2]

ggplot(data, aes(LDA1, fill = morph)) + geom_density(alpha = 0.5, color = NA) + theme(legend.position = "top")

ggplot(data, aes(LDA2, fill = morph)) + geom_density(alpha = 0.5, color = NA) + theme(legend.position = "top")

#plot
dev.off()
plot(m, abbrev = 1)

#manova follow-up
#tutorial here: https://my.ilstu.edu/~wjschne/444/MANOVA.html#(6) 
mm <- lm(cbind(B2,B3,S1U,S1Y,S1B,S1G,S1R,S6,S7,H1) ~ morph, data)
mm
mm %>% manova %>% summary(test = "Wilks")

#visualize
ggplot(data, aes(LDA1, LDA2, color = morph)) + 
  geom_point() + scale_color_manual(breaks = c("orange", "white", "yellow"),
                                    values=c("orange", "black", "yellow")) + stat_ellipse()

ggplot(data, aes(LDA1, LDA2, color = morph)) + 
  geom_point(pch=21, colour = "black", size = 2) + 
  scale_color_manual(breaks = c("orange", "white", "yellow"), values=c("orange", "black", "yellow")) + 
  stat_ellipse(size = 1)

myplot <- ggplot(data, aes(LDA1, LDA2, color = morph)) + 
  geom_point(pch=21, colour = "black", size = 2) + 
  scale_color_manual(breaks = c("orange", "white", "yellow"), values=c("orange", "black", "yellow")) + 
  stat_ellipse(size = 1)

myplot + theme(panel.background = element_blank(), 
               axis.line = element_line(colour = "black"), panel.grid.major = element_line(colour = "grey92"))

almostdone <- myplot + theme(panel.background = element_blank(), 
                             axis.line = element_line(colour = "black"))


almostdone + theme(axis.text.x = element_text(size = 15),
                   axis.text.y = element_text(size = 15),  
                   axis.title.x = element_text(size = 15),
                   axis.title.y = element_text(size = 15))

finalLDFA <- almostdone + theme(axis.text.x = element_text(size = 15),
                                axis.text.y = element_text(size = 15),  
                                axis.title.x = element_text(size = 15),
                                axis.title.y = element_text(size = 15))

finalLDFA +ggtitle("Linear Discriminant Analysis - 2017 Females") +
  theme(plot.title = element_text(hjust = 0.5))

# K means cluster analysis - FEMALE 2017 DATA - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Load the data:
data.col.f17 <- read.csv("dfacol_female2017vars.csv")

set.seed(788)

str(data.col.f17)
summary(data.col.f17)
any(is.na(data.col.f17))
# No NAs!

# Store the labels in a separate variable and exclude the "morph" column from your dataset in order to do clustering.
morph_label_f17 <- data.col.f17$morph
data.col.f17$morph <- NULL
str(data.col.f17)
# OK, morph is gone

# Also need to do this for first column of dataframe with names stored in "X"
ID_label_f17 <- data.col.f17$X
data.col.f17$X <- NULL
str(data.col.f17)

# Use R's scale() function to scale all your column values because we have 10 vars on different scales.
color_df_sc_f17 <- as.data.frame(scale(data.col.f17))
summary(color_df_sc_f17)
# OK, now everything has mean of 0 and STDEV of 1.

### Gap statistic
set.seed(123)
# Compute gap statistic for kmeans
# Recommended B value is 500
gap_stat_f <- clusGap(color_df_sc_f17, FUN = kmeans, nstart = 25,
                      K.max = 6, B = 500)
print(gap_stat_f, method = "firstmax") #3
print(gap_stat_f, method = "globalmax") #3
fviz_gap_stat(gap_stat_f, maxSE = list(method = "firstmax")) #3
fviz_gap_stat(gap_stat_f, maxSE = list(method = "globalmax")) #3

# Male Color Data Prep - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -#
library(pavo)
set.seed(1612217)

specs_OGmale <- getspec("~/your_path/2017_morph_male_spec", ext = "txt", decimal = ",", subdir = FALSE, subdir.names = FALSE)

specs.orange <- subset(specs_OGmale, "_OO")
head(specs.orange)[1:8]

specs.white <- subset(specs_OGmale, "_WW")
head(specs.white)[1:8]

specs.yellow <- subset(specs_OGmale, "_YY")
head(specs.yellow)[1:8]

mu.specs <- aggspec(specs_OGmale, by = 6, FUN = mean)
mu.specs[1:5, 1:5]
dim(mu.specs) 

morph <- gsub("\\d", '', names(mu.specs))[-1]
table(morph)

mu.morphspec <- aggspec(mu.specs, by = morph, FUN = mean)
round(mu.morphspec[1:5, ], 2)

morphspec.sm <- procspec(mu.morphspec, opt = "smooth", span = 0.2)

spec.sm <- procspec(mu.specs, opt = "smooth", span = 0.2)

library(dplyr)
OO.spec.sm <- select(spec.sm, ends_with("_OO"))
par(mfrow = c(1,1))
plot(OO.spec.sm, col = spec2rgb(spec.sm))

YY.spec.sm <- select(spec.sm, ends_with("_YY"))
par(mfrow = c(1,1))
plot(YY.spec.sm, col = spec2rgb(spec.sm))

WW.spec.sm <- select(spec.sm, ends_with("_WW"))
par(mfrow = c(1,1))
plot(WW.spec.sm, col = spec2rgb(spec.sm))

spec.sm.min <- procspec(spec.sm, opt = "min")

spec.bin <- procspec(mu.specs, opt = c("bin", "center"))

spec.bin <- t(spec.bin) # transpose so wavelength are variables for the PCA
colnames(spec.bin) <- spec.bin[1, ]  # names variables as wavelength bins
spec.bin <- spec.bin[-1, ]  # remove 'wl' column
pca1 <- prcomp(spec.bin, scale. = TRUE)

summary(pca1)

spec.sm.bin <- procspec(spec.sm, opt = c("bin", "center"))

spec.sm.bin <- t(spec.sm.bin) # transpose so wavelength are variables for the PCA
colnames(spec.sm.bin) <- spec.sm.bin[1, ]  # names variables as wavelength bins
spec.sm.bin <- spec.sm.bin[-1, ]  # remove 'wl' column
pca1.spec.sm <- prcomp(spec.sm.bin, scale. = TRUE)

summary(pca1.spec.sm)

colr.sm <- spec2rgb(spec.sm)
wls.sm <- as.numeric(colnames(spec.sm.bin))

aggplot(spec.sm, morph, 
        FUN.center = median,
        ylim = c(0, 100),
        alpha = 0.4, legend = TRUE, lcol = "black", lty = (values=c("longdash", "dotted", "solid")), shadecol = (values=c("darkorange", "gray", "yellow")))

summary(spec.sm, subset = FALSE, wlmin = NULL,
        wlmax = NULL)

write.csv(t(as.data.frame(summary(spec.sm, subset = FALSE, wlmin = NULL,
                                  wlmax = NULL))), file="dfacol_male2017.csv") 

#  Discriminant Function Analysis of Morph Types (Males 2017) - - - - - - - - - - - - - - - - - - - - - #
library(tidyverse)
library(caret)

data.col2017m <- read.csv("dfacol_male2017.csv") #smoothed data from Part 1, "morph" replaced w/ number 1,2,3
head(data.col2017m)

col.model <- lm(morph ~ B1 + B2 + B3 + S1U + S1V + S1B + S1G + S1Y + S1R + S2 + S3 + S4 + 
                  S5 + S6 + S7 + S8 + S9 + S10 + H1 + H2 + H3 + H4 + H5, data = data.col2017m)
summary(col.model)

X<-data.col2017m[,2:23]
library(GGally)
ggpairs(X)

data <- read.csv("dfa_male2017.csv")

set.seed(123)
training.samples <- data$morph %>%
  createDataPartition(p = 0.8, list = FALSE)
train.data <- data[training.samples, ]
test.data <- data[-training.samples, ]

train.data # N = 80 lizards
test.data # N = 16 lizards

# 2. Normalize the data. Categorical variables are automatically ignored:
# Estime preprocessing parameters:
preproc.param <- train.data %>% 
  preProcess(method = c("center", "scale"))
# Transform the data using the estimated parameters
train.transformed <- preproc.param %>% predict(train.data)
test.transformed <- preproc.param %>% predict(test.data)

# Linear discriminant function analysis - LDA
train.transformed[-1]

library(MASS)
# Fit the model (needed to add [-1] for this to work: 
#model1
model1 <- lda(morph~., data = train.transformed[-1]) # variables are colinear and violate assumptions of LDA
#model2 
model2 <- lda(morph ~ B2 + B3 + S1U + S1Y + S1B + S1G + S1R + S6 + S7 + H1 + H2, data = train.transformed[-1])
### Proceed w model 2 - uses all significant vars and eliminates collinear vars

predictions <- model2 %>% predict(test.transformed)

mean(predictions$class==test.transformed$morph)

model2

predictions <- model2 %>% predict(test.transformed)
names(predictions)

lda.data <- cbind(train.transformed, predict(model2)$x)
ggplot(lda.data, aes(LD1, LD2)) +
  geom_point(aes(color = morph))

library(pander)
library(magrittr)
library(dplyr)

data %>% head %>% pander

#data w descriptives
library(tables)
tabular( B2 + B3 + S1U + S1Y + S1B + S1G + S1R + S6 + S7 + H1 
         ~ Heading() * morph * (Heading(Mean) * mean + Heading(SD) *sd) * Format(sprintf("%.2f")), 
         data = data) %>% pander(style = "multiline")

#lda
library(MASS)
m <- lda(morph ~ B2 + B3 + S1U + S1Y + S1B + S1G + S1R + S6 + S7 + H1, data)
m

#observed vs. predicted morph categories
p <- predict(m)
freqtable <- table(p$class, data$morph) 
rownames(freqtable) <- paste0("Predicted ", data$morph %>% levels)
freqtable %>% addmargins %>% pander("Observed vs. Predicted Frequencies")

#proportions
prop.table(freqtable) %>% addmargins %>% pander("Proportions")

#discriminant functions
data.frame(p$posterior, 
           ObservedGroup = data$morph, 
           PredictedGroup = p$class, p$x) %>% 
  head(20) %>% pander

data$LDA1 <- p$x[,1]
data$LDA2 <- p$x[,2]

ggplot(data, aes(LDA1, fill = morph)) + geom_density(alpha = 0.5, color = NA) + theme(legend.position = "top")

ggplot(data, aes(LDA2, fill = morph)) + geom_density(alpha = 0.5, color = NA) + theme(legend.position = "top")

#plot
dev.off()
plot(m, abbrev = 1)

#manova follow-up
#tutorial here: https://my.ilstu.edu/~wjschne/444/MANOVA.html#(6) 
mm <- lm(cbind(B2,B3,S1U,S1Y,S1B,S1G,S1R,S6,S7,H1) ~ morph, data)
mm
mm %>% manova %>% summary(test = "Wilks")

#visualize
ggplot(data, aes(LDA1, LDA2, color = morph)) + 
  geom_point() + scale_color_manual(breaks = c("orange", "white", "yellow"),
                                    values=c("orange", "black", "yellow")) + stat_ellipse()

ggplot(data, aes(LDA1, LDA2, color = morph)) + 
  geom_point(pch=21, colour = "black", size = 2) + 
  scale_color_manual(breaks = c("orange", "white", "yellow"), values=c("orange", "black", "yellow")) + 
  stat_ellipse(size = 1)

myplot <- ggplot(data, aes(LDA1, LDA2, color = morph)) + 
  geom_point(pch=21, colour = "black", size = 2) + 
  scale_color_manual(breaks = c("orange", "white", "yellow"), values=c("orange", "black", "yellow")) + 
  stat_ellipse(size = 1)

myplot + theme(panel.background = element_blank(), 
               axis.line = element_line(colour = "black"), panel.grid.major = element_line(colour = "grey92"))

almostdone <- myplot + theme(panel.background = element_blank(), 
                             axis.line = element_line(colour = "black"))


almostdone + theme(axis.text.x = element_text(size = 15),
                   axis.text.y = element_text(size = 15),  
                   axis.title.x = element_text(size = 15),
                   axis.title.y = element_text(size = 15))

finalLDFA <- almostdone + theme(axis.text.x = element_text(size = 15),
                                axis.text.y = element_text(size = 15),  
                                axis.title.x = element_text(size = 15),
                                axis.title.y = element_text(size = 15))

finalLDFA +ggtitle("Linear Discriminant Analysis - 2017 Males") +
  theme(plot.title = element_text(hjust = 0.5))

# K means cluster analysis - MALE 2017 DATA - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Load the data:
data.col.m17 <- read.csv("dfacol_male2017vars.csv")

set.seed(788)

str(data.col.m17)
summary(data.col.m17)
any(is.na(data.col.m17))
# No NAs!

# Store the labels in a separate variable and exclude the "morph" column from your dataset in order to do clustering w/ no prior morph information.
morph_label_m17 <- data.col.m17$morph
data.col.m17$morph <- NULL
str(data.col.m17)
# OK, morph is gone

# Also need to do this for first column of dataframe with names stored in "X"
ID_label_m17 <- data.col.m17$X
data.col.m17$X <- NULL
str(data.col.m17)

# Use R's scale() function to scale all your column values because we have 10 vars on very different scales.
color_df_sc_m17 <- as.data.frame(scale(data.col.m17))
summary(color_df_sc_m17)
# OK, now everything has mean of 0 and STDEV of 1.

### Gap statistic
library(factoextra)
library(cluster)
set.seed(123)
# Compute gap statistic for kmeans
# Recommended value is B = 500
gap_stat <- clusGap(color_df_sc_m17, FUN = kmeans, nstart = 25,
                    K.max = 6, B = 500)
print(gap_stat, method = "firstmax") #3
print(gap_stat, method = "globalmax") #3
fviz_gap_stat(gap_stat) #3
fviz_gap_stat(gap_stat, maxSE = list(method = "globalSEmax")) #3
fviz_gap_stat(gap_stat, maxSE = list(method = "firstmax")) #3
fviz_gap_stat(gap_stat, maxSE = list(method = "globalmax")) #3
fviz_gap_stat(gap_stat, maxSE = list(method = "Tibs2001SEmax")) #3

# II. BITE FORCE #########################################################################################
library(effects)
library(MASS)
library(ape)
library(car)
library(ggplot2)
library(pander)
library(lsmeans)
library(lme4)
library(scales)
library(Hmisc)
library(vegan)
library(moments)
library(ISLR)
library(corrplot)

# MALE BITE FORCE ANALYSES - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
dat2017v2 <-read.csv("2017_male_biteforce_v2.csv")

# ANOVAs morphology
fit.svlmorph<- aov(svl~X3morph, data = dat2017v2_nonas)  
summary(fit.svlmorph)
thsd_svlmorph<-TukeyHSD(fit.svlmorph)
pander(thsd_svlmorph$`X3morph`)

fit.hlmorph<- aov(head_length~X3morph, data = dat2017v2_nonas)  
summary(fit.hlmorph)
thsd_hlmorph<-TukeyHSD(fit.hlmorph)
pander(thsd_hlmorph$`X3morph`)

fit.hwmorph<- aov(head_width~X3morph, data = dat2017v2_nonas)  
summary(fit.hwmorph)
thsd_hwmorph<-TukeyHSD(fit.hwmorph)
pander(thsd_hwmorph$`X3morph`)

fit.hhmorph<- aov(head_height~X3morph, data = dat2017v2_nonas)  
summary(fit.hhmorph)
thsd_hhmorph<-TukeyHSD(fit.hhmorph)
pander(thsd_hhmorph$`X3morph`)

#ANOVA morph biteforce
fit.bitemorph<- aov(maxbite~X3morph, data = dat2017v2_nonas)  
summary(fit.bitemorph)
thsd_bitemorph<-TukeyHSD(fit.bitemorph)
pander(thsd_bitemorph$`X3morph`)

#ANOVA morph biteforce residuals
fit.bitemorphsize<- aov(bf_resid~X3morph, data = dat2017v2_nonas)  
summary(fit.bitemorphsize)
thsd_bitemorphsize<-TukeyHSD(fit.bitemorphsize)
pander(thsd_bitemorphsize$`X3morph`)

lmer <- lmer(maxbite ~ svl + (1|X3morph), data = dat2017v2_nonas)
summary(lmer)
svlmaxbite.lm <- lm (maxbite ~ svl, data = dat2017v2_nonas)
summary(svlmaxbite.lm)
scatter <- ggplot(dat2017v2_nonas, aes(x=svl, y=maxbite, color=X3morph)) + 
  geom_point(aes(size=8)) + 
  labs(title="Relationship between SVL and Maximum Biteforce",x="SVL (mm)", y = "Max Bite Force")
scatter
scattermorphcolor <- scatter + scale_color_manual(values=c("#E69F00", "#FFFFFF", "#F0E442"))
scattermorphcolor
scattermorphcolor + geom_smooth(method=lm)

box.mbite <- ggplot(dat2017v2_nonas, aes(x = X3morph, y = maxbite, fill = X3morph)) +
  geom_boxplot()
box.mbite
boxpoints.mbite <- box.mbite + geom_jitter()
boxpoints.mbite
mk <- boxpoints.mbite + scale_fill_manual(values=c("#E69F00", "#FFFFFF", "#F0E442"))
mk
mk + theme(panel.background = element_blank(), 
           axis.line = element_line(colour = "black"))

mk1 <- mk + theme(panel.background = element_blank(), 
                  axis.line = element_line(colour = "black"))

mk1 + ggtitle("Male Maximum Bite Force by Morph") +
  theme(plot.title = element_text(hjust = 0.5))

mk2 <- mk1 +ggtitle("Male Maximum Bite Force by Morph") +
  theme(plot.title = element_text(hjust = 0.5))

mk2

mk2 + labs(x ="Color Morph", y = "Maximum Bite Force in Newtons (N)")

mk3 <- mk2 + labs(x ="Color Morph", y = "Maximum Bite Force in Newtons (N)")

mk3 + theme(axis.text.x = element_text(size = 15),
            axis.text.y = element_text(size = 15),  
            axis.title.x = element_text(size = 15),
            axis.title.y = element_text(size = 15))

box.msvl <- ggplot(dat2017v2_nonas, aes(x = X3morph, y = svl, fill = X3morph)) +
  geom_boxplot()
box.msvl
boxpoints.msvl <- box.msvl + geom_jitter()
boxpoints.msvl
msvl <- boxpoints.msvl + scale_fill_manual(values=c("#E69F00", "#FFFFFF", "#F0E442"))
msvl
msvl1 <- msvl + theme(panel.background = element_blank(), 
                      axis.line = element_line(colour = "black"))
msvl1
msvl1 + ggtitle ("Male SVL by Morph") +
  theme(plot.title = element_text(hjust = 0.5))
msvl2 <- msvl1 + ggtitle("Male SVL by Morph") +
  theme(plot.title = element_text(hjust = 0.5))

msvl2

msvl2 + labs(x ="Color Morph", y = "SVL (mm)")

msvl3 <- msvl2 + labs(x ="Color Morph", y = "SVL (mm)")

msvl3 + theme(axis.text.x = element_text(size = 15),
              axis.text.y = element_text(size = 15),  
              axis.title.x = element_text(size = 15),
              axis.title.y = element_text(size = 15))

# FEMALE BITE FORCE ANALYSES - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
femdat2017 <- read.csv("females_2017.csv")

lmfemfinal <- lm(maxbite ~ svl + X3morph, data = femdat2017)
summary(lmfemfinal) 

#ANOVA biteforce
fit.fbitemorph<- aov(maxbite~X3morph, data = femdat2017)  
summary(fit.fbitemorph)
thsd_fbitemorph<-TukeyHSD(fit.fbitemorph)
pander(thsd_fbitemorph$`X3morph`)

#ANOVA SVL
fit.fsvlmorph<- aov(svl~X3morph, data = femdat2017)  
summary(fit.fsvlmorph)
thsd_fsvlmorph<-TukeyHSD(fit.fsvlmorph)
pander(thsd_fsvlmorph$`X3morph`)

#ANOVA head length
fit.fhlmorph<- aov(head_length~X3morph, data = femdat2017)  
summary(fit.fhlmorph)
thsd_fhlmorph<-TukeyHSD(fit.fhlmorph)
pander(thsd_fhlmorph$`X3morph`)

#ANOVA head width
fit.fhwmorph<- aov(head_width~X3morph, data = femdat2017)  
summary(fit.fhwmorph)
thsd_fhwmorph<-TukeyHSD(fit.fhwmorph)
pander(thsd_fhwmorph$`X3morph`)

#ANOVA head height
fit.fhhmorph<- aov(head_height~X3morph, data = femdat2017)  
summary(fit.fhhmorph)
thsd_fhhmorph<-TukeyHSD(fit.fhhmorph)
pander(thsd_fhhmorph$`X3morph`)

#Is female max bite force positively correlated with body size (SVL)?
lm.fembitesize <- lm(maxbite ~ svl, dat = femdat2017)
summary(lm.fembitesize)
cor.svlbite <- cor.test(femdat2017$maxbite, femdat2017$svl, 
                        method = "pearson")
cor.svlbite
ggscatter(femdat2017, x = "svl", y = "maxbite", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Body size (SVL)", ylab = "Female maximum bite force")

box.fembite <- ggplot(femdat2017, aes(x = X3morph, y = maxbite, fill = X3morph)) +
  geom_boxplot()
box.fembite
boxpoints.fembite <- box.fembite + geom_jitter()
boxpoints.fembite
ok <- boxpoints.fembite + scale_fill_manual(values=c("#E69F00", "#FFFFFF", "#F0E442"))
ok

ok + theme(panel.background = element_blank(), 
           axis.line = element_line(colour = "black"))

ok1 <- ok + theme(panel.background = element_blank(), 
                  axis.line = element_line(colour = "black"))

ok1 + ggtitle("Female Maximum Bite Force by Morph") +
  theme(plot.title = element_text(hjust = 0.5))

ok2 <- ok1 +ggtitle("Female Maximum Bite Force by Morph") +
  theme(plot.title = element_text(hjust = 0.5))

ok2

ok2 + labs(x ="Color Morph", y = "Maximum Bite Force in Newtons (N)")

ok3 <- ok2 + labs(x ="Color Morph", y = "Maximum Bite Force in Newtons (N)")

ok3 + theme(axis.text.x = element_text(size = 15),
            axis.text.y = element_text(size = 15),  
            axis.title.x = element_text(size = 15),
            axis.title.y = element_text(size = 15))


### You can see orange still tend to bite a little harder. Perhaps could be different if had more data?

#boxplot svl ~ morphs:
box.femsvl <- ggplot(femdat2017, aes(x = X3morph, y = svl, fill = X3morph)) +
  geom_boxplot()
box.femsvl
boxpoints.femsvl <- box.femsvl + geom_jitter()
boxpoints.femsvl
boxpoints.femsvl + scale_fill_manual(values=c("#E69F00", "#FFFFFF", "#F0E442"))
fembox <- boxpoints.femsvl + scale_fill_manual(values=c("#E69F00", "#FFFFFF", "#F0E442"))
fembox
fembox + theme(panel.background = element_blank(), 
               axis.line = element_line(colour = "black"))

fembox1 <- fembox + theme(panel.background = element_blank(), 
                          axis.line = element_line(colour = "black"))

fembox1 + ggtitle("Female SVL by Morph") +
  theme(plot.title = element_text(hjust = 0.5))

fembox2 <- fembox1 +ggtitle("Female SVL by Morph") +
  theme(plot.title = element_text(hjust = 0.5))

fembox2 + labs(x ="Color Morph", y = "SVL (mm)")


fembox3 <- fembox2 + labs(x ="Color Morph", y = "SVL (mm)")

fembox3 + theme(axis.text.x = element_text(size = 15),
                axis.text.y = element_text(size = 15),  
                axis.title.x = element_text(size = 15),
                axis.title.y = element_text(size = 15))


#linear model and morph scatterplit maxbite ~ svl:
lmer.fembitesize <- lmer(maxbite ~ svl + (1|X3morph), data = femdat2017)
summary(lmer.fembitesize)
femsvlmaxbite.lm <- lm (maxbite ~ svl, data = femdat2017)
summary(femsvlmaxbite.lm)
fem.bitesizescatter <- ggplot(femdat2017, aes(x=svl, y=maxbite, color=X3morph)) + 
  geom_point(aes(size=8)) + 
  labs(title="Relationship between SVL and Maximum Biteforce",x="SVL (mm)", y = "Max Bite Force")
fem.bitesizescatter
fem.scattermorphcolor <- fem.bitesizescatter + scale_color_manual(values=c("#E69F00", "#FFFFFF", "#F0E442"))
fem.scattermorphcolor
fem.scattermorphcolor + geom_smooth(method=lm)


# FINAL FEMALE SCATTER:
femscatter <- ggplot(femdat2017, aes(x=svl, y=maxbite, color=X3morph)) + 
  geom_point(aes(size=8)) + 
  labs(title="Relationship between SVL and Maximum Bite Force", subtitle = "Females", x="SVL (mm)", y = "Maximum Bite Force in Newtons (N)") + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))
femscatter
femscattermorphcolor <- femscatter + scale_color_manual(values=c("#E69F00", "#000000", "#F0E442"))
femscattermorphcolor
fems1 <- femscattermorphcolor + geom_smooth(method=lm)
fems1
fems1 + theme(panel.background = element_blank(), 
              axis.line = element_line(colour = "black"))
fems2 <- fems1 + theme(panel.background = element_blank(), 
                       axis.line = element_line(colour = "black"))
fems2
fems2 + theme(axis.text.x = element_text(size = 15),
              axis.text.y = element_text(size = 15),  
              axis.title.x = element_text(size = 15),
              axis.title.y = element_text(size = 15))

# III. CHEMICAL ANALYSES #################################################################################
chemdata<-read.table("chemical_dataset.txt",header=T, sep = "\t")

pca<-prcomp(chemdata[, 4:84])
pca$rotation[, 1:4]
summary(pca)
scores<-(pca$x)
scores<-as.data.frame(scores[, 1:4])


lm1<-lm(PC1~log10(svl))
lm2<-lm(PC2~log10(svl))
lm3<-lm(PC3~log10(svl))
lm4<-lm(PC4~log10(svl))

lm5<-lm(PC1~morph)
lm5b<-LSD.test(aov(PC1~morph),trt="morph",p.adj="bonferroni",group=F) 
lm6<-lm(PC2~morph)
lm7<-lm(PC3~morph)
lm8<-lm(PC4~morph)