#### to not further mess up the NPvsMu_Gab file, this takes the slopes of PS and growth for temps ≤32 and plots them against each other. Then, we colour the dots by assayT to see where they lie. 
rm(list=ls())


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)
require(cowplot)
require(lme4)
require(visreg)
library(MuMIn)


setwd("~/Dropbox/Thally long fluc paper")
#### read in newest big file 
bf<-read.csv("20170327_giantbigtable.csv")
bf$growthT <- factor(bf$growthT, levels = c('Ancestor', '22', '26', '32', 'FS','FL'))
head(bf)
#make 20 22 and 25 26 ? if not already there 
bf$assayT[bf$assayT == 25] <- 26
bf$assayT[bf$assayT == 20] <- 22


#plot bivariate correlations of untransformed data per biorep per assayT==growthT OR growthT ==Ancestor 
bfsub<-subset(bf, assayT==growthT | growthT=="Ancestor")
#for plotting, keep only 

#growth and NP (net photo)
#growth and CUE
#growth and N:P 

#


# make summary per growthT 
bfav<-ddply(bf,c("growthT","assayT","generation") , function(df) return(c(ChlC.avg=mean(df$Chl.C),ChlC.sd=sd(df$Chl.C),CN.avg=mean(df$V),CN.sd=sd(df$V))))

#conversions
# NP over 12h
bf$NP_CperC <- ((bf$NP.umol.ml.min * 32 * bf$M * (12/44) * 60 * 24) / (bf$Ccell_umol * bf$cellsml)) * 0.5


bf$GP_CperC <- ((bf$GP.umol.ml.min * 32 * bf$M * (12/44) * 60 * 24) / (bf$Ccell_umol * bf$cellsml)) 

bf$R_CperC <- ((bf$R.umol.ml.min * 32 * bf$M * (12/44) * 60 * 24) / (bf$Ccell_umol * bf$cellsml)) 

# NP calculated the propper way
bf$NetC <- bf$GP_CperC * 0.5 - bf$R_CperC 

subset(bf,growthT=="FS"& assayT==26)$NetC/12

head(bf)
#make averages , this time for GP, R, NP
#make a new column that splits the FL into the proper fractions ...or use the none-averages 


bfav02<-ddply(bf,c("growthT","assayT","generation") , function(df) return(c(GP.avg=mean(df$Gpscaled.mumol.min.cell),GP.se=sd(df$Gpscaled.mumol.min.cell)/100,Resp.avg=mean(df$RESPscaled.mumol.min.cell),Resp.se=sd(df$RESPscaled.mumol.min.cell)/(100))))


bf<-within(bf, id.plot<-as.character(as.factor(growthT):as.factor(biorep)))

Gpplasplot02<-qplot(as.factor(assayT), log(Gpscaled.mumol.min.cell), data=subset(bf, assayT>20&assayT<33 &assayT!=30), col=growthT,shape=biorep)+
  geom_smooth(aes(group=id.plot), method="lm", se=F)+
  facet_wrap(~growthT)+
  scale_colour_manual(values=c("darkgrey","darkgreen","blue","red","purple","black"))+
  theme_classic(base_size = 14, base_family = 'Helvetica')
Gpplasplot02

#
slopefn<-function(d) { 
  m<-lm (log(Gpscaled.mumol.min.cell)~assayT,data=d) #can force through origin by setting offset to 0 
  sum.<-summary(m) 
  r2<-sum.$r.squared
  intercept<-m$coefficients[1]
  slope<-m$coefficients[2]
  output<-data.frame(slope,intercept,r2)
  output}
slopesub<-subset(bf, assayT>20&assayT<33 &assayT!=30)
summary(slopesub)
slopedata<-ddply(slopesub,.(slopesub$growthT,slopesub$biorep),slopefn)
#slopedata[1,3]
slopedata

slopefn02<-function(d) { 
  m<-lm (log(growthrate)~assayT,data=d) #can force through origin by setting offset to 0 
  sum.<-summary(m) 
  r2<-sum.$r.squared
  intercept<-m$coefficients[1]
  slope<-m$coefficients[2]
  output<-data.frame(slope,intercept,r2)
  output}
slopesub02<-subset(bf, assayT>20&assayT<33 &assayT!=30)
summary(slopesub02)
slopedata02<-ddply(slopesub02,.(slopesub02$growthT,slopesub02$biorep),slopefn02)
#slopedata[1,3]
slopedata02

colnames(slopedata02)[1]<-"growthT"
colnames(slopedata02)[2]<-"biorep"


Gpplasplot05<-qplot(as.factor(growthT), slope, data=slopedata02, geom="boxplot", fill=growthT, ylab="Evolved change in growth (Slopes of reaction norms for photosynthesis)",xlab="Selection environment")+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgrey","darkgreen","blue","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgrey","lightgreen","cornflowerblue","red","violet","darkgrey"))+
  theme_classic(base_size = 14, base_family = 'Helvetica')
Gpplasplot05


