llibrary(readxl)
Thermalsingle1 <- read_excel("~/Desktop/Thermal gametes/doi_10.5061_dryad.wstqjq2s9__v2/Thermalsingle1.xlsx")

View(Thermalsingle1)
View(Thermalsingle1)


Thermalsinglefin <- read_excel("~/Desktop/Thermal gametes/doi_10.5061_dryad.wstqjq2s9__v2/Thermalsinglefin.xlsx")
View(Thermalsinglefin)

EGG <- factor(c("W", "C"))
SPERM <- factor(c("W", "C"))
DEV <- factor(c("W", "C"))

str(Thermalsinglefin)

Thermalsinglefin$MALE <- as.factor(Thermalsinglefin$MALE)
Thermalsinglefin$FAMILY <- as.factor(Thermalsinglefin$FAMILY)
Thermalsinglefin$TRAY <- as.factor(Thermalsinglefin$TRAY)
Thermalsinglefin$FEMALE <- as.factor(Thermalsinglefin$FEMALE)
Thermalsinglefin$FRYD <- as.numeric (Thermalsinglefin$FRYD)
Thermalsinglefin$TPF <- as.numeric (Thermalsinglefin$TPF)
Thermalsinglefin$EGG <- as.factor(Thermalsinglefin$EGG)
Thermalsinglefin$SPERM <- as.factor(Thermalsinglefin$SPERM)
Thermalsinglefin$DEV <- as.factor(Thermalsinglefin$DEV)
Thermalsinglefin$eggcount <- as.numeric(Thermalsinglefin$eggcount)
Thermalsinglefin$FRYTOTPROP <- as.numeric(Thermalsinglefin$FRYTOTPROP)
Thermalsinglefin$DDAY <- as.numeric(Thermalsinglefin$DDAY)


Thermalsingle1$MALE <- as.factor(Thermalsingle1$MALE)
Thermalsingle1$FAMILY <- as.factor(Thermalsingle1$FAMILY)
Thermalsingle1$TRAY <- as.factor(Thermalsingle1$TRAY)
Thermalsingle1$FEMALE <- as.factor(Thermalsingle1$FEMALE)
Thermalsingle1$FRYD <- as.numeric (Thermalsingle1$FRYD)
Thermalsingle1$TPF <- as.numeric (Thermalsingle1$TPF)
Thermalsingle1$EGG <- as.factor(Thermalsingle1$EGG)
Thermalsingle1$SPERM <- as.factor(Thermalsingle1$SPERM)
Thermalsingle1$DEV <- as.factor(Thermalsingle1$DEV)
Thermalsingle1$eggcount <- as.numeric(Thermalsingle1$eggcount)
Thermalsingle1$FRYTOTPROP <- as.numeric(Thermalsingle1$FRYTOTPROP)





library (ggplot2) 
library(Rmisc)    
library(tidyverse) 
library(ggpubr)
library(ggsignif)
library(MASS)
library(DHARMa)
library(lme4)
library(car)
library(lmtest)
library(emmeans)


# modelling embryonic age at hatch (Degree days) across treatments, families as random factor

M1DDAY<-lmer(DDAY ~ EGG*SPERM*DEV +(1|FAMILY), data = Thermalsinglefin)
M1bDDAY<-lmer(DDAY ~ 1 +(1|FAMILY), data = Thermalsinglefin)
lrtest(M1DDAY)
anova(M1DDAY, M1bDDAY)
summary(M1DDAY)
coef(M1DDAY)
joint_tests(M1DDAY)
confint(M1DDAY)


Anova(M1DDAY, type="3")

s1M1DDAY <- simulateResiduals(M1DDAY)
plot(s1M1DDAY)
testDispersion(s1M1DDAY)

#Estimate the modelled means...also backtransformed to the response variable units from log scale
# by using type=response
library(emmeans)

m1dday.emm <- emmeans(M1DDAY, specs="EGG", type="response")
m1dday.emm

m1dday.emm_dt<-data.frame(m1dday.emm)

