#This code is for figure 01 (at time of submission) in "Rapid evolution of thermal tolerance and phenotypic plasticity in variable environments " by Schaum et al.  

#track keeping selexperiment 
library(ggplot2)
library(plyr)
library(plotrix)
library(nlme)
library(lme4)
library(cowplot)
library(MuMIn)
rm(list=ls())

#setwd("/Users/elisaschaum/Dropbox/Post doc ESI/Thally actual selection experiment/transfers and keeping track")
 setwd("~/Dropbox/Thally long fluc paper/Thally actual selection experiment/transfers and keeping track")

tr<-read.csv("track_keeper.csv") 
tr$gr<-(tr$mue*7)+log(100)
#subset for just FL (long flucs) 
subFL<-subset(tr,generations.passed<130 & nutrients=="FULL" & treatment=="FL"& fl.at==32)

#plot all BUT FL 
trp<-qplot(as.factor(week), mue, ylab=expression(Growth~rate~µ~day^{-1}),xlab="Week of experiment",data=subset(tr, nutrients=="FULL"&generations.passed<=300 &treatment!="FL"),facets=.~treatment, fill=as.factor(treatment), geom="boxplot")+scale_fill_manual(values=c('darkgreen', 'blue','darkred','purple','black')) +theme_classic(base_size = 14) +geom_smooth(aes(group=treatment), colour='black',method='lm', formula=y~poly(x,20)) + scale_x_discrete(name="Week of experiment", breaks=seq(0, 71,10)) +facet_wrap(~treatment,ncol=1) +theme(legend.position='none')
trp

trpnew<-qplot(as.factor(round(generations.passed,-1)), mue, ylab="Growth rate µ",xlab="Week of experiment",data=subset(tr, nutrients=="FULL"&generations.passed<=300 &treatment!="FL"), fill=as.factor(treatment), geom="boxplot")+scale_fill_manual(values=c('darkgreen', 'blue','purple','darkred','black')) +theme_classic(base_size = 14) +geom_smooth(aes(group=treatment), colour='black',method='lm', formula=y~poly(x,15)) + scale_x_discrete(name="Generations", breaks=seq(0, 300,20))  +theme(legend.position='none')
trpnew

trp2<-qplot(as.factor(week), gr,data=subset(tr, nutrients=="FULL"&generations.passed<=300 &treatment!="FL"),facets=.~treatment, fill=as.factor(treatment), geom="point")+scale_fill_manual(values=c('darkgreen', 'blue','darkred','purple','black')) +theme_classic(base_size = 14) +geom_smooth(aes(group=treatment), colour='black',method='lm', formula=y~poly(x,20)) + scale_x_discrete(name="Week of experiment", breaks=seq(0, 71,10)) +facet_wrap(~treatment,ncol=1) +theme(legend.position='none')
trp2

tr<- within(tr, tr_br<- as.character(factor(treatment):factor(biorep)))


trp3<-qplot(as.factor(week), mue,ylab=expression(Growtrate~µ~(d^{-1})),data=subset(tr, nutrients=="FULL"&generations.passed<=300 &treatment!="FL"),facets=.~treatment, fill=as.factor(treatment), geom="point", colour=as.factor(treatment))+scale_fill_manual(values=c('darkgreen', 'blue','darkred','purple','black')) +geom_smooth(aes(group=tr_br), colour='black',method='lm', formula=y~poly(x,10))+theme_classic(base_size = 14) + scale_x_discrete(name="Week of experiment", breaks=seq(0, 71,10)) +facet_wrap(~treatment,ncol=2) +theme(legend.position='top')+scale_colour_manual(values=c('darkgreen', 'blue','darkred','purple','black'))
trp3

#make geometric mean of all 
gm_mean <- function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

# GAMM output median log pop size 
trgm<-ddply(subset(tr1.2, generations.passed<300 & which=="end"),c("treatment","biorep") , function(df) return(c(wmean=gm_mean((df$gr)), grsd=sd(df$gr), var=var(df$gr), mmean=median(df$gr)))) 
mmp<-qplot(as.factor(treatment), mmean, ylab=expression(GAMM~output~median~log~population~size~(ln~(cells~ml^{-1}))),xlab="Selection environment",data=trgm,fill=as.factor(treatment), geom="boxplot")+scale_fill_manual(values=c('darkgreen', 'blue','darkred','purple','black')) +theme_classic(base_size = 14) +theme(legend.position='none')
mmp

