R Notebook for Fig.4 in second manuscript

This is an R Markdown Notebook. In this file, we will sort the burrowing depth data and make the plots by heatwave phases.

1. loading packages and function

library(tidyr)

Attaching package: ‘tidyr’

The following objects are masked from ‘package:Matrix’:

    expand, pack, unpack

library(gtools)
library(ggplot2)

Need help getting started? Try the R Graphics Cookbook: https://r-graphics.org

library(RColorBrewer)
library(ggpubr)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

library(lme4)
library(lmerTest)
library(car)

Loading required package: carData

Attaching package: ‘car’

The following object is masked from ‘package:gtools’:

    logit

source('summ_function.R')

2. loading data

depth.hf = read.csv('hf_deltaD.csv', header = T, stringsAsFactors=FALSE)
depth.hf$newtime = depth.hf$time
depth.hf$newtime = as.character(depth.hf$newtime)
depth.hf$newtime[depth.hf$newtime == 'T1'] = '9:00'
depth.hf$newtime[depth.hf$newtime == 'T2'] = '11:00'
depth.hf$newtime[depth.hf$newtime == 'T3'] = '12:30'
depth.hf$newtime[depth.hf$newtime == 'T4'] = '14:00'
depth.hf$newtime[depth.hf$newtime == 'T5'] = '16:00'
depth.hf$datetime = paste0(depth.hf$date, " ", depth.hf$newtime)
depth.hf = subset(depth.hf, phases != 'lumino')

depth.lf = read.csv('lf_deltaD.csv', header = T, stringsAsFactors=FALSE)
depth.lf$newtime = depth.lf$time
depth.lf$newtime = as.character(depth.lf$newtime)
depth.lf$newtime[depth.lf$newtime == 'T1'] = '9:00'
depth.lf$newtime[depth.lf$newtime == 'T2'] = '11:00'
depth.lf$newtime[depth.lf$newtime == 'T3'] = '12:30'
depth.lf$newtime[depth.lf$newtime == 'T4'] = '14:00'
depth.lf$newtime[depth.lf$newtime == 'T5'] = '16:00'
depth.lf$datetime = paste0(depth.lf$date, " ", depth.lf$newtime)
depth.lf = subset(depth.lf, phases != 'lumino')

3. subset data by heatwave types and phases

hf.calm = depth.hf[depth.hf$phases == 'cooling1' | 
                     depth.hf$phases == 'cooling2' | 
                     depth.hf$phases == 'cooling3' | 
                     depth.hf$phases == 'cooling4', ]
hf.hws = depth.hf[depth.hf$phases == 'heating1' | 
                     depth.hf$phases == 'heating2' | 
                     depth.hf$phases == 'heating3' | 
                     depth.hf$phases == 'heating4' |
                     depth.hf$phases == 'heating5', ]

lf.calm = depth.lf[depth.lf$phases == 'cooling1' | 
                     depth.lf$phases == 'cooling2', ]
lf.hws = depth.lf[depth.lf$phases == 'heating1' | 
                     depth.lf$phases == 'heating2' | 
                     depth.lf$phases == 'heating3', ]

4. summarize data (check)

hf.calm.c = hf.calm[hf.calm$treatment == 'control',]
hf.calm.h = hf.calm[hf.calm$treatment == 'heating',]

hf.hws.c = hf.hws[hf.hws$treatment == 'control',]
hf.hws.h = hf.hws[hf.hws$treatment == 'heating',]

lf.calm.c = lf.calm[lf.calm$treatment == 'control',]
lf.calm.h = lf.calm[lf.calm$treatment == 'heating',]

lf.hws.c = lf.hws[lf.hws$treatment == 'control',]
lf.hws.h = lf.hws[lf.hws$treatment == 'heating',]

abt.c = rbind(hf.calm.c, hf.hws.c, lf.calm.c, lf.hws.c)
abt.c$group = rep('ambient', length(abt.c$deltaL))

hf.calm.h$group = rep('calm_hf', length(hf.calm.h$deltaL))
lf.calm.h$group = rep('calm_lf', length(lf.calm.h$deltaL))

hf.hws.h$group = rep('hws_hf', length(hf.hws.h$deltaL))
lf.hws.h$group = rep('hws_lf', length(lf.hws.h$deltaL))

df.calm = rbind(abt.c, hf.calm.h, lf.calm.h)
df.hws = rbind(abt.c, hf.hws.h, lf.hws.h)

###########################################################
df.calm.summ = summarySE(df.calm, measurevar='deltaL', groupvars=c('group','time'), na.rm=T)

Attaching package: ‘plyr’