#contrasts pairwise differences
m1dday.emm.pairs <- contrast(m1dday.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m1dday.emm.pairs


m1dday.emm <- emmeans(M1DDAY, specs="SPERM", type="response")
m1dday.emm


m1dday.emm_dt<-data.frame(m1dday.emm)

#contrasts pairwise differences
m1dday.emm.pairs <- contrast(m1dday.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m1dday.emm.pairs

m1dday.emm <- emmeans(M1DDAY, specs="DEV", type="response")
m1dday.emm

m1dday.emm_dt<-data.frame(m1dday.emm)

#contrasts pairwise differences
m1dday.emm.pairs <- contrast(m1dday.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m1dday.emm.pairs

#Contrast test for trends due to presence of significant interactions

ESDDAY.emm <- emmeans(M1DDAY, ~ EGG * SPERM * DEV)
ESDDAY.emm

pairs(ESDDAY.emm, simple = list("EGG", c("SPERM", "DEV")))

pairs(ESDDAY.emm, simple = "each")

emmDDAY_s.t <- emmeans(M1DDAY, pairwise ~ EGG | DEV)

## NOTE: Results may be misleading due to involvement in interactions
emmDDAY_s.t
contrast(emmDDAY_s.t[[1]], "poly")   
ICDDAY_st <- contrast(regrid(emmDDAY_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
ICDDAY_st
coef(ICDDAY_st)

test(ICDDAY_st, joint = TRUE)


emm_s.t <- emmeans(M1, pairwise ~ SPERM | DEV)
## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)

#The three-way interaction may be explored via interaction contrasts too:

contrast(emmeans(M1, ~ EGG*SPERM*DEV),
         interaction = c("poly", "consec", "consec"))







# modelling embryonic age at hatch across treatments, families as random factor

M1<-lmer(TPF ~ EGG*SPERM*DEV +(1|FAMILY), data = Thermalsinglefin)
M1b<-lmer(TPF ~ 1 +(1|FAMILY), data = Thermalsinglefin)
lrtest(M1, M1b)
anova(M1, M1b)
summary(M1)
coef(M1)
joint_tests(M1)
confint(M1)


Anova(M1, type="3")

s1 <- simulateResiduals(M1)
plot(s1)
testDispersion(s1)







Anova(lm1 <-lmer(TPF ~ EGG*SPERM*DEV +(1|FAMILY), data = Thermalsinglefin, 
                 contrasts=list(EGG=contr.poly)), type="III")

Anova(lm1, type="3")
summary(lm1)
coef(summary(lm1))
joint_tests(lm1)

s1 <- simulateResiduals(lm1)
plot(s1)
testDispersion(s1)

#Estimate the modelled means...also backtransformed to the response variable units from log scale
# by using type=response
library(emmeans)

m1.emm <- emmeans(M1, specs="EGG", type="response")
m1.emm

m1.emm_dt<-data.frame(m1.emm)

#contrasts pairwise differences
m1.pairs <- contrast(m1.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m1.pairs


m1.emm <- emmeans(M1, specs="SPERM", type="response")
m1.emm


m1.emm_dt<-data.frame(m1.emm)

#contrasts pairwise differences
m1.pairs <- contrast(m1.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m1.pairs

m1.emm <- emmeans(M1, specs="DEV", type="response")
m1.emm

m1.emm_dt<-data.frame(m1.emm)

#contrasts pairwise differences
m1.pairs <- contrast(m1.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m1.pairs

#Contrast test for trends due to presence of significant interactions

ESD.emm <- emmeans(M1, ~ EGG * SPERM * DEV)
ESD.emm

pairs(ESD.emm, simple = list("EGG", c("SPERM", "DEV")))

pairs(ESD.emm, simple = "each")

emm_s.t <- emmeans(M1, pairwise ~ EGG | DEV)

## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)


emm_s.t <- emmeans(M1, pairwise ~ SPERM | DEV)
## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)

#The three-way interaction may be explored via interaction contrasts too:

contrast(emmeans(M1, ~ EGG*SPERM*DEV),
         interaction = c("poly", "consec", "consec"))





#modelling the proportion of fry hatched per day as response variable to egg, sperm and developmental temperature (C or W)

Thermal_switch_R_NEWtime01 <- read_excel("~/Desktop/Thermal gametes/doi_10.5061_dryad.wstqjq2s9__v2/Thermal_switch_R_NEWtime01.xlsx")
View(Thermal_switch_R_NEWtime01)
#just to double check that the matrix contains the values I need
y<-(cbind(Thermal_switch_R_NEWtime01$FRYD, Thermal_switch_R_NEWtime01$eggcount))

M2<-glmer(y ~ EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"), data = Thermal_switch_R_NEWtime01)
summary(M2)

M2<-glmer(cbind(FRYD, eggcount) ~ EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"), data = Thermal_switch_R_NEWtime01)

summary(M2)
coef(summary(M2))

M2TPF<-glmer(cbind(FRYD, eggcount) ~ EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"), data = Thermal_switch_R_NEWtime01)

summary(M2)
coef(summary(M2))
confint(M2)

joint_tests(M2)
#Estimate the modelled means for M2...also backtransformed to the response variable units from log scale
# by using type=response
library(emmeans)
m2.emm <- emmeans(M2, specs="EGG", type="response")
m2.emm

m2.emm_dt<-data.frame(m2.emm)

#contrasts pairwise differences
m2.pairs <- contrast(m2.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m2.pairs

m2.emm <- emmeans(M2, specs="SPERM", type="response")
m2.emm

m2.emm_dt<-data.frame(m2.emm)

#contrasts pairwise differences
m2.pairs <- contrast(m2.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m2.pairs


m2.emm <- emmeans(M2, specs="DEV", type="response")
m2.emm

m2.emm_dt<-data.frame(m2.emm)

#contrasts pairwise differences
m2.pairs <- contrast(m2.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m2.pairs

#Contrast test for trends due to presence of significant interactions

ESD.emm <- emmeans(M2, ~ EGG * SPERM * DEV)
ESD.emm

pairs(ESD.emm, simple = list("EGG", c("SPERM", "DEV")))

pairs(ESD.emm, simple = "each")

emm_s.t <- emmeans(M2, pairwise ~ EGG | DEV)

## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)


emm_s.t <- emmeans(M2, pairwise ~ SPERM | DEV)
## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)

#The three-way interaction may be explored via interaction contrasts too:

contrast(emmeans(M2, ~ EGG*SPERM*DEV),
         interaction = c("poly", "consec", "consec"))

#modelling the same thing as a coulmn with proportion of fry hatching per day 
M2<-glmer(FRYPOP ~ EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"), data = Thermal_switch_R_NEWtime01)

summary(M2)
coef(summary(M2))

#just try with two decimal points


s2 <- simulateResiduals(M2)
plot(s2)
testDispersion(s2)

plot(resid(M2, type = "pearson") ~ fitted(M2))
qqnorm(resid(M2, type = "pearson"))
qqline(resid(M2, type = "pearson"))


Anova(lm2<-glmer(cbind(eggcount/FRYD) ~ EGG*SPERM*DEV + (1|FAMILY), family=binomial(link=logit), data = Thermal_switch_R_NEWtime01, 
                 contrasts=list(EGG=contr.poly)), type="III")


Anova(lm2, type="III")

m <- lmer(cbind(FRYD, eggcount)~ EGG*SPERM*DEV + (1|FAMILY), data = Thermal_switch_R_NEWtime01)

Anova(m, type="III")


#Estimate the modelled means for M2...also backtransformed to the response variable units from log scale
# by using type=response
library(emmeans)
m2.emm <- emmeans(M2, specs="EGG", type="response")
m2.emm

m2.emm_dt<-data.frame(m2.emm)

#contrasts pairwise differences
m2.pairs <- contrast(m2.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m2.pairs

m2.emm <- emmeans(M2, specs="SPERM", type="response")
m2.emm

m2.emm_dt<-data.frame(m2.emm)

#contrasts pairwise differences
m2.pairs <- contrast(m2.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m2.pairs


m2.emm <- emmeans(M2, specs="DEV", type="response")
m2.emm

m2.emm_dt<-data.frame(m2.emm)

#contrasts pairwise differences
m2.pairs <- contrast(m2.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m2.pairs

#Contrast test for trends due to presence of significant interactions

ESD.emm <- emmeans(M2, ~ EGG * SPERM * DEV)
ESD.emm

pairs(ESD.emm, simple = list("EGG", c("SPERM", "DEV")))

pairs(ESD.emm, simple = "each")

emm_s.t <- emmeans(M2, pairwise ~ EGG | DEV)

## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)


emm_s.t <- emmeans(M2, pairwise ~ SPERM | DEV)
## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)

#The three-way interaction may be explored via interaction contrasts too:

contrast(emmeans(M2, ~ EGG*SPERM*DEV),
         interaction = c("poly", "consec", "consec"))

library(Rmisc)

summarizedhatch<-summarySE(Thermal_switch_R_NEWtime01, measurevar = "FRYPERTOT", groupvars=c("TRAY"), na.rm=TRUE)

library(writexl)

write_xlsx(summarizedhatch, "OneDrive - University of East Anglia/R stats/R stats\\summarizedhatchtreat1.xlsx")

library(readxl)

summarizedhatchtreat2 <- read_excel("~/OneDrive - University of East Anglia/R stats/summarizedhatchtreat2.xlsx", 
                                    col_types = c("text", "numeric", "numeric", 
                                                  "numeric", "numeric", "numeric", 
                                                  "text", "text", "text", "text"))
View(summarizedhatchtreat2)


#graphing

library(ggplot2)
require(pscl)
require(MASS)
require(boot)


#creating cumulative data for graphing
require(dplyr)

#adding a column with cumulative hatch for graphing purposes (from a different dataset because I need to average the hatchings per day on the same day to have cumulative averages, otherwise having eliminated the zeroes different samples have different lenghts)




library(readxl)
Thermal_switch_R_single_MAY <- read_excel("Thermal_switch_R_single_MAY.xlsx")
View(Thermal_switch_R_single_MAY)
library(readxl)
Thermal_switch_R_NEWtime02 <- read_excel("Thermal_switch_R_NEWtime02.xlsx")
View(Thermal_switch_R_NEWtime02)
library(magrittr)
Thermal_switch_R_single_MAY %<>%
  mutate(DEVELOPMENTALT = case_when(
    TREATMENT == "WWW" ~ "WARM",
    TREATMENT == "CWW" ~ "WARM",
    TREATMENT == "WCW" ~ "WARM",
    TREATMENT == "CCW" ~ "WARM",
    TREATMENT == "WWC" ~ "COLD",
    TREATMENT == "WCC" ~ "COLD",
    TREATMENT == "CWC" ~ "COLD",
    TREATMENT == "CCC" ~ "COLD"))

Thermal_switch_R_single_MAY$MALE <- as.factor(Thermal_switch_R_single_MAY$MALE)
Thermal_switch_R_single_MAY$TREATMENT <- as.factor(Thermal_switch_R_single_MAY$TREATMENT)
Thermal_switch_R_single_MAY$FAMILY <- as.factor(Thermal_switch_R_single_MAY$FAMILY)
Thermal_switch_R_single_MAY$TRAY <- as.factor(Thermal_switch_R_single_MAY$TRAY)
Thermal_switch_R_single_MAY$FEMALE <- as.factor(Thermal_switch_R_single_MAY$FEMALE)
Thermal_switch_R_single_MAY$FRYD <- as.numeric (Thermal_switch_R_single_MAY$FRYD)
Thermal_switch_R_single_MAY$TPF <- as.numeric (Thermal_switch_R_single_MAY$TPF)
Thermal_switch_R_single_MAY$EGG <- as.factor(Thermal_switch_R_single_MAY$EGG)
Thermal_switch_R_single_MAY$SPERM <- as.factor(Thermal_switch_R_single_MAY$SPERM)
Thermal_switch_R_single_MAY$DEV <- as.factor(Thermal_switch_R_single_MAY$DEV)
Thermal_switch_R_single_MAY$eggcount <- as.numeric(Thermal_switch_R_single_MAY$eggcount)
Thermal_switch_R_single_MAY$FRYPROP <- as.numeric(Thermal_switch_R_single_MAY$FRYPROP)


Thermal_switch_R_single_MAY$csum <- ave(Thermal_switch_R_single_MAY$FRYPROP, Thermal_switch_R_single_MAY$FAMILY,Thermal_switch_R_single_MAY$TREATMENT, FUN=cumsum)

View(Thermal_switch_R_single_MAY)
as.data.frame(Thermal_switch_R_single_MAY)
str(Thermal_switch_R_single_MAY)
Thermal_switch_R_single_MAY$TREATMENT <- factor(Thermal_switch_R_single_MAY$TREATMENT, levels = c("CCC", "WCC", "CWC", "WWC", "WWW", "CWW", "WCW", "CCW"))


Thermal_switch_R_single_MAY$csum <- as.numeric(Thermal_switch_R_single_MAY$csum)

library(ggplot2)

Thermal_switch_R_single_MAY%>% 
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT)) + scale_color_manual(values=c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59")) +
  theme_classic() 


install.packages("writexl")

write_xlsx(Thermal_switch_R_NEWtime02,"OneDrive - University of East Anglia/R stats/Thermal_switch_R_NEWtime02.xlsx")

na.omit(Thermal_switch_R_single_MAY)
#subset by  warm developed groups
myplotcsumwarm<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT),size = 2) + scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(52, 75) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 



myplotcsumwarm


#subset by cold developed groups
myplotcsumcold<-Thermal_switch_R_single_MAY %>%filter(DEVELOPMENTALT =="COLD")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT), size = 2) + scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(170, 180) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 



myplotcsumcold








#I left the following graphs here for further aesthetic trials later, however i think the previous two are the most explanatory ones.
myplot0<-ggplot(Thermal_switch_R_single_MAY, aes(TPF,csum,group = TREATMENT, colour=TREATMENT))+ geom_line() + geom_point()
myplot0

myplot1<- ggplot(Thermal_switch_R_single_MAY, aes(TPF,csum,color = TREATMENT))+
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width = 0.75)) + scale_fill_brewer(palette = "OrRd") + theme_bw()

myplot2<-myplot1  + scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(55, 66) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() + scale_color_manual(values=c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59"))

myplot2 + xlim(48, 180)
myplot2 + stat_summary(fun = mean,
                       fun.min = function(x) mean(x) - sd(x), 
                       fun.max = function(x) mean(x) + sd(x), 
                       geom = "SMOOTH") 

myplot2 + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x') + ylab("CUMULATIVE HATCH (%)") + xlim(48, 180) + ylim(0.1, 1) + theme(legend.position = "none")

myplot2 + facet_grid(DEVELOPMENTALT ~TREATMENT)

myplot2 + facet_grid(DEVELOPMENTALT ~ TREATMENT, scales='free')

myplot2 + facet_grid(.~ DEVELOPMENTALT)













library(magrittr)
Thermal_switch_R_single_MAY %<>%
  mutate(DEVELOPMENTALT = case_when(
    TREATMENT == "WWW" ~ "WARM",
    TREATMENT == "CWW" ~ "WARM",
    TREATMENT == "WCW" ~ "WARM",
    TREATMENT == "CCW" ~ "WARM",
    TREATMENT == "WWC" ~ "COLD",
    TREATMENT == "WCC" ~ "COLD",
    TREATMENT == "CWC" ~ "COLD",
    TREATMENT == "CCC" ~ "COLD"))


library(ggplot2)
library(dplyr)
library(Rmisc)
library(ggplot2)

#subset by  warm developed groups
myplotcsumwarm<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(55, 66) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 

myplotcsumwarm
myplotcsumwarm<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + stat_summary(aes(group = TREATMENT, fill = TREATMENT), 
                                                            fun.min = function(y) mean(y) - sd(y), 
                                                            fun.max = function(y) mean(y) + sd(y), 
                                                            geom = "ribbon", alpha = 0.1, na.rm = T, xmin = 0, xmax = Inf)+scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")

myplotcsumwarm + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",
                                                                                                                     fun.args = list(
                                                                                                                       mult = 1
                                                                                                                     ))+ xlab("Time post fertilization (days)") + xlim(57, 67) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 



myplotcsumwarm<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) +scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")

myplotcsumwarm + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + xlim(55, 65) + ylim(0, 1) + theme(legend.position = "none") +theme_classic()

myplotcsumwarmex<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) +scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")

myplotcsumwarmex1<-myplotcsumwarmex + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + xlim(55, 65) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() +geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
                                                                                                                                                                                                                                                                                                                                                                                                                 position = position_jitter(0.3), alpha = (0.4),
                                                                                                                                                                                                                                                                                                                                                                                                                 size = 1.5) +stat_summary(
                                                                                                                                                                                                                                                                                                                                                                                                                   aes(color = TREATMENT),
                                                                                                                                                                                                                                                                                                                                                                                                                   fun.data="mean_sdl",  fun.args = list(mult=1), 
                                                                                                                                                                                                                                                                                                                                                                                                                   geom = "pointrange",  size = 1
                                                                                                                                                                                                                                                                                                                                                                                                                 )  +   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) 


myplotcsumwarmex1



myplotcsumwarm + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1),alpha = 0.5)+ xlab("Time post fertilization (days)") + xlim(55, 65) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumwarm<-Thermal_switch_R_single_MAY%>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT),size = 1) + scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(52, 75) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 



myplotcsumwarm




myplotcsumcold<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="COLD")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) +scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)")

