## This code can be used for getting data ready to be plotted to recreate Figure 02 (at time of submission) in "Rapid evolution of thermal tolerance and phenotypic plasticity in variable environments" by Schaum et al 
## The code is not the cleanest, neither is it spectacularly elegant. Use at your own risk. For questions please contact elisa.schaum@uni-hamburg.de 

## the code looks in detail at data from t0, t100, and t300. The first half treats the FL (long fluctuations) populations as one. The second half splits them up according to their evolved temperature preference. Not all of this made it into the final version of the manuscript. 

# read libraries
library(data.table)
library(plyr)
library(zoo)
library(lattice)
library(ggplot2)
library(tools)
require(nlme)
require(minpack.lm)
library(reshape2)
#install.packages("devtools")
library(devtools)

#install_github('padpadpadpad/TeamPhytoplankton')
#install_github('padpadpadpad/nlsTools')
library(TeamPhytoplankton)
library(nlsLoop)
library(AICcmodavg)
library(cowplot)
library(MuMIn)
rm(list=ls())

#### read in and combine data ######

#read in logistic growth curve data t100
setwd("~/Desktop/Post doc ESI/Thally actual selection experiment/acclimation and translocation")
res_logis<-read.csv('20161120_res_logis_t000.csv') #t000
res_logis$generation<-rep("t000", times=length(res_logis$pa))
setwd("~/Desktop/Post doc ESI/Thally actual selection experiment/t100 growth assays updated")
#res_logis2<-read.csv("20161123_resloglint100.csv") 
#res_logis2$generation<-rep("t100", times=length(res_logis2$pa))

setwd("~/Desktop/Post doc ESI/Thally actual selection experiment/t300 growth assays 2")
res_logis3<-read.csv("20161123_resloglint300.csv")
res_logis3$generation<-rep("t300", times=length(res_logis3$pa))

head(res_logis)
head(res_logis2)
head(res_logis3)

res_logis<-res_logis[,c(2,3,7:10)]
res_logis2<-res_logis2[,c(2:5,7,8)]
res_logis3<-res_logis3[,c(2:5,7,8)]

res_logis_all<-rbind(res_logis,res_logis2,res_logis3)  
head(res_logis_all)
#drop t100
res_logis_all_st<-subset(res_logis_all, generation!="t100")
#make log column
res_logis_all_st$logr<-log(res_logis_all_st$r) # makes some NANs, we can get rid of them later 
str(res_logis_all_st)

res_logis_all_st$assayT<-as.numeric(res_logis_all_st$assayT)
res_logis_all_st$K<-res_logis_all_st$assayT+273.15

head(res_logis_all_st)


##### get rid of no growth samples #####
# mainly ok
#make an ID column with time as a factor, too
res_logis_all_st <- within(res_logis_all_st, id.gen <- as.character(factor(biorep):factor(assayT):factor(selected):factor(generation)))# create a pond:assayT identifier 

drop_logis <-c('b2:25:22:t100','b3:25:26:t100','b2:20:FS:t100','b6:20:FS:t100','b6:35:22:t300','b1:15:22:t000','b2:15:22:t000') #also drop the ones that -Inf or Nan
res_logis_clean<- res_logis_all_st[! res_logis_all_st$id.gen %in% drop_logis, ]
head(res_logis_clean)
#loop for individual bioreps -- need for stats 
#remove NAs
res_logis_clean<-subset(res_logis_clean[complete.cases(res_logis_clean[,7]),])
res_logis_clean<-within(res_logis_clean, id.plot<-as.factor(biorep):as.factor(selected):as.factor(generation))

###### run as mixed effects model ####
##need to make new id column that has gen, sel, biorep, but not assayT

resl_grp <- groupedData(logr ~ K | id.plot,subset( res_logis_clean, generation=="t300"))
resl_grp$growthT <- as.factor(resl_grp$selected)
resl_grp$generation <- as.factor(resl_grp$gen) 

## full model  - FL won't run?? if included with only FS and 26 -- so do full model and drop the ones that are not necessary later for plotting 
resl.mix <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                 data = resl_grp,
                 fixed = list(ln.c + Ea + Eh + Th ~ 1 + growthT),
                 random = pdDiag(ln.c + Ea + Eh +Th ~ 1),
                 start = c(-.9, rep(0,4), 1, rep(0,4), 8, rep(0,4), 306, rep(0,4)), # starting parameters from nlsLoop
                 na.action = na.omit,
                 method = 'ML',
                 control = nlmeControl(pnlsTol = 0.5))
summary(resl.mix)



## simplify treatment effects - singles
# remove treatment on effect on Ea 
resl.mix1<- update(resl.mix, fixed = list(Ea ~ 1, ln.c + Eh + Th ~ 1 + growthT), start = c(-0.8, rep(0,4), 1, 4, rep(0,4), 305, rep(0,4)))
# remove treatment on effect on Eh 
resl.mix2<- update(resl.mix, fixed = list(Eh ~ 1, ln.c + Ea + Th ~ 1 + growthT), start = c(-0.8, rep(0,4), 1, rep(0,4), 4, 305, rep(0,4))) 
# remove treatment on effect on Th 
resl.mix3<- update(resl.mix, fixed = list(ln.c + Ea + Eh ~ 1 + growthT, Th ~ 1), start = c(-0.8, rep(0,4), 1, rep(0,4), 4, rep(0,4), 305))
# remove treatment on effect on ln.c 
resl.mix4<- update(resl.mix, fixed = list(ln.c ~ 1, Ea + Eh + Th ~ 1 + growthT), start = c(-0.8, 1, rep(0,4), 4, rep(0,4), 305, rep(0,4)))

## simplify treatment effects - pairs
# remove treatment on effect on Ea + Eh no 
resl.mix5<- update(resl.mix, fixed = list( Ea + Eh ~ 1,ln.c + Th ~ 1 + growthT), start = c(-0.8, rep(0,4), 1, 4, 305, rep(0,4)))
# remove treatment on effect on Ea + Th
resl.mix6<- update(resl.mix, fixed = list(ln.c + Eh ~ 1 + growthT, Ea + Th ~ 1), start = c(-0.8, rep(0,4), 1, 4, rep(0,4), 305))
# remove treatment on effect on Eh + Th
resl.mix7<- update(resl.mix, fixed = list(ln.c + Ea ~ 1 + growthT, Eh + Th ~ 1), start = c(-0.8, rep(0,4), 1, rep(0,4), 4, 305))
# remove treatment on effect on ln.c + Ea

resl.mix8<- update(resl.mix, fixed = list(ln.c + Ea ~ 1, Eh + Th ~ 1 + growthT), start = c(-0.8, 1, 4, rep(0,4), 305, rep(0,4)))
# remove treatment on effect on ln.c + Eh no
resl.mix9<- update(resl.mix, fixed = list(ln.c + Eh ~ 1, Ea + Th ~ 1 + growthT), start = c(-0.8, 1, rep(0,4), 4, 305, rep(0,4)))
# remove treatment on effect on ln.c + Th

resl.mix10<- update(resl.mix, fixed = list(Ea + Eh ~ 1 + growthT, ln.c + Th ~ 1), start = c(-0.8, 1, rep(0,4), 4, rep(0,4), 305)) 

## simplify treatment effects - Three
# remove treatment on effect on Ea + Eh + Th
resl.mix11<- update(resl.mix, fixed = list(ln.c ~ 1 + growthT, Ea + Eh + Th ~ 1), start = c(-0.8, rep(0,4), 1, 4, 305))
# remove treatment on effect on ln.c + Ea + Eh
resl.mix12<- update(resl.mix, fixed = list(ln.c + Ea + Eh ~ 1, Th ~ 1 + growthT), start = c(-0.8, 1, 4, 305, rep(0,4)))
# remove treatment on effect on ln.c + Ea + Th
resl.mix13<- update(resl.mix, fixed = list(Eh ~ 1 + growthT, ln.c + Ea + Th ~ 1), start = c(-0.8, 1, 4, rep(0,4), 305))
# remove treatment on effect on ln.c + Eh + Th
resl.mix14<- update(resl.mix, fixed = list(Ea ~ 1 + growthT, ln.c + Eh + Th ~ 1), start = c(-0.8, 1, rep(0,4), 4, 305))
# remove treatment on effect on all -no
resl.mix15<- update(resl.mix, fixed = list(ln.c + Ea + Eh + Th ~ 1), start = c(-0.8, 1, 4, 305))

# multi model comparison
cands <- list(resl.mix,resl.mix1,resl.mix2,resl.mix3, resl.mix4 ,resl.mix6,resl.mix7,resl.mix10,resl.mix12 ,resl.mix13, resl.mix14,resl.mix15 )
namemod<-c("resl.mix","resl.mix1","resl.mix2","resl.mix3","resl.mix4","resl.mix6","resl.mix7","resl.mix10","resl.mix12", "resl.mix13", "resl.mix14","resl.mix15")
tab <- aictab(cand.set=cands, modnames=namemod,second.ord = TRUE)
tab


library(multcomp)

library(MuMIn)


#refit the best model here if only one model. else, make averages. now with the clean data frame, the full model is best 
resl.mix_final<-update(resl.mix, method="REML") #this is just fitting the full model 