gmp<-qplot(as.factor(treatment), (wmean), ylab=expression(Geometri~mean~population~size~(ln~(cells~ml^{-1}))),xlab="Selection environment",data=trgm,fill=as.factor(treatment), geom="boxplot")+scale_fill_manual(values=c('darkgreen', 'blue','darkred','purple','black')) +theme_classic(base_size = 14) +theme(legend.position='none')
gmp


 
#make averages just stable 

trav<-ddply(subset(tr, nutrients=="FULL"&generations.passed<=300 &treatment!="FL"&treatment!="FS"),c("treatment","week") , function(df) return(c(r.avg=mean(df$mue),r.sd=sd(df$mue))))
trpnewav<-qplot(as.factor(week), r.avg, ylab="Growth rate µ",xlab="Week of experiment",data=trav, colour=as.factor(treatment), fill=as.factor(treatment))+scale_colour_manual(values=c('darkgreen', 'blue','purple')) +theme_classic(base_size = 14)+scale_fill_manual(values=c('darkgreen', 'blue','purple')) 
trpnewav +geom_smooth(aes(group=treatment), colour='black',method='lm', formula=y~poly(x,20)) + scale_x_discrete(name="Week of experiment", breaks=seq(0, 71,10)) +geom_errorbar(aes(ymin=r.avg-r.sd/sqrt(6), ymax=r.avg+r.sd/sqrt(6))) +theme(legend.position='top')

last_plot() + theme(panel.grid.major = element_blank()) 
last_plot() + theme(panel.grid.minor = element_blank())


#also read all the correlated responses (i.e reciprocal transplants) - was in the end delegated to SI rather than main manuscript figure 
corr<-read.csv("corr_Response growth .csv") #read.csv("t100t300 reciprocal transplants.csv") #may need to delete space? 
head(corr)
#make selected assay column
corr <- within(corr, sel_assay  <- as.character(factor(selected):factor(assay)))

qplot(sel_assay, rel.to.22, ylab="Growth rate µ relative to growth of control at 22ºC",data=corr,facets=gen~., fill=as.factor(selected), geom="boxplot")+scale_fill_manual(values=c('darkgreen', 'blue','purple','darkred','black')) +theme_cowplot()+
  geom_segment(mapping=aes(x =1, y = 1.0, xend = 20, yend = 1.0), linetype=2)

#just mue, and with bioreplicates
qplot(selected, mue, ylab="Growth rate µ",data=subset(corr,gen=="t300"),facets=assay~., colour=as.factor(selected), geom="boxplot")+scale_colour_manual(values=c('darkgreen', 'blue','red','black','purple')) +theme_classic()+facet_wrap(~assay)+
 geom_jitter(aes(shape = factor(biorep)))+scale_shape_manual(values=c(1:7))
#+  geom_segment(mapping=aes(x =1, y = .645, xend = 20, yend = .645), linetype=2)

subcorr01<-subset(corr, gen=="t300")
MM4<-lme(mue~selected*assay, random=~1|selected/biorep, method="ML", data=subcorr01)
dd_MM4<-dredge(MM4)
dd_MM4

MM4_fin<-lme(mue~selected*assay, random=~1|selected/biorep, method="REML", data=subcorr01)
summary(MM4_fin)
library(lsmeans)
ll<-as.matrix(lsmeans(MM4_fin, pairwise ~ selected*assay, adjust = "tukey") )
#just mue, and just for FL at 32, 26, 22
qplot(sel_assay, mue, ylab="Growth rate µ",data=subset(corr,gen=="t300" & selected=="FL" & assay!="FS"), colour=as.factor(biorep), shape=biorep, size=0.2) +theme_classic() +scale_shape_manual(values=c(1:7))


#now per generation - we also need a growth rate compared to the ancestor or 22 at 22 plot 