myplotcsumcold + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + xlim(170, 180) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumcold + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1),alpha = 0.5)+ xlab("Time post fertilization (days)") +xlim(170, 180) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumcold1<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="COLD")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) +scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)")

myplotcsumcold1 + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + xlim(170, 180) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumcold2<-myplotcsumcold1 + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1),alpha = 0.5)+ xlab("Time post fertilization (days)") +xlim(170, 180) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() +geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
             position = position_jitter(0.3), alpha = (0.4),
             size = 1.5)+   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) +stat_summary(
               aes(color = TREATMENT),
               fun.data="mean_sdl",  fun.args = list(mult=1), 
               geom = "pointrange",  size = 1
             ) 

myplotcsumcold2
myplotcsumcold<-Thermal_switch_R_single_MAY%>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="COLD")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT),size = 2) + scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") +xlim(170, 180) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 



myplotcsumcold

Thermal_switch_R_single_MAY$TREATMENT <- factor(Thermal_switch_R_single_MAY$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))


myplotcsumgen<-Thermal_switch_R_single_MAY %>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) +scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + ylab("Cumulative hatch (proportion)")

myplotcsumgen + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + xlim(57, 180) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumgen + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1),alpha = 0.5)+ xlab("Time post fertilization (days)") +xlim(57, 180) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumcold<-Thermal_switch_R_single_MAY%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT),size = 2) + scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") +xlim(170, 180) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 

myplotcsumgen


myplotcsumgen + facet_grid(. ~fct_rev(DEVELOPMENTALT) , scales = 'free_x')
library(dplyr)



#subset by cold developed groups
myplotcsumcold<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="COLD")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT), size = 2) + scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(170, 180) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 



myplotcsumcold

myplotcsumwarm<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + stat_summary(aes(group = TREATMENT, fill = TREATMENT), 
                                                            fun.min = function(y) mean(y) - sd(y), 
                                                            fun.max = function(y) mean(y) + sd(y), 
                                                            geom = "ribbon", alpha = 0.1, na.rm = T, xmin = 0, xmax = Inf)+scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")

myplotcsumcoldcorr<-Thermal_switch_R_single_MAY%>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="COLD")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + stat_summary(aes(group = TREATMENT, fill = TREATMENT), 
                                                            fun.min = function(y) mean(y) - sd(y), 
                                                            fun.max = function(y) mean(y) + sd(y), 
                                                            geom = "ribbon", alpha = 0.1, na.rm = T, xmin = 0, xmax = Inf)+scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)")

myplotcsumcoldcorr


myplotcsumcoldcorr + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1),alpha = 2)+ xlab("Time post fertilization (days)") + xlim(160, 190) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() + scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(165, 180) +  ylim(0, 1) + theme(legend.position = "none") +theme_classic() 





myplot1corr<- Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>% ggplot(aes(TPF,csum,color = TREATMENT))+
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width = 5)) + scale_fill_brewer(palette = "OrRd") + theme_bw()






myplot3corr<- Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="COLD")%>% ggplot(aes(TPF,csum,color = TREATMENT))+ 
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width =0)) +  scale_fill_brewer(palette = "OrRd") + theme_bw()
myplot3corr

myplot4corr<-myplot3corr + scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(160, 180) +  ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplot4corr 


myplot7corr<- Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>% ggplot(aes(TPF,csum,color = TREATMENT))+ 
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width =0)) +  scale_fill_brewer(palette = "OrRd") + theme_bw()
myplot7corr

myplot8corr<-myplot7corr + scale_color_manual(values=c("#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(40, 100) +  ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplot8corr 


myplotTPFcorr<- Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="COLD")%>% ggplot(aes(TPF,csum,color = TREATMENT))+ 
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width =0)) +  scale_fill_brewer(palette = "OrRd") + theme_bw()
myplot3corr

