library(ggplot2)
library(nlme)
library(car)
library(factoextra)
library(vegan)
library(ppcor)
library(basicTrendline)
library(psych)
library(dplyr)
library(FSA)
library(emmeans)
library(flextable)
library(multcompView)
library(reshape2)     
   
###### Repeatability of variables
## Using intraclass correlation coefficient

# Tpref
RepeatabilityTpref <- read.csv("RepeatabilityTpref.csv", row.names = 1)
View(RepeatabilityTpref)
ggplot(RepeatabilityTpref, aes(x=Tpref.1, y=Tpref.2)) + geom_point()
ICC(RepeatabilityTpref)
# ICC3 is the best here. Treats the "raters" as fixed effects rather than random, 
# while subjects are random effects
ylab1 <- "Tpref Trial 2 (°C)"
xlab1 <- "Tpref Trial 1 (°C)"
ggplot(RepeatabilityTpref, aes(Tpref.1, Tpref.2)) + geom_point()+theme_classic()+
  labs(x=xlab1, y=ylab1) +
  theme(axis.title = element_text(size=18))
cor.test(RepeatabilityTpref$Tpref.1, RepeatabilityTpref$Tpref.2)


# TprefCV
RepeatabilityTprefCV <- read.csv("~/Desktop/U. ornatus 2019/Boldness:POLS Excel Sheets/RepeatabilityTprefCV.csv", row.names=1)
View(RepeatabilityTprefCV)
ggplot(RepeatabilityTprefCV, aes(x=TprefCV.1, y=TprefCV.2)) + geom_point()
ICC(RepeatabilityTprefCV)


# UTset
RepeatabilityUTset <- read.csv("RepeatabilityUTset.csv", row.names = 1)
View(RepeatabilityUTset)
ggplot(RepeatabilityUTset, aes(x=UTset.1, y=UTset.2)) + geom_point()
ICC(RepeatabilityUTset)

# LTset
RepeatabilityLTset <- read.csv("RepeatabilityLTset.csv",  row.names = 1)
View(RepeatabilityLTset)
ggplot(RepeatabilityLTset, aes(x=LTset.1, y=LTset.2)) + geom_point()
ICC(RepeatabilityLTset)
# Repeatable; p value is significant but ICC value is still pretty low

# Tset
RepeatabilityTset <- read.csv("RepeatabilityTset.csv", row.names = 1)
View(RepeatabilityTset)
ggplot(RepeatabilityTset, aes(x=Tset.1, y=Tset.2)) + geom_point()
ICC(RepeatabilityTset)
# Repeatable

# Shuttle
RepeatabilityShuttle <- read.csv("RepeatabilityShuttle.csv", row.names=1)
View(RepeatabilityShuttle)
ggplot(RepeatabilityShuttle, aes(x=shuttle.1, y=shuttle.2)) + geom_point()
ICC(RepeatabilityShuttle)

# Transitions
RepeatabilityTransitions <- read.csv("RepeatabilityTransitions.csv", row.names = 1)
View(RepeatabilityTransitions)
ggplot(RepeatabilityTransitions, aes(x=transitions.1, y=transitions.2)) + geom_point()
ICC(RepeatabilityTransitions)
# Repeatable

# Basking
RepeatabilityBasking <- read.csv("RepeatabilityBasking.csv", row.names = 1)
View(RepeatabilityBasking)
ggplot(RepeatabilityBasking, aes(x=basking.1, y=basking.2)) + geom_point()
ICC(RepeatabilityBasking)
# Repeatable

# Perching
RepeatabilityPerching <- read.csv("RepeatabilityPerching.csv", row.names = 1)
View(RepeatabilityPerching)
ggplot(RepeatabilityPerching, aes(x=perching.1, y=perching.2)) + geom_point()
ICC(RepeatabilityPerching)
# Repeatable


##### Tset, LTset and all 3 measures of boldness are repeatable
## Tpref is fine... it's almost statistically significant, but visually it looks good.




## POLS traits...
POLS <- read.csv("~/Desktop/U. ornatus 2019/Boldness:POLS Excel Sheets/POLS_stamina.csv")
View(POLS)
str(POLS)
POLS$sex <- as.factor(POLS$sex)
POLS$morph <- as.factor(POLS$morph)
POLS$morph.group <- as.factor(POLS$morph.group)
str(POLS)
# 79 total lizards 

Male <- POLS[POLS$sex=="male",]
View(Male)
str(Male)
# 57 males

Female <- POLS[POLS$sex=="female",]
View(Female)
str(Female)
# 22 females

##### Summary stats for main variables by sex, t tests

# Tpref
POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(Tpref),sd=sd(Tpref),)
plot(x=POLS$sex, y=POLS$Tpref)
t.test(Tpref~sex, var.equal=FALSE, POLS)

# Tpref CV
POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(TprefCV),sd=sd(TprefCV))
plot(x=POLS$sex, y=POLS$TprefCV)
t.test(TprefCV~sex, var.equal=FALSE, POLS)

# Shuttle
POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(shuttle),sd=sd(shuttle))
t.test(shuttle~sex, var.equal=FALSE,data=POLS)

# stamina time
POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(staminaTime),sd=sd(staminaTime))
t.test(staminaTime~sex, var.equal=FALSE,data=POLS)

# stamina distance
POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(staminaDistance),sd=sd(staminaDistance))
t.test(staminaDistance~sex, var.equal=FALSE,data=POLS)

# stamina velocity
POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(staminaVelocity),sd=sd(staminaVelocity))
t.test(staminaVelocity~sex, var.equal=FALSE,data=POLS)

# transitions
POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(transition),sd=sd(transition))
t.test(transition~sex, var.equal=FALSE,data=POLS)

# basking
basktib<-POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(basking),sd=sd(basking))
flextable(basktib)
t.test(basking~sex, var.equal=FALSE,data=POLS)

# perching
perchtib<-POLS %>%
  group_by(sex) %>%
  summarise(mean=mean(perching),sd=sd(perching))
flextable(perchtib)
t.test(perching~sex, var.equal=FALSE,data=POLS)



########## Do SVL and mass correlate with any PCA variables?
#SVL males

# Tpref
plot(Male$Tpref, Male$SVL)
m.svl.Tpref.lm <- lm(Tpref~SVL, data=Male)
m.svl.Tpref.res <- resid(m.svl.Tpref.lm)
plot(Male$SVL, m.svl.Tpref.res) + abline(0,0)
cor.test(Male$Tpref, Male$SVL, method="pearson")

# CV
plot(Male$TprefCV, Male$SVL)
m.svl.cv.lm <- lm(TprefCV~SVL, data=Male)
m.svl.cv.res <- resid(m.svl.cv.lm)
plot(Male$SVL, m.svl.cv.res) + abline(0,0)
cor.test(Male$TprefCV, Male$SVL, method="pearson")

# shuttle
plot(Male$shuttle, Male$SVL)
m.svl.shut.lm <- lm(shuttle~SVL, data=Male)
m.svl.shut.res <- resid(m.svl.shut.lm)
plot(Male$SVL, m.svl.shut.res) + abline(0,0)
cor.test(Male$shuttle, Male$SVL, method="pearson")

# transition
plot(Male$transition, Male$SVL)
m.svl.trans.lm <- lm(transition~SVL, data=Male)
m.svl.trans.res <- resid(m.svl.trans.lm)
plot(Male$SVL, m.svl.trans.res) + abline(0,0)
cor.test(Male$transition, Male$SVL, method="pearson")

# basking
plot(Male$basking, Male$SVL)
m.svl.bask.lm <- lm(basking~SVL, data=Male)
m.svl.bask.res <- resid(m.svl.bask.lm)
plot(Male$SVL, m.svl.bask.res) + abline(0,0)
cor.test(Male$basking, Male$SVL, method="pearson")

# perching
plot(Male$perching, Male$SVL)
m.svl.perch.lm <- lm(perching~SVL, data=Male)
m.svl.perch.res <- resid(m.svl.perch.lm)
plot(Male$SVL, m.svl.perch.res) + abline(0,0)
cor.test(Male$perching, Male$SVL, method="pearson")

# stamina time
plot(Male$staminaTime, Male$SVL)
m.svl.time.lm <- lm(staminaTime~SVL, data=Male)
m.svl.time.res <- resid(m.svl.time.lm)
plot(Male$SVL, m.svl.time.res) + abline(0,0)
cor.test(Male$staminaTime, Male$SVL, method="pearson")

# stamina distance
plot(Male$staminaDistance, Male$SVL)
m.svl.dist.lm <- lm(staminaDistance~SVL, data=Male)
m.svl.dist.res <- resid(m.svl.dist.lm)
plot(Male$SVL, m.svl.dist.res) + abline(0,0)
cor.test(Male$staminaDistance, Male$SVL, method="pearson")