trp2<-qplot(generations.passed, mue, ylab="growth rate µ",data=subset(tr,nutrients=="FULL"),facets=treatment~., colour=as.factor(actual.temp), shape=as.factor(biorep))
trp2+theme_bw()+scale_colour_manual(values = c('darkblue', 'orange', 'red', 'pink'))

trp2.d<-qplot(day.of.exp, mue, ylab="growth rate µ",data=subset(tr,nutrients=="FULL"),facets=treatment~., colour=as.factor(treatment), shape=as.factor(biorep))
trp2.d+theme_bw()


#per gen per biorep with smoother in temp facet - level ought to be 0.95 but obscures plot, so reduced to 0.4 

trp2.ds<-qplot(generations.passed, mue, ylab="growth rate µ",xlab="generations",data=subset(tr, nutrients=="FULL"&generations.passed<=320),facets=treatment~., colour=as.factor(treatment), fill=as.factor(treatment))+scale_colour_manual(values=c('darkgreen', 'blue','purple','darkred','black'))+scale_fill_manual(values=c('darkgreen', 'blue','purple','darkred','black')) +theme_cowplot()

trp2.ds +geom_smooth(aes(group=biorep),size=0.2, method='gam', formula=y~s(x, k=28), alpha=0.25, level=0.4) + scale_x_continuous(name="Generations", breaks=seq(0, 320,20))
#per day of experiment 


trp3.ds<-qplot(day.of.exp, mue, ylab="growth rate µ",data=subset(tr, nutrients=="FULL"&generations.passed<=320),facets=treatment~., colour=as.factor(treatment), fill=as.factor(treatment))+scale_colour_manual(values=c('darkgreen', 'blue','purple','darkred','black'))+scale_fill_manual(values=c('darkgreen', 'blue','purple','darkred','black')) +theme_cowplot()
trp3.ds +geom_smooth(aes(group=biorep),size=0.2, method='gam', formula=y~s(x, k=30), alpha=0.25, level=0.4) + scale_x_continuous(name="Day of Experiment", breaks=seq(0, 550,50))
#now just FL we know 1,2,5 are the one set and 3,4,6, the other set. we 
trp4.ds<-qplot(day.of.exp, mue, ylab="growth rate µ",
               data=subset(tr, nutrients=="FULL"&generations.passed<=320&treatment=="FL"),
                colour=as.factor(fl.at), 
               fill=as.factor(treatment), shape=as.factor(biorep))+geom_point(aes(size=0.02))+
  scale_colour_manual(values=c('black', 'darkred'))+scale_shape_manual(values=c(1,1,2,2,1,2))+
 scale_fill_manual(values=c('darkred', 'darkred')) +theme_cowplot()

trp4.ds +geom_smooth(aes(group=biorep),size=0.02, method='gam', formula=y~s(x, k=30), 
                     alpha=0.25, level=0.4) + scale_x_continuous(name="Day of Experiment", 
                                                                 breaks=seq(0, 550,50))





trp2.s<-qplot(generations.passed, mue, ylab="growth rate µ",data=subset(tr, nutrients=="FULL"), facets=biorep~., colour=treatment, shape=as.factor(biorep))
trp2.s+theme_bw()#+geom_line(aes(group=treatment), se=FALSE, fullrange=TRUE)

#only FL, coloured by selection environment temperature

trp2.2<-qplot(generations.passed, mue, ylab="growth rate µ",data=subset(tr, nutrients=="FULL" & treatment=="FL"), colour=as.factor(fl.at), facets=biorep~.)
trp2.2+theme_bw()+scale_colour_manual(values=c("blue","red"))

trp2.new<-qplot(as.factor(week), mue, ylab="growth rate µ",data=subset(tr, nutrients=="FULL" & treatment=="FL"), colour=as.factor(fl.at), facets=biorep~.)
trp2.new+theme_bw() +scale_colour_manual(values=c("black","purple"))#look at >week 33! - divergence??? biorep 1,2,5 much happier in 22C

# now with geometric mean, also 

gm_mean <- function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}