myplotTPF1corr<-myplotTPFcorr + scale_color_manual(values=c( "#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(160, 180) +  ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplot4corr 


myplot0<-ggplot(Thermal_switch_R_single_MAY, aes(TPF,csum,group = TREATMENT, colour=TREATMENT))+ geom_line() + geom_point()
myplot0

myplot1<- ggplot(Thermal_switch_R_single_MAY, aes(TPF,csum,color = TREATMENT))+
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width = 10.75)) + scale_fill_brewer(palette = "OrRd") + theme_bw()

myplot2<-myplot1 + scale_color_manual(values=c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59"))

myplot2 + xlim(48, 180)
myplot2 + stat_summary(fun = mean,
                       fun.min = function(x) mean(x) - sd(x), 
                       fun.max = function(x) mean(x) + sd(x), 
                       geom = "SMOOTH") 

myplot2 + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x') + ylab("CUMULATIVE HATCH (%)") + xlim(48, 180) + ylim(0, 1) + theme(legend.position = "none")

myplot2 + facet_grid(DEVELOPMENTALT ~TREATMENT)

myplot2 + facet_grid(DEVELOPMENTALT ~ TREATMENT, scales='free')

myplot2 + facet_grid(.~ DEVELOPMENTALT)



library(readxl)
Thermal_switch_R_NEWtime01 <- read_excel("Thermal_switch_R_NEWtime01.xlsx")
View(Thermal_switch_R_NEWtime01)









library(magrittr)
Thermal_switch_R_NEWtime01 %<>%
  mutate(DEVELOPMENTALT = case_when(
    TREATMENT == "WWW" ~ "WARM",
    TREATMENT == "CWW" ~ "WARM",
    TREATMENT == "WCW" ~ "WARM",
    TREATMENT == "CCW" ~ "WARM",
    TREATMENT == "WWC" ~ "COLD",
    TREATMENT == "WCC" ~ "COLD",
    TREATMENT == "CWC" ~ "COLD",
    TREATMENT == "CCC" ~ "COLD"))

#check the treatment corresponding colors!!!!!
#Plot cumulative events
library(RColorBrewer)

Thermal_switch_R_NEWtime0_reordered <- Thermal_switch_R_NEWtime01      #Replicate data
Thermal_switch_R_NEWtime0_reordered$TREATMENT <- factor(Thermal_switch_R_NEWtime0_reordered$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))

Thermal_switch_R_single_MAY$TREATMENT <- factor(Thermal_switch_R_single_MAY$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))

myplotcsumwarm<-Thermal_switch_R_single_MAY %>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) +scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")  +scale_fill_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)")

myplotcsumwarm + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + xlim(55, 65) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumwarm + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1),alpha = 0.5)+ xlab("Time post fertilization (days)") + xlim(55, 65) + ylim(0, 1) + theme(legend.position = "none") +theme_classic() 

myplotcsumwarm<-Thermal_switch_R_single_MAY%>%filter(Thermal_switch_R_single_MAY$DEVELOPMENTALT =="WARM")%>%
  ggplot(aes(x=TPF,y=csum,colour=TREATMENT)) + 
  geom_smooth(aes(group=TREATMENT),size = 1) + scale_color_manual(values=c( "#990000","#d7301f", "#ef6548", "#fc8d59")) + ylab("Cumulative hatch (proportion)") + xlab("Time post fertilization (days)") + xlim(52, 75) + ylim(0, 0.65) + theme(legend.position = "none") +theme_classic() 



myplotcsumwarm

myplotfryper<-ggplot(Thermal_switch_R_NEWtime0_reordered, aes(TPF,FRYPOP,group = TREATMENT, colour=TREATMENT))
myplotfryper + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + ylim(0, 0.5) + theme(legend.position = "none")  +theme_bw()

myplot1fryper<- ggplot(Thermal_switch_R_NEWtime0_reordered, aes(TPF,FRYPOP,color = TREATMENT)) + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + ylim(0, 0.7) + theme(legend.position = "none") +theme_bw() + geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                position = position_jitter(0.3), alpha = (0.4),
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                size = 1.5) +stat_summary(
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  aes(color = TREATMENT),
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  fun.data="mean_sdl",  fun.args = list(mult=1), 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  geom = "pointrange",  size = 1
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                )  +   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) 

myplot2fryper<-myplot1fryper + scale_color_manual(values=c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e" ))+  ylab("12 Hrs Hatching success (proportion)") +xlab("Time post fertilization (Days)")
myplot2fryper

##############
myplotfryper1<-ggplot(Thermal_switch_R_NEWtime0_reordered, aes(TPF,FRYPOP,group = TREATMENT, colour=TREATMENT))
myplotfryper2<-myplotfryper1 + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + ylim(0, 0.5) + theme(legend.position = "none")  +theme_bw()

myplot1fryper1<- ggplot(Thermal_switch_R_NEWtime0_reordered, aes(TPF,FRYPOP,color = TREATMENT)) + stat_summary(fun = mean, geom="line", size=1, position = position_dodge(width = 0))+   stat_summary(fun.data = "mean_sdl",                                                                                                     fun.args = list(mult = 1))+ xlab("Time post fertilization (days)") + ylim(0, 0.7) + theme(legend.position = "none") +theme_bw() 

myplot1fryper1 

myplot2fryper<-myplot1fryper1 + scale_color_manual(values=c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e" ))+ ylab("12 Hrs Hatching success (proportion)") +xlab("Time post fertilization (Days)")
myplot3fryper<-myplot2fryper + geom_jitter(aes( color = TREATMENT, shape= factor(FAMILY)), 
                          position = position_jitter(0.3), alpha = (0.3),
                          size = 2)+   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10))

myplot3fryper

############################

Thermal_switch_R_NEWtime0_reordered$EGG<-factor(Thermal_switch_R_NEWtime0_reordered$EGG, levels= c("W", "C"))
Thermal_switch_R_NEWtime0_reordered$FAMILY<-factor(Thermal_switch_R_NEWtime0_reordered$FAMILY, levels= c("1", "2","3", "4","5", "6","7", "8","9", "10"))
Thermal_switch_R_NEWtime0_reordered$SPERM<-factor(Thermal_switch_R_NEWtime0_reordered$SPERM, levels= c("W", "C"))
Thermal_switch_R_NEWtime0_reordered$DEV<-factor(Thermal_switch_R_NEWtime0_reordered$DEV, levels= c("W", "C"))
myplot2fryper + facet_grid(. ~DEVELOPMENTALT) + ylab("12 Hrs Hatching success (Proportion)") +xlab("Time post fertilization (Days)")

cmpd <- list(c("WWW", "CWW"), c("CWW", "WCW"), c("WCW", "CCW"), c("CCW", "WWC"), c("WWC", "WCC"), c("WCC", "CWC"),c("CWC", "CCC"), c("CWC", "WWC"), c("WWC", "WCW"), c("WCC", "WCW"),  c("WCW", "WWW"),c("CCC", "WWW") )

myplot3fryper + facet_grid(DEVELOPMENTALT ~TREATMENT)

myplot3fryper + facet_grid(DEVELOPMENTALT ~ TREATMENT, scales='free_x')


myplot3fryper + facet_grid(EGG ~ SPERM  ~ DEV, scales='free_x')
myplot3fryper + facet_grid(.~ DEVELOPMENTALT)

myplot3fryper + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x') + ylab("12 Hrs Hatching success (proportion)") +xlab("Time post fertilization (Days)")

# Split in vertical direction
myplot2 + facet_grid(DEVELOPMENTALT ~ .) + ylab("12 Hrs Hatching success (proportion)")
myplot2fryper + facet_grid(DEVELOPMENTALT ~ .)+ ylab("12 Hrs Hatching success (proportion)")

cmpd <- list(c("WWW", "CWW"), c("CWW", "WCW"), c("WCW", "CCW"), c("CCW", "WWC"), c("WWC", "WCC"), c("WCC", "CWC"),c("CWC", "CCC"), c("CWC", "WWC"), c("WWC", "WCW"), c("WCC", "WCW"),  c("WCW", "WWW"),c("CCC", "WWW") )
myplot1 + stat_compare_means(comparisons = cmpd, tip.length=0.01,
                             label = "p.signif", 
                             symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), 
                                                symbols = c("****", "***", "**", "*", "ns")))



myplotfrypercold<-subset(Thermal_switch_R_NEWtime0_reordered,TREATMENT=="CCC"|TREATMENT=="WCC"|TREATMENT== "CWC"|TREATMENT== "WWC")