The following object is masked from ‘package:ggpubr’:

    mutate

df.hws.summ = summarySE(df.hws, measurevar='deltaL', groupvars=c('group','time'), na.rm=T)

df.calm.summ$group[df.calm.summ$group == 'calm_hf'] = 'heatwaves: 3-day cycle'
df.calm.summ$group[df.calm.summ$group == 'calm_lf'] = 'heatwaves: 6-day cycle'

df.hws.summ$group[df.hws.summ$group == 'hws_hf'] = 'heatwaves: 3-day cycle'
df.hws.summ$group[df.hws.summ$group == 'hws_lf'] = 'heatwaves: 6-day cycle'

4. summarize data (optionnal)

# hf.calm.sum = summarySE(hf.calm, measurevar='deltaL', groupvars=c('treatment','time'), na.rm=T)
# hf.calm.sum$treatment[hf.calm.sum$treatment == 'control'] = 'control: 3-day cycle'
# hf.calm.sum$treatment[hf.calm.sum$treatment == 'heating'] = 'heatwaves: 3-day cycle'
# 
# hf.hws.sum = summarySE(hf.hws, measurevar='deltaL', groupvars=c('treatment','time'), na.rm=T)
# hf.hws.sum$treatment[hf.hws.sum$treatment == 'control'] = 'control: 3-day cycle'
# hf.hws.sum$treatment[hf.hws.sum$treatment == 'heating'] = 'heatwaves: 3-day cycle'
# 
# lf.calm.sum = summarySE(lf.calm, measurevar="deltaL", groupvars=c("treatment",'time'),na.rm=T)
# lf.calm.sum$treatment[lf.calm.sum$treatment == 'control'] = 'control: 6-day cycle'
# lf.calm.sum$treatment[lf.calm.sum$treatment == 'heating'] = 'heatwaves: 6-day cycle'
# 
# lf.hws.sum = summarySE(lf.hws, measurevar="deltaL", groupvars=c("treatment",'time'),na.rm=T)
# lf.hws.sum$treatment[lf.hws.sum$treatment == 'control'] = 'control: 6-day cycle'
# lf.hws.sum$treatment[lf.hws.sum$treatment == 'heating'] = 'heatwaves: 6-day cycle'
# 
# calm.sum = rbind(lf.calm.sum, hf.calm.sum)
# hws.sum = rbind(lf.hws.sum, hf.hws.sum)

5. change time format

df.calm.summ$time[df.calm.summ$time == 'T1'] = '9:00'
df.calm.summ$time[df.calm.summ$time == 'T2'] = '11:00'
df.calm.summ$time[df.calm.summ$time == 'T3'] = '12:30'
df.calm.summ$time[df.calm.summ$time == 'T4'] = '14:00'
df.calm.summ$time[df.calm.summ$time == 'T5'] = '16:00'

df.hws.summ$time[df.hws.summ$time == 'T1'] = '9:00'
df.hws.summ$time[df.hws.summ$time == 'T2'] = '11:00'
df.hws.summ$time[df.hws.summ$time == 'T3'] = '12:30'
df.hws.summ$time[df.hws.summ$time == 'T4'] = '14:00'
df.hws.summ$time[df.hws.summ$time == 'T5'] = '16:00'

6. save data

write.csv(df.calm.summ, 'df_calm_summ.csv', row.names = F)
write.csv(df.hws.summ, 'df_hws_summ.csv', row.names = F)

7. loading data

calm.profile = read.csv('df_calm_summ.csv', header = T)
calm.profile$time = factor(calm.profile$time, levels = c('9:00', '11:00', '12:30', '14:00', '16:00'))

hws.profile = read.csv('df_hws_summ.csv', header = T)
hws.profile$time = factor(hws.profile$time, levels = c('9:00', '11:00', '12:30', '14:00', '16:00'))

8. setting plot parameters

pd = position_dodge(0.05)

hw_frame = data.frame(start = 1.25, end = 4.75)

9. plotting for calm profile