# only FL, coloured by growth rate value 
# add a means per generation bin (weekly) to the column
sdw<-ddply(tr,c("week","treatment","biorep") , function(df) return(c(wmean=mean((df$mue)), Fsd=sd(df$mue), var=var(df$mue)))) # these are  for high chlorophyll data
tr<-merge(sdw, tr)
a<-"black"
b<-"red"
trp2.3<-qplot(as.factor(week), mue, ylab="growth rate µ",data=subset(tr, nutrients=="FULL" & treatment=="FL"), colour=ifelse(mue>wmean+Fsd,a,b), geom="point", shape=as.factor(fl.at), size=0.4, facets=biorep~.)
trp2.3+theme_bw()+scale_colour_manual(values=c("red","black"))+scale_shape_manual(values=c(1,17))
trp2.4<-qplot(as.factor(week), mue, ylab="growth rate µ",data=subset(tr, nutrients=="FULL" & treatment=="FL" & week>30), colour=as.factor(biorep), geom="point", shape=as.factor(fl.at), size=0.4)
trp2.4+theme_bw()+scale_colour_manual(values=c("red","black","purple","grey","blue","darkgrey"))+scale_shape_manual(values=c(1,17))

#or
trp2.5<-qplot(as.factor(week), mue, ylab="growth rate µ",data=subset(tr, nutrients=="FULL" & treatment=="FL" & week>31), colour=as.factor(biorep), geom="point", shape=as.factor(actual.temp), size=0.4)


#just for 22C

m22<-gm_mean(subset(tr,treatment=="22" & nutrients=="FULL"&generations.passed>250&generations.passed<350)$mue) # for each one indivually
manc22<-gm_mean(subset(tr,treatment=="22" & nutrients=="FULL"&week<4)$mue) # for each one indivually
m26<-gm_mean(subset(tr,treatment=="26" & nutrients=="FULL"&generations.passed>300&generations.passed<350)$mue) # for each one indivually
manc26<-gm_mean(subset(tr,treatment=="26" & nutrients=="FULL"&week<4)$mue) # for each one indivually
m32<-gm_mean(subset(tr,treatment=="32" & nutrients=="FULL"&generations.passed>250&generations.passed<350)$mue) # for each one indivually
manc32<-gm_mean(subset(tr,treatment=="32" & nutrients=="FULL"&week<4)$mue) # for each one indivually

mFS<-gm_mean(subset(tr,treatment=="FS" & nutrients=="FULL"&generations.passed>250&generations.passed<350)$mue) # for each one indivually
mancFS<-gm_mean(subset(tr,treatment=="FS" & nutrients=="FULL"&week<4)$mue) # for each one indivually


m22/manc22
m26/manc26
m32/manc32
mFS/mancFS

ma22<-gm_mean(subset(tr,treatment=="22" & nutrients=="FULL"&week>4)$mue) # for each one indivually
sd22<-sd(subset(tr,treatment=="22" & nutrients=="FULL"&week>4)$mue)/sqrt(6) # for each one indivually
ma26<-gm_mean(subset(tr,treatment=="26" & nutrients=="FULL"&week>4)$mue) # for each one indivually
sd26<-sd(subset(tr,treatment=="26" & nutrients=="FULL"&week>4)$mue)/sqrt(6) # for each one indivually
ma32<-gm_mean(subset(tr,treatment=="32" & nutrients=="FULL"&generations.passed<100)$mue) # for each one indivually -change bacl to week= 4 if needed
sd32<-sd(subset(tr,treatment=="32" & nutrients=="FULL"&week>4)$mue)/sqrt(6) # for each one indivually
maFS<-gm_mean(subset(tr,treatment=="FS" & nutrients=="FULL"&week>4)$mue) # for each one indivually

sFS<-sd(subset(tr,treatment=="FS" & nutrients=="FULL"&week>4)$mue)/sqrt(6) # for each one indivually
sFS<-sd(subset(tr,treatment=="FS" & nutrients=="FULL"&week>4)$mue)/sqrt(6) # for each one indivually

#FL! not necesasary gm, can just use mean
tr$fl.at<-as.numeric(tr$fl.at)
maFL22pref22<-gm_mean(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200  & fl.at== 22 & biorep=="1"|biorep=="2"|biorep=="5" )$mue)
maFL22pref22
sdFL22pref22<-sd(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200 & fl.at== 22 & biorep=="1"|biorep=="2"|biorep=="5" )$mue)
sdFL22pref22

