###### extracted parameters from FFRF for Thally02 paper. boxplots with bioreps super-imposed
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)
library(rgl)

# dir
setwd("~/Users/elisaschaum/Dropbox/Thally long fluc paper/Thally actual selection experiment/fastact by time and temp/t300 full curves")

#work computer wants this to be

("~/Dropbox/Thally long fluc paper/Thally actual selection experiment/fastact by time and temp/t300 full curves")


####read in rP data
rP_extract<-read.csv("rP_extracted_values.csv")

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


Iopt_rP<-qplot(as.factor(growthT), Iopt, data=rP_extract, geom="boxplot", fill=growthT, ylab="Iopt (PAR)", xlab="Selection Environment",facets=~assayT)+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgreen","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgreen","red","violet","darkgrey"))+
  theme_classic(base_size = 12, base_family = 'Helvetica')+theme(legend.position = "none")
Iopt_rP

Pmax_rP<-qplot(as.factor(growthT), Pmax, data=rP_extract, geom="boxplot", fill=growthT, ylab="rPmax", xlab="Selection Environment",facets=~assayT)+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgreen","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgreen","red","violet","darkgrey"))+
  theme_classic(base_size = 12, base_family = 'Helvetica')+theme(legend.position = "none")
Pmax_rP

a_rP<-qplot(as.factor(growthT), a, data=rP_extract, geom="boxplot", fill=growthT, ylab="alpha", xlab="Selection Environment",facets=~assayT)+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgreen","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgreen","red","violet","darkgrey"))+
  theme_classic(base_size = 12, base_family = 'Helvetica')+theme(legend.position = "none")
a_rP

grid_rP<-plot_grid(Iopt_rP,Pmax_rP, ncol=1, labels=c("A","B"))


#Iopt stats, may need a mixed model 

MMIopt<-lme(fixed = Iopt ~   assayT*growthT, random = ~ 1|biorep, method = 'ML', rP_extract) 
dM1<-dredge(MMIopt)
dM1

MMIoptfin<-lme(fixed = Iopt ~   assayT*growthT, random = ~ 1|biorep/growthT, method = 'REML', rP_extract) 
summary(MMIoptfin)

MMmax<-lme(fixed = Pmax ~   assayT*growthT, random = ~ 1|biorep/growthT, method = 'ML', rP_extract) 
dMax<-dredge(MMmax)
dMax

MMmaxfin<-lme(fixed = Pmax ~   assayT*growthT, random = ~ 1|biorep/growthT, method = 'REML', rP_extract) 
summary(MMmaxfin)


###now the psiPSII

psi_rP<-read.csv('psiPSI_slope_intercept.csv')

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

Slope<-qplot(as.factor(growthT), Slope, data=psi_rP, geom="boxplot", fill=growthT, ylab=expression(Slope~psi), xlab="Selection Environment",facets=~assayT)+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgreen","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgreen","red","violet","darkgrey"))+
  theme_classic(base_size = 12, base_family = 'Helvetica')+theme(legend.position = "none")
Slope

Intercept<-qplot(as.factor(growthT), Decay_intercept, data=psi_rP, geom="boxplot", fill=growthT, ylab=expression(Intercept~psi), xlab="Selection Environment",facets=~assayT)+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgreen","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgreen","red","violet","darkgrey"))+
  theme_classic(base_size = 12, base_family = 'Helvetica')+theme(legend.position = "none")
Intercept

grid_psi<-plot_grid(Slope,Intercept, ncol=1, labels=c("C","D"))

MMSlope<-lme(fixed = Slope ~   assayT*growthT, random = ~ 1|growthT/biorep, method = 'ML', psi_rP) 
dSlope<-dredge(MMSlope)
dSlope

MMSlopefin<-lme(fixed = Slope ~   assayT*growthT, random = ~ 1|biorep, method = 'REML', psi_rP) 
summary(MMSlopefin)

MMSlope<-lme(fixed = Slope ~   assayT*growthT, random = ~ 1|growthT/biorep, method = 'ML', psi_rP) 
dSlope<-dredge(MMSlope)
dSlope

MMint<-lme(fixed = Decay_intercept ~   assayT*growthT, random = ~ 1|biorep/growthT, method = 'ML', psi_rP) 
dInt<-dredge(MMint)
dInt
MMintfin<-lme(fixed = Decay_intercept ~   assayT*growthT, random = ~ 1|biorep/growthT, method = 'REML', psi_rP) 
summary(MMintfin)

#### now the closed centres

Closed_rP<-read.csv("CvaluesThally02.csv")

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

Ymax<-qplot(as.factor(growthT), ymax, data=Closed_rP, geom="boxplot", fill=growthT, ylab="Ymax", xlab="Selection Environment",facets=~assayT)+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgreen","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgreen","red","violet","darkgrey"))+
  theme_classic(base_size = 12, base_family = 'Helvetica')+theme(legend.position = "none")
Ymax

Ks<-qplot(as.factor(growthT), KS, data=Closed_rP, geom="boxplot", fill=growthT, ylab="Ks", xlab="Selection Environment",facets=~assayT)+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT))+
  scale_fill_manual(values=c("darkgreen","darkred","purple","black"))+
  scale_colour_manual(values=c("lightgreen","red","violet","darkgrey"))+
  theme_classic(base_size = 12, base_family = 'Helvetica')+theme(legend.position = "none")
Ks

grid_C<-plot_grid(Ymax,Ks, ncol=1, labels=c("E","F"))

grid_all<-plot_grid(grid_rP,grid_psi,grid_C, ncol=3)


MMymax<-lme(fixed = ymax ~   assayT*growthT, random = ~ 1|growthT/biorep, method = 'ML', Closed_rP) 
dym<-dredge(MMymax)
dym
MMymaxfin<-lme(fixed = ymax ~   assayT*growthT, random = ~ 1|growthT/biorep, method = 'REML', Closed_rP) 
summary(MMymaxfin)

MMKS<-lme(fixed = KS ~   assayT*growthT, random = ~ 1|biorep/growthT, method = 'ML', Closed_rP) 
dKS<-dredge(MMKS)
dKS
MMKSfin<-lme(fixed = KS ~   growthT, random = ~ 1|biorep/growthT, method = 'REML', Closed_rP) 
summary(MMKSfin)
