# Color Overlap and Pattern Energy Code
rm(list=ls())
cam<-read.table('camo.txt', header = T, na.string="NA")
attach(cam)

#transformations
shapiro.test(color.overlap.dorsal)
shapiro.test(color.overlap.lateral)
shapiro.test(pattern.energy.dif.dorsal)
shapiro.test(pattern.energy.dif.lateral)
col.dor <- sqrt(color.overlap.dorsal)
col.lat <- sqrt(color.overlap.lateral)
library(MASS)
boxcox(pattern.energy.dif.dorsal~as.factor(treatment)+svl, lambda = seq(-0.4,-0.3, 0.01)) 
pat.dor <-(((PE.dorsal)^(-.38))-1)/(-.38)
pat.lat <- log(pattern.energy.dif.lateral)
cube.cond <-resid(lm((mass^(1/3))~(svl), data=cam, na.action = na.exclude))
cam3 <- cbind(cam, col.dor, col.lat, pat.dor, pat.lat, cube.cond)

library(emmeans)

#final models (original models also included sex, condition, and body temperature as predictors)
#dorsal color overlap
m1 <- lm(col.dor~as.factor(treatment)+svl, data=cam, na.action = na.exclude)
m1.notreat <- lm(col.dor~svl, data=cam, na.action = na.exclude)
m1.nosvl <- lm(col.dor~as.factor(treatment), data=cam, na.action = na.exclude)
summary(m1)
anova(m1, m1.notreat, test = "F")
anova(m1, m1.nosvl, test = "F")
mod <- emmeans(m1, pairwise ~ as.factor(treatment))
pairs(mod, simple = "each")
vif(m1)

#lateral color overlap
m2 <- lm(col.lat~as.factor(treatment)+svl, data=cam, na.action = na.exclude)
m2.notreat <- lm(col.lat~svl, data=cam, na.action = na.exclude)
m2.nosvl <- lm(col.lat~as.factor(treatment), data=cam, na.action = na.exclude)
summary(m2)
anova(m2, m12.notreat, test = "F")
anova(m2, m12.nosvl, test = "F")
mod2 <- emmeans(m2, pairwise ~ as.factor(treatment))
pairs(mod2, simple = "each")
vif(m2)

#dorsal pattern matching
m13 <- lm(pat.dor~as.factor(treatment)+svl, data=cam, na.action = na.exclude)
summary(m13)
vif(m13)

m14 <- lm(pat.lat~as.factor(treatment)+svl, data=cam, na.action = na.exclude)
summary(m14)
vif(m14)

# ----------
# Edge Disruption (Gabrat) Code
rm(list=ls())
gab<-read.table('gab.txt', header = T, na.string="NA")
attach(gab)

#transformations
cube.cond <-resid(lm((mass^(1/3))~(svl), data=gab, na.action = na.exclude))
library(MASS)
boxcox(gabrat~as.factor(treatment)+svl, lambda = seq(0.1,0.15, 0.01)) 
gabrat.bc <-(((gabrat)^(.12))-1)/(.12)
gab2 <- cbind(gab, cube.cond, gabrat.bc)

body.only <- subset(gab2, region == "lateral")
dors.only <- subset(gab2, region == "dorsal")

library(lmerTest)
library(emmeans)

#final models (original models also included sex, condition, and body temperature as predictors)
#dorsal gabrat analysis
mbx <- lmer(log(gabrat)~as.factor(treatment)+svl+(1|lizard.id), data= dors.only, na.action = na.exclude, REML = FALSE)
mbx.notreat <- lmer(log(gabrat)~svl+(1|lizard.id), data= dors.only, na.action = na.exclude, REML = FALSE)
mbx.nosvl <- lmer(log(gabrat)~as.factor(treatment)+(1|lizard.id), data= dors.only, na.action = na.exclude, REML = FALSE)
anova(mbx,mbx.notreat) 
anova(mbx,mbx.nosvl) 
mbxe <- emmeans(mbx, pairwise ~ as.factor(treatment), pbkrtest.limit = 6700)
pairs(mbxe, simple = "each")
vif(mbx)

#lateral gabrat analysis
mcx <- lmer(log(gabrat)~as.factor(treatment)+svl+(1|lizard.id), data= body.only, na.action = na.exclude, REML = FALSE)
mcx.notreat <- lmer(log(gabrat)~svl+(1|lizard.id), data= body.only, na.action = na.exclude, REML = FALSE)
mcx.nosvl <- lmer(log(gabrat)~as.factor(treatment)+(1|lizard.id), data= body.only, na.action = na.exclude, REML = FALSE)
anova(mcx,mcx.notreat) 
anova(mcx,mcx.nosvl) 
mcxe <- emmeans(mcx, pairwise ~ as.factor(treatment), pbkrtest.limit = 6800)
pairs(mcxe, simple = "each")
vif(mcx)

# ----------
# Lizard Body Color Analysis - Supplementary Analysis 1

rm(list=ls())
cam2<-read.table('camrev.txt', header = T, na.string="NA")
attach(cam2)

new1 <- lm((log(luminance.final))~as.factor(treatment)+(luminance.initial), data=cam2, na.action = na.exclude)
summary(new1)
library(car)
Anova(new1, test.statistic="LR", type = 2)
n1 <- emmeans(new1, pairwise ~ as.factor(treatment))
pairs(n1, simple = "each")

new2 <- lm((redness.final)~as.factor(treatment)+redness.initial, data=cam2, na.action = na.exclude)
summary(new2)
Anova(new2, test.statistic="LR", type = 2)

new3 <- lm((blueness.final)~as.factor(treatment)+blueness.initial, data=cam2, na.action = na.exclude)
summary(new3)
Anova(new3, test.statistic="LR", type = 2)
n3 <- emmeans(new5, pairwise ~ as.factor(treatment))
pairs(n3, simple = "each")

diff <- lm((luminance.final-luminance.initial)~as.factor(treatment)+luminance.initial, data=cam2, na.action = na.exclude)
summary(diff)
Anova(diff, test.statistic="LR", type = 2)
nd <- emmeans(diff, pairwise ~ as.factor(treatment))
pairs(nd, simple = "each")

diff2 <- lm((redness.final-redness.initial)~as.factor(treatment)+redness.initial, data=cam2, na.action = na.exclude)
summary(diff2)
Anova(diff2, test.statistic="LR", type = 2)

diff3 <- lm((blueness.final-blueness.initial)~as.factor(treatment)+blueness.initial, data=cam2, na.action = na.exclude)
summary(diff3)
Anova(diff3, test.statistic="LR", type = 2)
n3 <- emmeans(diff3, pairwise ~ as.factor(treatment))
pairs(n3, simple = "each")