maFL32pref22<-gm_mean(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200& fl.at== 32 & biorep=="1"|biorep=="2"|biorep=="5" )$mue)
maFL32pref22
sdFL32pref22<-sd(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200 & fl.at== 32 & biorep=="1"|biorep=="2"|biorep=="5" )$mue)
sdFL32pref22


maFL22pref32<-gm_mean(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200  & fl.at== 22 & biorep=="3"|biorep=="4"|biorep=="6" )$mue)
maFL22pref32
sdFL22pref32<-sd(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200 & fl.at== 22 & biorep=="3"|biorep=="4"|biorep=="6" )$mue)
sdFL22pref32

maFL32pref32<-gm_mean(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200& fl.at== 32 & biorep=="3"|biorep=="4"|biorep=="6" )$mue)
maFL32pref32
sdFL32pref32<-sd(subset(tr,treatment=="FL" & nutrients=="FULL"&generations.passed>200 & fl.at== 32 & biorep=="3"|biorep=="4"|biorep=="6" )$mue)
sdFL32pref32


#now -- AFTER plastic response,  - 
subgeom<-subset(tr, nutrients=="FULL" & week >4 &treatment=="FL") #delete treatment ==Fl to look at others also, to make sure the difference isn't just an artefact of time (done so, difference only prevalent in FL)
subgeom_unique<-as.data.frame(unique(subgeom$treatment))

geom_means <- ddply(subgeom, .(treatment, fl.at), summarise, geom_mean = gm_mean(mue), stderror=std.error(mue))
geom_means$what<-(as.factor(geom_means$treatment):as.factor(geom_means$fl.at))
geo.mean.plot<-qplot(what,geom_mean , ylab="geometrical mean",data=geom_means, colour=treatment, size=0.3)

geo.mean.plot+theme_classic()+geom_errorbar(aes(ymin=geom_mean-stderror, ymax=geom_mean+stderror), size=0.1, width=0.1)



#now against plastic responses etc
#against growth at t0 (at 22C)
Str<-subset(tr, treatment==22)
trp4<-qplot(as.factor(week),muetot0 , ylab="µ compared to control at t0",data=subset(tr, nutrients=='FULL'), geom="boxplot", fill=treatment, facets=.~treatment)
trp4+theme_bw()+geom_segment(mapping=aes(x =1, y = 1.0, xend = 60, yend = 1.0), linetype=2)+geom_segment(mapping=aes(x =1, y = max(Str$muetot0), xend = 62, yend = max(Str$muetot0)), linetype=2)
#against growth at ca t20 - t30 (t in days not generations) at selection temperature (i.e. compared to acclimation in selection environment) - what the...
trp5<-qplot(as.factor(week),muetoplast0 , ylab="µ compared to growth in selection environment after  ca 20 generations",data=subset(tr, nutrients=='FULL'), geom="boxplot", fill=treatment, facets=treatment~.)

trp5+theme_bw()+geom_segment(mapping=aes(x =1, y = 1.0, xend = 69, yend = 1.0), linetype=2)+geom_segment(mapping=aes(x =1, y = max(Str$muetoplast0), xend = 69, yend = max(Str$muetoplast0)), linetype=2)

#this is one of the good ones 
trp6<-qplot(as.factor(treatment), muetoplast0, data=subset(tr, generations.passed >=300 & generations.passed<=320), fill=treatment,geom='boxplot') # change to week>4 for pooled data, or use generations passed 
trp6+theme_bw()+geom_segment(mapping=aes(x =0.2, y = 1.0, xend = 6, yend = 1.0), colour='black',linetype=2)+scale_fill_manual(values=c('darkgreen', 'blue','red','black','purple')) +theme_cowplot()
#Now the same again but without 26C and with the colour 