## t000
resl_t000gr <- groupedData(logr ~ K | id.plot,subset( res_logis_clean, generation=="t000"))
resl_t000gr$growthT <- as.factor(resl_t000gr$selected)
resl_t000gr$generation <- as.factor(resl_t000gr$gen) 

## full model
resl_t000 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                  data = resl_t000gr,
                  fixed = list(ln.c + Ea + Eh + Th ~ 1),
                  random = ln.c ~ 1,
                  start = c(-0.8, 1, 4, 305), 
                  na.action = na.omit,
                  method = 'REML',
                  control = nlmeControl(pnlsTol = 0.2))

summary(resl_t000)

resl_t0001 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                  data = resl_t000gr,
                  fixed = list(ln.c + Ea + Eh + Th ~ 1),
                  random = Ea ~ 1,
                  start = c(-0.8, 1, 4, 305), 
                  na.action = na.omit,
                  method = 'REML',
                  control = nlmeControl(pnlsTol = 0.2))

summary(resl_t0001)



resl_t0002 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                   data = resl_t000gr,
                   fixed = list(ln.c + Ea + Eh + Th ~ 1),
                   random = Eh ~ 1,
                   start = c(-0.8, 1, 4, 305), 
                   na.action = na.omit,
                   method = 'REML',
                   control = nlmeControl(pnlsTol = 0.2))

summary(resl_t0002)



resl_t0003 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                   data = resl_t000gr,
                   fixed = list(ln.c + Ea + Eh + Th ~ 1),
                   random = Th ~ 1,
                   start = c(-0.46, 0.49, 4.79, 304), #model super sensitive to params  - won't run with above ones
                   na.action = na.omit,
                   method = 'REML',
                   control = nlmeControl(pnlsTol = 0.2))

summary(resl_t0003) # will still not run 

cand.list.new<-list(resl_t000,resl_t0001,resl_t0002)
aictab(cand.list.new)

confint.mm <- data.frame(intervals(resl_t000, 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]
}

confint.mm

##### predicted data for plotting
# generate predicted curves from fixed effects for GP
predK <- seq(min(resl_grp $K), max(resl_grp $K), by = 0.1)
predTemp <- predK - 273.15

predgr_anc <- schoolfield.high(
  ln.c = fixef(resl_t000)[1], 
  Ea = fixef(resl_t000)[2], 
  Eh = fixef(resl_t000)[3], 
  Th = fixef(resl_t000)[4], 
  temp = predK, Tc = 18)

#replace with averaged model or updated REML - unhashtag as needed
predgr_22 <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[1], 
  Ea = fixef(resl.mix_final)[6], 
  Eh = fixef(resl.mix_final)[11], 
  Th = fixef(resl.mix_final)[16],  
  # ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2] , 
 #Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12], 
  #Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22],
 #Th = (model.avg(mm, subset = delta < 2.0))[[2]][32], 
  temp = predK, Tc = 18) 

predgr_26 <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[1]+fixef(resl.mix_final)[2], 
  Ea = fixef(resl.mix_final)[6] +fixef(resl.mix_final)[7], 
  Eh = fixef(resl.mix_final)[11]+fixef(resl.mix_final)[12], 
  Th = fixef(resl.mix_final)[16]+fixef(resl.mix_final)[17],
  #ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4], 
  #Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14], 
  #Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24],
  #Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34], 
  temp = predK, Tc = 18) 

predgr_32 <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[1]+fixef(resl.mix_final)[3], 
  Ea = fixef(resl.mix_final)[6]+fixef(resl.mix_final)[8],  
 Eh = fixef(resl.mix_final)[11]+fixef(resl.mix_final)[13], 
  Th = fixef(resl.mix_final)[16]+fixef(resl.mix_final)[18],
#  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6], 
#  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16], 
#  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26],
 # Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36], 
  temp = predK, Tc = 18)   

predgr_FL <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[1]+fixef(resl.mix_final)[4], 
 Ea = fixef(resl.mix_final)[6]+fixef(resl.mix_final)[9],  
  Eh = fixef(resl.mix_final)[11]+fixef(resl.mix_final)[14], 
  Th = fixef(resl.mix_final)[16]+fixef(resl.mix_final)[19],
#  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8], 
#  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18], 
#  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28],
#  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38], 
  temp = predK, Tc = 18)  

predgr_FS <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[1]+fixef(resl.mix_final)[5], 
  Ea = fixef(resl.mix_final)[6]+fixef(resl.mix_final)[9],  
  Eh = fixef(resl.mix_final)[11]+fixef(resl.mix_final)[15], 
  Th = fixef(resl.mix_final)[16]+fixef(resl.mix_final)[20],
#  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10], 
#  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20], 
#  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30],
#  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40], 
  
  temp = predK, Tc = 18)  

confint.mm <- data.frame(intervals(resl.mix_final, 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]
}

confint.mm


predgr <- data.frame(temp = rep(predTemp, times = 6), 
                     growthT = c(rep('Ancestor', times = length(predK)), rep('22', times = length(predK)), rep('26', times = length(predK)), rep('32', times = length(predK)),rep('FL', times = length(predK)),rep('FS', times = length(predK))),  
                     lngr = c(predgr_anc, predgr_22, predgr_26, predgr_32, predgr_FL,predgr_FS))



#### average across replicates for plotting ####
res_logis_clean$growthT <- as.character(res_logis_clean$selected)
res_logis_clean$growthT[res_logis_clean$generation == 't000'] <- 'Ancestor'

av_gr <- ddply(res_logis_clean, .(assayT, growthT), summarise, avgr =  mean(logr, na.rm=T),avgr02=mean(r, na.rm=T),segr02= sd(r, na.rm=T)/sqrt(length(r)),segr = sd(logr, na.rm=T)/sqrt(length(logr)))

av_gr$growthT <- factor(av_gr$growthT, c('Ancestor', '22', '26', '32','FL','FS'))
av_gr<-subset(av_gr, growthT=='Ancestor'| growthT=='22' | growthT=='26'| growthT=='32'|growthT=='FL'|growthT=='FS'&segr>0)