myplotfrypercold1<-ggplot(myplotfrypercold, aes(TPF,FRYPOP,group = TREATMENT, colour=TREATMENT))+ geom_line() + geom_point()
myplotfrypercold1

myplotfrypercold1<- ggplot(subset(Thermal_switch_R_NEWtime0_reordered,TREATMENT=="CCC"|TREATMENT=="WCC"|TREATMENT== "CWC"|TREATMENT== "WWC"), aes(TPF,FRYPOP,color = TREATMENT))+
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width = 0)) + scale_fill_brewer(palette = "OrRd") + theme_bw()

myplotfrypercold2<-myplotfrypercold1 + scale_color_manual(values=c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59"))
myplotfrypercold3<-myplotfrypercold2 + ylab("12 Hrs Hatching success (proportion)")+ xlab("Time post fertilization (Days)") + xlim(160, 180) + ylim(0, 1) + theme(legend.position = "none")
myplotfrypercold3

myplotfryperwarm<-subset(Thermal_switch_R_NEWtime0_reordered,TREATMENT=="WWW"|TREATMENT=="CWW"|TREATMENT== "WCW"|TREATMENT== "CCW")

myplotfryperwarm1<-ggplot(myplotfryperwarm, aes(TPF,FRYPOP,group = TREATMENT, colour=TREATMENT))+ geom_line() + geom_point()
myplotfryperwarm1

myplotfryperwarm1<- ggplot(subset(Thermal_switch_R_NEWtime0_reordered,TREATMENT=="WWW"|TREATMENT=="CWW"|TREATMENT== "WCW"|TREATMENT== "CCW"), aes(TPF,FRYPOP,color = TREATMENT))+
  stat_summary(fun= "mean", geom="line", size=1, position = position_dodge(width = 0)) + scale_fill_brewer(palette = "OrRd") + theme_bw()

myplotfryperwarm2<-myplotfryperwarm1 + scale_color_manual(values=c("#990000","#d7301f", "#ef6548", "#fc8d59"))
myplotfryperwarm3<-myplotfryperwarm2 + ylab("12 Hrs Hatching success (proportion)")+ xlab("Time post fertilization (Days)") + xlim(55, 65) + ylim(0, 1) + theme(legend.position = "none")
myplotfryperwarm3

blankPlot <- ggplot()+geom_blank(aes(1,1))+
  theme(
    plot.background = element_blank(), 
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(), 
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(), 
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.line = element_blank()
  )


library("gridExtra")
grid.arrange(myplot2fryper, blankPlot  , myplotfryperwarm3, myplotfrypercold3, 
             ncol=2, nrow=2, widths=c(4, 4), heights=c(1.4, 4))

grid.arrange(myplot2fryper, arrangeGrob(myplotfryperwarm3,myplotfrypercold3, ncol=2), heights=c(1.5/4, 2.5/4), ncol=1)


library(patchwork)
library("gridExtra")

library(tidyverse)
library(Rmisc)
library(dplyr)
library(ggplot2)

NEWThermal<- 
  Thermal_switch_R_NEWtime01 %>%
  group_by(TREATMENT, TPF) %>%
  dplyr::summarise(mean_csum = mean(csum), 
                   uci_csum = CI(csum)[1], 
                   lci_csum = CI(csum)[3]) %>%
  mutate(Time = TPF %>% as.factor())

#change color codes to these ones cause it might be important to describe differences in total hatch
z<- ggplot(NEWThermal, aes(TPF,mean_csum,color = TREATMENT))+ geom_line(position = position_dodge(width = 10)) + theme_bw() 
stat_summary(fun= "mean", geom="line", size=4.5) 

z + scale_color_manual(values=c("#634a25","#ae8b41", "#fdd49e","#fdbb84", "#990000","#d7301f", "#ef6548", "#fc8d59"))

l<-z+geom_ribbon(aes(ymin=NEWThermal$lci_csum, ymax=NEWThermal$uci_csum), linetype=2, alpha=0.05) + geom_line()
l  + scale_fill_brewer(palette= "OrRd")



#modelling total percentage of fry (FRYPERTOT) across time and treatment + interactions
#create new dataframe with summarised FRYPERTOT info
install.packages("dplyr") # Install dplyr package
library("dplyr")   # Load dplyr package
str(Thermalsingle)

Thermal_frypertotmean <- Thermalsingle %>% group_by(FRYTOTPROP, TREATMENT)
view(Thermal_frypertotmean)
Thermal_frypertotmean_TREAT_FRYPERTOT <- Thermal_frypertotmean %>% summarise(n = n())
view(Thermal_frypertotmean_TREAT_FRYPERTOT)

Thermal_frypertotmean1<-Thermalsingle%>%
  group_by(FAMILY, TREATMENT)%>%
  slice(1)%>%
  ungroup()

view(Thermal_frypertotmean1)

#cbind()


Thermal_frypertotmean<-Thermalsingle %>%         # Specify data frame
  group_by(EGG,SPERM, DEV, FAMILY) %>%                                                 # Specify group indicator
  summarise_at(vars(FRYPERTOT),              # Specify column
               list(FRYPERTOT = mean))               # Specify function

# linear mixed effect model for Total number of fry hatched per group.
library(lme4)
library(glmmTMB)
M6<-glmmTMB(cbind(TOTFRY, eggcount) ~ EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"), data = Thermal_frypertotmean1)
summary(M6)
es<-summary(M6)
coef(summary(M6))
Anova(M6,type="III")
joint_tests(M6)

s6 <- simulateResiduals(M6)
plot(s6)
testDispersion(s6)


M6<-glmer(cbind(TOTFRY, eggcount) ~ EGG*SPERM*DEV + (1|MALE), family=binomial(link="logit"), data = Thermal_frypertotmean1)
summary(M6)
es<-summary(M6)
coef(summary(M6))
anova(M6)
joint_tests(M6)

s6 <- simulateResiduals(M6)
plot(s6)
testDispersion(s6)


#Estimate the modelled means...also backtransformed to the response variable units from log scale
# by using type=response
library(emmeans)
m6.emm <- emmeans(M6, specs="EGG", type="response")
m6.emm

m6.emm_dt<-data.frame(m6.emm)

#contrasts pairwise differences
m6.pairs <- contrast(m6.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m6.pairs


m6.emm <- emmeans(M6, specs="SPERM", type="response")
m6.emm

m6.emm_dt<-data.frame(m6.emm)

#contrasts pairwise differences
m6.pairs <- contrast(m6.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m6.pairs


m6.emm <- emmeans(M6, specs="DEV", type="response")
m6.emm

m6.emm_dt<-data.frame(m6.emm)

#contrasts pairwise differences
m6.pairs <- contrast(m6.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m6.pairs



write.table(m6.emm, "emmeans.txt")
write.table(m6.pairs, "pairs.txt")
write.table(es, "summarymodel.txt")

#contrasts pairwise differences
m6.pairs <- contrast(m6.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m6.pairs

#Contrast test for trends due to presence of significant interactions

ESD.emm <- emmeans(M6, ~ EGG * SPERM * DEV)
ESD.emm

pairs(ESD.emm, simple = list("EGG", c("SPERM", "DEV")))

pairs(ESD.emm, simple = "each")

emm_s.t <- emmeans(M6, pairwise ~ EGG | DEV)

## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)


emm_s.t <- emmeans(M6, pairwise ~ SPERM | DEV)
## NOTE: Results may be misleading due to involvement in interactions
emm_s.t
contrast(emm_s.t[[1]], "poly")   
IC_st <- contrast(regrid(emm_s.t[[1]]), interaction = c("poly", "consec"), by = NULL)
IC_st
coef(IC_st)

test(IC_st, joint = TRUE)

#The three-way interaction may be explored via interaction contrasts too:

contrast(emmeans(M6, ~ EGG*SPERM*DEV),
         interaction = c("poly", "consec", "consec"))










M6<-lmer(FRYTOTPROP ~ EGG*SPERM*DEV + (1|FAMILY), data = Thermal_frypertotmean1)
summary(M6)
es<-summary(M6)
coef(summary(M6))

s6 <- simulateResiduals(M6)
plot(s6)
testDispersion(s6)

anova(M6)
#Estimate the modelled means...also backtransformed to the response variable units from log scale
# by using type=response
library(emmeans)
m6.emm <- emmeans(M6, specs="EGG", type="response")
m6.emm

m6.emm_dt<-data.frame(m6.emm)

#contrasts pairwise differences
m6.pairs <- contrast(m6.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m6.pairs


m6.emm <- emmeans(M6, specs="SPERM", type="response")
m6.emm

m6.emm_dt<-data.frame(m6.emm)

#contrasts pairwise differences
m6.pairs <- contrast(m6.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m6.pairs

#simulate the residuals
library(DHARMa)
output<-simulateResiduals(M6)
plot(M6)

par(mfrow = c(2, 2))
plot(M6)

anova(M6)

joint_tests(M6)

write.table(m6.emm, "emmeans.txt")
write.table(m6.pairs, "pairs.txt")
write.table(es, "summarymodel.txt")



m7<-lmer(FRYTOTPROP ~ EGG*SPERM + (1|FAMILY), data = Thermal_frypertotmean1%>%filter(DEV=="W"))
summary(m7)

coef(summary(m7))

s7 <- simulateResiduals(m7)
plot(s7)
testDispersion(s7)


#Estimate the modelled means...also backtransformed to the response variable units from log scale
# by using type=response
library(emmeans)
m7.emm <- emmeans(m7, specs="EGG", type="response")
m7.emm

m7.emm_dt<-data.frame(m7.emm)

#contrasts pairwise differences
m7.pairs <- contrast(m7.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m7.pairs


m7.emm <- emmeans(m7, specs="SPERM", type="response")
m7.emm

m7.emm_dt<-data.frame(m7.emm)

#contrasts pairwise differences
m7.pairs <- contrast(m7.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m7.pairs

anova(m7)
joint_tests(m7)

#bello proprio
m9<-lmer(FRYTOTPROP ~ EGG*SPERM + (1|FAMILY), data = Thermal_frypertotmean%>%filter(DEV=="C"))

summary(m9)

m9.emm <- emmeans(m9, specs="EGG", type="response")
m9.emm

m9.emm_dt<-data.frame(m7.emm)

m9.pairs <- contrast(m9.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m9.pairs

m9.emm <- emmeans(m9, specs="SPERM", type="response")
m9.emm

m9.emm_dt<-data.frame(m7.emm)

m9.pairs <- contrast(m9.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m9.pairs

anova(m9)
joint_tests(m9)


#boxplot differences in hatching success
# Basic boxplot

boxplot(Thermal_frypertotmean1$FRYPERTOT ~ Thermal_frypertotmean1$TREATMENT , col=terrain.colors(8) )

Thermal_frypertotmean1 %<>%
  mutate(DEVELOPMENTALT = case_when(
    TREATMENT == "WWW" ~ "WARM",
    TREATMENT == "CWW" ~ "WARM",
    TREATMENT == "WCW" ~ "WARM",
    TREATMENT == "CCW" ~ "WARM",
    TREATMENT == "WWC" ~ "COLD",
    TREATMENT == "WCC" ~ "COLD",
    TREATMENT == "CWC" ~ "COLD",
    TREATMENT == "CCC" ~ "COLD"))


library(ggplot2)
library(ggpubr)

Thermal_frypertotmean_reordered <- Thermalfrypertotmean_1      #Replicate data
Thermal_frypertotmean_reordered$TREATMENT <- factor(Thermal_frypertotmean_reordered$TREATMENT, levels = c("CCC", "WCC", "CWC", "WWC", "WWW", "CWW", "WCW", "CCW"))


myplot <- ggplot(Thermal_frypertotmean_reordered, aes(x =(TREATMENT), y = FRYTOTPROP, fill=TREATMENT)) + geom_boxplot() + ggtitle("Hatching success") + theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) + scale_y_continuous(name = bquote ("Hatching success (Proportion)")) + scale_x_discrete(name = "", limits = c("CCC", "WCC", "CWC", "WWC", "WWW", "CWW", "WCW", "CCW")) + scale_fill_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0) + coord_cartesian(y=c(0,1))

Thermalsingle1 %<>%
  mutate(DEVELOPMENTALT = case_when(
    TREATMENT == "WWW" ~ "WARM",
    TREATMENT == "CWW" ~ "WARM",
    TREATMENT == "WCW" ~ "WARM",
    TREATMENT == "CCW" ~ "WARM",
    TREATMENT == "WWC" ~ "COLD",
    TREATMENT == "WCC" ~ "COLD",
    TREATMENT == "CWC" ~ "COLD",
    TREATMENT == "CCC" ~ "COLD"))

Thermalsingle1$TREATMENT <- factor(Thermalsingle1$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))

violini<-ggplot(Thermalsinglefin, aes(x = TREATMENT, y = TPF, fill = TREATMENT)) +
  geom_violin(trim = FALSE) + 
  geom_boxplot(width = 0.07) +
  guides(fill = guide_legend(title = "Title")) + scale_y_continuous(name = bquote ("Age at hatching")) + scale_x_discrete(name = "", limits = c("CCC", "WCC", "CWC", "WWC", "WWW", "CWW", "WCW", "CCW")) + scale_fill_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0)
violini

library(dplyr)

Thermalsingle1$TREATMENT <- factor(Thermalsingle$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))


violini2<-Thermalsinglefin %>%filter(DEVELOPMENTALT =="WARM")%>%ggplot(aes(x = TREATMENT, y = TPF, fill = TREATMENT)) +
  geom_violin(trim = FALSE) +
  guides(fill = guide_legend(title = "Title")) + scale_y_continuous(name = bquote ("Age at hatching")) + scale_x_discrete(name = "", limits = c( "WWW", "CWW", "WCW", "CCW")) + scale_fill_manual(values = c( "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0)  + ylim(55, 65)

#violini2 + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x')
violini2 



violini3<-Thermalsinglefin%>%filter(DEVELOPMENTALT =="COLD")%>%ggplot(aes(x = TREATMENT, y = TPF, fill = TREATMENT)) +
  geom_violin(trim = FALSE)  +
  guides(fill = guide_legend(title = "Title")) + scale_y_continuous(name = bquote ("Age at hatching")) + scale_x_discrete(name = "", limits = c( "CCC", "WCC", "CWC", "WWC")) + scale_fill_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0)  + ylim(150,180)

#violini2 + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x')
violini3



violini4<-ggplot(Thermalsinglefin, aes(x = TREATMENT, y = DDAY, fill = TREATMENT)) +
  geom_violin(trim = FALSE) + 
  geom_boxplot(width = 0.07) +
  guides(fill = guide_legend(title = "Title")) + scale_y_continuous(name = bquote ("Age at hatching")) + scale_x_discrete(name = "", limits = c("CCC", "WCC", "CWC", "WWC", "WWW", "CWW", "WCW", "CCW")) + scale_fill_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0)
violini4

library(dplyr)

Thermalsingle1$TREATMENT <- factor(Thermalsingle$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))

Thermalsinglefin$TREATMENT <- factor(Thermalsinglefin$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))



violini5<-Thermalsinglefin %>%filter(DEVELOPMENTALT =="WARM")%>%ggplot(aes(x = TREATMENT, y = DDAY, fill = TREATMENT)) +
  geom_violin(trim = FALSE) +
  guides(fill = guide_legend(title = "Title")) + scale_y_continuous(name = bquote ("Age at hatching")) + scale_x_discrete(name = "", limits = c( "WWW", "CWW", "WCW", "CCW")) + scale_fill_manual(values = c( "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0)  + ylim(450, 500)

#violini2 + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x')
violini5



violini6<-Thermalsinglefin%>%filter(DEVELOPMENTALT =="COLD")%>%ggplot(aes(x = TREATMENT, y = DDAY, fill = TREATMENT)) +
  geom_violin(trim = FALSE)  +
  guides(fill = guide_legend(title = "Title")) + scale_y_continuous(name = bquote ("Age at hatching")) + scale_x_discrete(name = "", limits = c( "CCC", "WCC", "CWC", "WWC")) + scale_fill_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0)  + ylim(300,360)

#violini2 + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x')
violini6


basicggplot<- Thermalsingle1 %>%filter(DEVELOPMENTALT =="WARM")%>%ggplot(aes(x=TREATMENT, y=TPF, fill=TREATMENT)) 


basicggplot1<-basicggplot +  scale_x_discrete(name = "", limits = c("WWW", "CWW", "WCW", "CCW")) +theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) +  ggtitle("")  + scale_colour_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59", "#d7301f" )) + coord_cartesian(y=c(57,62))+  theme(text = element_text(size = 15)) + scale_y_continuous(name = bquote ("Age at hatch (days)")) 



+geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
             position = position_jitter(0.3), alpha = (0.3),
             size = 2)+   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10))

basicggplot2<-basicggplot1 +grids(linetype = "dashed")+ geom_hline(yintercept=59) + geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
             position = position_jitter(0.3), alpha = (0.3),
             size = 2)+   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) +stat_summary(
                                                                                                aes(color = TREATMENT),
                                                                                                fun.data="mean_sdl",  fun.args = list(mult=1), 
                                                                                                geom = "pointrange",  size = 1.5
                                                                                              )  