# stamina velocity
plot(Male$staminaVelocity, Male$SVL)
m.svl.vel.lm <- lm(staminaVelocity~SVL, data=Male)
m.svl.vel.res <- resid(m.svl.vel.lm)
plot(Male$SVL, m.svl.vel.res) + abline(0,0)
cor.test(Male$staminaVelocity, Male$SVL, method="pearson")

# SVL doesn't relate to anything in males


## mass in males
# Tpref
plot(Male$Tpref, Male$Mass)
m.mass.Tpref.lm <- lm(Tpref~Mass, data=Male)
m.mass.Tpref.res <- resid(m.mass.Tpref.lm)
plot(Male$Mass, m.mass.Tpref.res) + abline(0,0)
cor.test(Male$Tpref, Male$Mass, method="pearson")

# CV
plot(Male$TprefCV, Male$Mass)
m.mass.cv.lm <- lm(TprefCV~Mass, data=Male)
m.mass.cv.res <- resid(m.mass.cv.lm)
plot(Male$Mass, m.mass.cv.res) + abline(0,0)
cor.test(Male$TprefCV, Male$Mass, method="spearman")

# shuttle
plot(Male$shuttle, Male$Mass)
m.mass.shut.lm <- lm(shuttle~Mass, data=Male)
m.mass.shut.res <- resid(m.mass.shut.lm)
plot(Male$Mass, m.mass.shut.res) + abline(0,0)
cor.test(Male$shuttle, Male$Mass, method="pearson")

# transition
plot(Male$transition, Male$Mass)
m.mass.trans.lm <- lm(transition~Mass, data=Male)
m.mass.trans.res <- resid(m.mass.trans.lm)
plot(Male$Mass, m.mass.trans.res) + abline(0,0)
cor.test(Male$transition, Male$Mass, method="spearman")

# basking
plot(Male$basking, Male$Mass)
m.mass.bask.lm <- lm(basking~Mass, data=Male)
m.mass.bask.res <- resid(m.mass.bask.lm)
plot(Male$Mass, m.mass.bask.res) + abline(0,0)
cor.test(Male$basking, Male$Mass, method="spearman")

# perching
plot(Male$perching, Male$Mass)
m.mass.perch.lm <- lm(perching~Mass, data=Male)
m.mass.perch.res <- resid(m.mass.perch.lm)
plot(Male$Mass, m.mass.perch.res) + abline(0,0)
cor.test(Male$perching, Male$Mass, method="pearson")

# stamina time
plot(Male$staminaTime, Male$Mass)
m.mass.time.lm <- lm(staminaTime~Mass, data=Male)
m.mass.time.res <- resid(m.mass.time.lm)
plot(Male$Mass, m.mass.time.res) + abline(0,0)
cor.test(Male$staminaTime, Male$Mass, method="pearson")

# stamina distance
plot(Male$staminaDistance, Male$Mass)
m.mass.dist.lm <- lm(staminaDistance~Mass, data=Male)
m.mass.dist.res <- resid(m.mass.dist.lm)
plot(Male$Mass, m.mass.dist.res) + abline(0,0)
cor.test(Male$staminaDistance, Male$Mass, method="pearson")

# stamina velocity
plot(Male$staminaVelocity, Male$Mass)
m.mass.vel.lm <- lm(staminaVelocity~Mass, data=Male)
m.mass.vel.res <- resid(m.mass.vel.lm)
plot(Male$Mass, m.mass.vel.res) + abline(0,0)
cor.test(Male$staminaVelocity, Male$Mass, method="pearson")

### SVL females

# Tpref
plot(Female$Tpref, Female$SVL)
f.svl.Tpref.lm <- lm(Tpref~SVL, data=Female)
f.svl.Tpref.res <- resid(f.svl.Tpref.lm)
plot(Female$SVL, f.svl.Tpref.res) + abline(0,0)
cor.test(Female$Tpref, Female$SVL, method="pearson")
# significant

# CV
plot(Female$TprefCV, Female$SVL)
f.svl.cv.lm <- lm(TprefCV~SVL, data=Female)
f.svl.cv.res <- resid(f.svl.cv.lm)
plot(Female$SVL, f.svl.cv.res) + abline(0,0)
cor.test(Female$TprefCV, Female$SVL, method="pearson")

# shuttle
plot(Female$shuttle, Female$SVL)
f.svl.shut.lm <- lm(shuttle~SVL, data=Female)
f.svl.shut.res <- resid(f.svl.shut.lm)
plot(Female$SVL, f.svl.shut.res) + abline(0,0)
cor.test(Female$shuttle, Female$SVL, method="pearson")

# transition
plot(Female$transition, Female$SVL)
f.svl.trans.lm <- lm(transition~SVL, data=Female)
f.svl.trans.res <- resid(f.svl.trans.lm)
plot(Female$SVL, f.svl.trans.res) + abline(0,0)
cor.test(Female$transition, Female$SVL, method="pearson")

# basking
plot(Female$basking, Female$SVL)
f.svl.bask.lm <- lm(basking~SVL, data=Female)
f.svl.bask.res <- resid(f.svl.bask.lm)
plot(Female$SVL, f.svl.bask.res) + abline(0,0)
cor.test(Female$basking, Female$SVL, method="pearson")

# perching
plot(Female$perching, Female$SVL)
f.svl.perch.lm <- lm(perching~SVL, data=Female)
f.svl.perch.res <- resid(f.svl.perch.lm)
plot(Female$SVL, f.svl.perch.res) + abline(0,0)
cor.test(Female$perching, Female$SVL, method="pearson")

# stamina time
plot(Female$staminaTime, Female$SVL)
f.svl.time.lm <- lm(staminaTime~SVL, data=Female)
f.svl.time.res <- resid(f.svl.time.lm)
plot(Female$SVL, f.svl.time.res) + abline(0,0)
cor.test(Female$staminaTime, Female$SVL, method="pearson")

# stamina distance
plot(Female$staminaDistance, Female$SVL)
f.svl.dist.lm <- lm(staminaDistance~SVL, data=Female)
f.svl.dist.res <- resid(f.svl.dist.lm)
plot(Female$SVL, f.svl.dist.res) + abline(0,0)
cor.test(Female$staminaDistance, Female$SVL, method="pearson")

# stamina velocity
plot(Female$staminaVelocity, Female$SVL)
f.svl.vel.lm <- lm(staminaVelocity~SVL, data=Female)
f.svl.vel.res <- resid(f.svl.vel.lm)
plot(Female$SVL, f.svl.vel.res) + abline(0,0)
cor.test(Female$staminaVelocity, Female$SVL, method="pearson")

# SVL correlates with Tpref in females

## mass in females
# Tpref
plot(Female$Tpref, Female$Mass)
f.mass.Tpref.lm <- lm(Tpref~Mass, data=Female)
f.mass.Tpref.res <- resid(f.mass.Tpref.lm)
plot(Female$Mass, f.mass.Tpref.res) + abline(0,0)
cor.test(Female$Tpref, Female$Mass, method="pearson")
# significant

# CV
plot(Female$TprefCV, Female$Mass)
f.mass.cv.lm <- lm(TprefCV~Mass, data=Female)
f.mass.cv.res <- resid(f.mass.cv.lm)
plot(Female$Mass, f.mass.cv.res) + abline(0,0)
cor.test(Female$TprefCV, Female$Mass, method="pearson")

# shuttle
plot(Female$shuttle, Female$Mass)
f.mass.shut.lm <- lm(shuttle~Mass, data=Female)
f.mass.shut.res <- resid(f.mass.shut.lm)
plot(Female$Mass, f.mass.shut.res) + abline(0,0)
cor.test(Female$shuttle, Female$Mass, method="pearson")

# transition
plot(Female$transition, Female$Mass)
f.mass.trans.lm <- lm(transition~Mass, data=Female)
f.mass.trans.res <- resid(f.mass.trans.lm)
plot(Female$Mass, f.mass.trans.res) + abline(0,0)
cor.test(Female$transition, Female$Mass, method="pearson")

# basking
plot(Female$basking, Female$Mass)
f.mass.bask.lm <- lm(basking~Mass, data=Female)
f.mass.bask.res <- resid(f.mass.bask.lm)
plot(Female$Mass, f.mass.bask.res) + abline(0,0)
cor.test(Female$basking, Female$Mass, method="spearman")

# perching
plot(Female$perching, Female$Mass)
f.mass.perch.lm <- lm(perching~Mass, data=Female)
f.mass.perch.res <- resid(f.mass.perch.lm)
plot(Female$Mass, f.mass.perch.res) + abline(0,0)
cor.test(Female$perching, Female$Mass, method="spearman")

# stamina time
plot(Female$staminaTime, Female$Mass)
f.mass.time.lm <- lm(staminaTime~Mass, data=Female)
f.mass.time.res <- resid(f.mass.time.lm)
plot(Female$Mass, f.mass.time.res) + abline(0,0)
cor.test(Female$staminaTime, Female$Mass, method="pearson")

