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)))
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=30, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
lmPlot
ggsave('ldm_neo.png', lmPlot, units = 'in', width = 12, height = 9,dpi=600)
View(freq_summ)
View(lumino_coresum)
# statistics ANOVA
res.aov <- aov(bioturbated_area ~ Setting , data = lumino_coresum)
summary(res.aov)
summary(glht(res.aov, linfct = mcp(group = "Tukey")))
library(multcomp)
summary(glht(res.aov, linfct = mcp(group = "Tukey")))
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")))
TukeyHSD(res.aov)
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")), test = adjusted("holm"))
summary(glht(res.aov, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
library(multcomp)
freq_summ = summarySE(lumino_coresum, measurevar="bioturbated_area", groupvars=c('Setting'), na.rm=T)
summary(glht(res.aov, linfct = mcp(group = "Tukey")), test = adjusted("holm"))
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")), test = adjusted("holm"))
TukeyHSD(res.aov)
summary(glht(res.aov, linfct = mcp(group = "Tukey")))
View(res.aov)
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")))
# statistics ANOVA
lumino_coresum$bioturbated_area = as.factor(lumino_coresum$bioturbated_area)
res.aov <- aov(bioturbated_area ~ Setting , data = lumino_coresum)
summary(res.aov)
lumino_freq$bioturbated_area = (lumino_freq$area_adjusted_perc)/100 * 35.23865
lumino_coresum = aggregate(bioturbated_area ~ Day + Setting + Tank + Replicat, data = lumino_freq, FUN = sum)
# statistics ANOVA
lumino_coresum$Setting = as.factor(lumino_coresum$Setting)
res.aov <- aov(bioturbated_area ~ Setting , data = lumino_coresum)
summary(res.aov)
TukeyHSD(res.aov)
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")))
my_comparisons <- list( c("low", "high"), c("low", "control"), c("high", "control") )
my_comparisons <- list( c("low", "high"), c("low", "control"), c("high", "control") )
ggboxplot(lumino_coresum, x = "Setting", y = "bioturbated_area",
color = "Setting", palette = "jco")+
stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
stat_compare_means(label.y = 50)
library(ggplot2)
ggboxplot(lumino_coresum, x = "Setting", y = "bioturbated_area",
color = "Setting", palette = "jco")+
stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
stat_compare_means(label.y = 50)
library(ggpubr)
ggboxplot(lumino_coresum, x = "Setting", y = "bioturbated_area",
color = "Setting", palette = "jco")+
stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
stat_compare_means(label.y = 50)
hist(lumino_freq$bioturbated_area)
hist(lumino_freq$bioturbated_area, breaks = 50)
hist(lumino_freq$bioturbated_area, breaks = 100)
shapiro.test(lumino_freq$bioturbated_area)
res.aov1 = kruskal.test(bioturbated_area ~ Setting , data = lumino_coresum)
summary(res.aov1)
res.aov1
pairwise.wilcox.test(lumino_coresum$bioturbated_area, lumino_coresum$Setting,
p.adjust.method = "BH")
res.aov1 = kruskal.test(bioturbated_area ~ Setting , data = lumino_coresum)
res.aov1
pairwise.wilcox.test(lumino_coresum$bioturbated_area, lumino_coresum$Setting,
p.adjust.method = "BH")
pwt = pairwise.wilcox.test(lumino_coresum$bioturbated_area, lumino_coresum$Setting,
p.adjust.method = "BH")
pwt$p.value
pwt$method
View(pwt)
summary(res.aov)
TukeyHSD(res.aov)
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")))
res.aov1 = kruskal.test(bioturbated_area ~ Setting , data = lumino_coresum)
summary(res.aov1)
pwt = pairwise.wilcox.test(lumino_coresum$bioturbated_area, lumino_coresum$Setting,
p.adjust.method = "BH")
pwt
wilcox.test(bioturbated_area ~ Setting, data = lumino_coresum)
View(lumino_coresum)
wilcox_effsize(bioturbated_area ~ Setting, data = lumino_coresum)
library(tidyverse)
library(rstatix)
library(ggpubr)
wilcox_effsize(bioturbated_area ~ Setting, data = lumino_coresum)
install.packages('coin')
library(coin)
wilcox_effsize(bioturbated_area ~ Setting, data = lumino_coresum)
View(lumino_coresum)
lumino_coresum$bioturbated_areaDA = lumino_coresum$bioturbated_area/5
lumino_freq$bioturbated_area = (lumino_freq$area_adjusted_perc)/100 * 35.23865
lumino_coresum = aggregate(bioturbated_area ~ Day + Setting + Tank + Replicat, data = lumino_freq, FUN = sum)
lumino_coresum$bioturbated_areaDA = lumino_coresum$bioturbated_area/5
# statistics ANOVA
bioturbated_areaDA$Setting = as.factor(bioturbated_areaDA$Setting)
res.aov <- aov(bioturbated_area ~ Setting , data = bioturbated_areaDA)
lumino_freq$bioturbated_area = (lumino_freq$area_adjusted_perc)/100 * 35.23865
lumino_coresum = aggregate(bioturbated_area ~ Day + Setting + Tank + Replicat, data = lumino_freq, FUN = sum)
lumino_coresum$bioturbated_areaDA = lumino_coresum$bioturbated_area/5
# statistics ANOVA
bioturbated_areaDA$Setting = as.factor(bioturbated_areaDA$Setting)
library(ggplot2)
library(multcomp)
# 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)
}
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')
library(ggplot2)
library(multcomp)
# 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)
}
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')
library(ggplot2)
library(multcomp)
# 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)
}
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')
lumino_freq$bioturbated_area = (lumino_freq$area_adjusted_perc)/100 * 35.23865
lumino_coresum = aggregate(bioturbated_area ~ Day + Setting + Tank + Replicat, data = lumino_freq, FUN = sum)
lumino_coresum$bioturbated_areaDA = lumino_coresum$bioturbated_area/5
View(lumino_coresum)
# statistics ANOVA
lumino_coresum$Setting = as.factor(lumino_coresum$Setting)
lumino_freq$bioturbated_area = (lumino_freq$area_adjusted_perc)/100 * 35.23865
lumino_coresum = aggregate(bioturbated_area ~ Day + Setting + Tank + Replicat, data = lumino_freq, FUN = sum)
lumino_coresum$bioturbated_areaDA = lumino_coresum$bioturbated_area/5
# statistics ANOVA
lumino_coresum$Setting = as.factor(lumino_coresum$Setting)
res.aov <- aov(bioturbated_areaDA ~ Setting , data = lumino_coresum)
summary(res.aov)
TukeyHSD(res.aov)
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")))
res.aov1 = kruskal.test(bioturbated_areaDA ~ Setting , data = lumino_coresum)
summary(res.aov1)
pwt = pairwise.wilcox.test(lumino_coresum$bioturbated_areaDA, lumino_coresum$Setting,
p.adjust.method = "BH")
wilcox_effsize(bioturbated_areaDA ~ Setting, data = lumino_coresum)
library(rstatix)
pwt = pairwise.wilcox.test(lumino_coresum$bioturbated_areaDA, lumino_coresum$Setting,
p.adjust.method = "BH")
wilcox_effsize(bioturbated_areaDA ~ Setting, data = lumino_coresum)
freq_summ = summarySE(lumino_coresum, measurevar="bioturbated_areaDA", groupvars=c('Setting'), na.rm=T)
freq_summ$Setting = factor(freq_summ$Setting, levels = c("control", "high", 'low'))
pd = position_dodge(0.1)
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('Total bioturbbation area'~(cm^2))))+
theme(
plot.title = element_blank(),
panel.background = element_blank(),
legend.position="none",
axis.title.x = element_text(size = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold", colour=c('#3399FF', '#FF9900', '#FF0000')),
axis.text.y = element_text(size = 28, 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)))
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=30, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
lmPlot
ggsave('ldm_neo.png', lmPlot, units = 'in', width = 12, height = 9,dpi=600)
lmPlot = lmPlot +
geom_segment(aes(x=1, xend=3, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=10, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
lmPlot
ggsave('ldm_neo.png', lmPlot, units = 'in', width = 12, height = 9,dpi=600)
lmPlot = lmPlot +
geom_segment(aes(x=1, xend=3, y=-Inf, yend=-Inf), size = 2.5, color = 'black') +
geom_segment(aes(y=0, yend=10, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
lmPlot
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=10, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
lmPlot
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('Total bioturbbation area'~(cm^2))))+
theme(
plot.title = element_blank(),
panel.background = element_blank(),
legend.position="none",
axis.title.x = element_text(size = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold", colour=c('#3399FF', '#FF9900', '#FF0000')),
axis.text.y = element_text(size = 28, 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)))
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=8, x=-Inf, xend=-Inf), size = 2.5, color = 'black')
lmPlot
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('Total bioturbbation area'~(cm^2))))+
theme(
plot.title = element_blank(),
panel.background = element_blank(),
legend.position="none",
axis.title.x = element_text(size = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold", colour=c('#3399FF', '#FF9900', '#FF0000')),
axis.text.y = element_text(size = 28, 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)))
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)
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)
}
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')
lumino_freq$bioturbated_area = (lumino_freq$area_adjusted_perc)/100 * 35.23865
lumino_coresum = aggregate(bioturbated_area ~ Day + Setting + Tank + Replicat, data = lumino_freq, FUN = sum)
lumino_coresum$bioturbated_areaDA = lumino_coresum$bioturbated_area/5
# statistics ANOVA
lumino_coresum$Setting = as.factor(lumino_coresum$Setting)
res.aov <- aov(bioturbated_areaDA ~ Setting , data = lumino_coresum)
summary(res.aov)
TukeyHSD(res.aov)
summary(glht(res.aov, linfct = mcp(Setting = "Tukey")))
res.aov1 = kruskal.test(bioturbated_areaDA ~ Setting , data = lumino_coresum)
summary(res.aov1)
pwt = pairwise.wilcox.test(lumino_coresum$bioturbated_areaDA, lumino_coresum$Setting,
p.adjust.method = "BH")
wilcox_effsize(bioturbated_areaDA ~ Setting, data = lumino_coresum)
freq_summ = summarySE(lumino_coresum, measurevar="bioturbated_areaDA", groupvars=c('Setting'), na.rm=T)
freq_summ$Setting = factor(freq_summ$Setting, levels = c("control", "high", 'low'))
pd = position_dodge(0.1)
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 = 30, face="bold"),
axis.title.y = element_text(size = 30, face="bold"),
axis.text.x = element_text(size =28, face="bold", colour=c('#3399FF', '#FF9900', '#FF0000')),
axis.text.y = element_text(size = 28, 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)))
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)
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 = 25, face="bold"),
axis.title.y = element_text(size = 25, face="bold"),
axis.text.x = element_text(size =28, face="bold", colour=c('#3399FF', '#FF9900', '#FF0000')),
axis.text.y = element_text(size = 28, 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)))
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)
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)))
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)
pwt = pairwise.wilcox.test(lumino_coresum$bioturbated_areaDA, lumino_coresum$Setting,
p.adjust.method = "BH")
wilcox_effsize(bioturbated_areaDA ~ Setting, data = lumino_coresum)
res.aov1 = kruskal.test(bioturbated_areaDA ~ Setting , data = lumino_coresum)
summary(res.aov1)
res.aov1
wilcox_effsize(bioturbated_areaDA ~ Setting, data = lumino_coresum)
wilcox_effsize(bioturbated_areaDA ~ Setting, data = lumino_coresum)
View(pwt)
