---
title: "Heatwave_Figures code"
author: "Matthew Clements"
date: "`r Sys.Date()`"
output: html_document
---

#Load packages
```{r}
#Righting and Respiration
library(car)
library(readr)
library(ggplot2)
library(ggpubr)
library(car)
library(nlme)
library(Rmisc)
library(emmeans)
library(readxl)
library(cowplot)

#Survival
library(readr)
library(readxl)
library(survival)
library(survminer)
install.packages("RColorBrewer")
library(RColorBrewer)
```

#Righting (experiment 1)
##Load data
```{r}
RTo<-read_excel("/Users/matthewclements/Desktop/GCB/Experiment 1-Respiration & Righting.xls", sheet = "Righting", range = "A2:D36", col_names = c("Temperature",	"Righting",	"Diameter",	"Standardised"))
RT<-RTo[-c(28:35), ]

RTsum1<-summarySE(RT, measurevar="Standardised", group=c("Temperature"), na.rm=TRUE) # averaging between same ID 
```

##Figure (Format for Global Change Biology)
```{r}
GCB_right<-ggplot(RTsum1, aes(x=as.factor(Temperature), y=Standardised))+
  geom_bar(position="dodge", stat="identity",colour="black",fill="lightgrey")+
  ggtitle("","")+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        axis.text=element_text(size=14, family = "Times New Roman"),
        strip.background = element_blank(),
        strip.text.x = element_blank(),
        axis.title.y = element_text(size=18, family = "Times New Roman", margin = margin(t = 0, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(size=18, family = "Times New Roman", margin = margin(t = 10, r = 10, b = 0, l = 0)))+
  labs(x= "Temperature", y= "Standardised righting time (s)")+
  annotate("text", x = 1, y = 61, label = "ab", size = 3.8) +
  annotate("text", x = 2, y = 75, label = "b", size = 3.8) +
  annotate("text", x = 3, y = 61, label = "ab", size = 3.8) +
  annotate("text", x = 4, y = 32, label = "a", size = 3.8) +
  #annotate("text", x = 5, y = 35, label = "a", size = 3.7) +
  theme(text=element_text(family="Times New Roman", face="plain", size=14))+
  geom_errorbar(aes(ymin=Standardised-se, ymax=Standardised+se), width=.2, colour="black", position=position_dodge(.9))+
  scale_y_continuous(expand=c(0,0),limits = c(0, 80))
GCB_right
```

#Respiration (experiment 1)
##Load data
```{r}
resp<-read_excel("/Users/matthewclements/Desktop/GCB/Experiment 1-Respiration & Righting.xls", sheet = "Respiration", range = "A2:B37", col_names = c("Temperature",	"MO2"))

resp$Temperature<-as.factor(resp$`Temperature`)

respsum<-summarySE(resp, measurevar="MO2", group=c("Temperature"), na.rm=TRUE) # averaging between same ID 
```


##Figure (Format for Global Change Biology)
```{r}
#respiration
GCB_res<-ggplot(respsum, aes(x=Temperature, y=MO2))+
  geom_bar(position="dodge", stat="identity", colour="black",fill="lightgrey")+
  ggtitle("","")+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        axis.text=element_text(size=14, family = "Times New Roman"),
        strip.background = element_blank(),
        strip.text.x = element_blank(),
        axis.title.y = element_text(size=18, family = "Times New Roman", margin = margin(t = 0, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(size=18, family = "Times New Roman", margin = margin(t = 10, r = 0, b = 0, l = 0)))+
  labs(x= "Temperature (˚C)", y= expression(M0[2]~(mg~O[2]~ g^-1 ~ h^-1)), colour="")+
  annotate("text", x = 1, y = 0.07, label = "a", size = 3.7) +
  annotate("text", x = 2, y = 0.053, label = "a", size = 3.7) +
  annotate("text", x = 3, y = 0.075, label = "ab", size = 3.7) +
  annotate("text", x = 4, y = 0.102, label = "b", size = 3.7) +
  #annotate("text", x = 5, y = 0.052, label = "ab", size = 3.7) +
  theme(text=element_text(family="Times New Roman", face="plain", size=14))+
  scale_fill_manual(name="pH", labels = c("Ambient", "Low"), values=c("cornflowerblue", "brown3"))+
  geom_errorbar(aes(ymin=MO2-se, ymax=MO2+se), width=.2, colour="black", position=position_dodge(.9))+
  coord_cartesian(ylim = c(0, 0.12))+ 
  scale_y_continuous(expand=c(0,0), breaks=c(0,0.02,0.04,0.06, 0.08,0.10, 0.12))
GCB_res

```

