---
title: "R Notebook for Fig.4 in second manuscript"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. In this file, we will sort the burrowing depth data and make the plots by heatwave phases.

## 1. loading packages and function
```{r}
library(tidyr)
library(gtools)
library(ggplot2)
library(RColorBrewer)
library(ggpubr)
library(lme4)
library(lmerTest)
library(car)

source('summ_function.R')

```

## 2. loading data
```{r}
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
```{r}
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)
```{r}
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)
```{r}
# 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
```{r}
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
```{r}
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
```{r}
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
```{r}
m2 <- lmer(deltaL ~ group + day + (1|individual), df.calm)
summary(m2)
anova(m2)

m3 <- lmer(deltaL ~ group + day + (1|individual), df.hws)
summary(m3)
anova(m3)

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

#check data
leveneTest(deltaL ~ group, data = df.calm) 
plot(m2)
aov_residuals <- residuals(object = m2)
shapiro.test(x = aov_residuals)
ks.test(df.calm$deltaL, "pnorm")


```