### plot means as always
tdepGR <- ggplot(av_gr, aes(x = assayT, y = avgr, colour = growthT)) +
  geom_errorbar(aes(ymin = avgr - segr, ymax = avgr+ segr, width = 0.2), size=.9) +
  geom_point(size=2.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('darkgrey','darkgreen', 'blue','darkred','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +scale_y_continuous(limits = c(-4,0.8))+
  geom_line(data = predgr, aes(x = temp, y = lngr, colour = growthT), size=1.0) 
tdepGR 

#### fit with an order 2 polynomial, too - note this is for visualisation, not necessarily analysis 

tdepGR_poly <- ggplot(av_gr, aes(x = assayT, y = avgr, colour = growthT)) +
  geom_errorbar(aes(ymin = avgr - segr, ymax = avgr+ segr, width = 0.2), size=.9) +
  geom_point(size=2.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression((Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('darkgrey','darkgreen', 'blue','darkred','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +scale_y_continuous(limits = c(-4,0.8))+
  geom_smooth ( aes(group=growthT, colour=growthT),method="lm",se=FALSE,formula=y~poly(x,2)) 
tdepGR_poly


#or a gamm - note this is for visualisation, not necessarily analysis 


tdepGR_gam <- ggplot(av_gr, aes(x = assayT, y = avgr, colour = growthT)) +
  geom_errorbar(aes(ymin = avgr - segr, ymax = avgr+ segr, width = 0.2), size=.9) +
  geom_point(size=2.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('darkgrey','darkgreen', 'blue','darkred','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +scale_y_continuous(limits = c(-4,0.8))+
  geom_smooth ( aes(group=growthT, colour=growthT),method="gam",se=FALSE,formula= y ~ s(x, k = 3))
tdepGR_gam


pd <- position_dodge(0.5)

####POLYNOMIAL FITTED TO non ln######

#the poly on the non ln data

tdepGR_poly02 <- ggplot(av_gr, aes(x = assayT, y = avgr02, colour = growthT)) +
  geom_errorbar(aes(ymin = avgr02 - segr02, ymax = avgr02+ segr02, width = 0.2), size=.9) +
  geom_point(size=2.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression((Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('darkgreen','blue', 'darkred','darkgrey','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +
  geom_smooth ( aes(group=growthT, colour=growthT),method="lm",se=FALSE,formula=y~poly(x,4)) 
tdepGR_poly02

####GAMM FITTED TO non ln######

tdepGR_gam02 <- ggplot(av_gr, aes(x = assayT, y = avgr02, colour = growthT, alpga=growthT)) +
  geom_errorbar(aes(ymin = avgr02 - segr02, ymax = avgr02+ segr02, width = 0.2), size=.9) +
  geom_point(size=2.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression(Growth~rate~mu~day^{-1})) +
  scale_colour_manual(values = c('darkgrey','darkgreen', 'blue','red','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +
  geom_smooth ( aes(group=growthT, colour=growthT),method="gam",se=FALSE,formula= y ~ s(x, k = 5))
tdepGR_gam02

res_logis_all_st <- within(res_logis_all_st, id.plotnew <- as.character(factor(biorep):factor(selected)))# create a pond:assayT identifier 

tdepGR_gam03 <- ggplot(subset(res_logis_all_st,selected=="FS"|selected=="FL"|selected=="32"), aes(x = assayT, y = r, colour = selected)) +
  geom_point(size=2.5, position=pd) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression((Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('red','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +
  geom_smooth ( aes(group=id.plotnew, colour=selected),method="gam",se=FALSE,formula= y ~ s(x, k = 6))
tdepGR_gam03



#### AREAS under curve per growthT and biorep####


library(VGAM)
library(AUC)
library(MESS)

subtest<-subset(res_logis_clean, growthT=="FL")
MESS::auc(subtest$assayT,subtest$r, type = 'spline')


areafun<-function(d) { 
  auc_out<-auc(d$assayT,d$r,type='spline')
  out<-auc_out[1]
  output<-data.frame(out)
  output}
areadata01<-ddply(res_logis_clean,.(growthT),areafun)
areadata01 #this is per treatment 

areafun2<-function(d) { 
  auc_out<-auc(d$assayT,d$r,type='spline') 
  out<-auc_out[1]
  output<-data.frame(out)
  output}
areadata02<-ddply(res_logis_clean,.(biorep,growthT),areafun2)
areadata02 #this is per biorep per treatment 

#re-order the factors

areadata02$growthT<-factor(areadata02$growthT,  levels = c('Ancestor','22','26','32','FS','FL'))


Areaplot<-qplot(growthT,out,data=areadata02, geom="boxplot", fill=growthT)+scale_fill_manual(values=c("darkgrey","darkgreen",'blue','red','purple','black'))+theme(legend.position='none') 
Areaplot

#plot without 26 and with geom_points in there

Areaplot02<-qplot(growthT,out,data=subset(areadata02,growthT!="26"), geom="boxplot", fill=growthT)+scale_fill_manual(values=c("darkgrey","darkgreen",'red','purple','cornflowerblue'))+theme(legend.position='none') +geom_jitter(aes(colour=growthT, shape=biorep),size=2)+scale_colour_manual(values=c("darkgrey","darkgreen",'red','purple','cornflowerblue'))
Areaplot02

#### Some less fancy, but more straight forward stats ####
# can just do a mixed effects model or typeIII ANOVA on that? 

Mod0 <- lme(fixed = out ~ 1 +  growthT, random = ~ 1|biorep/growthT, method = 'ML', areadata02)

S1 <- lme(fixed = out ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', areadata02)

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

Sfinal <- lme(fixed = out ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'REML', areadata02)

krsk_S0<- kruskal.test(out ~ factor(growthT), areadata02)
krsk_S0 #yes, this finds that, too. 
library(lsmeans)
lsmeans(Sfinal, pairwise ~ growthT, adjust = "tukey") 


###### can also do this using the predicted slopes from the gam
require(gamm4)
require(MuMIn)
library(gam)
library(ROCR)


str(res_logis_clean)

res_logis_clean$growthT<-as.factor(res_logis_clean$growthT)

ts_all <- uGamm(r ~  growthT+s(assayT, by = growthT, bs = 'cs', k = 4) , random = ~ (1|biorep/growthT), REML=1, data =res_logis_clean)

# find best model
dd_ts <- dredge(ts_all, rank = 'AICc')
dd_ts
res_logis_clean$fits <- predict(ts_all$gam, type = "response")

#plot with slightly better and more sophisticated gamm
ts_plot_fits <- ggplot(res_logis_clean) +
  geom_smooth(aes(x = assayT, y = fits, colour = growthT)) +
  scale_colour_manual(values=c('darkgreen', 'blue', 'red', 'grey', 'black','purple')) +
  geom_point(aes(x = assayT, y = r, colour = growthT), alpha = 0.6,position=pd) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +
  theme(strip.background = element_blank()) +
  xlab('Assay Temperature (ºC)') +
  ylab(expression(µ[max]~(d^{-1}))) 
ts_plot_fits

#now use these for the area under curve 


areafun3<-function(d) { 
  auc_out<-auc(d$assayT,d$fits,type='spline') 
  out<-auc_out[1]
  output<-data.frame(out)
  output}
areadata03<-ddply(res_logis_clean,.(biorep,growthT),areafun3)

#re-order the factors

areadata03$growthT<-factor(areadata02$growthT,  levels = c('Ancestor','22','26','32','FS','FL'))

qplot(growthT,out,data=areadata03, fill=growthT)+scale_colour_manual(values=c("darkgrey","darkgreen",'blue','red','purple','black'))


#### curve analysis again - solving for rate at Topt and finding peaks #####
#can use pk rate? and also Topt 
#yes, do, and then calculate back to actual rates (or just extract max rate from values measured if in a complete and utter pickle and need values fast. Don't use this latter approach for proper analysis)


res_rmix <- as.data.frame(coef(resl.mix_final))
colnames(res_rmix) <- c('ln.c', 'ln.c26','ln.c32','lnc.FL','ln.cFS','Ea22','Ea26','Ea32','EaFL','EaFS','Eh22','Eh26','Eh32','EhFL','EhFS','Th22','Th26','Th32','ThFL','ThFS')
res_rmix$id <- rownames(res_rmix)
res_rmix$growthT <- substr(res_rmix$id,4,5)
res_rmix$biorep <- substr(res_rmix$id,1,2)
res_rmix$generation <- substr(res_rmix$id,7,10)
head(res_rmix)

# derive Topt - may need to melt dataframe... or just add up the random effects in the correct way
res_rmix <- as.data.frame(coef(resl.mix_final))
colnames(res_rmix) <- c('ln.c', 'ln.c26','ln.c32','lnc.FL','ln.cFS','Ea22','Ea26','Ea32','EaFL','EaFS','Eh22','Eh26','Eh32','EhFL','EhFS','Th22','Th26','Th32','ThFL','ThFS')
res_rmix$id <- rownames(res_rmix)
res_rmix$growthT <- substr(res_rmix$id,4,5)
res_rmix$biorep <- substr(res_rmix$id,1,2)
res_rmix$generation <- substr(res_rmix$id,7,10)

##### correct random coefficents for treatment effects
res_rmix$ln.c_fin <- ifelse(res_rmix$growthT == '22',  res_rmix$ln.c+fixef(resl.mix_final)[1],
                            ifelse (res_rmix$growthT=='26', res_rmix$ln.c+res_rmix$ln.c26+fixef(resl.mix_final)[2],
                                    ifelse(res_rmix$growthT=='32', res_rmix$ln.c+res_rmix$ln.c32+fixef(resl.mix_final)[3],
                                           ifelse(res_rmix$growthT=='FL',res_rmix$ln.c+res_rmix$lnc.FL+fixef(resl.mix_final)[4],res_rmix$ln.c+res_rmix$ln.cFS+fixef(resl.mix_final)[5])))) 

res_rmix$Ea_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[6],
                          ifelse (res_rmix$growthT=='26',res_rmix$Ea22+res_rmix$Ea26+ fixef(resl.mix_final)[7],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Ea22+res_rmix$Ea32  +fixef(resl.mix_final)[8],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Ea22+ fixef(resl.mix_final)[9],res_rmix$Ea22+ fixef(resl.mix_final)[10]))))


res_rmix$Eh_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[11],
                          ifelse (res_rmix$growthT=='26',res_rmix$Eh22+res_rmix$Eh26+ fixef(resl.mix_final)[12],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Eh22+res_rmix$Eh32  +fixef(resl.mix_final)[13],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Eh22+ fixef(resl.mix_final)[14],res_rmix$Eh22+ fixef(resl.mix_final)[15]))))

res_rmix$Th_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[16],
                          ifelse (res_rmix$growthT=='26',res_rmix$Th22+res_rmix$Th26+ fixef(resl.mix_final)[17],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Th22+res_rmix$Th32  +fixef(resl.mix_final)[18],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Th22+ fixef(resl.mix_final)[19],res_rmix$Th22+ fixef(resl.mix_final)[20]))))

res_rmix$Topt <- Topt(res_rmix$Eh_fin, res_rmix$Th_fin, res_rmix$Ea_fin, K = 'Y')

#now also the ancestor

res_rmt00 <- as.data.frame(coef(resl_t000))
colnames(res_rmt00) <- c('ln.c','Ea','Eh','Th')
res_rmt00$id <- rownames(res_rmt00)
res_rmt00$growthT <- substr(res_rmt00$id,4,5)
res_rmt00$biorep <- substr(res_rmt00$id,1,2)
res_rmt00$generation <- substr(res_rmt00$id,7,10)

# correct random coefficents for treatment effects
res_rmt00$ln.c_fin <- res_rmt00$ln.c+fixef(resl_t000)[1]

res_rmt00$Eh_fin <- rep(fixef(resl_t000)[3], length(res_rmt00$id))

res_rmt00$Th_fin <- rep( fixef(resl_t000)[4],length(res_rmt00$id))
res_rmt00$Ea_fin <- rep( fixef(resl_t000)[2],length(res_rmt00$id))



#### ToptC of eolved is here ####
res_rmix$ToptC<-res_rmix$Topt-273.15
res_rmix<-subset(res_rmix,Topt>0)

#### of ancestor####
# derive Topt
res_rmt00$ToptC <- Topt(res_rmt00$Eh_fin, res_rmt00$Th_fin, res_rmt00$Ea_fin, K = 'N')
res_rmt00<-subset(res_rmt00,Topt>0)

#########derive rates at Topt for peak rates/curve height is here###
res_rmix$pk.rate <- schoolfield.high(ln.c = res_rmix$ln.c_fin , Ea = res_rmix$Ea_fin  , Eh = res_rmix$Eh_fin , Th = res_rmix$Th_fin , Tc = 18, temp = res_rmix$Topt)
res_rmix$pk.rate_exp<-exp(res_rmix$pk.rate )

res_rmt00$pk.rate <- schoolfield.high(ln.c = res_rmt00$ln.c_fin , Ea = res_rmt00$Ea_fin  , Eh = res_rmt00$Eh_fin , Th = res_rmt00$Th_fin , Tc = 18, temp = res_rmt00$Topt)
res_rmt00$pk.rate_exp<-exp(res_rmt00$pk.rate )

#now make the frames more manageable
res_rmix<-res_rmix[,c(22,23,30:32)]
res_rmt00<-res_rmt00[,c(6,7,13:15)]
res_rmt00$growthT[res_rmt00$growthT == "22"] <- "Ancestor"

res_both<-rbind(res_rmix,res_rmt00)

#mixed model to solve it (top is for ToptC, pk is for peak rate)
#ToptC
Mod0top <- lme(fixed = ToptC ~ 1 +  growthT, random = ~ 1|biorep/growthT, method = 'ML', res_both)

S1top <- lme(fixed = ToptC ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', res_both)

S2top <- lme(fixed = ToptC ~ 1 , random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', res_both)
anova(S1top,S2top) # growthT is siginficant effect!!!
#need to think about some kind of post hoc too, either lsmeans or TukeyHSD? 
S1topfinal <- lme(fixed = ToptC ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'REML', res_both)

lsmeans(S1topfinal, pairwise ~ growthT, adjust = "tukey") 


krsk_S0top<- kruskal.test(ToptC ~ factor(growthT), res_both)
krsk_S0top #yes, this finds that, too. 

#peak rates 

Mod0peak <- lme(fixed = pk.rate_exp ~ 1 +  growthT, random = ~ 1|biorep/growthT, method = 'ML', res_both)

S1peak <- lme(fixed = pk.rate_exp ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', res_both)

S2peak <- lme(fixed = pk.rate_exp ~ 1 , random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', res_both)
anova(S1peak,S2peak) # growthT is siginficant effect!!!
#need to think about some kind of post hoc too, either lsmeans or TukeyHSD? 

krsk_S0peak<- kruskal.test(pk.rate_exp ~ factor(growthT), res_both)
krsk_S0peak #yes, this finds that, too. 

S1peakfin <- lme(fixed = pk.rate_exp ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'REML', res_both)

lsmeans(S1peakfin, pairwise ~ growthT, adjust = "tukey") 

#make a nice wee plot of both Topt and the peak rates

res_both$growthT<-factor(res_both$growthT,  levels = c('Ancestor','22','26','32','FS','FL'))


Toptplot02<-qplot(growthT,ToptC,data=res_both, geom="boxplot", fill=growthT)+scale_fill_manual(values=c("darkgrey","darkgreen",'blue','red','purple','black'))+ theme(legend.position='none') 
Toptplot02

#Topt without 26, but with geom_point()

Toptplot02<-qplot(growthT,ToptC,data=subset(res_both, growthT!=26), geom="boxplot", fill=growthT)+scale_fill_manual(values=c("darkgrey","darkgreen",'red','purple','cornflowerblue'))+ geom_jitter(aes(colour=growthT, shape=biorep), size=2)+theme(legend.position='none') +scale_colour_manual(values=c("darkgrey","darkgreen",'red','purple','cornflowerblue'))
Toptplot02

peakplot<-qplot(growthT,pk.rate_exp,data=res_both, geom="boxplot", fill=growthT)+scale_fill_manual(values=c("darkgrey","darkgreen",'blue','red','purple','black'))+theme(legend.position='none') 
peakplot

#peak without 26, but with geom_point()
peakplot02<-qplot(growthT,pk.rate_exp,data=subset(res_both, growthT!=26), geom="boxplot", fill=growthT)+
  scale_fill_manual(values=c("darkgrey","darkgreen",'red','purple','cornflowerblue'))+
  geom_jitter(aes(colour=growthT, shape=biorep), size=2)+
  theme(legend.position='none') +
  scale_colour_manual(values=c("darkgrey","darkgreen",'red','purple','cornflowerblue'))
peakplot02

####think about how to do width of curve #####

#find percantage of growthrate/temp pairs where growth is higher than 
#calculate growth at lowest temp/growth at lowest when growthT==22C
#then growth at highest temp/growth at highest when growthT==22C  - could just devide all of them by mean22? 
#no, we just look at the range of temperatures where growth treat > growth 22

##### Preliminary composite figure 02 ####

plot_grid(tdepGR_gam02,Areaplot,Toptplot02,peakplot, labels=c("A","B","C","D"),ncol=2)


### back to the old stuff #####
av_gr$flucs<-ifelse(av_gr$growthT=="FS"|av_gr$growthT=="FL","fluc","stable")
predgr$flucs<-ifelse(predgr$growthT=="FS"|predgr$growthT=="FL","fluc","stable")
pd <- position_dodge(0.5)



#plot same, but make alpha on stable ones, and drop the ancestor 
tdepGR2 <- ggplot(subset(av_gr, growthT!="Ancestor"), aes(x = assayT, y = avgr, colour = growthT, alpha=flucs)) +
  geom_errorbar(aes(ymin = avgr - segr, ymax = avgr+ segr, width = 0.2), size=.69) +
  geom_point(size=1.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('darkgreen', 'blue','darkred','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +scale_alpha_manual(values=c(0.9,0.36))+
  theme(legend.position='none') +
  geom_line(data = subset(predgr,growthT!='Ancestor'), aes(x = temp, y = lngr, colour = growthT), size=1.0) 
tdepGR2

res_logis_clean$flucs<-ifelse(res_logis_clean$growthT=="FS"|res_logis_clean$growthT=="FL", "fluc","stable")

#plot with a simple gam , include 26 and the flucs . a polynomial works better here -- finds the curves in b3 better, but have left in gam for now. 
tdepgr2.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="32"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('blue','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica') +geom_line(aes(group=growthT, alpha=flucs),stat="smooth",method="gam",se=FALSE,formula= y ~ s(x, k = 3), size=1)
  theme(legend.position='none') 
tdepgr2.1


#plot with a simple gam , include 32 and the flucs 

tdepgr2.2 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="26"), aes(x = assayT, y = r, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('red','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica') +geom_line(aes(group=growthT, alpha=flucs),stat="smooth",method="gam",se=FALSE,formula= y ~ s(x, k = 3), size=1)+
  theme(legend.position='none') 
tdepgr2.2

#now adjust for the random effect for bioreplicate to plot thermal tolerance curves for each parameter, then 
#look at randoms
randses<-data.frame(random.effects(resl.mix_final))
randses$id<-rownames(randses)

#now make predicted data frames... per biorep!! fun!! :P blergh ..there must be a smart way to do this? 

predgr_22_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b1:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b1:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b1:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b1:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_22_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b2:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b2:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b2:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b2:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_22_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b3:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b3:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b3:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b3:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_22_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b4:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b4:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b4:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b4:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_22_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b5:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b5:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b5:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b5:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_22_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b6:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b6:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b6:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b6:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b1:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b1:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b1:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b1:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b2:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b2:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b2:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b2:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_26_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b3:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b3:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b3:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b3:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b4:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b4:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b4:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b4:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b5:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b5:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b5:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b5:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b6:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b6:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b6:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b6:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b1:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b1:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b1:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b1:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_32_b2 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b2:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b2:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b2:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b2:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b3 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b3:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b3:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b3:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b3:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b4 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b4:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b4:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b4:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b4:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b5 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b5:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b5:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b5:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b5:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b6 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b6:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b6:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b6:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b6:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_FL_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b1:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b1:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b1:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b1:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b2:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b2:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b2:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b2:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b3:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b3:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b3:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b3:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b4:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b4:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b4:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b4:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b5:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b5:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b5:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b5:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b6:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b6:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b6:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b6:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  


predgr_FS_b1 <- schoolfield.high(

  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b1:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b1:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b1:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b1:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_FS_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b2:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b2:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b2:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b2:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b3:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b3:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b3:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b3:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b4:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b4:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b4:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b4:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b5:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b5:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b5:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b5:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b6:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b6:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b6:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b6:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  


#make new predgr data frame. am probably over-thinking this, and there is a more straight forward way. 
predgr_biorep <- data.frame(temp = rep(predTemp, times = 24), 
                     growthTid = c( rep('26_b1', times = length(predK)),rep('26_b2', times = length(predK)),rep('26_b3', times = length(predK)),rep('26_b4', times = length(predK)),rep('26_b5', times = length(predK)),rep('26_b6', times = length(predK)),
                                    rep('32_b1', times = length(predK)), rep('32_b2', times = length(predK)), rep('32_b3', times = length(predK)), rep('32_b4', times = length(predK)), rep('32_b5', times = length(predK)), rep('32_b6', times = length(predK)),
                                    rep('FL_b1', times = length(predK)),rep('FL_b2', times = length(predK)),rep('FL_b3', times = length(predK)),rep('FL_b4', times = length(predK)),rep('FL_b5', times = length(predK)),rep('FL_b6', times = length(predK)),
                                    rep('FS_b1', times = length(predK)),rep('FS_b2', times = length(predK)),rep('FS_b3', times = length(predK)),rep('FS_b4', times = length(predK)),rep('FS_b5', times = length(predK)),rep('FS_b6', times = length(predK))),  
                     lngr = c(predgr_26_b1,predgr_26_b2,predgr_26_b3,predgr_26_b4,predgr_26_b5,predgr_26_b6, 
                              predgr_32_b1,predgr_32_b2,predgr_32_b3,predgr_32_b4,predgr_32_b5,predgr_32_b6,
                              predgr_FL_b1,predgr_FL_b2,predgr_FL_b3,predgr_FL_b4,predgr_FL_b5,predgr_FL_b6,
                              predgr_FS_b1,predgr_FS_b2,predgr_FS_b3,predgr_FS_b4,predgr_FS_b5,predgr_FS_b6))
head(predgr_biorep)
predgr_biorep$growthT<-substr(predgr_biorep$growthTid,1,2)
predgr_biorep$biorep<-substr(predgr_biorep$growthTid,4,5)


pd <- position_dodge(0.5)

res_logis_clean$flucs<-ifelse(res_logis_clean$growthT=="FS"|res_logis_clean$growthT=="FL", "fluc","stable")
predgr_biorep$flucs<-ifelse(predgr_biorep$growthT=="FS"| predgr_biorep$growthT=="FL","fluc","stable")

#plot with fitted lines per biorep and facet per biorep, using 26º 
tdepgr3.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="32"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('blue','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="32"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.1

#plot with fitted lines per biorep and facet per biorep, using 32º 
tdepgr3.2 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="26"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('red','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="26"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.2


#facet by growthT rather than biorep for 26ºC

tdepgr3.1.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="32"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('blue','black','purple')) + facet_grid(~growthT)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="32"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.1.1 # or could make FS and FL individual, then 26 overlay 


# and now facet per growthT for 32ºC

tdepgr3.2.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="26"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('red','black','purple')) + facet_grid(~growthT)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="26"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.2.1


#get Topt and peak rates for growth, too

res_rmix <- as.data.frame(coef(resl.mix_final))
colnames(res_rmix) <- c('ln.c', 'ln.c26','ln.c32','lnc.FL','ln.cFS','Ea22','Ea26','Ea32','EaFL','EaFS','Eh22','Eh26','Eh32','EhFL','EhFS','Th22','Th26','Th32','ThFL','ThFS')
res_rmix$id <- rownames(res_rmix)
res_rmix$growthT <- substr(res_rmix$id,4,5)
res_rmix$biorep <- substr(res_rmix$id,1,2)
res_rmix$generation <- substr(res_rmix$id,7,10)

# correct random coefficents for treatment effects
res_rmix$ln.c_fin <- ifelse(res_rmix$growthT == '22',  res_rmix$ln.c+fixef(resl.mix_final)[1],
                            ifelse (res_rmix$growthT=='26', res_rmix$ln.c+res_rmix$ln.c26+fixef(resl.mix_final)[2],
                                    ifelse(res_rmix$growthT=='32', res_rmix$ln.c+res_rmix$ln.c32+fixef(resl.mix_final)[3],
                                           ifelse(res_rmix$growthT=='FL',res_rmix$ln.c+res_rmix$lnc.FL+fixef(resl.mix_final)[4],res_rmix$ln.c+res_rmix$ln.cFS+fixef(resl.mix_final)[5])))) 

res_rmix$Ea_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[6],
                          ifelse (res_rmix$growthT=='26',res_rmix$Ea22+res_rmix$Ea26+ fixef(resl.mix_final)[7],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Ea22+res_rmix$Ea32  +fixef(resl.mix_final)[8],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Ea22+ fixef(resl.mix_final)[9],res_rmix$Ea22+ fixef(resl.mix_final)[10]))))


res_rmix$Eh_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[11],
                          ifelse (res_rmix$growthT=='26',res_rmix$Eh22+res_rmix$Eh26+ fixef(resl.mix_final)[12],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Eh22+res_rmix$Eh32  +fixef(resl.mix_final)[13],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Eh22+ fixef(resl.mix_final)[14],res_rmix$Eh22+ fixef(resl.mix_final)[15]))))

res_rmix$Th_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[16],
                          ifelse (res_rmix$growthT=='26',res_rmix$Th22+res_rmix$Th26+ fixef(resl.mix_final)[17],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Th22+res_rmix$Th32  +fixef(resl.mix_final)[18],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Th22+ fixef(resl.mix_final)[19],res_rmix$Th22+ fixef(resl.mix_final)[20]))))



# derive Topt
res_rmix$Topt <- Topt(res_rmix$Eh_fin, res_rmix$Th_fin, res_rmix$Ea_fin, K = 'Y')

res_rmix$ToptC<-res_rmix$Topt-273.15
res_rmix<-subset(res_rmix,Topt>0)
#derive rates at Topt
res_rmix$pk.rate <- schoolfield.high(ln.c = res_rmix$ln.c_fin , Ea = res_rmix$Ea_fin  , Eh = res_rmix$Eh_fin , Th = res_rmix$Th_fin , Tc = 18, temp = res_rmix$Topt)
res_rmix$pk.rate_exp<-exp(res_rmix$pk.rate )

### now we make the means
res_rav<-ddply(res_rmix,.(growthT),summarise,meanTopt=mean(ToptC), seTopt=sd(ToptC)/sqrt(6),meanpk=mean(pk.rate_exp),sdpk=sd(pk.rate_exp))
res_rav
res_rav$logback<-log(res_rav$meanpk)

#### now the same for the ancestor 
res_rmt00 <- as.data.frame(coef(resl_t000))
colnames(res_rmt00) <- c('ln.c','Ea','Eh','Th')
res_rmt00$id <- rownames(res_rmt00)
res_rmt00$growthT <- substr(res_rmt00$id,4,5)
res_rmt00$biorep <- substr(res_rmt00$id,1,2)
res_rmt00$generation <- substr(res_rmt00$id,7,10)

# correct random coefficents for treatment effects
res_rmt00$ln.c_fin <- res_rmt00$ln.c+fixef(resl_t000)[1]

res_rmt00$Eh_fin <- rep(fixef(resl_t000)[3], length(res_rmt00$id))

res_rmt00$Th_fin <- rep( fixef(resl_t000)[4],length(res_rmt00$id))
res_rmt00$Ea_fin <- rep( fixef(resl_t000)[2],length(res_rmt00$id))

# derive Topt
res_rmt00$Topt <- Topt(res_rmt00$Eh_fin, res_rmt00$Th_fin, res_rmt00$Ea_fin, K = 'N')
res_rmt00<-subset(res_rmt00,Topt>0)
#derive rates at Topt
res_rmt00$pk.rate <- schoolfield.high(ln.c = res_rmt00$ln.c_fin , Ea = res_rmt00$Ea_fin  , Eh = res_rmt00$Eh_fin , Th = res_rmt00$Th_fin , Tc = 18, temp = res_rmt00$Topt+273.15)
res_rmt00$pk.rate_exp<-exp(res_rmt00$pk.rate )

#we keep Topt (and peak rate if needed for later) pf each, along with growthT, biorep, and generation. 
res_evo<-res_rmix[,c(21:24,30:32)]
res_anc<-res_rmt00[,c(5:8,13:15)]
colnames(res_anc)<-c("id","growthT"  ,   "biorep"    ,  "generation" ,  "ToptC",        "pk.rate"  ,   "pk.rate_exp")

Topt_allR<-rbind(res_evo,res_anc)

head(Topt_allR)
Topt_allR$growthT<-ifelse(Topt_allR$generation=="t000","Ancestor",Topt_allR$growthT)
Topt_allR$specialrep<-ifelse(Topt_allR$biorep=="b1"|Topt_allR$biorep=="b2"|Topt_allR$biorep=="b5","pref22","pref32")
Topt_allR$growthTnew<-ifelse(Topt_allR$growthT=="FL",as.character(factor(Topt_allR$growthT):factor(Topt_allR$specialrep)),Topt_allR$growthT) # could have done in one step. SIGH.

#now we can plot this 
#sort out order of Temps 
Topt_allR$growthTnew<-factor(Topt_allR$growthTnew,  levels = c('Ancestor','22','26','32','FS','FL:pref22','FL:pref32'))

#???looks like 32 and FS are swapped?? how?? 
qplot( as.factor(growthTnew), ToptC, data = subset(Topt_allR, growthT!="22" & growthT!="26"&growthT!="Ancestor"), geom="boxplot", fill=growthT)+scale_fill_manual(values=c("red","black","purple"))


######now the same but we want only FS, FL and 32ºC######

###### run as mixed effects model - need to make new id column that has gen, sel, biorep, but not assayT

resl_grp <- groupedData(logr ~ K | id.plot,subset( res_logis_clean, generation=="t300" & selected=="32" | selected=="FS" | selected=="FL"))
resl_grp$growthT <- as.factor(resl_grp$selected)
resl_grp$generation <- as.factor(resl_grp$gen) 

## full model  - FL won't run?? if included with only FS and 26 -- so do full model and drop the ones that are not necessary later for plotting 
resl.mix <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                 data = resl_grp,
                 fixed = list(ln.c + Ea + Eh + Th ~ 1 + growthT),
                 random = pdDiag(ln.c + Ea + Eh +Th ~ 1),
                 start = c(-.9, rep(0,2), 1, rep(0,2), 8, rep(0,2), 306, rep(0,2)), # starting parameters from nlsLoop
                 na.action = na.omit,
                 method = 'ML',
                 control = nlmeControl(pnlsTol = 0.5))
summary(resl.mix)



## simplify treatment effects - singles
# remove treatment on effect on Ea 
resl.mix1<- update(resl.mix, fixed = list(Ea ~ 1, ln.c + Eh + Th ~ 1 + growthT), start = c(-0.8, rep(0,2), 1, 4, rep(0,2), 305, rep(0,2)))
# remove treatment on effect on Eh 
resl.mix2<- update(resl.mix, fixed = list(Eh ~ 1, ln.c + Ea + Th ~ 1 + growthT), start = c(-0.8, rep(0,2), 1, rep(0,2), 4, 305, rep(0,2))) 
# remove treatment on effect on Th 
resl.mix3<- update(resl.mix, fixed = list(ln.c + Ea + Eh ~ 1 + growthT, Th ~ 1), start = c(-0.8, rep(0,2), 1, rep(0,2), 4, rep(0,2), 305))
# remove treatment on effect on ln.c 
resl.mix4<- update(resl.mix, fixed = list(ln.c ~ 1, Ea + Eh + Th ~ 1 + growthT), start = c(-0.8, 1, rep(0,2), 4, rep(0,2), 305, rep(0,2)))

## simplify treatment effects - pairs
# remove treatment on effect on Ea + Eh no 
resl.mix5<- update(resl.mix, fixed = list( Ea + Eh ~ 1,ln.c + Th ~ 1 + growthT), start = c(-0.8, rep(0,2), 1, 4, 305, rep(0,2)))
# remove treatment on effect on Ea + Th
resl.mix6<- update(resl.mix, fixed = list(ln.c + Eh ~ 1 + growthT, Ea + Th ~ 1), start = c(-0.8, rep(0,2), 1, 4, rep(0,2), 305))
# remove treatment on effect on Eh + Th
resl.mix7<- update(resl.mix, fixed = list(ln.c + Ea ~ 1 + growthT, Eh + Th ~ 1), start = c(-0.8, rep(0,2), 1, rep(0,2), 4, 305))
# remove treatment on effect on ln.c + Ea

resl.mix8<- update(resl.mix, fixed = list(ln.c + Ea ~ 1, Eh + Th ~ 1 + growthT), start = c(-0.8, 1, 4, rep(0,2), 305, rep(0,2)))
# remove treatment on effect on ln.c + Eh no
resl.mix9<- update(resl.mix, fixed = list(ln.c + Eh ~ 1, Ea + Th ~ 1 + growthT), start = c(-0.8, 1, rep(0,2), 4, 305, rep(0,2)))
# remove treatment on effect on ln.c + Th

resl.mix10<- update(resl.mix, fixed = list(Ea + Eh ~ 1 + growthT, ln.c + Th ~ 1), start = c(-0.8, 1, rep(0,2), 4, rep(0,2), 305)) 

## simplify treatment effects - Three
# remove treatment on effect on Ea + Eh + Th
resl.mix11<- update(resl.mix, fixed = list(ln.c ~ 1 + growthT, Ea + Eh + Th ~ 1), start = c(-0.8, rep(0,2), 1, 4, 305))
# remove treatment on effect on ln.c + Ea + Eh
resl.mix12<- update(resl.mix, fixed = list(ln.c + Ea + Eh ~ 1, Th ~ 1 + growthT), start = c(-0.8, 1, 4, 305, rep(0,2)))
# remove treatment on effect on ln.c + Ea + Th
resl.mix13<- update(resl.mix, fixed = list(Eh ~ 1 + growthT, ln.c + Ea + Th ~ 1), start = c(-0.8, 1, 4, rep(0,2), 305))
# remove treatment on effect on ln.c + Eh + Th
resl.mix14<- update(resl.mix, fixed = list(Ea ~ 1 + growthT, ln.c + Eh + Th ~ 1), start = c(-0.8, 1, rep(0,2), 4, 305))
# remove treatment on effect on all -no
resl.mix15<- update(resl.mix, fixed = list(ln.c + Ea + Eh + Th ~ 1), start = c(-0.8, 1, 4, 305))

# multi model comparison
cands<- list(resl.mix,resl.mix1,resl.mix2,resl.mix3 ,resl.mix6,resl.mix7,resl.mix10,resl.mix11,resl.mix12 ,resl.mix13, resl.mix14,resl.mix15 )
namemod<-c("resl.mix","resl.mix1","resl.mix2","resl.mix3","resl.mix6","resl.mix7","resl.mix10","resl.mix11","resl.mix12", "resl.mix13", "resl.mix14","resl.mix15")
tab <- aictab(cand.set=cands, modnames=namemod,second.ord = TRUE)
tab 

library(MuMIn)
library(multcomp)



# models 3, and full,  need to be averaged full is just below 2 lower, so all three needed 
mm<-model.sel(cands)
model.avg(mm)
model.avg(mm, subset = delta < 2.0)
#models with delta.aicc < 2
summary(model.avg(mm, subset = delta < 2.0)) # finds that ln.c, Eh, are different not Ea and Th 


#refit the best model here if only one model. else, make averages. now with the clean data frame, the full model is best 
resl.mix_final<-update(resl.mix, method="REML") 

## t000
resl_t000gr <- groupedData(logr ~ K | id.plot,subset( res_logis_clean, generation=="t000"))
resl_t000gr$growthT <- as.factor(resl_t000gr$selected)
resl_t000gr$generation <- as.factor(resl_t000gr$gen) 

## full model
resl_t000 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                  data = resl_t000gr,
                  fixed = list(ln.c + Ea + Eh + Th ~ 1),
                  random = ln.c ~ 1,
                  start = c(-0.8, 1, 4, 305), 
                  na.action = na.omit,
                  method = 'REML',
                  control = nlmeControl(pnlsTol = 0.2))

summary(resl_t000)

resl_t0001 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                   data = resl_t000gr,
                   fixed = list(ln.c + Ea + Eh + Th ~ 1),
                   random = Ea ~ 1,
                   start = c(-0.8, 1, 4, 305), 
                   na.action = na.omit,
                   method = 'REML',
                   control = nlmeControl(pnlsTol = 0.2))

summary(resl_t0001)



resl_t0002 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                   data = resl_t000gr,
                   fixed = list(ln.c + Ea + Eh + Th ~ 1),
                   random = Eh ~ 1,
                   start = c(-0.8, 1, 4, 305), 
                   na.action = na.omit,
                   method = 'REML',
                   control = nlmeControl(pnlsTol = 0.2))

summary(resl_t0002)



resl_t0003 <- nlme(logr ~ schoolfield.high(ln.c,Ea,Eh,Th,temp=K,Tc=18),
                   data = resl_t000gr,
                   fixed = list(ln.c + Ea + Eh + Th ~ 1),
                   random = Th ~ 1,
                   start = c(-0.46, 0.49, 4.79, 304), #model super sensitive to params  - won't run with above ones
                   na.action = na.omit,
                   method = 'REML',
                   control = nlmeControl(pnlsTol = 0.2))

summary(resl_t0003) # will still not run 

cand.list.new<-list(resl_t000,resl_t0001,resl_t0002)
aictab(cand.list.new)

confint.mm <- data.frame(intervals(resl_t000, 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]
}

confint.mm

##### predicted data for plotting
# generate predicted curves from fixed effects for GP
predK <- seq(min(resl_grp $K), max(resl_grp $K), by = 0.1)
predTemp <- predK - 273.15

predgr_anc <- schoolfield.high(
  ln.c = fixef(resl_t000)[1], 
  Ea = fixef(resl_t000)[2], 
  Eh = fixef(resl_t000)[3], 
  Th = fixef(resl_t000)[4], 
  temp = predK, Tc = 18)

#replace with averaged model or updated REML - unhashtag as needed
predgr_22 <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[2], 
  Ea = fixef(resl.mix_final)[1], 
  Eh = fixef(resl.mix_final)[7], 
  Th = fixef(resl.mix_final)[12],  
 # ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2] , 
#  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12], 
#  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22],
#  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32], 
  temp = predK, Tc = 18) 

predgr_26 <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[2]+fixef(resl.mix_final)[3], 
  Ea = fixef(resl.mix_final)[1], 
  Eh = fixef(resl.mix_final)[7]+fixef(resl.mix_final)[8], 
  Th = fixef(resl.mix_final)[12]+fixef(resl.mix_final)[13],
  #ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4], 
  #Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14], 
  #Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24],
  #Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34], 
  temp = predK, Tc = 18) 

predgr_32 <- schoolfield.high(
  ln.c = fixef(resl.mix_final)[2]+fixef(resl.mix_final)[4], 
    Ea = fixef(resl.mix_final)[1], 
    Eh = fixef(resl.mix_final)[7]+fixef(resl.mix_final)[9], 
    Th = fixef(resl.mix_final)[12]+fixef(resl.mix_final)[14],
  #ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6], 
  #Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16], 
  #Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26],
  #Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36], 
  temp = predK, Tc = 18)   

predgr_FL <- schoolfield.high(
  #  ln.c = fixef(resl.mix_final)[2]+fixef(resl.mix_final)[5], 
  #  Ea = fixef(resl.mix_final)[1], 
  #  Eh = fixef(resl.mix_final)[7]+fixef(resl.mix_final)[10], 
  #  Th = fixef(resl.mix_final)[12]+fixef(resl.mix_final)[15],
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8], 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18], 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28],
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38], 
  temp = predK, Tc = 18)  