# stamina distance
plot(Female$staminaDistance, Female$Mass)
f.mass.dist.lm <- lm(staminaDistance~Mass, data=Female)
f.mass.dist.res <- resid(f.mass.dist.lm)
plot(Female$Mass, f.mass.dist.res) + abline(0,0)
cor.test(Female$staminaDistance, Female$Mass, method="pearson")

# stamina velocity
plot(Female$staminaVelocity, Female$Mass)
f.mass.vel.lm <- lm(staminaVelocity~Mass, data=Female)
f.mass.vel.res <- resid(f.mass.vel.lm)
plot(Female$Mass, f.mass.vel.res) + abline(0,0)
cor.test(Female$staminaVelocity, Female$Mass, method="pearson")

# mass correlates to Tpref in females


####### Morph differences in traits?

# Males first, all morphs

# 12 B
# 6 OB
# 17 YB
# 2 O
# 16 Y
# 4 YO

# Tpref
ggplot(Male, aes(x=morph, y=Tpref))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphpref <- lm(Tpref~morph, Male)
m.morphpref.res <- resid(m.lm.morphpref)
shapiro.test(m.morphpref.res)
leveneTest(Tpref~morph, Male)
kruskal.test(Tpref~morph, Male)

# Tpref CV
ggplot(Male, aes(x=morph, y=TprefCV))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphcv <- lm(TprefCV~morph, Male)
m.morphcv.res <- resid(m.lm.morphcv)
shapiro.test(m.morphcv.res)
leveneTest(TprefCV~morph, Male)
kruskal.test(TprefCV~morph, Male)

# shuttle
ggplot(Male, aes(x=morph, y=shuttle))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphshut <- lm(shuttle~morph, Male)
m.morphshut.res <- resid(m.lm.morphshut)
shapiro.test(m.morphshut.res)
leveneTest(shuttle~morph, Male)
Anova(m.lm.morphshut)

# transitions
ggplot(Male, aes(x=morph, y=transition))+geom_boxplot()+theme_classic()+
  stat_summary(fun.y=mean, geom="point") + xlab("Morph") + ylab("Exploration (number of transitions)") + scale_x_discrete(labels=c('B','O/B','Y/B','O','Y',"O/Y"))+theme(axis.title = element_text(size=18))
m.lm.morphtrans <- lm(transition~morph, Male)
m.morphtrans.res <- resid(m.lm.morphtrans)
shapiro.test(m.morphtrans.res)
leveneTest(transition~morph, Male)
kruskal.test(transition~morph, Male)
dunnTest(transition~morph, Male)
# I read that if I have a priori comparisions, I don't adjust p value. Ask Miles
# blue and yellow/blue males are more exploratory than yellow


# basking
ggplot(Male, aes(x=morph, y=basking))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphbask <- lm(basking~morph, Male)
m.morphbask.res <- resid(m.lm.morphbask)
shapiro.test(m.morphbask.res)
leveneTest(basking~morph, Male)
kruskal.test(basking~morph, Male)

# perching
ggplot(Male, aes(x=morph, y=perching))+geom_boxplot()+theme_classic()+
  stat_summary(fun.y=mean, geom="point") + xlab("Morph") + ylab("Perching (seconds)") + scale_x_discrete(labels=c('B','O/B','Y/B','O','Y',"O/Y"))+theme(axis.title = element_text(size=18))
m.lm.morphperch <- lm(perching~morph, Male)
m.morphperch.res <- resid(m.lm.morphperch)
shapiro.test(m.morphperch.res)
leveneTest(perching~morph, Male)
kruskal.test(perching~morph, Male)
dunnTest(perching~morph, Male)
# yellow/blue males bolder than yellow males and yellow/orange males

# stamina time
ggplot(Male, aes(x=morph, y=staminaTime))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphtime <- lm(staminaTime~morph, Male)
m.morphtime.res <- resid(m.lm.morphtime)
shapiro.test(m.morphtime.res)
leveneTest(staminaTime~morph, Male)
Anova(m.lm.morphtime)

# stamina distance
ggplot(Male, aes(x=morph, y=staminaDistance))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphdist <- lm(staminaDistance~morph, Male)
m.morphdist.res <- resid(m.lm.morphdist)
shapiro.test(m.morphdist.res)
leveneTest(staminaDistance~morph, Male)
Anova(m.lm.morphdist)

# stamina velocity
ggplot(Male, aes(x=morph, y=staminaVelocity))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphvel <- lm(staminaVelocity~morph, Male)
m.morphvel.res <- resid(m.lm.morphvel)
shapiro.test(m.morphvel.res)
leveneTest(staminaVelocity~morph, Male)
Anova(m.lm.morphvel)

## In males, get morph differences in transitions and perching


##### What about morph groups?

# Tpref
ggplot(Male, aes(x=morph.group, y=Tpref))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGpref <- lm(Tpref~morph.group, Male)
m.morphGpref.res <- resid(m.lm.morphGpref)
shapiro.test(m.morphGpref.res)
leveneTest(Tpref~morph.group, Male)
kruskal.test(Tpref~morph.group, Male)

# Tpref CV
ggplot(Male, aes(x=morph.group, y=TprefCV))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGcv <- lm(TprefCV~morph.group, Male)
m.morphGcv.res <- resid(m.lm.morphGcv)
shapiro.test(m.morphGcv.res)
leveneTest(TprefCV~morph.group, Male)
kruskal.test(TprefCV~morph.group, Male)

# shuttle
ggplot(Male, aes(x=morph.group, y=shuttle))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGshut <- lm(shuttle~morph.group, Male)
m.morphGshut.res <- resid(m.lm.morphGshut)
shapiro.test(m.morphGshut.res)
leveneTest(shuttle~morph.group, Male)
Anova(m.lm.morphGshut)

# transitions
ggplot(Male, aes(x=morph.group, y=transition))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGtrans <- lm(transition~morph.group, Male)
m.morphGtrans.res <- resid(m.lm.morphGtrans)
shapiro.test(m.morphGtrans.res)
leveneTest(transition~morph.group, Male)
kruskal.test(transition~morph.group, Male)
dunnTest(transition~morph.group, Male)
# blue more exploratory than Y,O,OY and OB,OY more exploratory than Y,O,OY

# basking
ggplot(Male, aes(x=morph.group, y=basking))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGbask <- lm(basking~morph.group, Male)
m.morphGbask.res <- resid(m.lm.morphGbask)
shapiro.test(m.morphGbask.res)
leveneTest(basking~morph.group, Male)
kruskal.test(basking~morph.group, Male)

# perching
ggplot(Male, aes(x=morph.group, y=perching))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGperch <- lm(perching~morph.group, Male)
m.morphGperch.res <- resid(m.lm.morphGperch)
shapiro.test(m.morphGperch.res)
leveneTest(perching~morph.group, Male)
kruskal.test(perching~morph.group, Male)
dunnTest(perching~morph.group, Male)
# OB,YB bolder than O,Y,OY

# stamina time
ggplot(Male, aes(x=morph.group, y=staminaTime))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGtime <- lm(staminaTime~morph.group, Male)
m.morphGtime.res <- resid(m.lm.morphGtime)
shapiro.test(m.morphGtime.res)
leveneTest(staminaTime~morph.group, Male)
Anova(m.lm.morphGtime)

# stamina distance
ggplot(Male, aes(x=morph.group, y=staminaDistance))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGdist <- lm(staminaDistance~morph.group, Male)
m.morphGdist.res <- resid(m.lm.morphGdist)
shapiro.test(m.morphGdist.res)
leveneTest(staminaDistance~morph.group, Male)
Anova(m.lm.morphGdist)

# stamina velocity
ggplot(Male, aes(x=morph.group, y=staminaVelocity))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
m.lm.morphGvel <- lm(staminaVelocity~morph.group, Male)
m.morphGvel.res <- resid(m.lm.morphGvel)
shapiro.test(m.morphGvel.res)
leveneTest(staminaVelocity~morph.group, Male)
Anova(m.lm.morphGvel)

# Same as in all morphs, obviously a little different in terms of exact pairwise differences


##### Females

# 10 O
# 10 Y
# 2 W

# Tpref
ggplot(Female, aes(x=morph, y=Tpref))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphpref <- lm(Tpref~morph, Female)
f.morphpref.res <- resid(f.lm.morphpref)
shapiro.test(f.morphpref.res)
leveneTest(Tpref~morph, Female)
Anova(f.lm.morphpref)

# Tpref CV
ggplot(Female, aes(x=morph, y=TprefCV))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphcv <- lm(TprefCV~morph, Female)
f.morphcv.res <- resid(f.lm.morphcv)
shapiro.test(f.morphcv.res)
leveneTest(TprefCV~morph, Female)
Anova(f.lm.morphcv)

