---
title: "R Notebook"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook for Figure 4a.

## 1. loading packages and function
```{r}
library(ggplot2)
library(multcomp)
library(rstatix)

# Define a function to summarizes data.
summarySE = function(data=NULL, measurevar, 
                     groupvars=NULL, na.rm=FALSE,
                     conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # if na.rm==T, don't count them
  length2 = function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else length(x)}
  
  # for each group's data frame, return a vector with N, mean, and sd
  datac = ddply(data, groupvars, 
                .drop = .drop,.fun = function(xx, col) {
                  c(N = length2(xx[[col]], na.rm = na.rm),
                    mean = mean(xx[[col]], na.rm = na.rm),
                    sd = sd(xx[[col]], na.rm = na.rm))},
                measurevar)
  
  # rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  # Calculate standard error of the mean
  datac$se <- datac$sd / sqrt(datac$N)  
  
  # calculate t-statistic for confidence interval: 
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}

```

## 2. loading data
```{r}
lumino = read.csv('alldata.csv', header = T)
lumino = subset(lumino, Depth != '0' & Depth != '0.5' & Depth != '4'& Depth != '5')
lumino$Depth = as.numeric(lumino$Depth)

lumino_freq = subset(lumino, Setting != 'monitoring')
lumino_mntr = subset(lumino, Setting == 'monitoring')

```

## 3. summarize data
```{r}
lumino_freq$bioturbated_area = (lumino_freq$area_adjusted_perc)/100 * 35.23865

lumino_coresum = aggregate(bioturbated_area ~ Day + Setting + Tank + Replicat, data = lumino_freq, FUN = sum)
lumino_coresum$bioturbated_areaDA = lumino_coresum$bioturbated_area/5

# statistics ANOVA
lumino_coresum$Setting = as.factor(lumino_coresum$Setting)
res.aov <- aov(bioturbated_areaDA ~ Setting , data = lumino_coresum)
summary(res.aov)

TukeyHSD(res.aov)
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")))

res.aov1 = kruskal.test(bioturbated_areaDA ~ Setting , data = lumino_coresum)
summary(res.aov1)

pwt = pairwise.wilcox.test(lumino_coresum$bioturbated_areaDA, lumino_coresum$Setting,
                 p.adjust.method = "BH")
wilcox_effsize(bioturbated_areaDA ~ Setting, data = lumino_coresum)

freq_summ = summarySE(lumino_coresum, measurevar="bioturbated_areaDA", groupvars=c('Setting'), na.rm=T)

freq_summ$Setting = factor(freq_summ$Setting, levels = c("control", "high", 'low'))

```


## 4. setting parameters
```{r}
pd = position_dodge(0.1)

```

## 5. plotting
```{r}
lmPlot = ggplot(freq_summ, aes(x=Setting, y=bioturbated_areaDA, colour=Setting)) + 
  geom_errorbar(aes(ymin=bioturbated_areaDA - se, ymax=bioturbated_areaDA + se), 
                colour="black", width=.1, position=pd) +
  #geom_line(aes(group = 1, color = factor(setting))) +
  geom_point(position=pd, size=6, shape=16) +
  scale_colour_manual(name="Treatments", 
                      breaks=c("control", "high", 'low'),
                      labels = c("Ambient temperature", 
                      "3-day cycle heatwaves", 
                      "6-day cycle heatwaves"), 
                      values=c('#3399FF', '#FF9900', '#FF0000')) +
  scale_x_discrete(name = 'Heatwave scenarios', 
                     breaks=c("control", "high", 'low'),
                     labels = c("Ambient\ntemperature", 
                                 "3-day cycle\nheatwaves", 
                                 "6-day cycle\nheatwaves"))+
  ylab(expression(bold('Depth-averaged bioturbated surface area '~(cm^2))))+
  theme(
    plot.title = element_blank(),
    panel.background = element_blank(), 
    legend.position="none",
    axis.title.x = element_text(size = 24, face="bold"),
    axis.title.y = element_text(size = 24, face="bold"),
    axis.text.x = element_text(size =22, face="bold", colour=c('#3399FF', '#FF9900', '#FF0000')),
    axis.text.y = element_text(size = 22, face="bold"),
    axis.ticks = element_line(size = 2),
    axis.ticks.length = unit(0.08, "inch"),
    plot.margin = unit(c(0.1,0.1,0.1,0.1), "in"))+
  guides(color = guide_legend(override.aes = list(size = 5)))
lmPlot

lmPlot = lmPlot + 
  geom_segment(aes(x=1, xend=3, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
  geom_segment(aes(y=0, yend=6, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
lmPlot

ggsave('ldm_neo.png', lmPlot, units = 'in', width = 12, height = 9,dpi=600)
```