predgr_FS <- schoolfield.high(
  # ln.c = fixef(resl.mix_final)[2]+fixef(resl.mix_final)[6], 
  #  Ea = fixef(resl.mix_final)[1], 
  #  Eh = fixef(resl.mix_final)[7]+fixef(resl.mix_final)[11], 
  #  Th = fixef(resl.mix_final)[12]+fixef(resl.mix_final)[16],
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10], 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20], 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30],
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40], 
  
  temp = predK, Tc = 18)  

confint.mm <- data.frame(intervals(resl.mix_final, 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]
}

confint.mm


predgr <- data.frame(temp = rep(predTemp, times = 6), 
                     growthT = c(rep('Ancestor', times = length(predK)), rep('22', times = length(predK)), rep('26', times = length(predK)), rep('32', times = length(predK)),rep('FL', times = length(predK)),rep('FS', times = length(predK))),  
                     lngr = c(predgr_anc, predgr_22, predgr_26, predgr_32, predgr_FL,predgr_FS))



### average across replicates for plotting
res_logis_clean$growthT <- as.character(res_logis_clean$selected)
res_logis_clean$growthT[res_logis_clean$generation == 't000'] <- 'Ancestor'

av_gr <- ddply(res_logis_clean, .(assayT, growthT), summarise, avgr =  mean(logr, na.rm=T), segr = sd(logr, na.rm=T)/sqrt(length(logr)))