#make biorep id column for both growth and PS slope data 
slopedata$id <- as.character(factor(slopedata$growthT):factor(slopedata$biorep))
slopedata02$id <- as.character(factor(slopedata02$growthT):factor(slopedata02$biorep))
slopedata02<-slopedata02[,c(3,6)]
colnames(slopedata02)[1]<-"slope_growth"

plasmerge<-merge(slopedata,slopedata02)
#get rid of NAs

plasmerge <- plasmerge[! is.na(plasmerge$slope_growth),]


Gpplasplot06<-qplot(slope_growth, slope, data=plasmerge, colour=growthT,ylab="Plastic response",xlab="Growth response")+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  geom_smooth(aes(group=growthT), method="lm",se=FALSE)+
  scale_colour_manual(values=c("lightgrey","darkgreen","blue","red","purple","black"))+
  theme_classic(base_size = 14, base_family = 'Helvetica')
Gpplasplot06



MM2 <- lme(fixed=slope_growth~slope*growthT, random = ~  1|biorep, method="REML",    na.action = na.omit,subset(plasmerge, growthT=="22")) 
predict(MM2, type = "response")
visreg2d(MM2, "slope", "growthT", plot.type="persp",type="conditional",col="purple") # not sure how to interpret 
visreg2d(MM2, "slope", "slope_growth", plot.type="gg") # not sure how to interpret 




bf$CVol_total<-(bf$Ccell_umol/bf$V)*bf$cellsml*bf$V

bf$GP_V <- ((bf$GP.umol.ml.min * 32 * bf$M * (12/44) * 60 * 24) / (bf$CVol_total)) 

bf$R_V <- ((bf$R.umol.ml.min * 32 * bf$M * (12/44) * 60 * 24) / (bf$CVol_total)) 


# id 
bf$id <- as.character(factor(bf$growthT):factor(bf$assayT):factor(bf$biorep))

# drop values that are negative or outliers
drop <- c("22:20:b1", "22:20:b4",       "22:25:b1",       "22:25:b2",      
          "22:25:b3",       "22:25:b4",       "22:30:b1",       "22:30:b3",      
          "22:30:b4",       "22:30:b6",       "22:32:b1",       "22:32:b2",      
          "22:32:b3",       "22:32:b4",       "22:32:b6",       "26:25:b5",      
          "26:30:b5",       "26:32:b5",       "26:32:b6",       "26:35:b3",      
          "26:35:b5",       "32:20:b3",       "32:25:b3",       "32:25:b5",      
          "32:30:b1",       "32:30:b2",       "32:30:b3",       "32:32:b3",      
          "32:32:b5",       "32:35:b2",       "32:40:b2",       "Ancestor:20:b3",
          "Ancestor:20:b4", "Ancestor:20:b6", "Ancestor:25:b3", "Ancestor:25:b4",
          "Ancestor:25:b5", "Ancestor:25:b6", "Ancestor:30:b1", "Ancestor:30:b3",
          "Ancestor:30:b5", "Ancestor:30:b6", "Ancestor:32:b3", "Ancestor:32:b6",
          "Ancestor:35:b6", "FS:25:b6",       "FS:40:b2",       "FS:40:b3")  

bf <- bf[! bf$id %in% drop, ]

# id from plotting
bf$id_ran <- as.character(factor(bf$growthT):factor(bf$biorep))

# drop dodgy biorep
bf <- bf[! bf$id_ran %in% c('Ancestor:b6', 'Ancestor:b5', '26:b5', '26:b6'),]

# sort out growthT
bf$growthT <- factor(bf$growthT, levels = c('Ancestor', '22', '26', '32', 'FS','FL'))
bf$lnR <- log2(bf$R_CperC)
bf$K<-bf$assayT+273.15

bf<-subset(bf, NetC>0)#&NetC<4) #  & growthT!="Ancestor")

rgl_init <- function(new.device = FALSE, bg = "white", width = 640) { 
  if( new.device | rgl.cur() == 0 ) {
    rgl.open()
    par3d(windowRect = 50 + c( 0, 0, width, width ) )
    rgl.bg(color = bg )
  }
  rgl.clear(type = c("shapes", "bboxdeco"))
  rgl.viewpoint(theta = 15, phi = 20, zoom = 0.7)
}


bf$lnC<-log(bf$NetC)
bf$lngr<-log(bf$growthrate)
bf_av<-ddply(bf, .(assayT, growthT), summarise, avlnC =  mean(lnC, na.rm=T), seGP = sd(lnC, na.rm=T)/sqrt(length(lnC)),avgr=mean(lngr, na.rm=T))

rgl_init() #to open a new window to actually plot the 3d planes... 
cols<-list("green","blue","red","purple","black")


MM2 <- lme(fixed=exp(lngr)~exp(lnC)*+assayT, random = ~  1|biorep, method="REML", bf) 
predict(MM2, type = "response")
visreg2d(MM2, "lnC", "assayT", plot.type="persp",type="conditional",col="purple") # not sure how to interpret 
visreg2d(MM2, "lnC", "assayT", plot.type="rgl",col="purple") # not sure how to interpret 
