#####################################################
#####Final analysis file for Macquarie Springtails###
#####################################################
setwd("C:/RBook/")
library(plyr)
library(lme4)
library(caper)
library(ggplot2)
library(multcomp)
library(MASS)
library(lmodel2)

###########################
##Body mass (fresh) means##
###########################
FM<-read.csv("Macquarie_Collembola_Mass.csv",header=TRUE)
str(FM)

FMmeans<-ddply(FM, c("Species"), summarise,
               N = length(Mass),
               mean = mean(Mass),
               sd = sd(Mass),
               max = max(Mass),
               min = min(Mass))

#Use these means for mass in species-level comparisons for CTs
#Entered into species means files

###############################################
###Critical thermal limits#####################
###############################################
###Individual data for F0 and 10°C Acclimation###
#CTmax F0 10C Acclimation only###
MqCTmaxind<-read.csv("CTmax_F0_10.csv",header=TRUE)
model1<-lm(MqCTmaxind$CTmax~MqCTmaxind$Status)
summary(model1)

boxplot(MqCTmaxind$CTmax~MqCTmaxind$Status)

CTmax1<-ddply(MqCTmaxind, c("Status"), summarise,
              N = length(CTmax),
              mean = mean(CTmax),
              sd = sd(CTmax),
              max = max(CTmax),
              min = min(CTmax))

Mxden<-ggplot(MqCTmaxind, aes(x=CTmax, colour=Status, fill=Status))+geom_density()+facet_grid(cols=vars(Status), scales="fixed")+ylab("Density")
Mxden
Mxden<-Mxden+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()+scale_fill_colorblind()+theme(legend.positio='none')
Mxden
Mxden<-Mxden+scale_x_continuous(breaks = seq(20,45,by = 2)) 
Mxden

#CTmin F010C Acclimation only###
MqCTminind<-read.csv("CTmin_F0_10.csv",header=TRUE)
str(MqCTminind)
model2<-lm(MqCTminind$CTmin~MqCTminind$Status)
summary(model2)
boxplot(MqCTminind$CTmin~MqCTminind$Status)
CTmin1<-ddply(MqCTminind, c("Status"), summarise,
              N = length(CTmin),
              mean = mean(CTmin),
              sd = sd(CTmin),
              max = max(CTmin),
              min = min(CTmin))
Mnden<-ggplot(MqCTminind, aes(x=CTmin, colour=Status, fill=Status))+geom_density()+facet_grid(cols=vars(Status), scales="fixed")+ylab("Density")
Mnden
Mnden<-Mnden+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()+scale_fill_colorblind()+theme(legend.positio='none')
Mnden
Mnden<-Mnden+scale_x_continuous(breaks = seq(-8,2,by = 2)) 
Mnden

###PGLS approach###
#Phylogenetic Generalized Least Squares
#Tree built using (Grafen 1989 PTRSB)
#Tree exported to Newick format with branchlengths and imported here
setwd("C:/RBook/")

library(caper)
max.data<-read.csv("CTmax_means_F0_10.csv",header=TRUE)
max.tree<-read.tree("CTsF0.tre")
str(max.data)