av_gr$growthT <- factor(av_gr$growthT, c('Ancestor', '22', '26', '32','FL','FS'))
av_gr<-subset(av_gr, growthT=='Ancestor'| growthT=='22' | growthT=='26'| growthT=='32'|growthT=='FL'|growthT=='FS'&segr>0)


### plot means as always
tdepGR <- ggplot(av_gr, aes(x = assayT, y = avgr, colour = growthT)) +
  geom_errorbar(aes(ymin = avgr - segr, ymax = avgr+ segr, width = 0.2), size=.9) +
  geom_point(size=2.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('darkgrey','darkgreen', 'blue','darkred','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +
  theme(legend.position='top') +scale_y_continuous(limits = c(-4,0.8))+
  geom_line(data = predgr, aes(x = temp, y = lngr, colour = growthT), size=1.0) 
tdepGR 


av_gr$flucs<-ifelse(av_gr$growthT=="FS"|av_gr$growthT=="FL","fluc","stable")
predgr$flucs<-ifelse(predgr$growthT=="FS"|predgr$growthT=="FL","fluc","stable")
pd <- position_dodge(0.5)

#plot same, but make alpha on stable ones, and drop the ancestor 
tdepGR2 <- ggplot(subset(av_gr, growthT!="Ancestor"), aes(x = assayT, y = avgr, colour = growthT, alpha=flucs)) +
  geom_errorbar(aes(ymin = avgr - segr, ymax = avgr+ segr, width = 0.2), size=.69) +
  geom_point(size=1.5) +     labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(Growth~rate~mu~day^{-1}))) +
  scale_colour_manual(values = c('darkgreen', 'blue','darkred','black','purple')) +
  theme_cowplot(font_size = 14, font_family = 'Helvetica') +scale_alpha_manual(values=c(0.9,0.36))+
  theme(legend.position='none') +
  geom_line(data = subset(predgr,growthT!='Ancestor'), aes(x = temp, y = lngr, colour = growthT), size=1.0) 