# shuttle
ggplot(Female, aes(x=morph, y=shuttle))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphshut <- lm(shuttle~morph, Female)
f.morphshut.res <- resid(f.lm.morphshut)
shapiro.test(f.morphshut.res)
leveneTest(shuttle~morph, Female)
Anova(f.lm.morphshut)

# transitions
ggplot(Female, aes(x=morph, y=transition))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphtrans <- lm(transition~morph, Female)
f.morphtrans.res <- resid(f.lm.morphtrans)
shapiro.test(f.morphtrans.res)
leveneTest(transition~morph, Female)
kruskal.test(transition~morph, Female)

# basking
ggplot(Female, aes(x=morph, y=basking))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphbask <- lm(basking~morph, Female)
f.morphbask.res <- resid(f.lm.morphbask)
shapiro.test(f.morphbask.res)
leveneTest(basking~morph, Female)
kruskal.test(basking~morph, Female)

# perching
ggplot(Female, aes(x=morph, y=perching))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphperch <- lm(perching~morph, Female)
f.morphperch.res <- resid(f.lm.morphperch)
shapiro.test(f.morphperch.res)
leveneTest(perching~morph, Female)
kruskal.test(perching~morph, Female)

# stamina time
ggplot(Female, aes(x=morph, y=staminaTime))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphtime <- lm(staminaTime~morph, Female)
f.morphtime.res <- resid(f.lm.morphtime)
shapiro.test(f.morphtime.res)
leveneTest(staminaTime~morph, Female)
Anova(f.lm.morphtime)

# stamina distance
ggplot(Female, aes(x=morph, y=staminaDistance))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphdist <- lm(staminaDistance~morph, Female)
f.morphdist.res <- resid(f.lm.morphdist)
shapiro.test(f.morphdist.res)
leveneTest(staminaDistance~morph, Female)
Anova(f.lm.morphdist)

# stamina velocity
ggplot(Female, aes(x=morph, y=staminaVelocity))+geom_boxplot()+theme_bw()+
  stat_summary(fun.y=mean, geom="point")
f.lm.morphvel <- lm(staminaVelocity~morph, Female)
f.morphvel.res <- resid(f.lm.morphvel)
shapiro.test(f.morphvel.res)
leveneTest(staminaVelocity~morph, Female)
Anova(f.lm.morphvel)




################### Now, looking for significant correlations among variables

#### For each relationship: overall, males, females

####### Tpref

## CV
plot(POLS$Tpref, POLS$TprefCV)
tpref.cv.lm <- lm(Tpref~TprefCV, POLS)
tpref.cv.res <- resid(tpref.cv.lm)
plot(POLS$TprefCV, tpref.cv.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$TprefCV, method="pearson")
# negative correlation

plot(Male$Tpref, Male$TprefCV)
m.tpref.cv.lm <- lm(Tpref~TprefCV, Male)
m.tpref.cv.res <- resid(m.tpref.cv.lm)
plot(Male$TprefCV, m.tpref.cv.res) + abline(0,0)
cor.test(Male$Tpref, Male$TprefCV, method="pearson")
# negative in males

plot(Female$Tpref, Female$TprefCV)
f.tpref.cv.lm <- lm(Tpref~TprefCV, Female)
f.tpref.cv.res <- resid(f.tpref.cv.lm)
plot(Female$TprefCV, f.tpref.cv.res) + abline(0,0)
cor.test(Female$Tpref, Female$TprefCV, method="pearson")
# nothing in females


## shuttle
plot(POLS$Tpref, POLS$shuttle)
tpref.shut.lm <- lm(Tpref~shuttle, POLS)
tpref.shut.res <- resid(tpref.shut.lm)
plot(POLS$shuttle, tpref.shut.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$shuttle, method="pearson")
# negative correlation

plot(Male$Tpref, Male$shuttle)
m.tpref.shut.lm <- lm(Tpref~shuttle, Male)
m.tpref.shut.res <- resid(m.tpref.shut.lm)
plot(Male$shuttle, m.tpref.shut.res) + abline(0,0)
cor.test(Male$Tpref, Male$shuttle, method="pearson")
# negative in males

plot(Female$Tpref, Female$shuttle)
f.tpref.shut.lm <- lm(Tpref~shuttle, Female)
f.tpref.shut.res <- resid(f.tpref.shut.lm)
plot(Female$shuttle, f.tpref.shut.res) + abline(0,0)
cor.test(Female$Tpref, Female$shuttle, method="pearson")
# nothing in females


## transitions
plot(POLS$Tpref, POLS$transition)
tpref.trans.lm <- lm(Tpref~transition, POLS)
tpref.trans.res <- resid(tpref.trans.lm)
plot(POLS$transition, tpref.trans.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$transition, method="spearman")

plot(Male$Tpref, Male$transition)
m.tpref.trans.lm <- lm(Tpref~transition, Male)
m.tpref.trans.res <- resid(m.tpref.trans.lm)
plot(Male$transition, m.tpref.trans.res) + abline(0,0)
cor.test(Male$Tpref, Male$transition, method="spearman")

plot(Female$Tpref, Female$transition)
f.tpref.trans.lm <- lm(Tpref~transition, Female)
f.tpref.trans.res <- resid(f.tpref.trans.lm)
plot(Female$transition, f.tpref.trans.res) + abline(0,0)
cor.test(Female$Tpref, Female$transition, method="pearson")
# nothing 


## basking
plot(POLS$Tpref, POLS$basking)
tpref.bask.lm <- lm(Tpref~basking, POLS)
tpref.bask.res <- resid(tpref.bask.lm)
plot(POLS$basking, tpref.bask.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$basking, method="spearman")
xlab2 <- "Tpref (°C)"
ggplot(POLS, aes(Tpref, basking)) + geom_point()+theme_classic()+
  labs(x=xlab2, y="Basking (s)") +
  theme(axis.title = element_text(size=18))+stat_smooth(method="lm")
# Negative, doesn't look great

plot(Male$Tpref, Male$basking)
m.tpref.bask.lm <- lm(Tpref~basking, Male)
m.tpref.bask.res <- resid(m.tpref.bask.lm)
plot(Male$basking, m.tpref.bask.res) + abline(0,0)
cor.test(Male$Tpref, Male$basking, method="spearman")
# negative

plot(Female$Tpref, Female$basking)
f.tpref.bask.lm <- lm(Tpref~basking, Female)
f.tpref.bask.res <- resid(f.tpref.bask.lm)
plot(Female$basking, f.tpref.bask.res) + abline(0,0)
cor.test(Female$Tpref, Female$basking, method="pearson")
# nothing


## perching
plot(POLS$Tpref, POLS$perching)
tpref.perch.lm <- lm(Tpref~perching, POLS)
tpref.perch.res <- resid(tpref.perch.lm)
plot(POLS$perching, tpref.perch.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$perching, method="spearman")
ggplot(POLS, aes(Tpref, perching)) + geom_point()+theme_classic()+
  labs(x=xlab2, y="Perching (s)") +
  theme(axis.title = element_text(size=18))+stat_smooth(method="lm")
# Same as basking

plot(Male$Tpref, Male$perching)
m.tpref.perch.lm <- lm(Tpref~perching, Male)
m.tpref.perch.res <- resid(m.tpref.perch.lm)
plot(Male$perching, m.tpref.perch.res) + abline(0,0)
cor.test(Male$Tpref, Male$perching, method="spearman")
# negative

plot(Female$Tpref, Female$perching)
f.tpref.perch.lm <- lm(Tpref~perching, Female)
f.tpref.perch.res <- resid(f.tpref.perch.lm)
plot(Female$perching, f.tpref.perch.res) + abline(0,0)
cor.test(Female$Tpref, Female$perching, method="spearman")
# nothing


## Stamina Distance
plot(POLS$Tpref, POLS$staminaDistance)
tpref.dist.lm <- lm(Tpref~staminaDistance, POLS)
tpref.dist.res <- resid(tpref.dist.lm)
plot(POLS$staminaDistance, tpref.dist.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$staminaDistance, method="pearson")

plot(Male$Tpref, Male$staminaDistance)
m.tpref.dist.lm <- lm(Tpref~staminaDistance, Male)
m.tpref.dist.res <- resid(m.tpref.dist.lm)
plot(Male$staminaDistance, m.tpref.dist.res) + abline(0,0)
cor.test(Male$Tpref, Male$staminaDistance, method="pearson")

plot(Female$Tpref, Female$staminaDistance)
f.tpref.dist.lm <- lm(Tpref~staminaDistance, Female)
f.tpref.dist.res <- resid(f.tpref.dist.lm)
plot(Female$staminaDistance, f.tpref.dist.res) + abline(0,0)
cor.test(Female$Tpref, Female$staminaDistance, method="pearson")
# nothing 