plot.calm = ggplot(calm.profile, aes(x=time, y=deltaL, colour=group, group = group)) + 
  geom_errorbar(aes(ymin=deltaL-se, ymax=deltaL+se), colour="black", width=0.3, position=pd) +
  geom_line(position=pd, size = 2.5, aes(group = group)) +
  geom_point(position=pd, size=3, shape=16) + 
  geom_vline(xintercept = 1.5, color = "#9F6554", size=1.0) +
  geom_vline(xintercept = 4.5, linetype="dashed", color = "#9F6554", size=1.0) +
  xlab("measurement time") +
  ylab('depth change (mm)') +
  scale_y_continuous(limits = c(-3,4), breaks = c(-2,-1,0,1,2,3))+
  scale_colour_manual(name="Treatments", 
                      breaks = c('ambient', 'heatwaves: 3-day cycle', 'heatwaves: 6-day cycle'), 
                      values= c('#3399FF', '#FF9900', '#FF0000'), 
                      labels = c('Ambient temperature', 'Heatwaves: 3-day cycle','Heatwaves: 6-day cycle')) +
  theme(
    plot.title = element_blank(),
    panel.background = element_rect(fill = 'white', colour = 'white'),
    legend.title = element_text(size=28, face="bold", color = "black"), 
    legend.text = element_text(size=28, face="bold"),
    legend.position=c(0.3, 0.85),
    legend.background = element_rect(fill=alpha('white', 0.8),
                                     size=0.5, linetype="solid", 
                                     colour ="black"),
    legend.margin=margin(r = 1, unit='inch'),
    legend.key = element_blank(),
    legend.key.width = unit(0.5, 'inch'),
    legend.key.height = unit(0.5, 'inch'),
    panel.border = element_rect(colour = "black", fill=NA, size=2),
    axis.title.x = element_text(size = 32, face="bold"),
    axis.title.y = element_text(size = 32, face="bold"),
    axis.text.x = element_text(size = 27, face="bold"),
    axis.text.y = element_text(size = 28, face="bold"),
    plot.margin = unit(c(0.1,0.6,0.1,0.1), "in"))+
  guides(color = guide_legend(override.aes = list(size = 4)))
plot.calm

10. plotting for heatwave profiles

plot.hws = ggplot(hws.profile, aes(x=time, y=deltaL, colour=group, group = group)) +
  geom_errorbar(aes(ymin=deltaL-se, ymax=deltaL+se), colour="black", width=0.3, position=pd) +
  geom_line(position=pd, size = 2.5, aes(group = group)) +
  geom_point(position=pd, size=3, shape=16) + 
  geom_vline(xintercept = 1.5, color = "#9F6554", size=1.0) +
  geom_vline(xintercept = 4.5, linetype="dashed", color = "#9F6554", size=1.0) +
  xlab("measurement time") +
  ylab('depth change (mm)') +
  scale_y_continuous(limits = c(-3,4), breaks = c(-2,-1,0,1,2,3))+
  scale_colour_manual(name="Treatments", 
                      breaks = c('ambient', 'heatwaves: 3-day cycle', 'heatwaves: 6-day cycle'), 
                      values= c('#3399FF', '#FF9900', '#FF0000'), 
                      labels = c('Ambient temperature', 'Heatwaves: 3-day cycle','Heatwaves: 6-day cycle')) +
  theme(
    plot.title = element_blank(),
    panel.background = element_rect(fill = 'white', colour = 'white'),
    legend.title = element_text(size=28, face="bold", color = "black"), 
    legend.text = element_text(size=28, face="bold"),
    legend.position=c(0.3, 0.85),
    legend.background = element_rect(fill=alpha('white', 0.8),
                                     size=0.5, linetype="solid", 
                                     colour ="black"),
    legend.margin=margin(r = 1, unit='inch'),
    legend.key = element_blank(),
    legend.key.width = unit(0.5, 'inch'),
    legend.key.height = unit(0.5, 'inch'),
    panel.border = element_rect(colour = "black", fill=NA, size=2),
    axis.title.x = element_text(size = 32, face="bold"),
    axis.title.y = element_text(size = 32, face="bold"),
    axis.text.x = element_text(size = 27, face="bold"),
    axis.text.y = element_text(size = 28, face="bold"),
    plot.margin = unit(c(0.1,0.6,0.1,0.1), "in"))+
  guides(color = guide_legend(override.aes = list(size = 4)))
plot.hws = plot.hws + annotate("rect", xmin = 1.25, xmax = 4.75, ymin = -Inf, ymax = Inf, alpha = .15, fill = "#E3242B")

11. save images

ggsave('plot.calm.png', plot.calm, units = 'in', dpi = 300, width = 16, height = 12)
ggsave('plot.hws.png', plot.hws, units = 'in', dpi = 300, width = 16, height = 12)

12. merge plots

plot.merge = ggarrange(plot.calm, plot.hws, 
          labels = c("a", "b"),
          ncol = 1, nrow = 2, font.label = list(size = 35))
ggsave('plot_merge.png', plot.merge, units = 'in', dpi = 300, width = 16, height = 24)

13. statistics

ks.test(df.calm$deltaL, "pnorm")

Warning: ties should not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  df.calm$deltaL
D = 0.39201, p-value < 2.2e-16
alternative hypothesis: two-sided