R Notebook

This is an R Markdown Notebook for Figure 4a.

1. loading packages and function

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

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

4. setting parameters

pd = position_dodge(0.1)

5. plotting

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)))

Warning: Vectorized input to `element_text()` is not officially supported.
Results may be unexpected or may change in future versions of ggplot2.

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)