## Stamina Time
plot(POLS$Tpref, POLS$staminaTime)
tpref.time.lm <- lm(Tpref~staminaTime, POLS)
tpref.time.res <- resid(tpref.time.lm)
plot(POLS$staminaTime, tpref.time.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$staminaTime, method="pearson")

plot(Male$Tpref, Male$staminaTime)
m.tpref.time.lm <- lm(Tpref~staminaTime, Male)
m.tpref.time.res <- resid(m.tpref.time.lm)
plot(Male$staminaTime, m.tpref.time.res) + abline(0,0)
cor.test(Male$Tpref, Male$staminaTime, method="pearson")

plot(Female$Tpref, Female$staminaTime)
f.tpref.time.lm <- lm(Tpref~staminaTime, Female)
f.tpref.time.res <- resid(f.tpref.time.lm)
plot(Female$staminaTime, f.tpref.time.res) + abline(0,0)
cor.test(Female$Tpref, Female$staminaTime, method="pearson")
# nothing


## Stamina Velocity
plot(POLS$Tpref, POLS$staminaVelocity)
tpref.vel.lm <- lm(Tpref~staminaVelocity, POLS)
tpref.vel.res <- resid(tpref.vel.lm)
plot(POLS$staminaVelocity, tpref.vel.res) + abline(0,0)
cor.test(POLS$Tpref, POLS$staminaVelocity, method="pearson")

plot(Male$Tpref, Male$staminaVelocity)
m.tpref.vel.lm <- lm(Tpref~staminaVelocity, Male)
m.tpref.vel.res <- resid(m.tpref.vel.lm)
plot(Male$staminaVelocity, m.tpref.vel.res) + abline(0,0)
cor.test(Male$Tpref, Male$staminaVelocity, method="pearson")

plot(Female$Tpref, Female$staminaVelocity)
f.tpref.vel.lm <- lm(Tpref~staminaVelocity, Female)
f.tpref.vel.res <- resid(f.tpref.vel.lm)
plot(Female$staminaVelocity, f.tpref.vel.res) + abline(0,0)
cor.test(Female$Tpref, Female$staminaVelocity, method="pearson")
# nothing

pcor.test(POLS$Tpref, POLS$staminaTime, POLS$staminaTb)
pcor.test(POLS$Tpref, POLS$staminaDistance, POLS$staminaTb)
pcor.test(POLS$Tpref, POLS$staminaVelocity, POLS$staminaTb)
# Nothing

### Tpref correlated with CV and shuttling... lizards that prefer warmer 
# temperatures shuttle less and utilize a lower range of temperatures
# Additionally, Tpref and basking/perching time in males is correlated (negative)

######## CV

## shuttle
plot(POLS$TprefCV, POLS$shuttle)
cv.shut.lm <- lm(TprefCV~shuttle, POLS)
cv.shut.res <- resid(cv.shut.lm)
plot(POLS$shuttle, cv.shut.res) + abline(0,0)
cor.test(POLS$TprefCV, POLS$shuttle, method="pearson")
# correlation

plot(Male$TprefCV, Male$shuttle)
m.cv.shut.lm <- lm(TprefCV~shuttle, Male)
m.cv.shut.res <- resid(m.cv.shut.lm)
plot(Male$shuttle, m.cv.shut.res) + abline(0,0)
cor.test(Male$TprefCV, Male$shuttle, method="pearson")
# positive

plot(Female$TprefCV, Female$shuttle)
f.cv.shut.lm <- lm(TprefCV~shuttle, Female)
f.cv.shut.res <- resid(f.cv.shut.lm)
plot(Female$shuttle, f.cv.shut.res) + abline(0,0)
cor.test(Female$TprefCV, Female$shuttle, method="pearson")
# postive


## transition
plot(POLS$TprefCV, POLS$transition)
cv.trans.lm <- lm(TprefCV~transition, POLS)
cv.trans.res <- resid(cv.trans.lm)
plot(POLS$transition, cv.trans.res) + abline(0,0)
cor.test(POLS$TprefCV, POLS$transition, method="pearson")

plot(Male$TprefCV, Male$transition)
m.cv.trans.lm <- lm(TprefCV~transition, Male)
m.cv.trans.res <- resid(m.cv.trans.lm)
plot(Male$transition, m.cv.trans.res) + abline(0,0)
cor.test(Male$TprefCV, Male$transition, method="pearson")

plot(Female$TprefCV, Female$transition)
f.cv.trans.lm <- lm(TprefCV~transition, Female)
f.cv.trans.res <- resid(f.cv.trans.lm)
plot(Female$transition, f.cv.trans.res) + abline(0,0)
cor.test(Female$TprefCV, Female$transition, method="pearson")
# nothing


## basking
plot(POLS$TprefCV, POLS$basking)
cv.bask.lm <- lm(TprefCV~basking, POLS)
cv.bask.res <- resid(cv.bask.lm)
plot(POLS$basking, cv.bask.res) + abline(0,0)
cor.test(POLS$TprefCV, POLS$basking, method="spearman")

plot(Male$TprefCV, Male$basking)
m.cv.bask.lm <- lm(TprefCV~basking, Male)
m.cv.bask.res <- resid(m.cv.bask.lm)
plot(Male$basking, m.cv.bask.res) + abline(0,0)
cor.test(Male$TprefCV, Male$basking, method="pearson")

plot(Female$TprefCV, Female$basking)
f.cv.bask.lm <- lm(TprefCV~basking, Female)
f.cv.bask.res <- resid(f.cv.bask.lm)
plot(Female$basking, f.cv.bask.res) + abline(0,0)
cor.test(Female$TprefCV, Female$basking, method="pearson")
# nothing


## perching
plot(POLS$TprefCV, POLS$perching)
cv.perch.lm <- lm(TprefCV~perching, POLS)
cv.perch.res <- resid(cv.perch.lm)
plot(POLS$perching, cv.perch.res) + abline(0,0)
cor.test(POLS$TprefCV, POLS$perching, method="pearson")

plot(Male$TprefCV, Male$perching)
m.cv.perch.lm <- lm(TprefCV~perching, Male)
m.cv.perch.res <- resid(m.cv.perch.lm)
plot(Male$perching, m.cv.perch.res) + abline(0,0)
cor.test(Male$TprefCV, Male$perching, method="pearson")

plot(Female$TprefCV, Female$perching)
f.cv.perch.lm <- lm(TprefCV~perching, Female)
f.cv.perch.res <- resid(f.cv.perch.lm)
plot(Female$perching, f.cv.perch.res) + abline(0,0)
cor.test(Female$TprefCV, Female$perching, method="pearson")
# nothing


## stamina distance
plot(POLS$TprefCV, POLS$staminaDistance)
cv.dist.lm <- lm(TprefCV~staminaDistance, POLS)
cv.dist.res <- resid(cv.dist.lm)
plot(POLS$staminaDistance, cv.dist.res) + abline(0,0)
cor.test(POLS$TprefCV, POLS$staminaDistance, method="pearson")

plot(Male$TprefCV, Male$staminaDistance)
m.cv.dist.lm <- lm(TprefCV~staminaDistance, Male)
m.cv.dist.res <- resid(m.cv.dist.lm)
plot(Male$staminaDistance, m.cv.dist.res) + abline(0,0)
cor.test(Male$TprefCV, Male$staminaDistance, method="pearson")

plot(Female$TprefCV, Female$staminaDistance)
f.cv.dist.lm <- lm(TprefCV~staminaDistance, Female)
f.cv.dist.res <- resid(f.cv.dist.lm)
plot(Female$staminaDistance, f.cv.dist.res) + abline(0,0)
cor.test(Female$TprefCV, Female$staminaDistance, method="pearson")
# nothing


## stamina time
plot(POLS$TprefCV, POLS$staminaTime)
cv.time.lm <- lm(TprefCV~staminaTime, POLS)
cv.time.res <- resid(cv.time.lm)
plot(POLS$staminaTime, cv.time.res) + abline(0,0)
cor.test(POLS$TprefCV, POLS$staminaTime, method="pearson")

plot(Male$TprefCV, Male$staminaTime)
m.cv.time.lm <- lm(TprefCV~staminaTime, Male)
m.cv.time.res <- resid(m.cv.time.lm)
plot(Male$staminaTime, m.cv.time.res) + abline(0,0)
cor.test(Male$TprefCV, Male$staminaTime, method="pearson")

plot(Female$TprefCV, Female$staminaTime)
f.cv.time.lm <- lm(TprefCV~staminaTime, Female)
f.cv.time.res <- resid(f.cv.time.lm)
plot(Female$staminaTime, f.cv.time.res) + abline(0,0)
cor.test(Female$TprefCV, Female$staminaTime, method="pearson")
# nothing


## stamina velocity
plot(POLS$TprefCV, POLS$staminaVelocity)
cv.vel.lm <- lm(TprefCV~staminaVelocity, POLS)
cv.vel.res <- resid(cv.vel.lm)
plot(POLS$staminaVelocity, cv.vel.res) + abline(0,0)
cor.test(POLS$TprefCV, POLS$staminaVelocity, method="spearman")