tdepGR2

res_logis_clean$flucs<-ifelse(res_logis_clean$growthT=="FS"|res_logis_clean$growthT=="FL", "fluc","stable")

#plot with a simple gam , include 26 and the flucs . a polynomial works better here -- finds the curves in b3 better, but have left in gam for now. 
tdepgr2.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="32"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('blue','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica') +geom_line(aes(group=growthT, alpha=flucs),stat="smooth",method="gam",se=FALSE,formula= y ~ s(x, k = 3), size=1)+
  theme(legend.position='none') 
tdepgr2.1


#plot with a simple gam , include 32 and the flucs 

tdepgr2.2 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="26"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('red','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica') +geom_line(aes(group=growthT, alpha=flucs),stat="smooth",method="gam",se=FALSE,formula= y ~ s(x, k = 3), size=1)+
  theme(legend.position='none') 
tdepgr2.2

#now adjust for the random effect for bioreplicate to plot thermal tolerance curves for each parameter, then 
#look at randoms
randses<-data.frame(random.effects(resl.mix_final))
randses$id<-rownames(randses)

#now make predicted data frames... per biorep!! fun!! :P blergh ..there must be a smart way to do this? 

predgr_22_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b1:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b1:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b1:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b1:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_22_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b2:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b2:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b2:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b2:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_22_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b3:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b3:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b3:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b3:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_22_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b4:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b4:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b4:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b4:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_22_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b5:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b5:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b5:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b5:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_22_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+randses[randses$id =="b6:22:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +randses[randses$id =="b6:22:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+randses[randses$id =="b6:22:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+randses[randses$id =="b6:22:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b1:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b1:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b1:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b1:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b2:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b2:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b2:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b2:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_26_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b3:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b3:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b3:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b3:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b4:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b4:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b4:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b4:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b5:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b5:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b5:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b5:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_26_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][4]+randses[randses$id =="b6:26:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][14]+randses[randses$id =="b6:26:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][24]+randses[randses$id =="b6:26:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][34]+randses[randses$id =="b6:26:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b1:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b1:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b1:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b1:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_32_b2 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b2:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b2:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b2:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b2:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b3 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b3:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b3:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b3:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b3:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b4 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b4:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b4:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b4:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b4:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b5 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b5:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b5:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b5:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b5:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 

predgr_32_b6 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][6]+randses[randses$id =="b6:32:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][16]+randses[randses$id =="b6:32:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][26]+randses[randses$id =="b6:32:t300",]$Eh..Intercept., 
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][36]+randses[randses$id =="b6:32:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_FL_b1 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b1:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b1:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b1:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b1:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b2:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b2:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b2:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b2:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b3:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b3:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b3:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b3:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b4:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b4:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b4:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b4:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b5:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b5:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b5:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b5:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FL_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][8]+randses[randses$id =="b6:FL:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][18]+randses[randses$id =="b6:FL:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][28]+randses[randses$id =="b6:FL:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][38]+randses[randses$id =="b6:FL:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  


predgr_FS_b1 <- schoolfield.high(
  
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b1:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b1:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b1:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b1:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18) 


predgr_FS_b2 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b2:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b2:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b2:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b2:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b3 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b3:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b3:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b3:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b3:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b4 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b4:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b4:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b4:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b4:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b5 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b5:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b5:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b5:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b5:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  

predgr_FS_b6 <- schoolfield.high(
  ln.c =(model.avg(mm, subset = delta < 2.0))[[2]][2]+ (model.avg(mm, subset = delta < 2.0))[[2]][10]+randses[randses$id =="b6:FS:t300",]$ln.c..Intercept., 
  Ea = (model.avg(mm, subset = delta < 2.0))[[2]][12] +(model.avg(mm, subset = delta < 2.0))[[2]][20]+randses[randses$id =="b6:FS:t300",]$Ea..Intercept., 
  Eh = (model.avg(mm, subset = delta < 2.0))[[2]][22]+(model.avg(mm, subset = delta < 2.0))[[2]][30]+randses[randses$id =="b6:FS:t300",]$Eh..Intercept.,
  Th = (model.avg(mm, subset = delta < 2.0))[[2]][32]+(model.avg(mm, subset = delta < 2.0))[[2]][40]+randses[randses$id =="b6:FS:t300",]$Th..Intercept., 
  temp = predK, Tc = 18)  


#make new predgr data frame. am probably over-thinking this, and there is a more straight forward way. 
predgr_biorep <- data.frame(temp = rep(predTemp, times = 24), 
                            growthTid = c( rep('26_b1', times = length(predK)),rep('26_b2', times = length(predK)),rep('26_b3', times = length(predK)),rep('26_b4', times = length(predK)),rep('26_b5', times = length(predK)),rep('26_b6', times = length(predK)),
                                           rep('32_b1', times = length(predK)), rep('32_b2', times = length(predK)), rep('32_b3', times = length(predK)), rep('32_b4', times = length(predK)), rep('32_b5', times = length(predK)), rep('32_b6', times = length(predK)),
                                           rep('FL_b1', times = length(predK)),rep('FL_b2', times = length(predK)),rep('FL_b3', times = length(predK)),rep('FL_b4', times = length(predK)),rep('FL_b5', times = length(predK)),rep('FL_b6', times = length(predK)),
                                           rep('FS_b1', times = length(predK)),rep('FS_b2', times = length(predK)),rep('FS_b3', times = length(predK)),rep('FS_b4', times = length(predK)),rep('FS_b5', times = length(predK)),rep('FS_b6', times = length(predK))),  
                            lngr = c(predgr_26_b1,predgr_26_b2,predgr_26_b3,predgr_26_b4,predgr_26_b5,predgr_26_b6, 
                                     predgr_32_b1,predgr_32_b2,predgr_32_b3,predgr_32_b4,predgr_32_b5,predgr_32_b6,
                                     predgr_FL_b1,predgr_FL_b2,predgr_FL_b3,predgr_FL_b4,predgr_FL_b5,predgr_FL_b6,
                                     predgr_FS_b1,predgr_FS_b2,predgr_FS_b3,predgr_FS_b4,predgr_FS_b5,predgr_FS_b6))
head(predgr_biorep)
predgr_biorep$growthT<-substr(predgr_biorep$growthTid,1,2)
predgr_biorep$biorep<-substr(predgr_biorep$growthTid,4,5)


pd <- position_dodge(0.5)

res_logis_clean$flucs<-ifelse(res_logis_clean$growthT=="FS"|res_logis_clean$growthT=="FL", "fluc","stable")
predgr_biorep$flucs<-ifelse(predgr_biorep$growthT=="FS"| predgr_biorep$growthT=="FL","fluc","stable")

#plot with fitted lines per biorep and facet per biorep, using 26º 
tdepgr3.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="32"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('blue','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="32"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.1

#plot with fitted lines per biorep and facet per biorep, using 32º 
tdepgr3.2 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="26"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('red','black','purple')) + facet_grid(~biorep)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="26"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.2


#facet by growthT rather than biorep for 26ºC

tdepgr3.1.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="32"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('blue','black','purple')) + facet_grid(~growthT)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="32"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.1.1 # or could make FS and FL individual, then 26 overlay 


