m3 <- lmer(deltaL ~ group * day + (1 + individual), df.hws)
m3 <- lmer(deltaL ~ group * day + (1 | individual), df.hws)
summary(m3)
anova(m3)
m2 <- lmer(deltaL ~ group * day + (1|individual), df.calm)
summary(m2)
anova(m2)
library(multcomp)
summary(glht(m3, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(m2, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(m3, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(m2, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(m3, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
m3 <- lmer(deltaL ~ group * day + (1 + day|individual), df.hws)
m3 <- lmer(deltaL ~ group + day + group:day + (1 + day|individual), df.hws)
m3 <- lmer(deltaL ~ group * day + (1|individual), df.hws)
summary(m3)
anova(m3)
m2 <- lmer(deltaL ~ group * day + (1|individual), df.calm)
summary(m2)
anova(m2)
library(multcomp)
summary(glht(m3, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(m2, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
library(emmeans)
pairs(emmeans(m3, "Treatment"))
pairs(emmeans(m3, "group"))
pairs(emmeans(m2, "group"))
summary(glht(m3, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
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"))
m2 <- lmer(deltaL ~ group + day + (1|individual), df.calm)
m2 <- lmer(deltaL ~ group + day + (1|individual), df.calm)
summary(m2)
anova(m2)
m3 <- lmer(deltaL ~ group + day + (1|individual), df.hws)
m3 <- lmer(deltaL ~ group + day + (1|individual), df.hws)
summary(m3)
summary(m3)
anova(m3)
anova(m2)
anova(m3)
library(multcomp)
summary(glht(m3, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(m2, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(m3, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
fit0 = kruskal.test(area_adjusted_perc ~ Setting + Depth, data = lumino_freq)
library(tidyr)
library(gtools)
library(ggplot2)
library(RColorBrewer)
library(ggpubr)
library(lme4)
library(lmerTest)
source('summ_function.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')
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', ]
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'
# 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)
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'
write.csv(df.calm.summ, 'df_calm_summ.csv', row.names = F)
write.csv(df.hws.summ, 'df_hws_summ.csv', row.names = F)
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'))
pd = position_dodge(0.05)
hw_frame = data.frame(start = 1.25, end = 4.75)
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
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")
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)
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)
calm.aov <- aov(deltaL ~ group + time, data = df.calm)
summary(calm.aov)
hws.aov <- aov(deltaL ~ group + time, data = df.hws)
summary(hws.aov)
TukeyHSD(hws.aov)
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 * day, data = df.calm) #plot(fit2)
library(car)
#check data
leveneTest(deltaL ~ group * day, data = df.calm) #plot(fit2)
aov_residuals <- residuals(object = m2)
shapiro.test(x = aov_residuals)
#check data
leveneTest(deltaL, data = df.calm) #plot(fit2)
#check data
leveneTest(df.calm$deltaL) #plot(fit2)
#check data
leveneTest(deltaL ~ group, data = df.calm) #plot(fit2)
#check data
leveneTest(deltaL ~ group + day, data = df.calm) #plot(fit2)
#check data
leveneTest(deltaL ~ group * day + (1|individual), data = df.calm) #plot(fit2)
#check data
leveneTest(deltaL ~ group * day * individual, data = df.calm) #plot(fit2)
aov_residuals <- residuals(object = m2)
shapiro.test(x = aov_residuals)
plot(m2)
plot(m2)
#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")
plot(m2, 2)
plot(m2)
View(m2)
hist(df.calm$deltaL, breaks = 1000)
hist(df.calm$deltaL, breaks = 10)
hist(df.calm$deltaL, breaks = 100)
hist(df.hws$deltaL, breaks = 100)
hist(df.calm$deltaL, breaks = 100)
hist(df.calm$deltaL, breaks = 100, xlim = c(-10,10))
hist(df.calm$deltaL, breaks = 1000, xlim = c(-10,10))
hist(df.calm$deltaL, breaks = 100, xlim = c(-10,10))
hist(df.calm$deltaL, breaks = seq(-10,10,1), xlim = c(-10,10))
hist(df.calm$deltaL, breaks = seq(-100,100,1), xlim = c(-10,10))
df.calm$treatment
hist(df.calm$deltaL[df.calm$treatment == 'control'], breaks = seq(-100,100,1), xlim = c(-10,10))
df.calm$time
hist(df.calm$deltaL[df.calm$treatment == 'control' & df.calm$time == 'T3'], breaks = seq(-100,100,1), xlim = c(-10,10))
hist(df.calm$deltaL[df.calm$treatment == 'control' & df.calm$time == 'T3'], breaks = seq(-100,100,1))
shapiro.test(df.calm$deltaL[df.calm$treatment == 'control' & df.calm$time == 'T3'])
shapiro.test(df.calm$deltaL[df.calm$treatment == 'low' & df.calm$time == 'T3'])
df.calm$treatment
unique(df.calm$treatment)
shapiro.test(df.calm$deltaL[df.calm$treatment == 'heating' & df.calm$time == 'T3'])
hist(df.calm$deltaL[df.calm$treatment == 'heating' & df.calm$time == 'T3'], breaks = seq(-100,100,1))
df.calm$group
hist(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'], breaks = seq(-100,100,1))
shapiro.test(df.calm$deltaL[df.calm$treatment == 'ambient' & df.calm$time == 'T3'])
shapiro.test(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'])
unique(df.calm$treatment)
qqnorm(hist(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'], breaks = seq(-100,100,1)))
x = df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3']
qqnorm(x)
hist(x)
qqnorm(x[x > -10 & x < 10])
qqnorm(x[x > -7 & x < 7])
shapiro.test(x[x > -7 & x < 7])
shapiro.test(x[x > -10 & x < 10])
shapiro.test(x[x > -5 & x < 5])
shapiro.test(x[x > -1 & x < 1])
hist(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'], breaks = seq(-100,100,1))
hist(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'], breaks = seq(-100,100,1))
hist(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'], breaks = seq(-100,100,1))
hist(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'], breaks = seq(-100,100,1))
qqnorm(hist(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'], breaks = seq(-100,100,1)))
shapiro.test(df.calm$deltaL[df.calm$group == 'ambient' & df.calm$time == 'T3'])
library(tidyr)
library(gtools)
library(ggplot2)
library(RColorBrewer)
library(ggpubr)
library(lme4)
library(lmerTest)
library(car)
source('summ_function.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')
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', ]
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'
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'
write.csv(df.calm.summ, 'df_calm_summ.csv', row.names = F)
write.csv(df.hws.summ, 'df_hws_summ.csv', row.names = F)
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'))
pd = position_dodge(0.05)
hw_frame = data.frame(start = 1.25, end = 4.75)
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
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")
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)
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)
calm.aov <- aov(deltaL ~ group + time, data = df.calm)
summary(calm.aov)
hws.aov <- aov(deltaL ~ group + time, data = df.hws)
summary(hws.aov)
TukeyHSD(hws.aov)
# revised data analysis
m2 <- lmer(deltaL ~ group + day + (1|individual), df.calm)
summary(m2)
anova(m2)
m3 <- lmer(deltaL ~ group + day + (1|individual), df.hws)
m3 <- lmer(deltaL ~ group + day + (1|individual), df.hws)
summary(m3)
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"))
data(varespec)
data(varechem)
sol <- bioenv(wisconsin(varespec) ~ log(N) + P + K + Ca + pH + Al, varechem)
sol
summary(sol)
library(vegan)
data(varespec)
data(varechem)
sol <- bioenv(wisconsin(varespec) ~ log(N) + P + K + Ca + pH + Al, varechem)
sol
summary(sol)
View(varespec)
View(varechem)
View(varespec)
#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")