plot(Male$TprefCV, Male$staminaVelocity)
m.cv.vel.lm <- lm(TprefCV~staminaVelocity, Male)
m.cv.vel.res <- resid(m.cv.vel.lm)
plot(Male$staminaVelocity, m.cv.vel.res) + abline(0,0)
cor.test(Male$TprefCV, Male$staminaVelocity, method="spearman")

plot(Female$TprefCV, Female$staminaVelocity)
cor.test(Female$TprefCV, Female$staminaVelocity, method="pearson")


########### Shuttling

## transition
plot(POLS$shuttle, POLS$transition)
shuttle.trans.lm <- lm(shuttle~transition, POLS)
shuttle.trans.res <- resid(shuttle.trans.lm)
plot(POLS$transition, shuttle.trans.res) + abline(0,0)
cor.test(POLS$shuttle, POLS$transition, method="pearson")

plot(Male$shuttle, Male$transition)
m.shut.trans.lm <- lm(shuttle~transition, Male)
m.shut.trans.res <- resid(m.shut.trans.lm)
plot(Male$transition, m.shut.trans.res) + abline(0,0)
cor.test(Male$shuttle, Male$transition, method="pearson")

plot(Female$shuttle, Female$transition)
f.shut.trans.lm <- lm(shuttle~transition, Female)
f.shut.trans.res <- resid(f.shut.trans.lm)
plot(Female$transition, f.shut.trans.res) + abline(0,0)
cor.test(Female$shuttle, Female$transition, method="pearson")
# nothing


## basking
plot(POLS$shuttle, POLS$basking)
shut.bask.lm <- lm(shuttle~basking, POLS)
shut.bask.res <- resid(shut.bask.lm)
plot(POLS$basking, shut.bask.res) + abline(0,0)
cor.test(POLS$shuttle, POLS$basking, method="spearman")

plot(Male$shuttle, Male$basking)
m.shut.bask.lm <- lm(shuttle~basking, Male)
m.shut.bask.res <- resid(m.shut.bask.lm)
plot(Male$basking, m.shut.bask.res) + abline(0,0)
cor.test(Male$shuttle, Male$basking, method="spearman")

plot(Female$shuttle, Female$basking)
f.shut.bask.lm <- lm(shuttle~basking, Female)
f.shut.bask.res <- resid(f.shut.bask.lm)
plot(Female$basking, f.shut.bask.res) + abline(0,0)
cor.test(Female$shuttle, Female$basking, method="spearman")
# nothing


## perching
plot(POLS$shuttle, POLS$perching)
shut.perch.lm <- lm(shuttle~perching, POLS)
shut.perch.res <- resid(shut.perch.lm)
plot(POLS$perching, shut.perch.res) + abline(0,0)
cor.test(POLS$shuttle, POLS$perching, method="pearson")

plot(Male$shuttle, Male$perching)
m.shut.perch.lm <- lm(shuttle~perching, Male)
m.shut.perch.res <- resid(m.shut.perch.lm)
plot(Male$perching, m.shut.perch.res) + abline(0,0)
cor.test(Male$shuttle, Male$perching, method="pearson")

plot(Female$shuttle, Female$perching)
f.shut.perch.lm <- lm(shuttle~perching, Female)
f.shut.perch.res <- resid(f.shut.perch.lm)
plot(Female$perching, f.shut.perch.res) + abline(0,0)
cor.test(Female$shuttle, Female$perching, method="spearman")
# nothing


## stamina distance
plot(POLS$shuttle, POLS$staminaDistance)
cor.test(POLS$shuttle, POLS$staminaDistance, method="pearson")

plot(Male$shuttle, Male$staminaDistance)
m.shut.dist.lm <- lm(shuttle~staminaDistance, Male)
m.shut.dist.res <- resid(m.shut.dist.lm)
plot(Male$staminaDistance, m.shut.dist.res) + abline(0,0)
cor.test(Male$shuttle, Male$staminaDistance, method="pearson")

plot(Female$shuttle, Female$staminaDistance)
cor.test(Female$shuttle, Female$staminaDistance, method="pearson")
# nothing


## stamina time
plot(POLS$shuttle, POLS$staminaTime)
cor.test(POLS$shuttle, POLS$staminaTime, method="pearson")

plot(Male$shuttle, Male$staminaTime)
cor.test(Male$shuttle, Male$staminaTime, method="pearson")

plot(Female$shuttle, Female$staminaTime)
cor.test(Female$shuttle, Female$staminaTime, method="pearson")
# nothing


## stamina velocity
plot(POLS$shuttle, POLS$staminaVelocity)
shut.vel.lm <- lm(shuttle~staminaVelocity, POLS)
shut.vel.res <- resid(shut.vel.lm)
plot(POLS$staminaVelocity, shut.vel.res) + abline(0,0)
cor.test(POLS$shuttle, POLS$staminaVelocity, method="pearson")

plot(Male$shuttle, Male$staminaVelocity)
cor.test(Male$shuttle, Male$staminaVelocity, method="pearson")

plot(Female$shuttle, Female$staminaVelocity)
cor.test(Female$shuttle, Female$staminaVelocity, method="pearson")



############## Transitions

## basking
plot(POLS$transition, POLS$basking)
trans.bask.lm <- lm(transition~basking, POLS)
trans.bask.res <- resid(trans.bask.lm)
plot(POLS$basking, trans.bask.res) + abline(0,0)
cor.test(POLS$transition, POLS$basking, method="spearman")
# positive correlation

plot(Male$transition, Male$basking)
m.trans.bask.lm <- lm(transition~basking, Male)
m.trans.bask.res <- resid(m.trans.bask.lm)
plot(Male$basking, m.trans.bask.res) + abline(0,0)
cor.test(Male$transition, Male$basking, method="spearman")
# positive

plot(Female$transition, Female$basking)
cor.test(Female$transition, Female$basking, method="spearman")


## perching
plot(POLS$transition, POLS$perching)
trans.perch.lm <- lm(transition~perching, POLS)
trans.perch.res <- resid(trans.perch.lm)
plot(POLS$perching, trans.perch.res) + abline(0,0)
cor.test(POLS$transition, POLS$perching, method="spearman")
# positive correlation

plot(Male$transition, Male$perching)
cor.test(Male$transition, Male$perching, method="spearman")
# positive

plot(Female$transition, Female$perching)
cor.test(Female$transition, Female$perching, method="spearman")


## stamina distance
plot(POLS$transition, POLS$staminaDistance)
trans.dist.lm <- lm(transition~staminaDistance, POLS)
trans.dist.res <- resid(trans.dist.lm)
plot(POLS$staminaDistance, trans.dist.res) + abline(0,0)
cor.test(POLS$transition, POLS$staminaDistance, method="spearman")

plot(Male$transition, Male$staminaDistance)
m.trans.dist.lm <- lm(transition~staminaDistance, Male)
m.trans.dist.res <- resid(m.trans.dist.lm)
plot(Male$staminaDistance, m.trans.dist.res) + abline(0,0)
cor.test(Male$transition, Male$staminaDistance, method="spearman")

plot(Female$transition, Female$staminaDistance)
f.trans.dist.lm <- lm(transition~staminaDistance, Female)
f.trans.dist.res <- resid(f.trans.dist.lm)
plot(Female$staminaDistance, f.trans.dist.res) + abline(0,0)
cor.test(Female$transition, Female$staminaDistance, method="spearman")


## stamina time
plot(POLS$transition, POLS$staminaTime)
trans.time.lm <- lm(transition~staminaTime, POLS)
trans.time.res <- resid(trans.time.lm)
plot(POLS$staminaTime, trans.time.res) + abline(0,0)
cor.test(POLS$transition, POLS$staminaTime, method="spearman")

plot(Male$transition, Male$staminaTime)
cor.test(Male$transition, Male$staminaTime, method="spearman")

plot(Female$transition, Female$staminaTime)
cor.test(Female$transition, Female$staminaTime, method="spearman")
# nothing


## stamina velocity
plot(POLS$transition, POLS$staminaVelocity)
cor.test(POLS$transition, POLS$staminaVelocity, method="spearman")

plot(Male$transition, Male$staminaVelocity)
cor.test(Male$transition, Male$staminaVelocity, method="spearman")

plot(Female$transition, Female$staminaVeloicity)
cor.test(Female$transition, Female$staminaVelocity, method="spearman")
# nothing

################## Basking

## perching
plot(POLS$basking, POLS$perching)
bask.perch.lm <- lm(basking~perching, POLS)
bask.perch.res <- resid(bask.perch.lm)
plot(POLS$perching, bask.perch.res) + abline(0,0)
cor.test(POLS$basking, POLS$perching, method="spearman")
# positive correlation

plot(Male$basking, Male$perching)
cor.test(Male$basking, Male$perching, method="spearman")

