#loading necessary libraries - install any of these that you do not already have
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 4.4.0
library(plyr)
## Warning: package 'plyr' was built under R version 4.3.3
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'ggplot2' was built under R version 4.4.2
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'purrr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'stringr' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.4.2
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::arrange() masks plyr::arrange()
## ✖ purrr::compact() masks plyr::compact()
## ✖ dplyr::count() masks plyr::count()
## ✖ dplyr::desc() masks plyr::desc()
## ✖ dplyr::failwith() masks plyr::failwith()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::id() masks plyr::id()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::mutate() masks plyr::mutate()
## ✖ dplyr::rename() masks plyr::rename()
## ✖ dplyr::summarise() masks plyr::summarise()
## ✖ dplyr::summarize() masks plyr::summarize()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.4.2
##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:plyr':
##
## mutate
library(rstatix)
## Warning: package 'rstatix' was built under R version 4.3.3
##
## Attaching package: 'rstatix'
##
## The following objects are masked from 'package:plyr':
##
## desc, mutate
##
## The following object is masked from 'package:stats':
##
## filter
library(ggsignif)
## Warning: package 'ggsignif' was built under R version 4.4.2
library(glmmTMB)
## Warning: package 'glmmTMB' was built under R version 4.3.3
## Warning in checkMatrixPackageVersion(getOption("TMB.check.Matrix", TRUE)): Package version inconsistency detected.
## TMB was built with Matrix ABI version 1
## Current Matrix ABI version is 0
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## Warning in check_dep_version(): ABI version mismatch:
## lme4 was built with Matrix ABI version 1
## Current Matrix ABI version is 0
## Please re-install lme4 from source or restore original 'Matrix' package
library(lme4)
## Warning: package 'lme4' was built under R version 4.3.3
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(MASS)
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:rstatix':
##
## select
##
## The following object is masked from 'package:dplyr':
##
## select
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 4.3.3
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(plotrix)
## Warning: package 'plotrix' was built under R version 4.3.2
library(reprex)
## Warning: package 'reprex' was built under R version 4.3.3
library(rcartocolor)
## Warning: package 'rcartocolor' was built under R version 4.3.3
library(paletteer)
## Warning: package 'paletteer' was built under R version 4.3.3
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(tinytex)
## Warning: package 'tinytex' was built under R version 4.3.3
#loading in solo crab predator trial CSV - alter path to match the path for where you've downloaded the Solo Crab Pred csv on your own device
prf <- read.csv("C:/Users/Cassian/Downloads/SURF Control Proof - Solo Crab Experiments.csv", sep = ",", header = TRUE)
View(prf)
#Performing a Wilcoxon ranked-sum test to see if there are significant differences between males and females for the predator responses
compare_means(Time_to_Burrow ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow Control Predator 0.434 0.43 0.43 ns Wilcoxon
compare_means(Aggression ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aggression Control Predator 0.0433 0.043 0.043 * Wilcoxon
compare_means(Flee ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Flee Control Predator 0.000724 0.00072 0.00072 *** Wilcoxon
compare_means(Pace_Slow ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Pace_Slow Control Predator 0.0227 0.023 0.023 * Wilcoxon
compare_means(Burrow_Dur ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur Control Predator 0.452 0.45 0.45 ns Wilcoxon
compare_means(Freeze_Dur ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Freeze_Dur Control Predator 0.000000134 0.00000013 1.3e-07 **** Wilcoxon
compare_means(Water_Dur ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Water_Dur Control Predator 0.0000248 0.000025 2.5e-05 **** Wilcoxon
compare_means(Climb_Dur ~ Treatment, data = prf, paired = FALSE)
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Climb_Dur Control Predator 0.0000922 0.000092 9.2e-05 **** Wilcoxon
#loading in solo crab predator (SCP) trial CSV - alter path to match the path for where you've downloaded the Solo Crab Pred csv on your own device
scp <- read.csv("C:/Users/Cassian/OneDrive/Documents/Ocean Environment/SURF Data Analysis - Solo Crab Pred.csv", sep = ",", header = TRUE)
View(scp)
#creating object for SCP males
male <- scp %>% filter(Sex == "M")
View(male)
#creating object for SCP females
female <- scp %>% filter(Sex == "F")
View(female)
#counting how many crabs showed each response - solo trials
#looking at male crabs only
magg <- male %>% filter(Aggression == 0)
View(magg)
#count: 21/25 did not show agg
mflee <- male %>% filter(Flee == 0)
View(mflee)
#count: 14/25 did not flee
mpace <- male %>% filter(Pace_Slow == 0)
View(mpace)
#count: 19/25 did not slow
mfreeze <- male %>% filter(Freeze_Dur == 0)
View(mfreeze)
#count: 3/25 did not freeze
mwater <- male %>% filter(Water_Dur == 0)
View(mwater)
#count: 15/25 did not submerge
#looking at female crabs only
fflee <- female %>% filter(Flee == 0)
View(fflee)
#count: 20/25 did not flee
fpace <- female %>% filter(Pace_Slow == 0)
View(fpace)
#count: 22/25 did not slow
ffreeze <- female %>% filter(Freeze_Dur == 0)
View(ffreeze)
#count: 7/25 did not freeze
fwater <- female %>% filter(Water_Dur == 0)
View(fwater)
#count: 16/25 did not submerge
#Performing a Wilcoxon ranked-sum test to see if there are significant differences between males and females for the predator responses
compare_means(Time_to_Burrow ~ Sex, data = scp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow M F 0.0455 0.046 0.046 * Wilcoxon
compare_means(Aggression ~ Sex, data = scp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aggression M F 0.0410 0.041 0.041 * Wilcoxon
compare_means(Flee ~ Sex, data = scp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Flee M F 0.0453 0.045 0.045 * Wilcoxon
compare_means(Pace_Slow ~ Sex, data = scp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Pace_Slow M F 0.324 0.32 0.32 ns Wilcoxon
compare_means(Burrow_Dur ~ Sex, data = scp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur M F 0.0345 0.034 0.034 * Wilcoxon
compare_means(Freeze_Dur ~ Sex, data = scp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Freeze_Dur M F 0.0194 0.019 0.019 * Wilcoxon
compare_means(Water_Dur ~ Sex, data = scp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Water_Dur M F 0.491 0.49 0.49 ns Wilcoxon
#calculating mean/SD for freeze durations and creating a dataframe to use in graph creation
sd(male$Freeze_Dur)
## [1] 23.81577
mean(male$Freeze_Dur)
## [1] 27.452
sd(female$Freeze_Dur)
## [1] 14.06097
mean(female$Freeze_Dur)
## [1] 12.968
freeze_df <- scp %>%
group_by(Sex) %>%
summarise(mean= mean(Freeze_Dur),
sd = std.error(Freeze_Dur)) #use standard error for the SD in the dataframe
View(freeze_df) #check that dataframe was created correctly
#making a bar graph of mean freeze times
freezegraph <- ggplot(freeze_df, aes(x = Sex, y = mean, fill=factor(Sex))) +
geom_bar(stat ="identity", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=12), axis.title.y = element_text(size=12),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) + #set text sizes
labs (x="Sex", y="Average Freeze Duration (s)", fill="Sex") +
ggtitle("Freeze Duration") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#D95F02","#D691C1FF")) + #add colors, first is for F, last is for M
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.2, colour="black", alpha=0.9, size=1) + #add error bars
geom_segment(aes(x="F", xend="M", y=35, yend=35)) + #add line for significance label
annotate("text", x=c(1.5),y=c(36.5),label=c("P=0.019"),size=4) #inclusion of significance label
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
freezegraph
ggsave("freezegraph.tiff") #save the graph
## Saving 7 x 5 in image
#calculate mean/SD for fleeing, make dataframe for graph creation
mean(female$Flee)
## [1] 0.24
sd(female$Flee)
## [1] 0.5228129
mean(male$Flee)
## [1] 0.76
sd(male$Flee)
## [1] 1.011599
flee_df <- scp %>% #making dataframe
group_by(Sex) %>%
summarise(mean= mean(Flee),
sd = std.error(Flee)) #use standard error
View(flee_df) #check to ensure dataframe generated correctly
#create graph for fleeing
fleegraph <- ggplot(flee_df, aes(x = Sex, y = mean, fill=factor(Sex))) +
geom_bar(stat ="identity", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=12), axis.title.y = element_text(size=12),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
labs (x="Sex", y="Average # of Fleeing Occurrences", fill="Sex") +
ggtitle("Fleeing") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#D95F02","#D691C1FF")) +
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.2, colour="black", alpha=0.9, size=1) +
geom_segment(aes(x="F", xend="M", y=1, yend=1)) +
annotate("text", x=c(1.5),y=c(1.05),label=c("P=0.045"),size=4)
fleegraph
ggsave("fleeplot.tiff")
## Saving 7 x 5 in image
#combine first two figures to make a more compact figure
predplotfig2 <- ggarrange(freezegraph, fleegraph,
labels = c("A", "B"), hjust=-2,vjust = 2,
ncol = 2, nrow =1,
widths = c(10,10),
heights=c(10, 10))
predplotfig2
ggsave("predplotfig2.tiff")
## Saving 7 x 5 in image
#loading in CSV of all solo crab data to analyze burrowing
bcp <- read.csv("C:/Users/Cassian/Downloads/SURF Control Proof - Solo Crab Experiments.csv", sep = ",", header = TRUE)
View(bcp)
#creating object for BCP males
mbcp <- bcp %>% filter(Sex == "M")
View(mbcp)
#creating object for BCP females
fbcp <- bcp %>% filter(Sex == "F")
View(fbcp)
#creating object for crabs that burrowed to count how many crabs burrowed
bbcp <- bcp %>% filter(Burrowed == "Y")
View(bbcp)
#testing significance between the number of female and male crabs that burrowed
binom.test(11,38, 0.5)
##
## Exact binomial test
##
## data: 11 and 38
## number of successes = 11, number of trials = 38, p-value = 0.01385
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.1542463 0.4590321
## sample estimates:
## probability of success
## 0.2894737
compare_means(Time_to_Burrow ~ Sex, data = bcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow M F 0.00150 0.0015 0.0015 ** Wilcoxon
compare_means(Burrow_Dur ~ Sex, data = bcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur M F 0.00160 0.0016 0.0016 ** Wilcoxon
#calculating averages for time to burrow and making a dataframe for graph creation
mean(fbcp$Time_to_Burrow)
## [1] 427.1
sd(fbcp$Time_to_Burrow)
## [1] 207.5107
mean(mbcp$Time_to_Burrow)
## [1] 531.16
sd(mbcp$Time_to_Burrow)
## [1] 150.5627
ttb_df <- bcp %>%
group_by(Sex) %>%
summarise(mean= mean(Time_to_Burrow),
sd = std.error(Time_to_Burrow))
View(ttb_df)
#creating graph for burrow latency
ttbgraph <- ggplot(ttb_df, aes(x = Sex, y = mean, fill=factor(Sex))) +
geom_bar(stat ="identity", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=12), axis.title.y = element_text(size=12),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
labs (x="Sex", y="Average Burrow Latency (s)", fill="Sex") +
ggtitle("Burrow Latency") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#D95F02","#D691C1FF")) +
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.2, colour="black", alpha=0.9, size=1) +
geom_segment(aes(x="F", xend="M", y=595, yend=595)) +
annotate("text", x=c(1.5),y=c(615),label=c("P=1.50e-3"),size=4)
ttbgraph
ggsave("ttbgraph.tiff")
## Saving 7 x 5 in image
#calculating averages for burrow duration and making dataframe for graph creation
mean(mbcp$Burrow_Dur)
## [1] 61.65714
sd(mbcp$Burrow_Dur)
## [1] 143.7935
mean(fbcp$Burrow_Dur)
## [1] 174.9667
sd(fbcp$Burrow_Dur)
## [1] 210.2323
burrowtime_df <- bcp %>%
group_by(Sex) %>%
summarise(mean= mean(Burrow_Dur),
sd = std.error(Burrow_Dur))
View(burrowtime_df)
#create graph for burrow duration
burrowgraph <- ggplot(burrowtime_df, aes(x = Sex, y = mean, fill=factor(Sex))) +
geom_bar(stat ="identity", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=12), axis.title.y = element_text(size=12),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
labs (x="Sex", y="Average Time in Burrow (s)", fill="Sex") +
ggtitle("Burrowing Duration") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#D95F02","#D691C1FF")) +
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.2, colour="black", alpha=0.9, size=1) +
geom_segment(aes(x="F", xend="M", y=215, yend=215)) +
annotate("text", x=c(1.5),y=c(225),label=c("P=1.60e-3"),size=4)
burrowgraph
ggsave("burrowgraph.tiff")
## Saving 7 x 5 in image
#combine first two figures to make a more compact figure
predplotfig3 <- ggarrange(ttbgraph,burrowgraph,
labels = c("A", "B"), hjust=-3,vjust = 2,
ncol = 2, nrow =1,
widths = c(15,15),
heights=c(5, 5))
predplotfig3
ggsave("predplotfig3pt2.tiff")
## Saving 7 x 5 in image
paird <- read.csv("C:/Users/Cassian/Downloads/Dual Crab Pred - Dual Crab Experiments.csv", sep = ",", header = TRUE)
View(paird)
#counting how many crabs showed each response - paired trials
fpair <- paird %>% filter(Sex == "F")
View(fpair)
mpair <- paird %>% filter(Sex == "M")
#looking at male crabs only
mpagg <- mpair %>% filter(Aggression == 0)
View(mpagg)
#count: 19/30 did not show agg
mpflee <- mpair %>% filter(Flee == 0)
View(mpflee)
#count: 23/30 did not flee
mppace <- mpair %>% filter(Pace_Slow == 0)
View(mppace)
#count: 15/30 did not slow
mpfreeze <- mpair %>% filter(Freeze_Dur == 0)
View(mpfreeze)
#count: 7/30 did not freeze
mpwater <- mpair %>% filter(Water_Dur == 0)
View(mpwater)
#count: 25/30 did not submerge
#looking at female crabs only
fpflee <- fpair %>% filter(Flee == 0)
View(fpflee)
#count: 28/30 did not flee
fppace <- fpair %>% filter(Pace_Slow == 0)
View(fppace)
#count: 19/30 did not slow
fpfreeze <- fpair %>% filter(Freeze_Dur == 0)
View(fpfreeze)
#count: 9/30 did not freeze
fpwater <- fpair %>% filter(Water_Dur == 0)
View(fpwater)
#count: 26/30 did not submerge
#load in CSV for paired AND solo crabs - replace path with path for wherever you saved the mixed crab pred trials CSV on your device
#note: this CSV includes all paired and solo crab predator trial data and allows comparison between paired and solo trials
pcp <- read.csv("C:/Users/Cassian/Downloads/Dual Crab Pred - Mixed Crab Pred Trials (1).csv", sep = ",", header = TRUE)
View(pcp)
#Performing Wilcoxon ranked-sum tests on the predator response variables, using a Benjamini-Hochburg adjustment. Note that some of these will appear significant, and this significance will disappear upon a second run with data filtered down to include only one sex
compare_means(Time_to_Burrow ~ Trial_Type, data = pcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 6 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow FF MF 0.617 0.74 0.62 ns Wilcoxon
## 2 Time_to_Burrow FF MM 0.200 0.6 0.20 ns Wilcoxon
## 3 Time_to_Burrow FF Solo 0.951 0.95 0.95 ns Wilcoxon
## 4 Time_to_Burrow MF MM 0.536 0.74 0.54 ns Wilcoxon
## 5 Time_to_Burrow MF Solo 0.570 0.74 0.57 ns Wilcoxon
## 6 Time_to_Burrow MM Solo 0.187 0.6 0.19 ns Wilcoxon
compare_means(Aggression ~ Trial_Type, data = pcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 6 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aggression FF MF 0.163 0.24 0.16259 ns Wilcoxon
## 2 Aggression FF MM 0.000926 0.0028 0.00093 *** Wilcoxon
## 3 Aggression FF Solo 0.202 0.24 0.20154 ns Wilcoxon
## 4 Aggression MF MM 0.0139 0.028 0.01392 * Wilcoxon
## 5 Aggression MF Solo 0.758 0.76 0.75798 ns Wilcoxon
## 6 Aggression MM Solo 0.000141 0.00084 0.00014 *** Wilcoxon
compare_means(Flee ~ Trial_Type, data = pcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 6 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Flee FF MF 1 1 1.000 ns Wilcoxon
## 2 Flee FF MM 0.225 0.34 0.225 ns Wilcoxon
## 3 Flee FF Solo 0.0479 0.14 0.048 * Wilcoxon
## 4 Flee MF MM 0.225 0.34 0.225 ns Wilcoxon
## 5 Flee MF Solo 0.0479 0.14 0.048 * Wilcoxon
## 6 Flee MM Solo 0.401 0.48 0.401 ns Wilcoxon
compare_means(Pace_Slow ~ Trial_Type, data = pcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 6 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Pace_Slow FF MF 0.331 0.43 0.331 ns Wilcoxon
## 2 Pace_Slow FF MM 0.859 0.86 0.859 ns Wilcoxon
## 3 Pace_Slow FF Solo 0.0209 0.063 0.021 * Wilcoxon
## 4 Pace_Slow MF MM 0.359 0.43 0.359 ns Wilcoxon
## 5 Pace_Slow MF Solo 0.184 0.37 0.184 ns Wilcoxon
## 6 Pace_Slow MM Solo 0.0146 0.063 0.015 * Wilcoxon
compare_means(Burrow_Dur ~ Trial_Type, data = pcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 6 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur FF MF 0.689 0.83 0.69 ns Wilcoxon
## 2 Burrow_Dur FF MM 0.227 0.68 0.23 ns Wilcoxon
## 3 Burrow_Dur FF Solo 0.927 0.93 0.93 ns Wilcoxon
## 4 Burrow_Dur MF MM 0.510 0.83 0.51 ns Wilcoxon
## 5 Burrow_Dur MF Solo 0.636 0.83 0.64 ns Wilcoxon
## 6 Burrow_Dur MM Solo 0.181 0.68 0.18 ns Wilcoxon
compare_means(Freeze_Dur ~ Trial_Type, data = pcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 6 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Freeze_Dur FF MF 0.869 0.87 0.87 ns Wilcoxon
## 2 Freeze_Dur FF MM 0.530 0.81 0.53 ns Wilcoxon
## 3 Freeze_Dur FF Solo 0.292 0.81 0.29 ns Wilcoxon
## 4 Freeze_Dur MF MM 0.540 0.81 0.54 ns Wilcoxon
## 5 Freeze_Dur MF Solo 0.334 0.81 0.33 ns Wilcoxon
## 6 Freeze_Dur MM Solo 0.855 0.87 0.86 ns Wilcoxon
compare_means(Water_Dur ~ Trial_Type, data = pcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 6 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Water_Dur FF MF 0.0895 0.21 0.0895 ns Wilcoxon
## 2 Water_Dur FF MM 0.299 0.36 0.2988 ns Wilcoxon
## 3 Water_Dur FF Solo 0.00908 0.054 0.0091 ** Wilcoxon
## 4 Water_Dur MF MM 0.510 0.51 0.5104 ns Wilcoxon
## 5 Water_Dur MF Solo 0.243 0.36 0.2430 ns Wilcoxon
## 6 Water_Dur MM Solo 0.104 0.21 0.1038 ns Wilcoxon
#creating objects for male paired crabs and female paired crabs
m_data <- pcp %>% filter(Sex== "M")
View(m_data)
f_data <- pcp %>% filter(Sex== "F")
View(f_data)
#Performing Wilcoxon ranked-sum tests on the predator response variables for males only, using a Benjamini-Hochburg adjustment
compare_means(Time_to_Burrow ~ Trial_Type, data = m_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow MF MM 0.489 0.73 0.49 ns Wilcoxon
## 2 Time_to_Burrow MF Solo 0.472 0.73 0.47 ns Wilcoxon
## 3 Time_to_Burrow MM Solo 0.797 0.8 0.80 ns Wilcoxon
compare_means(Aggression ~ Trial_Type, data = m_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aggression MF MM 0.177 0.27 0.177 ns Wilcoxon
## 2 Aggression MF Solo 0.717 0.72 0.717 ns Wilcoxon
## 3 Aggression MM Solo 0.0139 0.042 0.014 * Wilcoxon
compare_means(Flee ~ Trial_Type, data = m_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Flee MF MM 0.787 0.79 0.787 ns Wilcoxon
## 2 Flee MF Solo 0.133 0.2 0.133 ns Wilcoxon
## 3 Flee MM Solo 0.0953 0.2 0.095 ns Wilcoxon
compare_means(Pace_Slow ~ Trial_Type, data = m_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Pace_Slow MF MM 0.732 0.73 0.73 ns Wilcoxon
## 2 Pace_Slow MF Solo 0.273 0.41 0.27 ns Wilcoxon
## 3 Pace_Slow MM Solo 0.115 0.35 0.12 ns Wilcoxon
compare_means(Burrow_Dur ~ Trial_Type, data = m_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur MF MM 0.450 0.67 0.45 ns Wilcoxon
## 2 Burrow_Dur MF Solo 0.413 0.67 0.41 ns Wilcoxon
## 3 Burrow_Dur MM Solo 0.797 0.8 0.80 ns Wilcoxon
compare_means(Freeze_Dur ~ Trial_Type, data = m_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Freeze_Dur MF MM 0.493 0.49 0.493 ns Wilcoxon
## 2 Freeze_Dur MF Solo 0.0853 0.26 0.085 ns Wilcoxon
## 3 Freeze_Dur MM Solo 0.181 0.27 0.181 ns Wilcoxon
compare_means(Water_Dur ~ Trial_Type, data = m_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Water_Dur MF MM 0.760 0.76 0.760 ns Wilcoxon
## 2 Water_Dur MF Solo 0.253 0.38 0.253 ns Wilcoxon
## 3 Water_Dur MM Solo 0.0868 0.26 0.087 ns Wilcoxon
#Performing Wilcoxon ranked-sum tests on the predator response variables for females only, using a Benjamini-Hochburg adjustment
compare_means(Time_to_Burrow ~ Trial_Type, data = f_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow FF MF 0.533 0.53 0.53 ns Wilcoxon
## 2 Time_to_Burrow FF Solo 0.446 0.53 0.45 ns Wilcoxon
## 3 Time_to_Burrow MF Solo 0.180 0.53 0.18 ns Wilcoxon
compare_means(Flee ~ Trial_Type, data = f_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Flee FF MF 0.334 0.35 0.33 ns Wilcoxon
## 2 Flee FF Solo 0.354 0.35 0.35 ns Wilcoxon
## 3 Flee MF Solo 0.141 0.35 0.14 ns Wilcoxon
compare_means(Pace_Slow ~ Trial_Type, data = f_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Pace_Slow FF MF 0.267 0.4 0.27 ns Wilcoxon
## 2 Pace_Slow FF Solo 0.0200 0.06 0.02 * Wilcoxon
## 3 Pace_Slow MF Solo 0.568 0.57 0.57 ns Wilcoxon
compare_means(Burrow_Dur ~ Trial_Type, data = f_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur FF MF 0.533 0.53 0.53 ns Wilcoxon
## 2 Burrow_Dur FF Solo 0.402 0.53 0.40 ns Wilcoxon
## 3 Burrow_Dur MF Solo 0.194 0.53 0.19 ns Wilcoxon
compare_means(Freeze_Dur ~ Trial_Type, data = f_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Freeze_Dur FF MF 0.772 0.91 0.77 ns Wilcoxon
## 2 Freeze_Dur FF Solo 0.908 0.91 0.91 ns Wilcoxon
## 3 Freeze_Dur MF Solo 0.593 0.91 0.59 ns Wilcoxon
compare_means(Water_Dur ~ Trial_Type, data = f_data, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 3 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Water_Dur FF MF 0.0740 0.11 0.074 ns Wilcoxon
## 2 Water_Dur FF Solo 0.0211 0.063 0.021 * Wilcoxon
## 3 Water_Dur MF Solo 0.746 0.75 0.746 ns Wilcoxon
#creating dataframes for each trial type and sex
pmm = filter(pcp, Trial_Type == "MM")
pff = filter(pcp, Trial_Type == "FF")
phetm = filter(m_data, Trial_Type == "MF")
phetf = filter(f_data, Trial_Type == "MF")
pms = filter(m_data, Trial_Type == "Solo")
pfs = filter(f_data, Trial_Type == "Solo")
#calculating means for significantly different behaviors to compare them between trial types
mean(pmm$Aggression)
## [1] 1.55
mean(pms$Aggression)
## [1] 0.16
mean(phetm$Aggression)
## [1] 0.5
mean(pff$Pace_Slow)
## [1] 0.75
mean(pfs$Pace_Slow)
## [1] 0.52
mean(phetf$Pace_Slow)
## [1] 0.4
mean(pff$Water_Dur)
## [1] 0.81
mean(pfs$Water_Dur)
## [1] 3.602
mean(phetf$Water_Dur)
## [1] 2.5
#creating dataframe for female submergence duration, to use in graph creation
water_df <- f_data %>%
group_by(Trial_Type) %>%
summarise(mean= mean(Water_Dur),
sd = std.error(Water_Dur))
View(water_df)
#creating a graph for female submergence duration, comparing between trial types
watergraph <- ggplot(water_df, aes(x = Trial_Type, y = mean, fill=factor(Trial_Type))) +
geom_bar(stat ="identity", color="black", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=14), axis.title.y = element_text(size=14),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
labs (x="Trial Type", y="Average Submerge Duration (s)", fill="Trial Type") +
ggtitle("Average Water Duration - Females") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#6C568CFF","#BFCDD9FF", "#E394BBFF")) + #these colors are for FF pairs, MF pairs, and solo females respectively
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.1, colour="black", alpha=0.9, size=0.5) +
geom_segment(aes(x="FF", xend="Solo", y=6.2, yend=6.2)) +
annotate("text", x=c(2),y=c(6.5),label=c("P=0.021"),size=4) +
geom_segment(aes(x="FF", xend="MF", y=5.7, yend=5.7)) +
annotate("text", x=c(1.5),y=c(6),label=c("P=0.074"),size=4) +
geom_segment(aes(x="MF", xend="Solo", y=5.3, yend=5.3)) +
annotate("text", x=c(2.5),y=c(5.6),label=c("P=0.746"),size=4)
watergraph
ggsave("waterplot.tiff")
## Saving 7 x 5 in image
#creating dataframe for female pace slowing behavior to use for graph creation
pace_df <- f_data %>%
group_by(Trial_Type) %>%
summarise(mean= mean(Pace_Slow),
sd = std.error(Pace_Slow))
View(pace_df)
#creating graph displaying differences in female pace slowing behavior between different trial types
pacegraph <- ggplot(pace_df, aes(x = Trial_Type, y = mean, fill=factor(Trial_Type))) +
geom_bar(stat ="identity", color="black", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=14), axis.title.y = element_text(size=14),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
labs (x="Trial Type", y="Average # of Pace Slow Events", fill="Trial Type") +
ggtitle("Average Pace Slow Events - Females") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#6C568CFF","#BFCDD9FF", "#E394BBFF")) +
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.1, colour="black", alpha=0.9, size=0.5) +
geom_segment(aes(x="FF", xend="Solo", y=1.6, yend=1.6)) +
annotate("text", x=c(2),y=c(1.68),label=c("P=0.02"),size=4) +
geom_segment(aes(x="FF", xend="MF", y=1.45, yend=1.45)) +
annotate("text", x=c(1.53),y=c(1.53),label=c("P=0.267"),size=4) +
geom_segment(aes(x="MF", xend="Solo", y=1.3, yend=1.3)) +
annotate("text", x=c(2.5),y=c(1.38),label=c("P=0.568"),size=4)
pacegraph
ggsave("pace slow plot.tiff")
## Saving 7 x 5 in image
#creating dataframe for male aggressive behavior to use in graph creation
agg_df <- m_data %>%
group_by(Trial_Type) %>%
summarise(mean= mean(Aggression),
sd = std.error(Aggression))
View(agg_df)
#graph of male aggression to compare between trial types
agggraph <- ggplot(agg_df, aes(x = Trial_Type, y = mean, fill=factor(Trial_Type))) +
geom_bar(stat ="identity", color="black", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=14), axis.title.y = element_text(size=14),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
scale_x_discrete(limits = c("MM","MF","Solo")) +
labs (x="Trial Type", y="Average # of Aggression Events", fill="Trial Type") +
ggtitle("Average Aggression Events") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(breaks=c("MM","MF","Solo"),values=c("#5D761EFF","#BFCDD9FF","#E394BBFF")) + #these colors are for MM pairs, MF pairs, and solo males respectively
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.1, colour="black", alpha=0.9, size=0.5) +
geom_segment(aes(x="MM", xend="Solo", y=2.6, yend=2.6)) +
annotate("text", x=c(2),y=c(2.75),label=c("P=0.014"),size=4) +
geom_segment(aes(x="MM", xend="MF", y=2.2, yend=2.2)) +
annotate("text", x=c(1.5),y=c(2.35),label=c("P=0.177"),size=4) +
geom_segment(aes(x="MF", xend="Solo", y=2, yend=2)) +
annotate("text", x=c(2.5),y=c(2.15),label=c("P=0.717"),size=4)
agggraph
ggsave("agg male plot.tiff")
## Saving 7 x 5 in image
#combining the first two paired graphs
predplotfig4 <- ggarrange(watergraph, pacegraph,
labels = c("A", "B"), hjust=-1,vjust = 3,
ncol = 2, nrow =1,
widths = c(10,10),
heights=c(10, 10))
predplotfig4
ggsave("predplotfig4.tiff")
## Saving 7 x 5 in image
#annotating the last paired graph to make it a figure
predplotfig4pt2 <- ggarrange(agggraph,
labels = c("C"), hjust=-1,vjust = 3,
ncol = 1, nrow =1,
widths = c(10),
heights=c(10))
predplotfig4pt2
ggsave("predplotfig4pt2.tiff")
## Saving 7 x 5 in image
#loading in CSV of all crab data to analyze burrowing
pbcp <- read.csv("C:/Users/Cassian/Downloads/SURF Control Proof - Mixed Crab Proofs (2).csv", sep = ",", header = TRUE)
View(pbcp)
#creating object for all males
mbpp <- pbcp %>% filter(Sex == "M")
View(mbpp)
#creating object for all females
fbpp <- pbcp %>% filter(Sex == "F")
View(fbpp)
#creating object for paired crabs only
pairb <- pbcp %>% filter(Paired == "Y")
View(pairb)
#creating object for solo crabs only
solob <- pbcp %>% filter(Paired == "N")
View(solob)
#creating object for paired males
mpp <- pairb %>% filter(Sex == "M")
View(mpp)
#creating object for paired females
fpp <- pairb %>% filter(Sex == "F")
View(fpp)
#creating object for crabs that burrowed to count how many crabs burrowed
bpp <- pairb %>% filter(Burrowed == "Y")
View(bpp)
#testing significance between the number of female and male crabs that burrowed
binom.test(17,31, 0.5)
##
## Exact binomial test
##
## data: 17 and 31
## number of successes = 17, number of trials = 31, p-value = 0.7201
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.3603423 0.7268350
## sample estimates:
## probability of success
## 0.5483871
#testing significance between the number of solo and paired crabs that burrowed
binom.test(31,69, 0.5)
##
## Exact binomial test
##
## data: 31 and 69
## number of successes = 31, number of trials = 69, p-value = 0.4704
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.3292298 0.5738142
## sample estimates:
## probability of success
## 0.4492754
#testing if there are significant differences between paired and solo crabs in terms of burrowing behavior
compare_means(Time_to_Burrow ~ Paired, data = pbcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow N Y 0.0990 0.099 0.099 ns Wilcoxon
compare_means(Burrow_Dur ~ Paired, data = pbcp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur N Y 0.117 0.12 0.12 ns Wilcoxon
#testing if there are significant differences between male and female burrowing behavior for paired crabs
compare_means(Time_to_Burrow ~ Sex, data = pairb, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow M F 0.680 0.68 0.68 ns Wilcoxon
compare_means(Burrow_Dur ~ Sex, data = pairb, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur M F 0.472 0.47 0.47 ns Wilcoxon
#testing if there are significant differences between paired and solo crabs of each sex in terms of burrowing behavior
compare_means(Time_to_Burrow ~ Paired, data = fbpp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow N Y 0.0104 0.01 0.01 * Wilcoxon
compare_means(Burrow_Dur ~ Paired, data = fbpp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur N Y 0.0175 0.017 0.017 * Wilcoxon
compare_means(Time_to_Burrow ~ Paired, data = mbpp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Time_to_Burrow N Y 0.686 0.69 0.69 ns Wilcoxon
compare_means(Burrow_Dur ~ Paired, data = mbpp, method = "wilcox.test", paired = FALSE, p.adjust.method = "BH")
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Burrow_Dur N Y 0.829 0.83 0.83 ns Wilcoxon
#calculating averages for burrow latency
mean(fpp$Time_to_Burrow)
## [1] 508.1833
sd(fpp$Time_to_Burrow)
## [1] 175.2288
mean(fbcp$Time_to_Burrow)
## [1] 427.1
sd(fbcp$Time_to_Burrow)
## [1] 207.5107
mean(mpp$Time_to_Burrow)
## [1] 504.15
sd(mpp$Time_to_Burrow)
## [1] 184.5238
mean(mbcp$Time_to_Burrow)
## [1] 531.16
sd(mbcp$Time_to_Burrow)
## [1] 150.5627
mean(solob$Time_to_Burrow)
## [1] 479.13
sd(solob$Time_to_Burrow)
## [1] 187.7962
mean(pairb$Time_to_Burrow)
## [1] 506.1667
sd(pairb$Time_to_Burrow)
## [1] 179.1902
latent_df <- fbpp %>%
group_by(Paired) %>%
summarise(mean= mean(Time_to_Burrow),
sd = std.error(Time_to_Burrow))
View(latent_df)
#creating graph for burrow latency between solo and paired females
latentgraph <- ggplot(latent_df, aes(x = Paired, y = mean, fill=factor(Paired))) +
geom_bar(stat ="identity", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=12), axis.title.y = element_text(size=12),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
labs (x="Paired?", y="Average Burrow Latency (s)", fill="Paired") +
ggtitle("Burrow Latency - Females") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#E394BBFF","#EEAD0E")) +
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.2, colour="black", alpha=0.9, size=1) +
geom_segment(aes(x="N", xend="Y", y=595, yend=595)) +
annotate("text", x=c(1.5),y=c(615),label=c("P=1.04e-2"),size=4)
latentgraph
#ggsave("ttbgraph.png")
#calculating averages for burrow duration
mean(mpp$Burrow_Dur)
## [1] 90.83978
sd(mpp$Burrow_Dur)
## [1] 181.9832
mean(mbcp$Burrow_Dur)
## [1] 61.65714
sd(mbcp$Burrow_Dur)
## [1] 143.7935
mean(fpp$Burrow_Dur)
## [1] 93.0216
sd(fpp$Burrow_Dur)
## [1] 168.593
mean(fbcp$Burrow_Dur)
## [1] 174.9667
sd(fbcp$Burrow_Dur)
## [1] 210.2323
mean(solob$Burrow_Dur)
## [1] 118.3119
sd(solob$Burrow_Dur)
## [1] 188.0203
mean(pairb$Burrow_Dur)
## [1] 91.93069
sd(pairb$Burrow_Dur)
## [1] 174.6808
burr_df <- fbpp %>%
group_by(Paired) %>%
summarise(mean= mean(Burrow_Dur),
sd = std.error(Burrow_Dur))
View(burr_df)
#create graph for burrow duration
burrowgraphfem <- ggplot(burr_df, aes(x = Paired, y = mean, fill=factor(Paired))) +
geom_bar(stat ="identity", position=position_dodge()) +
theme_classic() +
theme(axis.title.x = element_text(size=12), axis.title.y = element_text(size=12),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
labs (x="Paired?", y="Average Time in Burrow (s)", fill="Paired") +
ggtitle("Burrowing Duration - Females") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values=c("#E394BBFF","#EEAD0E")) +
theme(legend.position = "right") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), position = position_dodge(width=0.5), width=0.2, colour="black", alpha=0.9, size=1) +
geom_segment(aes(x="N", xend="Y", y=215, yend=215)) +
annotate("text", x=c(1.5),y=c(225),label=c("P=1.75e-2"),size=4)
burrowgraphfem
#ggsave("burrowgraph.png")
#combine first two figures to make a more compact figure
predplotfigx <- ggarrange(latentgraph,burrowgraphfem,
labels = c("A", "B"), hjust=-3,vjust = 2,
ncol = 2, nrow =1,
widths = c(15,15),
heights=c(5, 5))
predplotfigx
ggsave("predplotfigx.tiff")
## Saving 7 x 5 in image
#importing the CSV for both paired and solo crab movement metrics. Change the path to match the path for this CSV on your own device
msf <- read.csv("C:/Users/Cassian/Downloads/SURF Dual Site Fidelity - Mixed Crab Fidelity (2).csv", sep = ",", header = TRUE)
View(msf)
#creating a model to see if treatment and trial type have a mixed effect on movement
msfmodel1 <- glm(Quad_Transitions ~ Treatment * Trial_Type + Crab_ID, family = poisson, data = msf)
msfmodel1
##
## Call: glm(formula = Quad_Transitions ~ Treatment * Trial_Type + Crab_ID,
## family = poisson, data = msf)
##
## Coefficients:
## (Intercept) TreatmentPredator
## 1.7187332 -0.5464848
## Trial_TypeMF Trial_TypeMM
## -0.1724724 -0.0478567
## Trial_TypeSolo Crab_ID
## 0.0179714 0.0001037
## TreatmentPredator:Trial_TypeMF TreatmentPredator:Trial_TypeMM
## 0.4981229 0.7178019
## TreatmentPredator:Trial_TypeSolo
## 0.6230780
##
## Degrees of Freedom: 215 Total (i.e. Null); 207 Residual
## Null Deviance: 819.5
## Residual Deviance: 784 AIC: 1460
summary(msfmodel1)
##
## Call:
## glm(formula = Quad_Transitions ~ Treatment * Trial_Type + Crab_ID,
## family = poisson, data = msf)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7187332 0.1129024 15.223 < 2e-16 ***
## TreatmentPredator -0.5464848 0.1590525 -3.436 0.000591 ***
## Trial_TypeMF -0.1724724 0.1471927 -1.172 0.241299
## Trial_TypeMM -0.0478567 0.1387447 -0.345 0.730150
## Trial_TypeSolo 0.0179714 0.1551726 0.116 0.907799
## Crab_ID 0.0001037 0.0009192 0.113 0.910146
## TreatmentPredator:Trial_TypeMF 0.4981229 0.2194669 2.270 0.023226 *
## TreatmentPredator:Trial_TypeMM 0.7178019 0.2061919 3.481 0.000499 ***
## TreatmentPredator:Trial_TypeSolo 0.6230780 0.1787850 3.485 0.000492 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 819.52 on 215 degrees of freedom
## Residual deviance: 783.96 on 207 degrees of freedom
## AIC: 1460
##
## Number of Fisher Scoring iterations: 5
anova(msfmodel1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Quad_Transitions
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 215 819.52
## Treatment 1 0.0813 214 819.44 0.7755203
## Trial_Type 3 20.6722 211 798.76 0.0001231 ***
## Crab_ID 1 0.0013 210 798.76 0.9708228
## Treatment:Trial_Type 3 14.7978 207 783.96 0.0019979 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#testing a second model for interaction between treatment and trial type affecting movement
msfmodel2 <- glm.nb(Quad_Transitions ~ Treatment * Trial_Type + Crab_ID, data = msf)
msfmodel2
##
## Call: glm.nb(formula = Quad_Transitions ~ Treatment * Trial_Type +
## Crab_ID, data = msf, init.theta = 1.83288897, link = log)
##
## Coefficients:
## (Intercept) TreatmentPredator
## 1.7375363 -0.5465410
## Trial_TypeMF Trial_TypeMM
## -0.1745820 -0.0461322
## Trial_TypeSolo Crab_ID
## 0.0559772 -0.0002222
## TreatmentPredator:Trial_TypeMF TreatmentPredator:Trial_TypeMM
## 0.5010418 0.7157007
## TreatmentPredator:Trial_TypeSolo
## 0.6219998
##
## Degrees of Freedom: 215 Total (i.e. Null); 207 Residual
## Null Deviance: 260.1
## Residual Deviance: 250.5 AIC: 1196
summary(msfmodel2)
##
## Call:
## glm.nb(formula = Quad_Transitions ~ Treatment * Trial_Type +
## Crab_ID, data = msf, init.theta = 1.83288897, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7375363 0.2266158 7.667 1.76e-14 ***
## TreatmentPredator -0.5465410 0.2879988 -1.898 0.0577 .
## Trial_TypeMF -0.1745820 0.2869034 -0.609 0.5429
## Trial_TypeMM -0.0461322 0.2772017 -0.166 0.8678
## Trial_TypeSolo 0.0559772 0.3110602 0.180 0.8572
## Crab_ID -0.0002222 0.0018233 -0.122 0.9030
## TreatmentPredator:Trial_TypeMF 0.5010418 0.4042312 1.239 0.2152
## TreatmentPredator:Trial_TypeMM 0.7157007 0.3933362 1.820 0.0688 .
## TreatmentPredator:Trial_TypeSolo 0.6219998 0.3338137 1.863 0.0624 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8329) family taken to be 1)
##
## Null deviance: 260.10 on 215 degrees of freedom
## Residual deviance: 250.48 on 207 degrees of freedom
## AIC: 1195.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.833
## Std. Err.: 0.257
##
## 2 x log-likelihood: -1175.908
anova(msfmodel2, test = "Chisq")
## Warning in anova.negbin(msfmodel2, test = "Chisq"): tests made without
## re-estimating 'theta'
## Analysis of Deviance Table
##
## Model: Negative Binomial(1.8329), link: log
##
## Response: Quad_Transitions
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 215 260.11
## Treatment 1 0.0205 214 260.08 0.8862
## Trial_Type 3 5.3965 211 254.69 0.1450
## Crab_ID 1 0.0350 210 254.65 0.8516
## Treatment:Trial_Type 3 4.1730 207 250.48 0.2434
#checking which model fits the data the best
AIC(msfmodel1, msfmodel2)
## df AIC
## msfmodel1 9 1460.000
## msfmodel2 10 1195.908
#lowest AIC - best model fit
anova(msfmodel1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Quad_Transitions
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 215 819.52
## Treatment 1 0.0813 214 819.44 0.7755203
## Trial_Type 3 20.6722 211 798.76 0.0001231 ***
## Crab_ID 1 0.0013 210 798.76 0.9708228
## Treatment:Trial_Type 3 14.7978 207 783.96 0.0019979 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(msfmodel1)
##
## Call:
## glm(formula = Quad_Transitions ~ Treatment * Trial_Type + Crab_ID,
## family = poisson, data = msf)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7187332 0.1129024 15.223 < 2e-16 ***
## TreatmentPredator -0.5464848 0.1590525 -3.436 0.000591 ***
## Trial_TypeMF -0.1724724 0.1471927 -1.172 0.241299
## Trial_TypeMM -0.0478567 0.1387447 -0.345 0.730150
## Trial_TypeSolo 0.0179714 0.1551726 0.116 0.907799
## Crab_ID 0.0001037 0.0009192 0.113 0.910146
## TreatmentPredator:Trial_TypeMF 0.4981229 0.2194669 2.270 0.023226 *
## TreatmentPredator:Trial_TypeMM 0.7178019 0.2061919 3.481 0.000499 ***
## TreatmentPredator:Trial_TypeSolo 0.6230780 0.1787850 3.485 0.000492 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 819.52 on 215 degrees of freedom
## Residual deviance: 783.96 on 207 degrees of freedom
## AIC: 1460
##
## Number of Fisher Scoring iterations: 5
#testing a model for no interaction between treatment and trial type when affecting movement
msfmodel_no_interaction <- glm(Quad_Transitions ~ Treatment + Trial_Type, family = poisson, data = msf)
msfmodel_no_interaction
##
## Call: glm(formula = Quad_Transitions ~ Treatment + Trial_Type, family = poisson,
## data = msf)
##
## Coefficients:
## (Intercept) TreatmentPredator Trial_TypeMF Trial_TypeMM
## 1.47972 -0.01012 0.05280 0.29177
## Trial_TypeSolo
## 0.31876
##
## Degrees of Freedom: 215 Total (i.e. Null); 211 Residual
## Null Deviance: 819.5
## Residual Deviance: 798.8 AIC: 1467
anova(msfmodel_no_interaction, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Quad_Transitions
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 215 819.52
## Treatment 1 0.0813 214 819.44 0.7755203
## Trial_Type 3 20.6722 211 798.76 0.0001231 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(msfmodel1,msfmodel_no_interaction) #testing AIC for new model vs previous best model
## df AIC
## msfmodel1 9 1460.000
## msfmodel_no_interaction 5 1466.799
#testing which pairing types have a significant effect on movement rate
msfmodeltest <- emmeans(msfmodel1, pairwise ~ Treatment|Trial_Type, data = msf, adjust = "tukey")
summary(msfmodeltest)
## $emmeans
## Trial_Type = FF:
## Treatment emmean SE df asymp.LCL asymp.UCL
## Control 1.73 0.1110 Inf 1.513 1.95
## Predator 1.18 0.1320 Inf 0.925 1.44
##
## Trial_Type = MF:
## Treatment emmean SE df asymp.LCL asymp.UCL
## Control 1.56 0.1190 Inf 1.325 1.79
## Predator 1.51 0.1150 Inf 1.283 1.74
##
## Trial_Type = MM:
## Treatment emmean SE df asymp.LCL asymp.UCL
## Control 1.68 0.1070 Inf 1.472 1.89
## Predator 1.85 0.1000 Inf 1.657 2.05
##
## Trial_Type = Solo:
## Treatment emmean SE df asymp.LCL asymp.UCL
## Control 1.75 0.0805 Inf 1.591 1.91
## Predator 1.82 0.0788 Inf 1.670 1.98
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## Trial_Type = FF:
## contrast estimate SE df z.ratio p.value
## Control - Predator 0.5465 0.1590 Inf 3.436 0.0006
##
## Trial_Type = MF:
## contrast estimate SE df z.ratio p.value
## Control - Predator 0.0484 0.1510 Inf 0.320 0.7492
##
## Trial_Type = MM:
## contrast estimate SE df z.ratio p.value
## Control - Predator -0.1713 0.1310 Inf -1.306 0.1917
##
## Trial_Type = Solo:
## contrast estimate SE df z.ratio p.value
## Control - Predator -0.0766 0.0816 Inf -0.938 0.3482
##
## Results are given on the log (not the response) scale.
#creating graph to visualize best-fit model for how movement changes with treatment and trial type
movement <- ggplot(msf, aes(x = Treatment, y = Quad_Transitions, color = Trial_Type)) +
stat_summary(fun = mean, geom = "point", position = position_dodge(0.2), size = 2.5) +
stat_summary(fun = mean, geom = "line", aes(group = Trial_Type), size = 2) +
labs(title = "Interaction Between Treatment and Trial Type",
y = "Mean Quadrant Transitions",
x = "Treatment", color = "Trial Type") +
theme_classic2() +
theme(axis.title.x = element_text(size=14), axis.title.y = element_text(size=14),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black")) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_color_manual(values=c("#6C568CFF","#BFCDD9FF", "#5D761EFF", "#E394BBFF")) + #these colors are for FF, MF, MM, and solo crab trials respectively
theme(legend.position = "right")
movement
ggsave("interaction graph.tiff")
## Saving 7 x 5 in image