basicggplot2

basicggplot3<- Thermalsingle1 %>%filter(DEVELOPMENTALT =="COLD")%>%ggplot(aes(x=TREATMENT, y=TPF, fill=TREATMENT)) 

basicggplot4<-basicggplot3 +  scale_x_discrete(name = "", limits = c("CCC", "WCC", "CWC", "WWC")) +theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) +  ggtitle("")  + scale_colour_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + coord_cartesian(y=c(165,180))+  theme(text = element_text(size = 15)) + scale_y_continuous(name = bquote ("Age at hatch (days)")) 

basicggplot5<-basicggplot4+grids(linetype = "dashed")+ geom_hline(yintercept=174) + geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
                                                                                                position = position_jitter(0.3), alpha = (0.3),
                                                                                                size = 2)+   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) +stat_summary(
                                                                                                aes(color = TREATMENT),
                                                                                                fun.data="mean_sdl",  fun.args = list(mult=1), 
                                                                                                geom = "pointrange",  size = 1.5
                                                                                              )  

basicggplot5



basicggplot6<- Thermalsinglefin%>%filter(DEVELOPMENTALT =="WARM")%>%ggplot(aes(x=TREATMENT, y=DDAY, fill=TREATMENT)) 


basicggplot7<-basicggplot6 +  scale_x_discrete(name = "", limits = c("WWW", "CWW", "WCW", "CCW")) +theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) +  ggtitle("")  + scale_colour_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59", "#d7301f" )) + coord_cartesian(y=c(460,500))+  theme(text = element_text(size = 15)) + scale_y_continuous (name = bquote ("Degree-days after fertilisation")) 