plot(Female$basking, Female$perching)
cor.test(Female$basking, Female$perching, method="spearman")


## stamina distance
plot(POLS$basking, POLS$staminaDistance)
bask.dist.lm <- lm(basking~staminaDistance, POLS)
bask.dist.res <- resid(bask.dist.lm)
plot(POLS$staminaDistance, bask.dist.res) + abline(0,0)
cor.test(POLS$basking, POLS$staminaDistance, method="spearman")

plot(Male$basking, Male$staminaDistance)
cor.test(Male$basking, Male$staminaDistance, method="spearman")

plot(Female$basking, Female$staminaDistance)
cor.test(Female$basking, Female$staminaDistance, method="spearman")


## stamina time
plot(POLS$basking, POLS$staminaTime)
bask.time.lm <- lm(basking~staminaTime, POLS)
bask.time.res <- resid(bask.time.lm)
plot(POLS$staminaTime, bask.time.res) + abline(0,0)
cor.test(POLS$basking, POLS$staminaTime, method="spearman")

plot(Male$basking, Male$staminaTime)
cor.test(Male$basking, Male$staminaTime, method="spearman")

plot(Female$basking, Female$staminaTime)
cor.test(Female$basking, Female$staminaTime, method="spearman")


## stamina velocity
plot(POLS$basking, POLS$staminaVelocity)
bask.vel.lm <- lm(basking~staminaVelocity, POLS)
bask.vel.res <- resid(bask.vel.lm)
plot(POLS$staminaVelocity, bask.vel.res) + abline(0,0)
cor.test(POLS$basking, POLS$staminaVelocity, method="spearman")

plot(Male$basking, Male$staminaVelocity)
cor.test(Male$basking, Male$staminaVelocity, method="spearman")

plot(Female$basking, Female$staminaVelocity)
cor.test(Female$basking, Female$staminaVelocity, method="spearman")


####################### perching

## stamina distance
plot(POLS$perching, POLS$staminaDistance)
cor.test(POLS$perching, POLS$staminaDistance, method="spearman")

plot(Male$perching, Male$staminaDistance)
cor.test(Male$perching, Male$staminaDistance, method="spearman")

plot(Female$perching, Female$staminaDistance)
cor.test(Female$perching, Female$staminaDistance, method="spearman")


## stamina time
plot(POLS$perching, POLS$staminaTime)
cor.test(POLS$perching, POLS$staminaTime, method="spearman")

plot(Male$perching, Male$staminaTime)
cor.test(Male$perching, Male$staminaTime, method="spearman")

plot(Female$perching, Female$staminaTime)
cor.test(Female$perching, Female$staminaTime, method="spearman")


## stamina velocity
plot(POLS$perching, POLS$staminaVelocity)
cor.test(POLS$perching, POLS$staminaVelocity, method="spearman")

plot(Male$perching, Male$staminaVelocity)
cor.test(Male$perching, Male$staminaVelocity, method="spearman")

plot(Female$perching, Female$staminaVelocity)
cor.test(Female$perching, Female$staminaTime, method="spearman")


################# stamina distance

## stamina time
plot(POLS$staminaDistance, POLS$staminaTime)
dist.time.lm <- lm(staminaDistance~staminaTime, POLS)
dist.time.res <- resid(dist.time.lm)
plot(POLS$staminaTime, dist.time.res) + abline(0,0)
cor.test(POLS$staminaDistance, POLS$staminaTime, method="pearson")
# positive correlation

plot(Male$staminaDistance, Male$staminaTime)
cor.test(Male$staminaDistance, Male$staminaTime, method="pearson")
# positive

plot(Female$staminaDistance, Female$staminaTime)
cor.test(Female$staminaDistance, Female$staminaTime, method="pearson")
# positive


## stamina velocity
plot(POLS$staminaDistance, POLS$staminaVelocity)
dist.vel.lm <- lm(staminaDistance~staminaVelocity, POLS)
dist.vel.res <- resid(dist.vel.lm)
plot(POLS$staminaVelocity, dist.vel.res) + abline(0,0)
cor.test(POLS$staminaDistance, POLS$staminaVelocity, method="pearson")

plot(Male$staminaDistance, Male$staminaVelocity)
cor.test(Male$staminaDistance, Male$staminaVelocity, method="pearson")

plot(Female$staminaDistance, Female$staminaVelocity)
cor.test(Female$staminaDistance, Female$staminaVelocity, method="pearson")


############ stamina time

## velocity
plot(POLS$staminaTime, POLS$staminaVelocity)
time.vel.lm <- lm(staminaTime~staminaVelocity, POLS)
time.vel.res <- resid(time.vel.lm)
plot(POLS$staminaVelocity, time.vel.res) + abline(0,0)
cor.test(POLS$staminaTime, POLS$staminaVelocity, method="pearson")

plot(Male$staminaTime, Male$staminaVelocity)
cor.test(Male$staminaTime, Male$staminaVelocity, method="pearson")

plot(Female$staminaTime, Female$staminaVelocity)
cor.test(Female$staminaTime, Female$staminaVelocity, method="pearson")



####### Final correlations:
# 1) Negative correlation between Tpref and shuttling: overall and in males
# 2) Negative correlation between Tpref and Tpref CV: overall and in males
# 3) Negative correlation between Tpref and basking/perching: overall and in males
# 3.5?) Tpref CV and velocity
# 4) Positive correlation between transitions and basking/perching: overall and in males
# 5) correlations between basking and perching, stamina distance and time, and stamina distance and velocity: overall and in males and in females


### So we have lots of correlations between thermal traits and correlations between our behavioral traits, but nothing with stamina and the only thermal/behavioral crossover is pretty weak



###### heat map
map <- POLS[,7:15]
View(map)
cormap <- round(cor(map),2)
head(cormap)

meltmap<-melt(cormap)
head(meltmap)

ggplot(meltmap, aes(x=Var1, y=Var2, fill=value)) + geom_title()

get_lower_tri<-function(cormatp){
  cormap[upper.tri(cormap)] <- NA
  return(cormap)
}

get_upper_tri <- function(cormap){
  cormap[lower.tri(cormap)]<- NA
  return(cormap)
}

upper_tri <- get_upper_tri(cormap)
upper_tri

melted_cormat <- melt(upper_tri, na.rm=TRUE)

ggplot(melted_cormat, aes(Var2, Var1, fill=value)) + geom_tile() + scale_fill_gradient2(low="blue", high="red", mid="white", midpoint=0, limit=c(-1,1), space="Lab", name="Correlation Coefficient") + theme_classic() + theme(axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text(angle=45, vjust=.6))



## PCA for males, morph groups
ornatus.male.pca.groups <- read.csv("~/Desktop/U. ornatus 2019/Boldness:POLS Excel Sheets/U.ornatus.male.pca.morphgroups.csv", row.names = "lizard")
View(ornatus.male.pca.groups)
pca.male.groups <- princomp(ornatus.male.pca.groups[,-10], cor=TRUE)
fviz_eig(pca.male.groups, addlabels=TRUE)
# 3 PCs
fviz_pca_var(pca.male.groups, repel = TRUE)

# individuals
fviz_pca_biplot(pca.male.groups, label ="var", repel = TRUE, habillage=ornatus.male.pca.groups$color, palette = c("brown", "blue", "gold"), legend.title = "Throat Morph", xlab="PC1 (28%)", ylab="PC2 (25%)", title="", addEllipses = FALSE, ellipse.type = "convex")


eig.val <- get_eigenvalue(pca.male.groups)
eig.val

res.var <- get_pca_var(pca.male.groups)
res.var$contrib 
loadings(pca.male.groups)


## Total variance of 3 axes: 72.6%
## PC1: 28%. Lizards with higher Tprefs have lower CV and shuttling
## PC2: 25%. Bolder lizards make more transitions and have lower Tpref
## PC3: 20%. Dominated by stamina


# What about if we look at differences based on PERMANOVA? (are centroids and dispersions equal between groups?)
all.individuals<- ornatus.male.pca.groups[,1:9]
View(all.individuals)
str(all.individuals)

treatment<-c("Y/B, O/B","Yellow/Orange","Yellow/Orange","Blue","Yellow/Orange","Blue","Yellow/Orange", "Yellow/Orange","Yellow/Orange","Y/B, O/B","Blue","Yellow/Orange","Y/B, O/B","Y/B, O/B","Yellow/Orange","Blue", "Blue", "Yellow/Orange","Blue","Yellow/Orange","Yellow/Orange", "Y/B, O/B", "Blue", "Yellow/Orange", "Yellow/Orange", "Y/B, O/B", "Y/B, O/B", "Yellow/Orange","Y/B, O/B", "Blue", "Y/B, O/B", "Y/B, O/B", "Y/B, O/B", "Y/B, O/B", "Yellow/Orange", "Yellow/Orange", "Yellow/Orange", "Y/B, O/B", "Yellow/Orange", "Y/B, O/B", "Y/B, O/B", "Blue", "Blue", "Y/B, O/B", "Blue", "Yellow/Orange", "Y/B, O/B", "Y/B, O/B", "Yellow/Orange", "Y/B, O/B", "Y/B, O/B", "Blue", "Yellow/Orange", "Yellow/Orange", "Y/B, O/B", "Y/B, O/B","Y/B, O/B")
str(treatment)
class(treatment)