# and now facet per growthT for 32ºC

tdepgr3.2.1 <- ggplot(subset(res_logis_clean, growthT!="Ancestor"&growthT!="22"&growthT!="26"), aes(x = assayT, y = logr, colour = growthT, alpha=flucs, shape=biorep)) +
  geom_point(size=1.5, position=pd)  + scale_alpha_manual(values=c(0.9,0.56))+
  scale_colour_manual(values = c('red','black','purple')) + facet_grid(~growthT)+
  labs(x=expression(Assay~Temperature~(degree*C)), y=expression(ln~(italic(µ)~day^{-1}))) +
  theme_cowplot(font_size = 12, font_family = 'Helvetica')+
  geom_line(data = subset(predgr_biorep,growthT!="26"), aes(x = temp, y = lngr, colour = growthT, alpha=flucs),size=1) +theme(legend.position='none') 
tdepgr3.2.1


#get Topt and peak rates for growth, too

res_rmix <- as.data.frame(coef(resl.mix_final))
colnames(res_rmix) <- c('ln.c', 'ln.c26','ln.c32','lnc.FL','ln.cFS','Ea22','Ea26','Ea32','EaFL','EaFS','Eh22','Eh26','Eh32','EhFL','EhFS','Th22','Th26','Th32','ThFL','ThFS')
res_rmix$id <- rownames(res_rmix)
res_rmix$growthT <- substr(res_rmix$id,4,5)
res_rmix$biorep <- substr(res_rmix$id,1,2)
res_rmix$generation <- substr(res_rmix$id,7,10)