#so we use this for our wee mixed model 
tr$treatment<-as.factor(tr$treatment)
trmod<-subset(tr, treatment!="26"&generations.passed >=300 & generations.passed<=320)
MM1 <- lme(muetoplast0~treatment, random = ~  1|treatment/biorep, method="ML", data=trmod) 
predict(MM1, type = "response")
plot(MM1) #is well behaved
hist(MM1$residuals) # yay! 

#now dredge this 
dd_all<- dredge(MM1, rank = 'AICc')
dd_all # best model finds significant differences
#now refit as REML 
MM1_fin <- lme(muetoplast0~treatment, random = ~  1|treatment/biorep, method="REML", data=trmod)
#make confindence intervals (or use st.errors)

confint.mm <- data.frame(intervals(MM1_fin, which='fixed')[1])
#for the rest, we actually want a loop because it does get lengthy otherwise 
for(i in 1:nrow(confint.mm)){
  confint.mm$stdev[i] <- confint.mm$fixed.upper[i]-confint.mm$fixed.est[i]
}

## now also split the groups 


### now split the FL in two for peak rate and Topt, too #####
trmod$treatment<-as.character(trmod$treatment)
trmod$biorep<-as.character(trmod$biorep)
trmod$Treatment2<-ifelse(trmod$treatment=="FL"&trmod$biorep=="1"|trmod$treatment=="FL"&trmod$biorep=="2"|trmod$treatment=="FL"&trmod$biorep=="5","FL2",
                         ifelse(trmod$treatment=="FL"&trmod$biorep=="3"|trmod$treatment=="FL"&trmod$biorep=="4"|trmod$treatment=="FL"&trmod$biorep=="6","FL1",as.character(trmod$treatment)))
trmod$treatment<-as.factor(trmod$treatment)
trmod$biorep<-as.factor(trmod$biorep)
trmod$Treatment2<-as.factor(trmod$Treatment2)


Mod0 <- lme(fixed = muetoplast0 ~ 1 +  Treatment2, random = ~ 1|biorep/Treatment2, method = 'ML', trmod)

S1 <- lme(fixed = muetoplast0 ~ 1 +  Treatment2, random = ~ 1|biorep/Treatment2, weights = varIdent( ~ 1|Treatment2), method = 'ML', trmod)

S2 <- lme(fixed = muetoplast0 ~ 1 , random = ~ 1|biorep/Treatment2, weights = varIdent( ~ 1|Treatment2), method = 'ML', trmod)
anova(S1,S2) # growthT is siginficant effect

Sfinal <- lme(fixed = muetoplast0 ~ 1 +  Treatment2, random = ~ 1|biorep/Treatment2, weights = varIdent( ~ 1|Treatment2), method = 'REML', trmod)
confint.mm <- data.frame(intervals(Sfinal, which='fixed')[1])
confint.mm


#now we look at the variance of specific environments at that point in time
data2<-subset(tr, generations.passed >=300 & generations.passed<=320 & treatment=="FS")
var(data2$muetoplast0)
data3<-subset(tr, generations.passed >=300 & generations.passed<=320 & treatment=="FL")
var(data3$muetoplast0)
#as is this 
trp7<-qplot(as.factor(treatment), muetoplast0, ylab="evolutionary response of fluctuating samples to warming",data=subset(tr, week>4), geom='boxplot', facets=.~fl.at)
trp7+theme_bw()+geom_segment(mapping=aes(x =0.2, y = 1.0, xend = 6, yend = 1.0), colour='black',linetype=2)#+facet_wrap(~week)

trp8<-qplot(as.factor(treatment), muetot0, data=subset(tr, week>4), geom='boxplot', facets=.~fl.at)
trp8+theme_bw()+geom_segment(mapping=aes(x =0.2, y = 1.0, xend = 6, yend = 1.0), colour='black',linetype=2)#+facet_wrap(~week)


# we need a boxplot of geommeans and then another of fold change compared to22c at the end 
# or is this the better geom mean function
gm_mean <- function(x, na.rm=TRUE, zero.propagate = FALSE){
  if(any(x < 0, na.rm = TRUE)){
    return(NaN)
  }
  if(zero.propagate){
    if(any(x == 0, na.rm = TRUE)){
      return(0)
    }
    exp(mean(log(x), na.rm = na.rm))
  } else {
    exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
  }
}