#CTmax
pg.model1<-pgls(CTmax~Status*Mass, comparative.data(max.tree, max.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model1)

pg.model1nint<-pgls(CTmax~Status, comparative.data(max.tree, max.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model1nint)

###PGLS approach###
min.data<-read.csv("CTmin_means_F0_10.csv",header=TRUE)
min.tree<-read.tree("CTsF0.tre")
str(min.data)

pg.model2<-pgls(CTmin~Status*Mass, comparative.data(min.tree, min.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model2)

pg.model2nint<-pgls(CTmin~Status+Mass, comparative.data(min.tree, min.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model2nint)

pg.model2final<-pgls(CTmin~Status, comparative.data(min.tree, min.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model2final)

###PGLS approach###
range.data<-read.csv("CTRange_means_F0_10.csv",header=TRUE)
range.tree<-read.tree("CTsF0.tre")
str(range.data)

pg.model3<-pgls(CTrange~Status*Mass, comparative.data(range.tree, range.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model2)

pg.model3nint<-pgls(CTrange~Status+Mass, comparative.data(range.tree, range.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model3nint)

pg.model3final<-pgls(CTrange~Status, comparative.data(range.tree, range.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.model3final)

#############################################
##Compare F0 and F2 CTs######################
##Can only use 10C Acclimation data##########
#############################################
###PGLS and Model II regression approaches###
#Phylogenetic Generalized Least Squares
#Tree built using (Grafen 1989 PTRSB)
#Tree exported to Newick format with branchlengths and imported here
setwd("C:/RBook/")

library(caper)
F0F2.data<-read.csv("F0F2mxmnmean.csv",header=TRUE)
F0F2.tree<-read.tree("CTsF2.tre")
str(F0F2.data)

#CTmax
#means
pg.F0F2x1<-pgls(F0CTmax~F2CTmax, comparative.data(F0F2.tree, F0F2.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.F0F2x1)
lm.F0F2x1<-lm(F0F2.data$F0CTmax~F0F2.data$F2CTmax)
summary(lm.F0F2x1)
confint(lm.F0F2x1)

library(lmodel2)
Model2r1<-lmodel2(F0CTmax~F2CTmax, data=F0F2.data, "relative", "relative", 99)
Model2r1

#CTmin
#means
pg.F0F2m1<-pgls(F0CTmin~F2CTmin, comparative.data(F0F2.tree, F0F2.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.F0F2m1)
lm.F0F2m1<-lm(F0F2.data$F0CTmin~F0F2.data$F2CTmin)
summary(lm.F0F2m1)
confint(lm.F0F2m1)

Model2r2<-lmodel2(F0CTmin~F2CTmin, data=F0F2.data, "interval", "interval", 99)
Model2r2

#CTrange
#means
pg.F0F2r1<-pgls(F0range~F2range, comparative.data(F0F2.tree, F0F2.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.F0F2r1)
lm.F0F2r1<-lm(F0F2.data$F0range~F0F2.data$F2range)
summary(lm.F0F2r1)
confint(lm.F0F2r1)
Model2r3<-lmodel2(F0range~F2range, data=F0F2.data, "relative", "relative", 99)
Model2r3

############################
##F2 Acclimation data#######
############################
##ARR by species for CTmax##
############################
F2max<-read.csv("CTmax_F2_LinAcc.csv",header=TRUE)
str(F2max)

#C.denticulata
F2cdent<-F2max[ which(F2max$Species=='C. denticulata'), ]
str(F2cdent)
model12<-lm(F2cdent$CTmax~F2cdent$Acclimation)
summary(model12)
plot(F2cdent$CTmax~F2cdent$Acclimation)

#H. purpurescens
F2hpurp<-F2max[ which(F2max$Species=='H. purpurescens'), ]
str(F2hpurp)
model13<-lm(F2hpurp$CTmax~F2hpurp$Acclimation)
summary(model13)
plot(F2hpurp$CTmax~F2hpurp$Acclimation)

#H. viatica
F2hviat<-F2max[ which(F2max$Species=='H. viatica'), ]
str(F2hviat)
model14<-lm(F2hviat$CTmax~F2hviat$Acclimation)
summary(model14)
plot(F2hviat$CTmax~F2hviat$Acclimation)

#L. violacea
F2lviol<-F2max[ which(F2max$Species=='L. violaceus'), ]
str(F2lviol)
model15<-lm(F2lviol$CTmax~F2lviol$Acclimation)
summary(model15)
plot(F2lviol$CTmax~F2lviol$Acclimation)

#M. caeca
F2mcaec<-F2max[ which(F2max$Species=='M. caeca'), ]
str(F2mcaec)
model16<-lm(F2mcaec$CTmax~F2mcaec$Acclimation)
summary(model16)
plot(F2mcaec$CTmax~F2mcaec$Acclimation)

#P. fimata
F2pfima<-F2max[ which(F2max$Species=='P. fimata'), ]
str(F2pfima)
model17<-lm(F2pfima$CTmax~F2pfima$Acclimation)
summary(model17)
plot(F2pfima$CTmax~F2pfima$Acclimation)

#P. insularis
F2pinsu<-F2max[ which(F2max$Species=='P. insularis'), ]
str(F2pinsu)
model18<-lm(F2pinsu$CTmax~F2pinsu$Acclimation)
summary(model18)
plot(F2pinsu$CTmax~F2pinsu$Acclimation)

#P.notablis
F2pnotab<-F2max[ which(F2max$Species=='P. notabilis'), ]
str(F2pnotab)
model19<-lm(F2pnotab$CTmax~F2pnotab$Acclimation)
summary(model19)
plot(F2pnotab$CTmax~F2pnotab$Acclimation)

#Proisotoma 
F2prois<-F2max[ which(F2max$Species=='Proisotoma sp.'), ]
str(F2prois)
model20<-lm(F2prois$CTmax~F2prois$Acclimation)
summary(model20)
plot(F2prois$CTmax~F2prois$Acclimation)

#T. bisetosa
F2tbise<-F2max[ which(F2max$Species=='T. bisetosa'), ]
str(F2tbise)
model21<-lm(F2tbise$CTmax~F2tbise$Acclimation)
summary(model21)
plot(F2tbise$CTmax~F2tbise$Acclimation)

###############################
####ARR by species for CTmin###
###############################
F2min<-read.csv("CTmin_F2_LinAcc.csv",header=TRUE)
str(F2min)

#C.denticulata
F2cdent<-F2min[ which(F2min$Species=='C. denticulata'), ]
str(F2cdent)
model22<-lm(F2cdent$CTmin~F2cdent$Acclimation)
summary(model22)
plot(F2cdent$CTmin~F2cdent$Acclimation)

#H. purpurescens
F2hpurp<-F2min[ which(F2min$Species=='H. purpurescens'), ]
str(F2hpurp)
model23<-lm(F2hpurp$CTmin~F2hpurp$Acclimation)
summary(model23)
plot(F2hpurp$CTmin~F2hpurp$Acclimation)

#H. viatica
F2hviat<-F2min[ which(F2min$Species=='H. viatica'), ]
str(F2hviat)
model24<-lm(F2hviat$CTmin~F2hviat$Acclimation)
summary(model24)
plot(F2hviat$CTmin~F2hviat$Acclimation)

#L. violacea
F2lviol<-F2min[ which(F2min$Species=='L. violaceus'), ]
str(F2lviol)
model25<-lm(F2lviol$CTmin~F2lviol$Acclimation)
summary(model25)
plot(F2lviol$CTmin~F2lviol$Acclimation)

#M. caeca
F2mcaec<-F2min[ which(F2min$Species=='M. caeca'), ]
str(F2mcaec)
model26<-lm(F2mcaec$CTmin~F2mcaec$Acclimation)
summary(model26)
plot(F2mcaec$CTmin~F2mcaec$Acclimation)

#P. fimata
F2pfima<-F2min[ which(F2min$Species=='P. fimata'), ]
str(F2pfima)
model27<-lm(F2pfima$CTmin~F2pfima$Acclimation)
summary(model27)
plot(F2pfima$CTmin~F2pfima$Acclimation)

#P. insularis
F2pinsu<-F2min[ which(F2min$Species=='P. insularis'), ]
str(F2pinsu)
model28<-lm(F2pinsu$CTmin~F2pinsu$Acclimation)
summary(model28)
plot(F2pinsu$CTmin~F2pinsu$Acclimation)

#P.notablis
F2pnotab<-F2min[ which(F2min$Species=='P. notabilis'), ]
str(F2pnotab)
model29<-lm(F2pnotab$CTmin~F2pnotab$Acclimation)
summary(model29)
plot(F2pnotab$CTmin~F2pnotab$Acclimation)

#Proisotoma 
F2prois<-F2min[ which(F2min$Species=='Proisotoma sp.'), ]
str(F2prois)
model30<-lm(F2prois$CTmin~F2prois$Acclimation)
summary(model30)
plot(F2prois$CTmin~F2prois$Acclimation)

#T. bisetosa
F2tbise<-F2min[ which(F2min$Species=='T. bisetosa'), ]
str(F2tbise)
model31<-lm(F2tbise$CTmin~F2tbise$Acclimation)
summary(model31)
plot(F2tbise$CTmin~F2tbise$Acclimation)

################################
##PGLS ARR comparisons###########
#################################
library(caper)
ARRF2.data<-read.csv("ARRdataF2.csv",header=TRUE)
ARRF2.tree<-read.tree("CTsF2.tre")
str(ARRF2.data)

#CTmax
pg.ARRx1<-pgls(CTmaxARR~Status, comparative.data(ARRF2.tree, ARRF2.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.ARRx1)

ARRCTmax<-ddply(ARRF2.data, c("Status"), summarise,
                N = length(CTmaxARR),
                mean = mean(CTmaxARR),
                sd = sd(CTmaxARR),
                max = max(CTmaxARR),
                min = min(CTmaxARR))



#CTmin
pg.ARRn1<-pgls(CTminARR~Status, comparative.data(ARRF2.tree, ARRF2.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.ARRn1)


ARRCTmin<-ddply(ARRF2.data, c("Status"), summarise,
                N = length(CTminARR),
                mean = mean(CTminARR),
                sd = sd(CTminARR),
                max = max(CTminARR),
                min = min(CTminARR))

#########################################
###Compare ARR for CTmin and for CTmax###
#########################################
ARRcomp<-read.csv("ARRdataF2mxmncomp.csv",header=TRUE)
str(ARRcomp)
Compmod1<-lm(ARRcomp$ARR~ARRcomp$Type)
summary(Compmod1)
ARRcomp2<-ddply(ARRcomp, c("Type"), summarise,
                N = length(ARR),
                mean = mean(ARR),
                sd = sd(ARR),
                max = max(ARR),
                min = min(ARR))

#########################################
###Extreme events acclimation effects####
#########################################
#########
##CTmax##
#########

F2allmax<-read.csv("CTmax_F2_AllAcc.csv",header=TRUE)
str(F2allmax)
F2allmax$Acclimation<-as.factor(F2allmax$Acclimation)
str(F2allmax)
library(multcomp)

#C.denticulata
F2cdent2<-F2allmax[ which(F2allmax$Species=='C. denticulata'), ]
str(F2cdent2)
model32<-lm(CTmax~Acclimation, data=F2cdent2)
summary(model32)
l32<-glht(model32, linfct = mcp(Acclimation = "Tukey"))
summary(l32)
boxplot(F2cdent2$CTmax~F2cdent2$Acclimation)
mAxF2 <- ggplot(F2cdent2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF2
mAxF2+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#H. purpurescens
F2hpurp2<-F2allmax[ which(F2allmax$Species=='H. purpurescens'), ]
str(F2hpurp2)
model33<-lm(CTmax~Acclimation, data=F2hpurp2)
summary(model33)
l33<-glht(model33, linfct = mcp(Acclimation = "Tukey"))
summary(l33)
boxplot(F2hpurp2$CTmax~F2hpurp2$Acclimation)
mAxF3 <- ggplot(F2hpurp2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF3
mAxF3+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#H. viatica
F2hviat2<-F2allmax[ which(F2allmax$Species=='H. viatica'), ]
str(F2hviat2)
model34<-lm(CTmax~Acclimation, data=F2hviat2)
summary(model34)
l34<-glht(model34, linfct = mcp(Acclimation = "Tukey"))
summary(l34)
boxplot(F2hviat2$CTmax~F2hviat2$Acclimation)
mAxF4 <- ggplot(F2hviat2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF4
mAxF4+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#L. violacea
F2lviol2<-F2allmax[ which(F2allmax$Species=='L. violaceus'), ]
str(F2lviol2)
model35<-lm(CTmax~Acclimation, data=F2lviol2)
summary(model35)
l35<-glht(model35, linfct = mcp(Acclimation = "Tukey"))
summary(l35)
boxplot(F2lviol2$CTmax~F2lviol2$Acclimation)
mAxF5 <- ggplot(F2lviol2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF5
mAxF5+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#M. caeca
F2mcaec2<-F2allmax[ which(F2allmax$Species=='M. caeca'), ]
str(F2mcaec2)
model36<-lm(CTmax~Acclimation, data=F2mcaec2)
summary(model36)
l36<-glht(model36, linfct = mcp(Acclimation = "Tukey"))
summary(l36)
boxplot(F2mcaec2$CTmax~F2mcaec2$Acclimation)
mAxF6 <- ggplot(F2mcaec2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF6
mAxF6+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#P. fimata
F2pfima2<-F2allmax[ which(F2allmax$Species=='P. fimata'), ]
str(F2pfima2)
model37<-lm(CTmax~Acclimation, data=F2pfima2)
summary(model37)
l37<-glht(model37, linfct = mcp(Acclimation = "Tukey"))
summary(l37)
boxplot(F2pfima2$CTmax~F2pfima2$Acclimation)
mAxF7 <- ggplot(F2pfima2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF7
mAxF7+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#P. insularis
F2pinsu2<-F2allmax[ which(F2allmax$Species=='P. insularis'), ]
str(F2pinsu2)
model38<-lm(CTmax~Acclimation, data=F2pinsu2)
summary(model38)
l38<-glht(model38, linfct = mcp(Acclimation = "Tukey"))
summary(l38)
boxplot(F2pinsu2$CTmax~F2pinsu2$Acclimation)
mAxF8 <- ggplot(F2pinsu2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF8
mAxF8+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#P.notablis
F2pnotab2<-F2allmax[ which(F2allmax$Species=='P. notabilis'), ]
str(F2pnotab2)
model39<-lm(CTmax~Acclimation, data=F2pnotab2)
summary(model39)
l39<-glht(model39, linfct = mcp(Acclimation = "Tukey"))
summary(l39)
boxplot(F2pnotab2$CTmax~F2pnotab2$Acclimation)
mAxF9 <- ggplot(F2pnotab2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF9
mAxF9+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#Proisotoma 
F2prois2<-F2allmax[ which(F2allmax$Species=='Proisotoma sp.'), ]
str(F2prois2)
model40<-lm(CTmax~Acclimation, data=F2prois2)
summary(model40)
l40<-glht(model40, linfct = mcp(Acclimation = "Tukey"))
summary(l40)
boxplot(F2prois2$CTmax~F2prois2$Acclimation)
mAxF10 <- ggplot(F2prois2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF10
mAxF10+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#T. bisetosa
F2tbise2<-F2allmax[ which(F2allmax$Species=='T. bisetosa'), ]
str(F2tbise2)
model41<-lm(CTmax~Acclimation, data=F2tbise2)
summary(model41)
l41<-glht(model41, linfct = mcp(Acclimation = "Tukey"))
summary(l41)
boxplot(F2tbise2$CTmax~F2tbise2$Acclimation)
mAxF11 <- ggplot(F2tbise2, aes(x=Acclimation, y=CTmax)) + 
  geom_boxplot()
mAxF11
mAxF11+ylab("CTmax (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#########
##CTmin##
#########
F2allmin<-read.csv("CTmin_F2_AllAcc.csv",header=TRUE)
str(F2allmin)
F2allmin$Acclimation<-as.factor(F2allmin$Acclimation)
str(F2allmin)

#C.denticulata
F2cdent3<-F2allmin[ which(F2allmin$Species=='C. denticulata'), ]
str(F2cdent3)
model42<-lm(CTmin~Acclimation, data=F2cdent3)
summary(model42)
l42<-glht(model42, linfct = mcp(Acclimation = "Tukey"))
summary(l42)
boxplot(F2cdent3$CTmin~F2cdent3$Acclimation)
mAxF12 <- ggplot(F2cdent3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF12
mAxF12+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#H. purpurescens
F2hpurp3<-F2allmin[ which(F2allmin$Species=='H. purpurescens'), ]
str(F2hpurp3)
model43<-lm(CTmin~Acclimation, data=F2hpurp3)
summary(model43)
l43<-glht(model43, linfct = mcp(Acclimation = "Tukey"))
summary(l43)
boxplot(F2hpurp3$CTmin~F2hpurp3$Acclimation)
mAxF13 <- ggplot(F2hpurp3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF13
mAxF13+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#H. viatica
F2hviat3<-F2allmin[ which(F2allmin$Species=='H. viatica'), ]
str(F2hviat3)
model44<-lm(CTmin~Acclimation, data=F2hviat3)
summary(model44)
l44<-glht(model44, linfct = mcp(Acclimation = "Tukey"))
summary(l44)
boxplot(F2hviat3$CTmin~F2hviat3$Acclimation)
mAxF14 <- ggplot(F2hviat3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF14
mAxF14+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#L. violacea
F2lviol3<-F2allmin[ which(F2allmin$Species=='L. violaceus'), ]
str(F2lviol3)
model45<-lm(CTmin~Acclimation, data=F2lviol3)
summary(model45)
l45<-glht(model45, linfct = mcp(Acclimation = "Tukey"))
summary(l45)
boxplot(F2lviol3$CTmin~F2lviol3$Acclimation)
mAxF15 <- ggplot(F2lviol3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF15
mAxF15+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#M. caeca
F2mcaec3<-F2allmin[ which(F2allmin$Species=='M. caeca'), ]
str(F2mcaec3)
model46<-lm(CTmin~Acclimation, data=F2mcaec3)
summary(model46)
l46<-glht(model46, linfct = mcp(Acclimation = "Tukey"))
summary(l46)
boxplot(F2mcaec3$CTmin~F2mcaec3$Acclimation)
mAxF16 <- ggplot(F2mcaec3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF16
mAxF16+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#P. fimata
F2pfima3<-F2allmin[ which(F2allmin$Species=='P. fimata'), ]
str(F2pfima3)
model47<-lm(CTmin~Acclimation, data=F2pfima3)
summary(model47)
l47<-glht(model47, linfct = mcp(Acclimation = "Tukey"))
summary(l47)
boxplot(F2pfima3$CTmin~F2pfima3$Acclimation)
mAxF17 <- ggplot(F2pfima3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF17
mAxF17+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#P. insularis
F2pinsu3<-F2allmin[ which(F2allmin$Species=='P. insularis'), ]
str(F2pinsu3)
model48<-lm(CTmin~Acclimation, data=F2pinsu3)
summary(model48)
l48<-glht(model48, linfct = mcp(Acclimation = "Tukey"))
summary(l48)
boxplot(F2pinsu3$CTmin~F2pinsu3$Acclimation)
mAxF18 <- ggplot(F2pinsu3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF18
mAxF18+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#P.notablis
F2pnotab3<-F2allmin[ which(F2allmin$Species=='P. notabilis'), ]
str(F2pnotab3)
model49<-lm(CTmin~Acclimation, data=F2pnotab3)
summary(model49)
l49<-glht(model49, linfct = mcp(Acclimation = "Tukey"))
summary(l49)
boxplot(F2pnotab3$CTmin~F2pnotab3$Acclimation)
mAxF19 <- ggplot(F2pnotab3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF19
mAxF19+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#Proisotoma 
F2prois3<-F2allmin[ which(F2allmin$Species=='Proisotoma sp.'), ]
str(F2prois3)
model50<-lm(CTmin~Acclimation, data=F2prois3)
summary(model50)
l50<-glht(model50, linfct = mcp(Acclimation = "Tukey"))
summary(l50)
boxplot(F2prois3$CTmin~F2prois3$Acclimation)
mAxF20 <- ggplot(F2prois3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF20
mAxF20+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

#T. bisetosa
F2tbise3<-F2allmin[ which(F2allmin$Species=='T. bisetosa'), ]
str(F2tbise3)
model51<-lm(CTmin~Acclimation, data=F2tbise3)
summary(model51)
l51<-glht(model51, linfct = mcp(Acclimation = "Tukey"))
summary(l51)
boxplot(F2tbise3$CTmin~F2tbise3$Acclimation)
mAxF21 <- ggplot(F2tbise3, aes(x=Acclimation, y=CTmin)) + 
  geom_boxplot()
mAxF21
mAxF21+ylab("CTmin (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')

###################################################
#############Dessication data analysis#############
###################################################
ssetwd("C:/RBook/")
library(plyr)

MqF0des<-read.csv("Macquarie_Desiccation_F0.csv",header=TRUE)
str(MqF0des)
MqF0des$Acclimation<-as.factor(MqF0des$Acclimation)
MqF0des$Test.Temp<-as.factor(MqF0des$Test.Temp)
str(MqF0des)
MqF0des$logtime<-log10(MqF0des$Time.to.death)
MqF0des$logmass<-log10(MqF0des$Individual.dry.mass)
str(MqF0des)

###Histogram
library(ggplot2)
Deshisto<-ggplot(MqF0des, aes(x=Time.to.death)) + 
  geom_histogram(binwidth=10, color="black", fill="white")
Deshisto
Deshisto+ylab("Number of individuals")+xlab("Time to death (minutes)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))+theme(legend.position='none')
Deshisto<-Deshisto+theme_few(base_size=14, base_family = "")
Deshisto
Deshisto2<-ggplot(MqF0des, aes(x=logtime)) + 
  geom_histogram(binwidth=0.1, colour="black", fill = "white")
Deshisto2<-Deshisto2+theme_few(base_size=14, base_family = "")
Deshisto2

#############################
##F0 versus F2 comparisons###
#############################
MqF0F2des<-read.csv("Macquarie_Desiccation_F0_F2.csv",header=TRUE)
str(MqF0F2des)

MqF0F2des$Acclimation<-as.factor(MqF0F2des$Acclimation)
MqF0F2des$Test.Temp<-as.factor(MqF0F2des$Test.Temp)
str(MqF0F2des)
MqF0F2des$logtime<-log10(MqF0F2des$Time.to.death)
MqF0F2des$logmass<-log10(MqF0F2des$Individual.dry.mass)
str(MqF0F2des)

######################
##Individual species##
######################
#M. caeca
F0F2dmcaec<-MqF0F2des[ which(MqF0F2des$Species=='M. caeca'), ]
str(F0F2dmcaec)
dfm9<-glm(F0F2dmcaec$Time.to.death~F0F2dmcaec$Generation, quasipoisson(log))
summary(dfm9)

#C. denticulata
F0F2dcdent<-MqF0F2des[ which(MqF0F2des$Species=='C. denticulata' & MqF0F2des$Acclimation==10), ]
str(F0F2dcdent)
dfm10<-glm(F0F2dcdent$Time.to.death~F0F2dcdent$Generation, quasipoisson(log))
summary(dfm10)

#P. fimata
F0F2dpfima<-MqF0F2des[ which(MqF0F2des$Species=='P. fimata' & MqF0F2des$Acclimation==10), ]
str(F0F2dpfima)
dfm11<-glm(F0F2dpfima$Time.to.death~F0F2dpfima$Generation, quasipoisson(log))
summary(dfm11)

#P. insularis
F0F2dpinsu<-MqF0F2des[ which(MqF0F2des$Species=='P. insularis' & MqF0F2des$Acclimation==10), ]
str(F0F2dpinsu)
dfm12<-glm(F0F2dpinsu$Time.to.death~F0F2dpinsu$Generation, quasipoisson(log))
summary(dfm12)

#Proistoma
F0F2dprois<-MqF0F2des[ which(MqF0F2des$Species=='Proisotoma sp.' & MqF0F2des$Acclimation==10), ]
str(F0F2dprois)
dfm13<-glm(F0F2dprois$Time.to.death~F0F2dprois$Generation, quasipoisson(log))
summary(dfm13)
plot(F0F2dprois$Time.to.death~F0F2dprois$Generation)
plot(F0F2dprois$logtime~F0F2dprois$Generation)

Proisomean<-ddply(F0F2dprois, c("Generation"), summarise,
                  N = length(Time.to.death),
                  mean = mean(Time.to.death),
                  sd = sd(Time.to.death),
                  median=median(Time.to.death),
                  max = max(Time.to.death),
                  min = min(Time.to.death))

#################
###F0 analysis###
#################

##########################################################
###F0 data on individuals, no species identity included###
##########################################################
MqF0des10<-MqF0des[ which(MqF0des$Acclimation=='10' & MqF0des$Test=='10'), ]
str(MqF0des10)

Dmodel1<-glm(MqF0des10$Time.to.death~MqF0des10$Status*MqF0des10$logmass, quasipoisson(log))
summary(Dmodel1)
Dmodel2<-glm(MqF0des10$Time.to.death~MqF0des10$Status+MqF0des10$logmass, quasipoisson(log))
summary(Dmodel2)

timeplot<-ggplot(MqF0des, aes(x=logmass, y=logtime, color=Status))+geom_point(shape=19, alpha=1/4, size=5)+
  geom_smooth(method=lm, fullrange=TRUE)+theme_bw()+theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank())
timeplot<-timeplot+ylab("Log time to death")+xlab("Log dry mass") + theme(axis.title.x=element_text(size=20), axis.title.y=element_text(size=20)) + theme(axis.text.x=element_text(size=20), axis.text.y=element_text(size=20))
timeplot
timeplot<-timeplot+theme_few(base_size=14, base_family = "")+scale_colour_colorblind()
timeplot

Timemean<-ddply(MqF0des10, c("Status"), summarise,
                N = length(Time.to.death),
                mean = mean(Time.to.death),
                sd = sd(Time.to.death),
                median = median(Time.to.death),
                min = min(Time.to.death),
                max = max(Time.to.death))

###############################################
############Means analysis#####################
###############################################
###############################################
###File includes F0 (10, 10) ##################
###############################################
###PGLS###
library(caper)
Desmeans.data<-read.csv("Macquarie_Des_Means_10_10_F0.csv",header=TRUE)
Desmeans.tree<-read.tree("DesF0.tre")
str(Desmeans.data)
pg.dmmodel1<-pgls(Logtime~Status*Logmass, comparative.data(Desmeans.tree, Desmeans.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.dmmodel1)

pg.dmmodel1a<-pgls(Logtime~Status+Logmass, comparative.data(Desmeans.tree, Desmeans.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.dmmodel1a)

pg.dmmodel1b<-pgls(Logtime~Status, comparative.data(Desmeans.tree, Desmeans.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.dmmodel1b)

###################################
#####F2 data analysis##############
###################################
MqF2des<-read.csv("Macquarie_Desiccation_F2_V2.csv",header=TRUE)
str(MqF2des)
MqF2des$Acclimation<-as.factor(MqF2des$Acclimation)
MqF2des$Test.Temp<-as.factor(MqF2des$Test.Temp)
str(MqF2des)
MqF2des$logtime<-log10(MqF2des$Time.to.death)
MqF2des$logmass<-log10(MqF2des$Individual.dry.mass)
str(MqF2des)

DestimeF2<-ddply(MqF2des, c("Species", "Acclimation", "Test.Temp"), summarise,
                 N = length(Time.to.death),
                 mean = mean(Time.to.death),
                 sd = sd(Time.to.death),
                 median=median(Time.to.death),
                 max = max(Time.to.death),
                 min = min(Time.to.death))

################################################
###F2 acclimation effects analysis per species##
################################################
#M. caeca
F2desmcaec<-MqF2des[ which(MqF2des$Species=='M. caeca'), ]
str(F2desmcaec)
dac1<-glm(F2desmcaec$Time.to.death~F2desmcaec$Acclimation*F2desmcaec$Test.Temp, quasipoisson(log))
summary(dac1)
dac1p <- ggplot(F2desmcaec, aes(x=Test.Temp, y=logtime, color=Acclimation)) + 
  geom_boxplot()
dac1p
dac1p<-dac1p+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
dac1p+ylab("Log time to death")+xlab("Test temperature (°C)")+theme(legend.positio='none')

#C. denticulata#
F2descdent<-MqF2des[ which(MqF2des$Species=='C. denticulata'), ]
str(F2descdent)
dac2<-glm(F2descdent$Time.to.death~F2descdent$Acclimation*F2descdent$Test.Temp, quasipoisson(log))
summary(dac2)
dac2p <- ggplot(F2descdent, aes(x=Test.Temp, y=logtime, color=Acclimation)) + 
  geom_boxplot()
dac2p
dac2p<-dac2p+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
dac2p+ylab("Log time to death")+xlab("Test temperature (°C)")+theme(legend.positio='none')

#H. purpurescens#
F2deshpurp<-MqF2des[ which(MqF2des$Species=='H. purpurescens'), ]
str(F2deshpurp)
dac3<-glm(F2deshpurp$Time.to.death~F2deshpurp$Acclimation*F2deshpurp$Test.Temp, quasipoisson(log))
summary(dac3)
dac3p <- ggplot(F2deshpurp, aes(x=Test.Temp, y=logtime, color=Acclimation)) + 
  geom_boxplot()
dac3p
dac3p<-dac3p+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
dac3p+ylab("Log time to death")+xlab("Test temperature (°C)")+theme(legend.positio='none')

#P. fimata#
F2despfima<-MqF2des[ which(MqF2des$Species=='P. fimata'), ]
str(F2despfima)
dac4<-glm(F2despfima$Time.to.death~F2despfima$Acclimation*F2despfima$Test.Temp, quasipoisson(log))
summary(dac4)
dac4p <- ggplot(F2despfima, aes(x=Test.Temp, y=logtime, color=Acclimation)) + 
  geom_boxplot()
dac4p
dac4p<-dac4p+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
dac4p+ylab("Log time to death")+xlab("Test temperature (°C)")+theme(legend.positio='none')

#P. insularis
F2despinsu<-MqF2des[ which(MqF2des$Species=='P. insularis'), ]
str(F2despinsu)
dac5<-glm(F2despinsu$Time.to.death~F2despinsu$Acclimation*F2despinsu$Test.Temp, quasipoisson(log))
summary(dac5)
dac5p <- ggplot(F2despinsu, aes(x=Test.Temp, y=logtime, color=Acclimation)) + 
  geom_boxplot()
dac5p
dac5p<-dac5p+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
dac5p+ylab("Log time to death")+xlab("Test temperature (°C)")+theme(legend.positio='none')

#P. notabilis#
F2despnota<-MqF2des[ which(MqF2des$Species=='P. notabilis'), ]
str(F2despnota)
dac6<-glm(F2despnota$Time.to.death~F2despnota$Acclimation*F2despnota$Test.Temp, quasipoisson(log))
summary(dac6)
dac6p <- ggplot(F2despnota, aes(x=Test.Temp, y=logtime, color=Acclimation)) + 
  geom_boxplot()
dac6p
dac6p<-dac6p+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
dac6p+ylab("Log time to death")+xlab("Test temperature (°C)")+theme(legend.positio='none')

#Proisotoma sp.
F2desprois<-MqF2des[ which(MqF2des$Species=='Proisotoma sp.'), ]
str(F2desprois)
dac7<-glm(F2desprois$Time.to.death~F2desprois$Acclimation*F2desprois$Test.Temp, quasipoisson(log))
summary(dac7)
dac7p <- ggplot(F2desprois, aes(x=Test.Temp, y=logtime, color=Acclimation)) + 
  geom_boxplot()
dac7p
dac7p<-dac7p+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
dac7p+ylab("Log time to death")+xlab("Test temperature (°C)")+theme(legend.positio='none')

##########################
###All species together###
##########################
dacall<-glm(MqF2des$Time.to.death~MqF2des$Acclimation*MqF2des$Test.Temp*MqF2des$Status, quasipoisson(log))
summary(dacall)

###############################################
###Rate of development analysis################
###############################################

################################################
##Means for each species at each temperature####
################################################
setwd("C:/RBook/")
Devrate<-read.csv("Macquarie_DevelopmentTime_F2_slc.csv",header=TRUE)
str(Devrate)
Devrate$Temperature<-as.factor(Devrate$Temperature)
Devrate$Days.to.hatching<-as.integer(Devrate$Days.to.hatching)
str(Devrate)

############################
#####Calculate the means####
############################
library(plyr)
library(ggplot2)
Devratemeans<-ddply(Devrate, c("Species", "Temperature"), summarise,
                    N = sum(!is.na(Rate)),
                    mean = mean(Rate, na.rm=TRUE),
                    median = median(Rate, na.rm=TRUE),
                    sd = sd(Rate, na.rm=TRUE),
                    max = max(Rate, na.rm=TRUE),
                    min = min(Rate, na.rm=TRUE))

##################
##Plot the lot####
##################
Devratemeansplot<- ggplot(Devratemeans, aes(x=Temperature, y=mean, group=Species, colour=Species)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=1,
                position=position_dodge(0.1))
Devratemeansplot
Devratemeansplot<-Devratemeansplot+theme_few(base_size=22, base_family = "")+scale_colour_colorblind()
Devratemeansplot+ylab("Log time to death")+xlab("Test temperature (°C)")

########################################################
###Identifying slopes, Topt and Umax for each species###
########################################################
#C.denticulata
DRCdent<-Devratemeans[ which(Devratemeans$Species=='C. denticulata'), ]
DRCdentp<- ggplot(DRCdent, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRCdentp
DRCdentp+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRCdent$Temperature<-as.numeric(levels(DRCdent$Temperature))[DRCdent$Temperature]
###Select the linear range
DRCdent2<-DRCdent[ which(DRCdent$Temperature >=5 & DRCdent$Temperature<=25), ]
##Apply a simple linear model
Dsmodel1<-lm(DRCdent2$mean~DRCdent2$Temperature)
summary(Dsmodel1)
plot(DRCdent2$mean~DRCdent2$Temperature)
abline(Dsmodel1)


#H. purpurescens
DRHpurp<-Devratemeans[ which(Devratemeans$Species=='H. purpurescens'), ]
DRHpurpp<- ggplot(DRHpurp, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRHpurpp
DRHpurpp+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRHpurp$Temperature<-as.numeric(levels(DRHpurp$Temperature))[DRHpurp$Temperature]
###Select the linear range
DRHpurp2<-DRHpurp[ which(DRHpurp$Temperature >=0 & DRHpurp$Temperature<=20), ]
##Apply a simple linear model
Dsmodel2<-lm(DRHpurp2$mean~DRHpurp2$Temperature)
summary(Dsmodel2)
plot(DRHpurp2$mean~DRHpurp2$Temperature)
abline(Dsmodel2)

#H. viatica
DRHviat<-Devratemeans[ which(Devratemeans$Species=='H. viatica'), ]
DRHviatp<- ggplot(DRHviat, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRHviatp
DRHviatp+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRHviat$Temperature<-as.numeric(levels(DRHviat$Temperature))[DRHviat$Temperature]
###Select the linear range
DRHviat2<-DRHviat[ which(DRHviat$Temperature >=10 & DRHviat$Temperature<=20), ]
##Apply a simple linear model
Dsmodel3<-lm(DRHviat2$mean~DRHviat2$Temperature)
summary(Dsmodel3)
plot(DRHviat2$mean~DRHviat2$Temperature)
abline(Dsmodel3)

#M. caeca (I)
DRMcaec<-Devratemeans[ which(Devratemeans$Species=='M. caeca'), ]
DRMcaecp<- ggplot(DRMcaec, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRMcaecp
DRMcaecp+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRMcaec$Temperature<-as.numeric(levels(DRMcaec$Temperature))[DRMcaec$Temperature]
###Select the linear range
DRMcaec2<-DRMcaec[ which(DRMcaec$Temperature >=5 & DRMcaec$Temperature<=20), ]
##Apply a simple linear model
Dsmodel4<-lm(DRMcaec2$mean~DRMcaec2$Temperature)
summary(Dsmodel4)
plot(DRMcaec2$mean~DRMcaec2$Temperature)
abline(Dsmodel4)

#P. fimata
DRPfima<-Devratemeans[ which(Devratemeans$Species=='P. fimata'), ]
DRPfimap<- ggplot(DRPfima, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRPfimap
DRPfimap+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRPfima$Temperature<-as.numeric(levels(DRPfima$Temperature))[DRPfima$Temperature]
###Select the linear range
DRPfima2<-DRPfima[ which(DRPfima$Temperature >=5 & DRPfima$Temperature<=25), ]
##Apply a simple linear model
Dsmodel5<-lm(DRPfima2$mean~DRPfima2$Temperature)
summary(Dsmodel5)
plot(DRPfima2$mean~DRPfima2$Temperature)
abline(Dsmodel5)

#P. insularis (I)
DRPinsu<-Devratemeans[ which(Devratemeans$Species=='P. insularis'), ]
DRPinsup<- ggplot(DRPinsu, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRPinsup
DRPinsup+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRPinsu$Temperature<-as.numeric(levels(DRPinsu$Temperature))[DRPinsu$Temperature]
###Select the linear range
DRPinsu2<-DRPinsu[ which(DRPinsu$Temperature >=5 & DRPinsu$Temperature<=15), ]
##Apply a simple linear model
Dsmodel6<-lm(DRPinsu2$mean~DRPinsu2$Temperature)
summary(Dsmodel6)
plot(DRPinsu2$mean~DRPinsu2$Temperature)
abline(Dsmodel6)

#P. notabilis
DRPnota<-Devratemeans[ which(Devratemeans$Species=='P. notabilis'), ]
DRPnotap<- ggplot(DRPnota, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRPnotap
DRPnotap+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRPnota$Temperature<-as.numeric(levels(DRPnota$Temperature))[DRPnota$Temperature]
###Select the linear range
DRPnota2<-DRPnota[ which(DRPnota$Temperature >=5 & DRPnota$Temperature<=25), ]
##Apply a simple linear model
Dsmodel7<-lm(DRPnota2$mean~DRPnota2$Temperature)
summary(Dsmodel7)
plot(DRPnota2$mean~DRPnota2$Temperature)
abline(Dsmodel7)

#Proisotoma sp.
DRProis<-Devratemeans[ which(Devratemeans$Species=='Proisotoma sp.'), ]
DRProisp<- ggplot(DRProis, aes(x=Temperature, y=mean)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)
DRProisp
DRProisp+ylab("Mean ± s.d. development rate")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))
###Get temperature back to a numeric value
DRProis$Temperature<-as.numeric(levels(DRProis$Temperature))[DRProis$Temperature]
###Select the linear range
DRProis2<-DRProis[ which(DRProis$Temperature >=5 & DRProis$Temperature<=20), ]
##Apply a simple linear model
Dsmodel8<-lm(DRProis2$mean~DRProis2$Temperature)
summary(Dsmodel8)
plot(DRProis2$mean~DRProis2$Temperature)
abline(Dsmodel8)

#################################################
######Electron Volt calculations##################
##################################################
setwd("C:/RBook/")
DeveV<-read.csv("Macquarie_DevelopmentTime_F2_slc_eV.csv",header=TRUE)
str(DeveV)
DeveV$Temperature<-as.factor(DeveV$Temperature)

DeveVmeans<-ddply(DeveV, c("Species", "Temperature"), summarise,
                  N = sum(!is.na(lnrate)),
                  mean = mean(lnrate, na.rm=TRUE),
                  median = median(lnrate, na.rm=TRUE),
                  sd = sd(lnrate, na.rm=TRUE),
                  max = max(lnrate, na.rm=TRUE),
                  min = min(lnrate, na.rm=TRUE))


##Slopes for eV
setwd("C:/RBook/")
Evfinal<-read.csv("RateDevineV.csv",header=TRUE)
str(Evfinal)

###Slopes per species###
#C.denticulata
DevCdent<-Evfinal[ which(Evfinal$Species=='C. denticulata'), ]
Devm1<-lm(DevCdent$lnRate~DevCdent$One_kT)
summary(Devm1)
plot(DevCdent$lnRate~DevCdent$One_kT)
abline(Devm1)

#H. purpurescens
DevHpurp<-Evfinal[ which(Evfinal$Species=='H. purpurescens'), ]
Devm2<-lm(DevHpurp$lnRate~DevHpurp$One_kT)
summary(Devm2)
plot(DevHpurp$lnRate~DevHpurp$One_kT)
abline(Devm2)

#H. viatica
DevHviat<-Evfinal[ which(Evfinal$Species=='H. viatica'), ]
Devm3<-lm(DevHviat$lnRate~DevHviat$One_kT)
summary(Devm3)
plot(DevHviat$lnRate~DevHviat$One_kT)
abline(Devm3)

#M. caeca (I)
DevMcaec<-Evfinal[ which(Evfinal$Species=='M. caeca'), ]
Devm4<-lm(DevMcaec$lnRate~DevMcaec$One_kT)
summary(Devm4)
plot(DevMcaec$lnRate~DevMcaec$One_kT)
abline(Devm4)

#P. fimata
DevPfima<-Evfinal[ which(Evfinal$Species=='P. fimata'), ]
Devm5<-lm(DevPfima$lnRate~DevPfima$One_kT)
summary(Devm5)
plot(DevPfima$lnRate~DevPfima$One_kT)
abline(Devm5)

#P. insularis (I)
DevPinsu<-Evfinal[ which(Evfinal$Species=='P. insularis'), ]
Devm6<-lm(DevPinsu$lnRate~DevPinsu$One_kT)
summary(Devm6)
plot(DevPinsu$lnRate~DevPinsu$One_kT)
abline(Devm6)

#P. notabilis
DevPnota<-Evfinal[ which(Evfinal$Species=='P. notabilis'), ]
Devm7<-lm(DevPnota$lnRate~DevPnota$One_kT)
summary(Devm7)
plot(DevPnota$lnRate~DevPnota$One_kT)
abline(Devm7)

#Proisotoma sp.
DevProis<-Evfinal[ which(Evfinal$Species=='Proisotoma sp.'), ]
Devm8<-lm(DevProis$lnRate~DevProis$One_kT)
summary(Devm8)
plot(DevProis$lnRate~DevProis$One_kT)
abline(Devm8)

###############################################################
###Alien versus indigenous development rate curve comparison###
###############################################################
###PGLS###
library(caper)
Devrate.data<-read.csv("Development_curve_stats.csv",header=TRUE)
Devrate.tree<-read.tree("Devrate.tre")
str(Devrate.data)

pg.devm1<-pgls(Slope~Status, comparative.data(Devrate.tree, Devrate.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.devm1)

pg.devm2<-pgls(ev~Status, comparative.data(Devrate.tree, Devrate.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.devm2)

pg.devm3<-pgls(Topt~Status, comparative.data(Devrate.tree, Devrate.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.devm3)

pg.devm4<-pgls(Umax~Status, comparative.data(Devrate.tree, Devrate.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.devm4)

#####################
###Linear models#####
#####################
devm1<-lm(Slope~Status, data=Devrate.data)
summary(devm1)

devm2<-lm(ev~Status, data=Devrate.data)
summary(devm2)

devm3<-lm(Topt~Status, data=Devrate.data)
summary(devm3)

devm4<-lm(Umax~Status, data=Devrate.data)
summary(devm4)

###Now excluding H. purpurescens###
Devrate2.data<-Devrate.data[ which(Devrate.data$Speciesfull!='H. purpurescens'), ]
str(Devrate2.data)

devm5<-lm(Slope~Status, data=Devrate2.data)
summary(devm5)

devm6<-lm(ev~Status, data=Devrate2.data)
summary(devm6)

devm7<-lm(Topt~Status, data=Devrate2.data)
summary(devm7)

devm8<-lm(Umax~Status, data=Devrate2.data)
summary(devm8)

##################################################
###Hatching success as a function of temperature##
##################################################
setwd("C:/RBook/")
Spsurv<-read.csv("DevSurv_F2_slc.csv",header=TRUE)
str(Spsurv)

Spsurvmeans<-ddply(Spsurv, c("Species", "Temperature"), summarise,
                   N = sum(!is.na(Proportion)),
                   mean = mean(Proportion, na.rm=TRUE),
                   median = median(Proportion, na.rm=TRUE),
                   sd = sd(Proportion, na.rm=TRUE),
                   max = max(Proportion, na.rm=TRUE),
                   min = min(Proportion, na.rm=TRUE))

##################
##Plot the lot####
##################
Spsurvmeansplot<- ggplot(Spsurvmeans, aes(x=Temperature, y=mean, group=Species, colour=Species)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=1,
                position=position_dodge(0.1))
Spsurvmeansplot
Spsurvmeansplot+ylab("Mean ± s.d. hatching success")+xlab("Temperature (°C)")+theme(axis.title.x=element_text(size=22), axis.title.y=element_text(size=22)) + 
  theme(axis.text.x=element_text(size=22), axis.text.y=element_text(size=22))

#####Subset the data by species
#C. denticulata
Spsurvcdent<-Spsurv[ which(Spsurv$Species=='C. denticulata'), ]
str(Spsurvcdent)
plot(Spsurvcdent$Temperature,Spsurvcdent$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#H. purpurescens
Spsurvhpurp<-Spsurv[ which(Spsurv$Species=='H. purpurescens'), ]
str(Spsurvhpurp)
plot(Spsurvhpurp$Temperature,Spsurvhpurp$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#H. viatica
Spsurvhviat<-Spsurv[ which(Spsurv$Species=='H. viatica'), ]
str(Spsurvhviat)
plot(Spsurvhviat$Temperature,Spsurvhviat$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#M. caeca (I)
Spsurvmcaec<-Spsurv[ which(Spsurv$Species=='M. caeca'), ]
str(Spsurvmcaec)
plot(Spsurvmcaec$Temperature,Spsurvmcaec$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#P. fimata
Spsurvpfima<-Spsurv[ which(Spsurv$Species=='P. fimata'), ]
str(Spsurvpfima)
plot(Spsurvpfima$Temperature,Spsurvpfima$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#P. insularis (I)
Spsurvpinsu<-Spsurv[ which(Spsurv$Species=='P. insularis'), ]
str(Spsurvpinsu)
plot(Spsurvpinsu$Temperature,Spsurvpinsu$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#P. notabilis
Spsurvpnota<-Spsurv[ which(Spsurv$Species=='P. notabilis'), ]
str(Spsurvpnota)
plot(Spsurvpnota$Temperature,Spsurvpnota$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#Proisotoma sp.
Spsurvprois<-Spsurv[ which(Spsurv$Species=='Proisotoma sp.'), ]
str(Spsurvprois)
plot(Spsurvprois$Temperature,Spsurvprois$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")


##Clear that LLT50 cannot be assessed because many species fail to hit zero hatching at low temperature

########################################
####Assess ULT50 for hatching success###
########################################
#C. denticulata
Spsurvcdentu<-Spsurvcdent[ which(Spsurvcdent$Temperature > 9), ]
str(Spsurvcdentu)
plot(Spsurvcdentu$Temperature,Spsurvcdentu$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#H. purpurescens
Spsurvhpurpu<-Spsurvhpurp[ which(Spsurvhpurp$Temperature > 9), ]
str(Spsurvhpurpu)
plot(Spsurvhpurpu$Temperature,Spsurvhpurpu$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#H. viatica
Spsurvhviatu<-Spsurvhviat[ which(Spsurvhviat$Temperature > 9), ]
str(Spsurvhviatu)
plot(Spsurvhviatu$Temperature,Spsurvhviatu$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#M. caeca (I)
Spsurvmcaecu<-Spsurvmcaec[ which(Spsurvmcaec$Temperature > 9), ]
str(Spsurvmcaecu)
plot(Spsurvmcaecu$Temperature,Spsurvmcaecu$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#P. fimata
Spsurvpfimau<-Spsurvpfima[ which(Spsurvpfima$Temperature > 9), ]
str(Spsurvpfimau)
plot(Spsurvpfimau$Temperature,Spsurvpfimau$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#P. insularis (I)
Spsurvpinsuu<-Spsurvpinsu[ which(Spsurvpinsu$Temperature > 9), ]
str(Spsurvpinsuu)
plot(Spsurvpinsuu$Temperature,Spsurvpinsuu$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#P. notabilis
Spsurvpnotau<-Spsurvpnota[ which(Spsurvpnota$Temperature > 9), ]
str(Spsurvpnotau)
plot(Spsurvpnotau$Temperature,Spsurvpnotau$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

#Proisotoma sp.
Spsurvproisu<-Spsurvprois[ which(Spsurvprois$Temperature > 9), ]
str(Spsurvproisu)
plot(Spsurvproisu$Temperature,Spsurvproisu$PercHatch,xlab="Temperature (°C)", ylab = "Hatching success (%)")

###Upper survival analyses
##Crawley's dose response method for estimating mean and s.e. of ULT50##
library(MASS)

#C.denticulata
survm1<-glm(cbind(Spsurvcdentu$Alive,Spsurvcdentu$Dead)~Spsurvcdentu$Temperature,binomial)
summary(survm1)
anova(survm1,test="Chi")
dose.p(survm1,p=c(0.5,0.9,0.95))

#H. purpurescens
survm2<-glm(cbind(Spsurvhpurpu$Alive,Spsurvhpurpu$Dead)~Spsurvhpurpu$Temperature,binomial)
summary(survm2)
anova(survm2,test="Chi")
dose.p(survm2,p=c(0.5,0.9,0.95))

#H. viatica
survm3<-glm(cbind(Spsurvhviatu$Alive,Spsurvhviatu$Dead)~Spsurvhviatu$Temperature,binomial)
summary(survm3)
anova(survm3,test="Chi")
dose.p(survm3,p=c(0.5,0.9,0.95))

#M. caeca (I)
survm4<-glm(cbind(Spsurvmcaecu$Alive,Spsurvmcaecu$Dead)~Spsurvmcaecu$Temperature,binomial)
summary(survm4)
anova(survm4,test="Chi")
dose.p(survm4,p=c(0.5,0.9,0.95))

#P. fimata
survm5<-glm(cbind(Spsurvpfimau$Alive,Spsurvpfimau$Dead)~Spsurvpfimau$Temperature,binomial)
summary(survm5)
anova(survm5,test="Chi")
dose.p(survm5,p=c(0.5,0.9,0.95))

#P. insularis (I)
survm6<-glm(cbind(Spsurvpinsuu$Alive,Spsurvpinsuu$Dead)~Spsurvpinsuu$Temperature,binomial)
summary(survm6)
anova(survm6,test="Chi")
dose.p(survm6,p=c(0.5,0.9,0.95))

#P. notabilis
survm7<-glm(cbind(Spsurvpnotau$Alive,Spsurvpnotau$Dead)~Spsurvpnotau$Temperature,binomial)
summary(survm7)
anova(survm7,test="Chi")
dose.p(survm7,p=c(0.5,0.9,0.95))

#Proisotoma sp.
survm8<-glm(cbind(Spsurvproisu$Alive,Spsurvproisu$Dead)~Spsurvproisu$Temperature,binomial)
summary(survm8)
anova(survm8,test="Chi")
dose.p(survm8,p=c(0.5,0.9,0.95))

##################################
###PGLS comparison among groups###
##################################
###PGLS###
library(caper)
Devrate.data<-read.csv("Development_curve_stats.csv",header=TRUE)
Devrate.tree<-read.tree("Devrate.tre")
str(Devrate.data)

pg.devm5<-pgls(HSULT50~Status, comparative.data(Devrate.tree, Devrate.data, Species, vcv=TRUE, vcv.dim=3), lambda="ML")
summary(pg.devm5)

####################################
###Here concludes the analyses######
####################################