basicggplot8<-basicggplot7 +grids(linetype = "dashed")+ geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
                                                                    position = position_jitter(0.3), alpha = (0.3),
                                                                    size = 2)+   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10))+ geom_hline(yintercept=472)+geom_jitter(aes( color = TREATMENT), 
                                                                                               position = position_jitter(0.3), alpha = (0.1),
                                                                                               size = 2) +stat_summary(
                                                                                                 aes(color = TREATMENT),
                                                                                                 fun.data="mean_sdl",  fun.args = list(mult=1), 
                                                                                                 geom = "pointrange",  size = 1.5
                                                                                               )  

basicggplot8

basicggplot9<- Thermalsinglefin %>%filter(DEVELOPMENTALT =="COLD")%>%ggplot(aes(x=TREATMENT, y=DDAY, fill=TREATMENT)) 

basicggplot10<-basicggplot9 +  scale_x_discrete(name = "", limits = c("CCC", "WCC", "CWC", "WWC")) +theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) +  ggtitle("")  + scale_colour_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + coord_cartesian(y=c(320,360))+  theme(text = element_text(size = 15)) + scale_y_continuous(name = bquote ("Degree-days after fertilisation")) 

basicggplot11<-basicggplot10+grids(linetype = "dashed")+ geom_hline(yintercept=347)+geom_jitter(aes( color = TREATMENT), 
                                                                                                position = position_jitter(0.3), alpha = (0.1),
                                                                                                size = 2) +stat_summary(
                                                                                                  aes(color = TREATMENT),
                                                                                                  fun.data="mean_sdl",  fun.args = list(mult=1), 
                                                                                                  geom = "pointrange",  size = 1.5
                                                                                                )  

basicggplot11



library(sjmisc)
library(sjPlot)

# simple forest plot
plot_model(M1)



hatchfinal2 + facet_grid(. ~DEVELOPMENTALT, scales = 'free_x')


#634a25","#ae8b41", "#D1BD94", "#fdd49e",
#scale_fill_manual
myplot


Thermal_frypertotmean_reordered$TREATMENT <- factor(Thermal_frypertotmean_reordered$TREATMENT, levels = c("WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC" ))
Thermal_frypertotmean_reordered$FAMILY <- factor(Thermal_frypertotmean_reordered$FAMILY, levels = c("1", "2", "3", "4","5", "6", "7", "8", "9", "10"))

Thermal_frypertotmean_reordered$TREATMENT<- as.factor(Thermal_frypertotmean_reordered$TREATMENT)
hatchfinal<-  ggplot(Thermal_frypertotmean_reordered, aes(x=TREATMENT, y=FRYTOTPROP, fill=TREATMENT)) 

hatchfinal1<-hatchfinal  +theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) +  ggtitle("Hatching success")  + scale_colour_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + coord_cartesian(y=c(0,1))+  theme(text = element_text(size = 15)) + scale_y_continuous(name = bquote ("Hatching success (Proportion)")) 

hatchfinal2<-hatchfinal1 +grids(linetype = "dashed")+ geom_hline(yintercept=0.5)+ geom_jitter(aes( color = TREATMENT, shape=FAMILY), 
                                                                                              position = position_jitter(0.3), alpha = (0.3),
                                                                                              size = 2)+   scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) + scale_fill_discrete(breaks=c('1','2','3','4', '5', '6', '7', '8', '9', '10'))+ stat_summary(
                                                                                               aes(color = TREATMENT),
                                                                                               fun.data="mean_sdl",  fun.args = list(mult=1), 
                                                                                               geom = "pointrange",  size = 1
                                                                                             )  

hatchfinal2

Thermal_frypertotmean_reordered
hatchfinal2 + facet_grid(. ~fct_rev(DEVELOPMENTALT) , scales = 'free_x')

compare_means(FRYTOTPROP ~ TREATMENT, data = Thermal_frypertotmean_reordered)
m

myplot4 <- ggplot(Thermal_frypertotmean_reordered, aes(x =(TREATMENT), y = FRYTOTPROP, fill=TREATMENT)) + geom_boxplot() + ggtitle("Hatching success  (Proportion)") + theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) + scale_y_continuous(name = bquote ("Hatching success (Proportion)"))  + scale_fill_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0) + coord_cartesian(y=c(0,1))
myplot
z

myplot4 <-  ggplot(Thermal_frypertotmean_reordered, aes(x=TREATMENT, y=FRYTOTPROP, fill=TREATMENT))
myplot4
myplot5<-myplot4+  scale_x_discrete(name = "", limits = c("CCC", "WCC", "CWC", "WWC", "WWW", "CWW", "WCW", "CCW")) +theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) +  ggtitle("Hatching success")  + scale_colour_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + coord_cartesian(y=c(0,1))+  theme(text = element_text(size = 15)) + scale_y_continuous(name = bquote ("Hatching success (Proportion)")) 

myplot5

myplot6<-myplot5  + facet_grid(. ~DEVELOPMENTALT, scale="free") 
myplot6

cmpr <- list( c("CCC", "CWC"), c("CCC", "WWC"), c("WWW", "CWW"), c("WWW", "WCW"), c("WWW", "CCW"))
myplot 

myplot6 <- ggplot(subset(Thermal_frypertotmean_reordered,TREATMENT=="CCC"|TREATMENT=="WCC"|TREATMENT== "CWC"|TREATMENT== "WWC"), aes(x =(TREATMENT), y = FRYTOTPROP, fill=TREATMENT)) + geom_boxplot() + ggtitle("Hatching success (Proportion)") + theme(text=element_text(size = 16), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(lineheight=.8, face="bold", hjust=0.5)) + scale_y_continuous(name = bquote ("Hatching success (Proportion)"))  + scale_fill_manual(values = c("#634a25","#ae8b41", "#D1BD94", "#fdd49e", "#990000","#d7301f", "#ef6548", "#fc8d59" )) + geom_jitter(width = 0) + coord_cartesian(y=c(0,1))

myplot7<-myplot6  + facet_grid(. ~DEVELOPMENTALT, scale="free") 
myplot6




library(readxl)
Thermalfrypertotmean_1 <- read_excel("Thermalfrypertotmean_1.xlsx")
View(Thermalfrypertotmean_1)
View(Thermalfrypertotmean_1)
Thermalfrypertotmean_2<-Thermalfrypertotmean_1


