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)
library(gtools)
library(ggplot2)
library(RColorBrewer)
library(ggpubr)
library(lme4)
library(lmerTest)
library(car)

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

summary(glht(m2, linfct = mcp(group = "Tukey")), test = adjusted("holm"))

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = deltaL ~ group + day + (1 | individual), data = df.calm)

Linear Hypotheses:
                        Estimate Std. Error z value Pr(>|z|)
calm_hf - ambient == 0  0.100469   0.110215   0.912    0.918
calm_lf - ambient == 0  0.094993   0.092813   1.023    0.918
calm_lf - calm_hf == 0 -0.005476   0.134321  -0.041    0.967
(Adjusted p values reported -- holm method)