#Combined righting and respiration plots
```{r}

GCB_res<-plot_grid(GCB_right, GCB_res, ncol = 2, labels = c("(a)", "(b)"),
               hjust = 0,
               label_fontfamily = "Times New Roman",
               label_size = 18)
GCB_res
```

##Export Figure
```{r}
ggsave(GCB_res, 
       filename = "GCB.png", 
       path = "/Users/matthewclements/Desktop/Heatwave coral and COTS juvie fig/GCB", 
       scale = 1, width = 15,
       dpi = 600)

ggsave(GCB_res, 
       filename = "GCB.jpg", 
       path = "/Users/matthewclements/Desktop/Heatwave coral and COTS juvie fig/GCB", 
       scale = 1, width = 15, 
       dpi = 600)

ggsave(GCB_res, 
       filename = "GCB.TIFF", 
       path = "/Users/matthewclements/Desktop/Heatwave coral and COTS juvie fig/GCB", 
       scale = 1, width = 15,
       dpi = 600)

GCB<-plot_grid(GCB_right, GCB_res, ncol = 2, labels = c("(a)", "(b)"),
               hjust = 0,
               label_fontfamily = "Times New Roman",
               label_size = 18)

# Closing the graphical device
dev.off() 
```

#Survival (experiment 2)

##Import data
```{r}
df<-read_excel("/Users/matthewclements/Desktop/GCB/Experiment 2-Survival.xls", sheet = "Survival", range = "A2:F41", col_names = c("Individual", 	"start",	"end",	"censor",	"Days",	"Temperature"))

df2 <- df[-c(27:40),]

# need to fix censoring 
# currently, 1 = censored, 0 = died  
# change to: FALSE = censored, TRUE = died
```

##Filter data by censor
```{r}
df_fix <- df %>% 
  # change 1 to something else first
  mutate(censor = replace(censor, censor == 0, 2)) %>%
  # change 0 to FALSE
  mutate(censor = replace(censor, censor == 1, FALSE)) %>%
  # change 2 to TRUE
  mutate(censor = replace(censor, censor == 2, TRUE))

df_fix2 <- df2 %>% 
  # change 1 to something else first
  mutate(censor = replace(censor, censor == 0, 2)) %>%
  # change 0 to FALSE
  mutate(censor = replace(censor, censor == 1, FALSE)) %>%
  # change 2 to TRUE
  mutate(censor = replace(censor, censor == 2, TRUE))

# check treatment levels:
unique(df_fix$Temperature)

unique(df_fix2$Temperature)
```

##Define variables
```{r}
# defininf the variables in df:
# Days: survival time in days
# censor: censored = FALSE; died = TRUE
# treatment: "CCA", "Amphiroa", "Biofilm"

# check survival object:
Surv(df_fix$Days, df_fix$censor)

Surv(df_fix2$Days, df_fix2$censor)

# create survival object:
smodel <- survfit(Surv(Days, censor) ~ Temperature, data = df_fix)
smodel

smodel2 <- survfit(Surv(Days, censor) ~ Temperature, data = df_fix2, conf.type="plain", conf.int=0.95)
smodel2

# plot
plot(smodel) # check
```

##Plots
```{r}
# publication-quality plot:
ggsurvplot(smodel2, 
           conf.int = T,
           # ggtheme = theme_bw(), # Change ggplot2 theme
           legend = "right",
           legend.title= c("Temperature"),
           legend.labs=c("27-30˚C","31˚C","32˚C","34˚C","36˚C"),
           palette = c("lightskyblue","tan1","lightcoral","red4","Black"))+
  xlab("Time (Days)")


# check results of logrank test
# p value shows differences between treatments
lgrank <- survdiff(Surv(Days, censor) ~ Temperature, data = df_fix)
lgrank
?survdiff
# post hoc pairwise using Benjamini-Hochberg adjusted p value:
pairwise_survdiff(Surv(Days, censor) ~ Temperature, data = df_fix)

```

##Export Figure
```{r}
ggsave(file = "ggsurv.png", 
       path = "/Users/matthewclements/Desktop/Work/SOLES/Projects/Asteroidea/COTS/Juvenile/Heatwave exp/Heatwave coral and COTS juvie fig/Base/Survival fig",
       scale = 1,
       height = 8,
       width = 10,
       dpi = 300)
```