ornatus.permpca <- princomp(all.individuals, cor=TRUE)

scores <- as.data.frame(cbind(ornatus.permpca$scores[,1], ornatus.permpca$scores[,2], ornatus.permpca$scores[,3]))
View(scores)
scores$site <- (rownames(scores))  
scores$group<-treatment
View(scores)


# Do morph groups vary by the centroid space between all 3 axes (3 response variables, manova)
ggplot(data=scores) + 
  stat_ellipse(aes(x=V1,y=V2,colour=treatment),level = 0.50) +
  geom_point(aes(x=V1,y=V2,shape=treatment,colour=treatment),size=4) +theme_bw()

adon.pcaresults<-adonis(all.individuals ~ treatment, method="euclidean",perm=999)
print(adon.pcaresults)

centroid.test<-cbind(scores$V1,scores$V2,scores$V3)
centroid.an<-lm(centroid.test~scores$group)
Anova(centroid.an)
# Yes... which axes are driving it?

# axis 1
group.v1<-lm(V1~group, data=scores)
Anova(group.v1)
emmeans(group.v1, specs = pairwise ~ group)
# overall difference but not pairwise...

# axis 2
group.v2<-lm(V2~group, data=scores)
Anova(group.v2)

# axis 3
group.v3<-lm(V3~group, data=scores)
Anova(group.v3)
emmeans(group.v3, specs = pairwise ~ group)
# same story as axis 1


######## PCA for males, all morphs
ornatus.male.pca.morph <- read.csv("~/Desktop/U. ornatus 2019/Boldness:POLS Excel Sheets/U.ornatus.male.pca.morph.csv", row.names="lizard")
View(ornatus.male.pca.morph)
pca.male <- princomp(ornatus.male.pca.morph[,-10], cor=TRUE)
fviz_eig(pca.male, addlabels=TRUE)
# 3 PCs
fviz_pca_var(pca.male, repel = TRUE)

# individuals
fviz_pca_biplot(pca.male, label ="var", repel = TRUE, habillage=ornatus.male.pca.morph$morph, palette = c("blue", "orange", "brown", "red", "gold", "dark green"), legend.title = "Throat Morph", xlab="PC1 (28%)", ylab="PC2 (25%)", title="", addEllipses = FALSE, ellipse.type = "convex")

scale_fill_discrete(labels=c('B','O','O/B','Y','Y/B','O/Y'))

eig.val.morph <- get_eigenvalue(pca.male)
eig.val.morph

res.var.morph <- get_pca_var(pca.male)
res.var.morph$contrib 
loadings(pca.male)

## Total variance of 3 axes: 72.6%
## PC1: 28%. Lizards with higher Tprefs have lower CV and shuttling
## PC2: 25%. Bolder lizards make more transitions and have lower Tpref
## PC3: 20%. Dominated by stamina


# What about if we look at differences based on PERMANOVA? (are centroids and dispersions equal between groups?)
all.individuals.morph<- ornatus.male.pca.morph[,1:9]
View(all.individuals.morph)
str(all.individuals.morph)

treatment.morph<-c("O/B","Y","Y","B","Y/O","B","Y", "Y","Y","Y/B","B","Y","Y/B","Y/B","Y","B", "B", "Y","B","Y","Y", "Y/B", "B", "Y", "Y", "Y/B", "Y/B", "Y/O","Y/B", "B", "O/B", "O/B", "Y/B", "Y/B", "O", "Y/O", "Y", "Y/B", "Y", "Y/B", "O/B", "B", "B", "O/B", "B", "Y", "Y/B", "Y/B", "Y/O", "Y/B", "Y/B", "B", "O", "Y", "Y/B", "O/B","Y/B")
str(treatment.morph)
class(treatment.morph)

ornatus.permpca.morph <- princomp(all.individuals.morph, cor=TRUE)

scores.morph <- as.data.frame(cbind(ornatus.permpca.morph$scores[,1], ornatus.permpca.morph$scores[,2], ornatus.permpca.morph$scores[,3]))
View(scores.morph)
scores.morph$site <- (rownames(scores.morph))  
scores.morph$group<-treatment.morph
View(scores.morph)


# Do morphs vary by the centroid space between all 3 axes (3 response variables, manova)
ggplot(data=scores.morph) + 
  stat_ellipse(aes(x=V1,y=V2,colour=treatment.morph),level = 0.50) +
  geom_point(aes(x=V1,y=V2,shape=treatment.morph,colour=treatment.morph),size=4) +theme_bw()

adon.pcaresults.morph<-adonis(all.individuals.morph ~ treatment.morph, method="euclidean",perm=999)
print(adon.pcaresults.morph)

centroid.test.morph<-cbind(scores.morph$V1,scores.morph$V2,scores.morph$V3)
centroid.an.morph<-lm(centroid.test.morph~scores.morph$group)
Anova(centroid.an.morph)
# No... let's look at axes anyway since close (curiosity)

# axis 1
group.v1.morph<-lm(V1~group, data=scores.morph)
Anova(group.v1.morph)

# axis 2
group.v2.morph<-lm(V2~group, data=scores.morph)
Anova(group.v2.morph)

# axis 3
group.v3.morph<-lm(V3~group, data=scores.morph)
Anova(group.v3.morph)
# Nothing



########## Let's do a PCA for females
female.pca <- read.csv("~/Desktop/U. ornatus 2019/Boldness:POLS Excel Sheets/U.ornatus.female.pca.csv", row.names="lizard")
View(female.pca)
pca.female <- princomp(female.pca[,-10], cor=TRUE)
fviz_eig(pca.female, addlabels=TRUE)
# 4 PCs
fviz_pca_var(pca.female, repel = TRUE)

# individuals
fviz_pca_biplot(pca.female, label ="var", repel = TRUE, habillage=female.pca$morph, palette = c("orange", "black", "gold"), legend.title = "Throat Morph", xlab="PC1 (29%)", ylab="PC2 (22%)", title="", addEllipses = FALSE, ellipse.type = "convex")

eig.val.female <- get_eigenvalue(pca.female)
eig.val.female

res.var.female <- get_pca_var(pca.female)
res.var.female$contrib 
loadings(pca.female)

## Total variance of 4 axes: 81.6%
## PC1: 29%. 
## PC2: 22%. 
## PC3: 16%. 
## PC4: 15%
############# But these all have loadings in the 0.2-0.4 range. Nothing really jumps out


# What about if we look at differences based on PERMANOVA? (are centroids and dispersions equal between groups?)
all.individuals.female<- female.pca[,1:9]
View(all.individuals.female)
str(all.individuals.female)

treatment.female<-c("O","W","O","W","O","Y","O", "Y","O","Y","Y","O","O","O","O","Y", "Y", "Y","Y","Y","Y", "O")
str(treatment.female)
class(treatment.female)

ornatus.permpca.female <- princomp(all.individuals.female, cor=TRUE)

scores.female <- as.data.frame(cbind(ornatus.permpca.female$scores[,1], ornatus.permpca.female$scores[,2], ornatus.permpca.female$scores[,3], ornatus.permpca.female$scores[,4]))
View(scores.female)
scores.female$site <- (rownames(scores.female))  
scores.female$group<-treatment.female
View(scores.female)


# Do morph groups vary by the centroid space between all 3 axes (3 response variables, manova)
ggplot(data=scores.female) + 
  stat_ellipse(aes(x=V1,y=V2,colour=treatment.female),level = 0.50) +
  geom_point(aes(x=V1,y=V2,shape=treatment.female,colour=treatment.female),size=4) +theme_bw()

adon.pcaresults.female<-adonis(all.individuals.female ~ treatment.female, method="euclidean",perm=999)
print(adon.pcaresults.female)

centroid.test.female<-cbind(scores.female$V1,scores.female$V2,scores.female$V3)
centroid.an.female<-lm(centroid.test.female~scores.female$group)
Anova(centroid.an.female)
# Nothing whatsoever

# axis 1
group.v1.morph.f<-lm(V1~group, data=scores.female)
Anova(group.v1.morph.f)

# axis 2
group.v2.morph.f<-lm(V2~group, data=scores.female)
Anova(group.v2.morph.f)

# axis 3
group.v3.morph.f<-lm(V3~group, data=scores.female)
Anova(group.v3.morph.f)
# Nothing

# axis 4
group.v4.morph.f<-lm(V4~group, data=scores.female)
Anova(group.v4.morph.f)