# correct random coefficents for treatment effects
res_rmix$ln.c_fin <- ifelse(res_rmix$growthT == '22',  res_rmix$ln.c+fixef(resl.mix_final)[1],
                            ifelse (res_rmix$growthT=='26', res_rmix$ln.c+res_rmix$ln.c26+fixef(resl.mix_final)[2],
                                    ifelse(res_rmix$growthT=='32', res_rmix$ln.c+res_rmix$ln.c32+fixef(resl.mix_final)[3],
                                           ifelse(res_rmix$growthT=='FL',res_rmix$ln.c+res_rmix$lnc.FL+fixef(resl.mix_final)[4],res_rmix$ln.c+res_rmix$ln.cFS+fixef(resl.mix_final)[5])))) 

res_rmix$Ea_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[6],
                          ifelse (res_rmix$growthT=='26',res_rmix$Ea22+res_rmix$Ea26+ fixef(resl.mix_final)[7],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Ea22+res_rmix$Ea32  +fixef(resl.mix_final)[8],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Ea22+ fixef(resl.mix_final)[9],res_rmix$Ea22+ fixef(resl.mix_final)[10]))))


res_rmix$Eh_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[11],
                          ifelse (res_rmix$growthT=='26',res_rmix$Eh22+res_rmix$Eh26+ fixef(resl.mix_final)[12],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Eh22+res_rmix$Eh32  +fixef(resl.mix_final)[13],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Eh22+ fixef(resl.mix_final)[14],res_rmix$Eh22+ fixef(resl.mix_final)[15]))))

res_rmix$Th_fin <- ifelse(res_rmix$growthT == '22', fixef(resl.mix_final)[16],
                          ifelse (res_rmix$growthT=='26',res_rmix$Th22+res_rmix$Th26+ fixef(resl.mix_final)[17],
                                  ifelse(res_rmix$growthT=='32',res_rmix$Th22+res_rmix$Th32  +fixef(resl.mix_final)[18],
                                         ifelse(res_rmix$growthT=='FL',res_rmix$Th22+ fixef(resl.mix_final)[19],res_rmix$Th22+ fixef(resl.mix_final)[20]))))



# derive Topt
res_rmix$Topt <- Topt(res_rmix$Eh_fin, res_rmix$Th_fin, res_rmix$Ea_fin, K = 'Y')

res_rmix$ToptC<-res_rmix$Topt-273.15
res_rmix<-subset(res_rmix,Topt>0)
#derive rates at Topt
res_rmix$pk.rate <- schoolfield.high(ln.c = res_rmix$ln.c_fin , Ea = res_rmix$Ea_fin  , Eh = res_rmix$Eh_fin , Th = res_rmix$Th_fin , Tc = 18, temp = res_rmix$Topt)
res_rmix$pk.rate_exp<-exp(res_rmix$pk.rate )

### now we make the means
res_rav<-ddply(res_rmix,.(growthT),summarise,meanTopt=mean(ToptC), seTopt=sd(ToptC)/sqrt(6),meanpk=mean(pk.rate_exp),sdpk=sd(pk.rate_exp))
res_rav
res_rav$logback<-log(res_rav$meanpk)

#### now the same for the ancestor 
res_rmt00 <- as.data.frame(coef(resl_t000))
colnames(res_rmt00) <- c('ln.c','Ea','Eh','Th')
res_rmt00$id <- rownames(res_rmt00)
res_rmt00$growthT <- substr(res_rmt00$id,4,5)
res_rmt00$biorep <- substr(res_rmt00$id,1,2)
res_rmt00$generation <- substr(res_rmt00$id,7,10)

# correct random coefficents for treatment effects
res_rmt00$ln.c_fin <- res_rmt00$ln.c+fixef(resl_t000)[1]

res_rmt00$Eh_fin <- rep(fixef(resl_t000)[3], length(res_rmt00$id))

res_rmt00$Th_fin <- rep( fixef(resl_t000)[4],length(res_rmt00$id))
res_rmt00$Ea_fin <- rep( fixef(resl_t000)[2],length(res_rmt00$id))

# derive Topt
res_rmt00$Topt <- Topt(res_rmt00$Eh_fin, res_rmt00$Th_fin, res_rmt00$Ea_fin, K = 'N')
res_rmt00<-subset(res_rmt00,Topt>0)
#derive rates at Topt
res_rmt00$pk.rate <- schoolfield.high(ln.c = res_rmt00$ln.c_fin , Ea = res_rmt00$Ea_fin  , Eh = res_rmt00$Eh_fin , Th = res_rmt00$Th_fin , Tc = 18, temp = res_rmt00$Topt+273.15)
res_rmt00$pk.rate_exp<-exp(res_rmt00$pk.rate )

#we keep Topt (and peak rate if needed for later) pf each, along with growthT, biorep, and generation. 
res_evo<-res_rmix[,c(21:24,30:32)]
res_anc<-res_rmt00[,c(5:8,13:15)]
colnames(res_anc)<-c("id","growthT"  ,   "biorep"    ,  "generation" ,  "ToptC",        "pk.rate"  ,   "pk.rate_exp")

Topt_allR<-rbind(res_evo,res_anc)

head(Topt_allR)
Topt_allR$growthT<-ifelse(Topt_allR$generation=="t000","Ancestor",Topt_allR$growthT)
Topt_allR$specialrep<-ifelse(Topt_allR$biorep=="b1"|Topt_allR$biorep=="b2"|Topt_allR$biorep=="b5","pref22","pref32")
Topt_allR$growthTnew<-ifelse(Topt_allR$growthT=="FL",as.character(factor(Topt_allR$growthT):factor(Topt_allR$specialrep)),Topt_allR$growthT) # could have done in one step. SIGH.

#now we can plot this 
#sort out order of Temps 
Topt_allR$growthTnew<-factor(Topt_allR$growthTnew,  levels = c('Ancestor','22','26','32','FS','FL:pref22','FL:pref32'))

#???looks like 32 and FS are swapped?? how?? 
qplot( as.factor(growthTnew), ToptC, data = subset(Topt_allR, growthT!="22" & growthT!="26"&growthT!="Ancestor"), geom="boxplot", fill=growthT)+scale_fill_manual(values=c("red","black","purple"))