library(tidyverse)
library(Rmisc)
library(dplyr)
library(magrittr)

Thermalfrypertotmean_2 %<>%
  mutate(DEVELOPMENTALT = case_when(
    TREATMENT == "WWW" ~ "WARM",
    TREATMENT == "CWW" ~ "WARM",
    TREATMENT == "WCW" ~ "WARM",
    TREATMENT == "CCW" ~ "WARM",
    TREATMENT == "WWC" ~ "COLD",
    TREATMENT == "WCC" ~ "COLD",
    TREATMENT == "CWC" ~ "COLD",
    TREATMENT == "CCC" ~ "COLD"))


library(ggplot2)
library(ggpubr)

Thermalfrypertotmean_3 <- Thermalfrypertotmean_2      #Replicate data
Thermalfrypertotmean_3$TREATMENT <- factor(Thermalfrypertotmean_3$TREATMENT, levels = c("CCC", "WCC", "CWC", "WWC", "WWW", "CWW", "WCW", "CCW"))

str(Thermalfrypertotmean_3)

Thermalfrypertotmean_3$EGG<-as.factor(Thermalfrypertotmean_3$EGG)
Thermalfrypertotmean_3$SPERM<-as.factor(Thermalfrypertotmean_3$SPERM)
Thermalfrypertotmean_3$DEV<-as.factor(Thermalfrypertotmean_3$DEV)



#Modelling proportion of malformed individuals (Cyclops) among the different thermal treatments

m12<-glmer(cbind(malfnoeye, eggcount) ~ EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"),control = glmerControl(optimizer ="Nelder_Mead"), data = Thermalfrypertotmean_3)

m12<-glmmTMB(cbind(malfnoeye, eggcount) ~ EGG*SPERM*DEV + (1|FAMILY),family=binomial(link="logit"), data = Thermalfrypertotmean_3)

m12<-lmer(malfnoeyeprop ~ EGG*SPERM*DEV + (1|FAMILY), data = Thermalfrypertotmean_3)


summary(m12)
anova(m12)
Anova(m12,type="III")
joint_tests(m12)


library(emmeans)

m12.emm <- emmeans(m12, specs="EGG", type="response")
m12.emm

m12.emm_dt<-data.frame(m12.emm)

m12.pairs <- contrast(m12.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m12.pairs




emm_s.t <- emmeans(m12, pairwise ~ DEV)

emm_s.t


m12.emm <- emmeans(m12, specs="DEV", type="response")
m12.emm

m12.emm_dt<-data.frame(m12.emm)

m12.pairs <- contrast(m12.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m12.pairs




emm_s.t <- emmeans(m12, pairwise ~ DEV)

emm_s.t




#Modelling proportion of individuals dead on hatch among the different thermal treatments

m13<-glmer(cbind(deadonhatch, eggcount) ~  EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"), data = Thermalfrypertotmean_3)


m13<-glmmTMB(deadonhatchprop ~ EGG*SPERM*DEV + (1|FAMILY), data = Thermalfrypertotmean_3)


summary(m13)
anova(m13)
Anova(m13,type="III")
joint_tests(m13)

m13.emm <- emmeans(m13, specs="EGG", type="response")
m13.emm

m13.emm_dt<-data.frame(m13.emm)

m13.pairs <- contrast(m13.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m13.pairs

m13.emm <- emmeans(m13, specs="SPERM", type="response")
m13.emm

m13.emm_dt<-data.frame(m13.emm)

m13.pairs <- contrast(m13.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m13.pairs


m13.emm <- emmeans(m13, specs="DEV", type="response")
m13.emm

m13.emm_dt<-data.frame(m13.emm)

m13.pairs <- contrast(m13.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m13.pairs




emm_s.t <- emmeans(m13, pairwise ~ TREATMENT)

emm_s.t


deaddev<-lmer(deadonhatchprop ~  EGG*SPERM*DEV +(1|FAMILY), data = Thermalfrypertotmean_3)


summary(deaddev)


deaddev.emm <- emmeans(deaddev, specs="DEV", type="response")
deaddev.emm

deaddev.emm_dt<-data.frame(deaddev)

deaddev.mpairs <- contrast(deaddev.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
deaddev.mpairs

anova(deaddev)
joint_tests(deaddev)


deaddevemm_s.t <- emmeans(deaddev, pairwise ~ DEV)

deaddevemm_s.t



#Modelling proportion of individuals dead after having reached the eyed stage among the different thermal treatments
m14<-glmer(cbind(eyeddead, eggcount) ~  EGG*SPERM*DEV + (1|FAMILY), family=binomial(link="logit"), data = Thermalfrypertotmean_3)


m14<-lmer(eyeddeadprop ~ EGG*SPERM*DEV + (1|FAMILY), data = Thermalfrypertotmean_3)


summary(m14)


m14.emm <- emmeans(m14, specs="DEV", type="response")
m14.emm

m14.emm_dt<-data.frame(m14.emm)

m14.pairs <- contrast(m14.emm, method="revpairwise", adjust="none") %>%
  summary(infer=c(TRUE, TRUE))
m14.pairs

anova(m14)
joint_tests(m14)


emm_s.t <- emmeans(m14, pairwise ~ DEV)

emm_s.t


#+++++++++++++++++++++++++
# Function to calculate the mean and the standard deviation
# for each group
#+++++++++++++++++++++++++
# data : a data frame
# varname : the name of a column containing the variable
#to be summariezed
# groupnames : vector of column names to be used as
# grouping variables


library(dplyr)
Thermalfrypertotmean_3mean<-Thermalfrypertotmean_3 %>% group_by(TREATMENT) %>% summarise_each(funs(mean, sd))

view(Thermalfrypertotmean_3mean)


Thermalfrypertotmean_3$TREATMENT <- factor(Thermalfrypertotmean_3$TREATMENT, levels = c( "WWW", "CWW", "WCW", "CCW","CCC", "WCC", "CWC", "WWC"))


# Barplot

malfBAR1<-ggplot(Thermalfrypertotmean_3, aes(x=TREATMENT, y=malfnoeyeprop, fill=TREATMENT)) 

malfBAR1 + scale_fill_manual(values = c( "#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + theme_classic()



#geom_bar(stat="summary", fun.y="mean_sdl")

malfBAR2<-malfBAR1 + stat_summary(fun.y = mean, geom="bar") + stat_summary(fun.data=mean_cl_boot,geom="errorbar")  + geom_jitter()
malfBAR3<-malfBAR2+ scale_fill_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e" ))

malfBAR3
malfBAR3  + facet_grid(. ~fct_rev(DEVELOPMENTALT), scale="free") + scale_y_continuous(name = bquote ("Malformed embryos (proportion)")) + theme_classic()


deadonhatchBAR1<-ggplot(Thermalfrypertotmean_3, aes(x=TREATMENT, y=deadonhatchprop, fill=TREATMENT)) 

deadonhatchBAR1 + scale_fill_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e" ))



#geom_bar(stat="summary", fun.y="mean_sdl")

deadonhatchBAR2<-deadonhatchBAR1 + stat_summary(fun.y = mean, geom="bar") + stat_summary(fun.data=mean_cl_boot,geom="errorbar")  + geom_jitter()
deadonhatchBAR3<-deadonhatchBAR2+ scale_fill_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e")) + theme_classic()
deadonhatchBAR3

deadonhatchBAR3+ facet_grid(. ~fct_rev(DEVELOPMENTALT), scale="free")+ scale_y_continuous(name = bquote ("Embryos dead on hatch (proportion)")) 




eyeddeadBAR1<-ggplot(Thermalfrypertotmean_3, aes(x=TREATMENT, y=eyeddeadprop, fill=TREATMENT)) 

eyeddeadBAR1 + scale_fill_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e" ))

#geom_bar(stat="summary", fun.y="mean_sdl")

eyeddeadBAR2<-eyeddeadBAR1 + stat_summary(fun.y = mean, geom="bar") + stat_summary(fun.data=mean_cl_boot,geom="errorbar")  + geom_jitter()
eyeddeadBAR3<-eyeddeadBAR2+ scale_fill_manual(values = c("#990000","#d7301f", "#ef6548", "#fc8d59","#634a25","#ae8b41", "#D1BD94", "#fdd49e" )) + theme_classic()
eyeddeadBAR3

eyeddeadBAR3 +facet_grid(. ~fct_rev(DEVELOPMENTALT), scale="free") + scale_y_continuous(name = bquote ("Embryos dead at eyed stage (proportion)"))