15_field_ANOVA-copepods
2023-08-24
This code runs a 2-way ANOVA test for Ditrichocorycaeus anglicus and medium calanoid copepods, comparing mean densities among months and location in relation to jellyfish aggregations. The data are square-root transformed and plotted. Only data from Sinclair Inlet and Quartermaster Harbor are included here.
- The RMarkdown file is here: https://rpubs.com/HailaSchultz/field_ANOVA-copepods
load packages
library(forcats)
library(dplyr)
library(ggplot2)
library(stringr)
library(scales)
library(ggpubr)
library(emmeans)Set up data
Read database into R. If current database changes, just change the path
Database <- read.csv("/Users/hailaschultz/Dropbox/Other studies/Aurelia project/Data Analysis/data/current_data/Final_Aurelia_Database_Jan11_2023.csv")subset data
#select only field data
Field_data<-subset(Database,Trial.Type=="Field")
#subset to QM and SC
Field_SCQM<-subset(Field_data, Site %in% c("QM","SC"))recode copepod names
Field_SCQM$Species_lifestage <- paste(Field_SCQM$Genus.species, Field_SCQM$Life.History.Stage, sep="_")
Field_SCQM <- Field_SCQM %>%
mutate(Species_lifestage_combined = fct_recode(Species_lifestage,
"CALANOIDA_Medium" = "ACARTIA_Copepodite",
"CALANOIDA_Medium" = "ACARTIA_Female, Adult",
"CALANUS PACIFICUS" = "CALANUS PACIFICUS_C5-adult",
"CALANUS PACIFICUS" = "CALANUS PACIFICUS_C5-Adult",
"CALANUS PACIFICUS" = "CALANUS PACIFICUS_Copepodite",
"CALANUS PACIFICUS" = "CALANUS PACIFICUS_Female, Adult",
"CALANUS PACIFICUS" = "CALANUS PACIFICUS_Male, Adult",
"DITRICHOCORYCAEUS ANGLICUS" = "DITRICHOCORYCAEUS ANGLICUS_Large",
"DITRICHOCORYCAEUS ANGLICUS" = "DITRICHOCORYCAEUS ANGLICUS_Small"))
unique(Field_SCQM$Species_lifestage_combined)## [1] CALANOIDA_Medium AETIDEUS_Copepodite
## [3] AETIDEUS_Female, Adult AURELIA LABIATA_Ephyra
## [5] BARNACLES_Cyprid larva BARNACLES_Nauplius
## [7] BIVALVIA_Veliger BRACHYURA_Unknown
## [9] BRYOZOA_Cyphonaut CALANOIDA_Copepodite
## [11] CALANOIDA_Large CALANOIDA_Small
## [13] CALANUS PACIFICUS CENTROPAGES_Female, Adult
## [15] Chaet/Euphaus Egg_Egg CHAETOGNATHA_Unknown
## [17] CLADOCERA_Unknown CLYTIA GREGARIA_Medusa
## [19] COPEPODA_Egg COPEPODA_Nauplius
## [21] CRABS_Megalopa CRABS_Zoea
## [23] Diatom-Centric_Unknown DITRICHOCORYCAEUS ANGLICUS
## [25] ECHINODERMATA_Pluteus larva EUPHAUSIIDAE_Calyptopis
## [27] EUPHAUSIIDAE_Furcilia EUPHAUSIIDAE_Nauplius
## [29] EUTONINA INDICANS_Medusa FISH_Egg
## [31] FISH_Larva GAETANUS_Adult
## [33] GAMMARIDEA_Unknown GASTROPODA_Veliger
## [35] HARPACTICOIDA_Unknown HYPERIIDEA_Unknown
## [37] Insecta_Unknown ISOPODA_Unknown
## [39] JELLYFISHES_Ephyra JELLYFISHES_Medusa
## [41] JELLYFISHES_Planula larvae LARVACEA_Unknown
## [43] LITTORINA_Egg METRIDIA_Female, Adult
## [45] MOLLUSCA_Trochophore larva Nemertea_Pilidium
## [47] NEOTRYPAEA CALIFORNIENSIS_Unknown NOCTILUCA_Unknown
## [49] OITHONA_Unknown POLYCHAETA_Unknown
## [51] SHRIMP_Unknown SIPHONOPHORA_Bract
## [53] SIPHONOPHORA_Gonophore SIPHONOPHORA_Unknown
## [55] UNKNOWN_Larva
## 55 Levels: CALANOIDA_Medium AETIDEUS_Copepodite ... UNKNOWN_LarvaDeal with duplicate station names
# combine station and date for unique stations
Field_SCQM$Station_unique <- paste(Field_SCQM$Station, Field_SCQM$Sample.Date)
unique(Field_SCQM$Station_unique)## [1] "QM1 08/28/2019" "SC4 09/25/2019" "SC2 09/25/2019" "SC3 09/25/2019"
## [5] "SC1 09/25/2019" "QM10 08/24/2021" "QM5b 07/29/2021" "QM6b 07/29/2021"
## [9] "QM2c 08/22/2021" "SC11 08/26/2021" "SC5c 08/25/2021" "SC6c 08/25/2021"
## [13] "SC1a 07/27/2021" "SC4a 07/27/2021" "SC16 08/26/2021" "SC10 08/26/2021"
## [17] "SC7 08/25/2021" "QM1b 07/29/2021" "SC2c 08/24/2021" "QM1c 08/22/2021"
## [21] "SC14 08/26/2021" "SC3a 07/27/2021" "QM9 08/23/2021" "SC5d 08/25/2021"
## [25] "QM2b 07/28/2021" "SC3c 08/24/2021" "QM6c 08/22/2021" "QM5c 08/22/2021"
## [29] "SC17 08/27/2021" "QM3c 08/22/2021" "SC5 09/25/2019" "QM6 08/28/2019"
## [33] "QM5 08/28/2019" "QM2 08/28/2019" "SC6 09/25/2019" "QM3 08/28/2019"
## [37] "QM4 08/28/2019" "QM3a 08/29/2020" "QM1a 08/28/2020" "SC2a 07/27/2021"
## [41] "QM5a 08/29/2020" "QM4a 08/29/2020" "SC8 08/25/2021" "SC12 08/26/2021"
## [45] "SC4c 08/25/2021" "SC6a 07/27/2021" "QM11 08/24/2021" "QM4b 07/29/2021"
## [49] "QM12 08/24/2021" "SC5a 07/27/2021" "QM3b 07/29/2021" "QM2 07/24/2019"
## [53] "QM6 07/24/2019" "QM1 07/24/2019" "QM2a 08/28/2020" "QM8a 08/29/2020"
## [57] "QM3 07/24/2019" "QM5 07/24/2019" "QM4 07/24/2019" "QM7a 08/29/2020"
## [61] "QM4c 08/22/2021"Field_SCQM$Station_unique <-recode(Field_SCQM$Station_unique,
'QM1 07/24/2019'='QM1s 07/24/2019',
'QM2 07/24/2019'='QM2s 07/24/2019',
'QM3 07/24/2019'='QM3s 07/24/2019',
'QM4 07/24/2019'='QM4s 07/24/2019',
'QM5 07/24/2019'='QM5s 07/24/2019',
'QM6 07/24/2019'='QM6s 07/24/2019',
)
#remove date
Field_SCQM$Station_unique<-str_sub(Field_SCQM$Station_unique, end=-12)
unique(Field_SCQM$Station_unique)## [1] "QM1" "SC4" "SC2" "SC3" "SC1" "QM10" "QM5b" "QM6b" "QM2c" "SC11"
## [11] "SC5c" "SC6c" "SC1a" "SC4a" "SC16" "SC10" "SC7" "QM1b" "SC2c" "QM1c"
## [21] "SC14" "SC3a" "QM9" "SC5d" "QM2b" "SC3c" "QM6c" "QM5c" "SC17" "QM3c"
## [31] "SC5" "QM6" "QM5" "QM2" "SC6" "QM3" "QM4" "QM3a" "QM1a" "SC2a"
## [41] "QM5a" "QM4a" "SC8" "SC12" "SC4c" "SC6a" "QM11" "QM4b" "QM12" "SC5a"
## [51] "QM3b" "QM2s" "QM6s" "QM1s" "QM2a" "QM8a" "QM3s" "QM5s" "QM4s" "QM7a"
## [61] "QM4c"Ditrichocorycaeus anglicus
subset
Field_DA<-subset(Field_SCQM, Species_lifestage_combined == "DITRICHOCORYCAEUS ANGLICUS")add up multiple lines per station
colnames(Field_DA)## [1] "BugSampleID" "Project"
## [3] "Sample.Code" "Sampling.Group"
## [5] "Station" "Site"
## [7] "Site.Name" "Basin"
## [9] "Sub.Basin" "Latitude"
## [11] "Longitude" "Sample.Date"
## [13] "Sample.Year" "Sample.Month"
## [15] "Sample.Time" "Tow.Type"
## [17] "Mesh.Size" "Station.Depth..m."
## [19] "Flow.meter..revs." "Broad.Group"
## [21] "Mid.Level.Group" "X1st.Word.Taxa"
## [23] "Genus.species" "Life.History.Stage"
## [25] "Total.Ct" "Density....m3."
## [27] "Vol.Filtered..m3." "Jelly.Mass..g."
## [29] "Number.of.Jellies" "Trial.Time"
## [31] "Trial.Type" "Jelly.Size"
## [33] "Jelly.Density....m3." "Jelly.Density..g.m3."
## [35] "Location" "Species_lifestage"
## [37] "Species_lifestage_combined" "Station_unique"Field_DA <- Field_DA %>%
group_by(Sample.Code,Station_unique,Sample.Year,Sample.Month,Site,Location) %>%
summarise(
copepod_density = sum(Density....m3.))add taxa column
Field_DA$taxa<-"Ditrichocorycaeus"Medium Calanoids
subset
Field_Cal<-subset(Field_SCQM, Species_lifestage_combined == "CALANOIDA_Medium")add up multiple lines per station
colnames(Field_Cal)## [1] "BugSampleID" "Project"
## [3] "Sample.Code" "Sampling.Group"
## [5] "Station" "Site"
## [7] "Site.Name" "Basin"
## [9] "Sub.Basin" "Latitude"
## [11] "Longitude" "Sample.Date"
## [13] "Sample.Year" "Sample.Month"
## [15] "Sample.Time" "Tow.Type"
## [17] "Mesh.Size" "Station.Depth..m."
## [19] "Flow.meter..revs." "Broad.Group"
## [21] "Mid.Level.Group" "X1st.Word.Taxa"
## [23] "Genus.species" "Life.History.Stage"
## [25] "Total.Ct" "Density....m3."
## [27] "Vol.Filtered..m3." "Jelly.Mass..g."
## [29] "Number.of.Jellies" "Trial.Time"
## [31] "Trial.Type" "Jelly.Size"
## [33] "Jelly.Density....m3." "Jelly.Density..g.m3."
## [35] "Location" "Species_lifestage"
## [37] "Species_lifestage_combined" "Station_unique"Field_Cal <- Field_Cal %>%
group_by(Sample.Code,Station_unique,Sample.Year,Sample.Month,Site,Location) %>%
summarise(
copepod_density = sum(Density....m3.))add taxa column
Field_Cal$taxa<-"Med_Calanoid"combine species
Copepod_QMSC<-rbind(Field_DA,Field_Cal)
Copepod_QMSC$Sample.Year<-as.character(Copepod_QMSC$Sample.Year)Location and Month
ANOVA
Cop_mod1 <- aov(copepod_density ~ taxa+Location+Sample.Month+taxa:Sample.Month+taxa:Location+Location:Sample.Month, data = Copepod_QMSC)
summary(Cop_mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## taxa 1 28701650 28701650 74.434 6.01e-14 ***
## Location 2 7707384 3853692 9.994 0.000104 ***
## Sample.Month 2 5665812 2832906 7.347 0.001020 **
## taxa:Sample.Month 2 2569623 1284812 3.332 0.039430 *
## taxa:Location 2 1161109 580554 1.506 0.226507
## Location:Sample.Month 4 762887 190722 0.495 0.739699
## Residuals 108 41644833 385600
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(Cop_mod1, which = "Location")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = copepod_density ~ taxa + Location + Sample.Month + taxa:Sample.Month + taxa:Location + Location:Sample.Month, data = Copepod_QMSC)
##
## $Location
## diff lwr upr p adj
## OB-IN 559.0947 210.8073 907.3822 0.0006612
## OS-IN 456.0708 148.0864 764.0551 0.0018244
## OS-OB -103.0240 -478.6108 272.5629 0.7917120check normality assumptions
hist(Field_DA$copepod_density)
qqnorm(Cop_mod1$residuals)
qqline(Cop_mod1$residuals)
shapiro.test(Cop_mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: Cop_mod1$residuals
## W = 0.95869, p-value = 0.0008857looks like copepods has a log distribution, the qqplot has a slight s-shape, and the shapiro-wilk test meaning that the data significantly deviate from normal distribution. A transformation is necessary.
##log tranformation
Cop_mod2 <- aov(log(copepod_density) ~ taxa+Location+Sample.Month+taxa:Sample.Month+taxa:Location+Location:Sample.Month, data = Copepod_QMSC)
summary(Cop_mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## taxa 1 95.58 95.58 87.273 1.46e-15 ***
## Location 2 39.62 19.81 18.088 1.68e-07 ***
## Sample.Month 2 25.63 12.82 11.703 2.51e-05 ***
## taxa:Sample.Month 2 1.50 0.75 0.683 0.507030
## taxa:Location 2 18.80 9.40 8.582 0.000347 ***
## Location:Sample.Month 4 1.96 0.49 0.448 0.773303
## Residuals 108 118.28 1.10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(Cop_mod2, which = "Location")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(copepod_density) ~ taxa + Location + Sample.Month + taxa:Sample.Month + taxa:Location + Location:Sample.Month, data = Copepod_QMSC)
##
## $Location
## diff lwr upr p adj
## OB-IN 1.3257368 0.7387679 1.9127056 0.0000014
## OS-IN 0.9656004 0.4465544 1.4846465 0.0000694
## OS-OB -0.3601363 -0.9931129 0.2728402 0.3697589check normality assumptions
qqnorm(Cop_mod2$residuals)
qqline(Cop_mod2$residuals)
shapiro.test(Cop_mod2$residuals)##
## Shapiro-Wilk normality test
##
## data: Cop_mod2$residuals
## W = 0.95266, p-value = 0.0003002The distribution looks slightly better, but it is right skewed. The qq plot looks weird to the left side. Shapiro-wilk is still significant, although the p-value has improved. Try square root transformation.
Square root transformation
ANOVA
Cop_mod3 <- aov(sqrt(copepod_density) ~ taxa+Location+Sample.Month+taxa:Sample.Month+taxa:Location+Location:Sample.Month+taxa:Sample.Month:Location, data = Copepod_QMSC)check normality assumptions
qqnorm(Cop_mod3$residuals)
qqline(Cop_mod3$residuals)
shapiro.test(Cop_mod3$residuals)##
## Shapiro-Wilk normality test
##
## data: Cop_mod3$residuals
## W = 0.99091, p-value = 0.6061These results look the best, and shapiro-wilk is not significant! Proceed with sqrt transformation. location is significant for bartlett, but since the sample sizes are more even for location, I think this is okay
summary(Cop_mod3)## Df Sum Sq Mean Sq F value Pr(>F)
## taxa 1 9582 9582 100.151 < 2e-16 ***
## Location 2 3279 1639 17.134 3.70e-07 ***
## Sample.Month 2 2133 1067 11.148 4.11e-05 ***
## taxa:Sample.Month 2 609 304 3.180 0.0457 *
## taxa:Location 2 732 366 3.826 0.0249 *
## Location:Sample.Month 4 129 32 0.336 0.8529
## taxa:Location:Sample.Month 4 489 122 1.279 0.2832
## Residuals 104 9950 96
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is an interaction effect between taxa and sample month and taxa and location. However, the three-way interaction is insignificant. I should then look into the two-way interactions and look at simple effects between those two pairs. Therefore, I will analyze the two taxa differently.
Ditrichocorycaeus
sqrt transformation
DA_mod1 <- aov(sqrt(copepod_density) ~ Location+Sample.Month+Sample.Month:Location, data = Field_DA)check normality assumptions
qqnorm(DA_mod1$residuals)
qqline(DA_mod1$residuals)
shapiro.test(DA_mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: DA_mod1$residuals
## W = 0.95721, p-value = 0.03216shapiro test isn’t quite above 0.05, but it is the best of the options
summary(DA_mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## Location 2 3315 1657.4 30.832 1.48e-09 ***
## Sample.Month 2 312 156.2 2.905 0.0637 .
## Location:Sample.Month 4 271 67.7 1.260 0.2977
## Residuals 52 2795 53.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(DA_mod1, which = "Location")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sqrt(copepod_density) ~ Location + Sample.Month + Sample.Month:Location, data = Field_DA)
##
## $Location
## diff lwr upr p adj
## OB-IN 18.38294 12.478877 24.287005 0.0000000
## OS-IN 10.32847 5.107613 15.549330 0.0000446
## OS-OB -8.05447 -14.421305 -1.687634 0.0098411get sample sizes
#count sample sizes per category
Field_DA %>%
group_by(Sample.Month,Location) %>%
summarize(count = n()) ## # A tibble: 9 × 3
## # Groups: Sample.Month [3]
## Sample.Month Location count
## <chr> <chr> <int>
## 1 Aug IN 21
## 2 Aug OB 5
## 3 Aug OS 11
## 4 Jul IN 6
## 5 Jul OB 6
## 6 Jul OS 6
## 7 Sep IN 2
## 8 Sep OB 2
## 9 Sep OS 2make plot
#reorder factors
Field_DA$Sample.Month <- factor(Field_DA$Sample.Month,
c("Jul", "Aug", "Sep"))
#recode stations
Field_DA <- Field_DA %>%
mutate(Location = fct_recode(Location,
"Inside Aggregation" = "IN",
"Outside Aggregation" = "OS",
"Outside Bay"="OB"))
Field_DA$Location <- factor(Field_DA$Location,
c("Inside Aggregation", "Outside Aggregation", "Outside Bay"))
DA_plot<- ggplot(Field_DA, aes(x=Sample.Month, y=sqrt(copepod_density), fill=Location)) +
geom_boxplot() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+theme_classic()+
xlab(bquote('Month'))+
ylab(expression(plain('sqrt (')~ italic('D. anglicus') ~ plain('Density )') ~ plain('( # /' ~ m^3 ~ ')')))+
scale_fill_manual(values = c("Inside Aggregation" = "#619CFF", "Outside Aggregation" = "#F8766D", "Outside Bay" = "#00BA38"))+
theme(axis.text.y=element_text(size=15,colour="black"),
legend.position=c(0.78, 0.91),
legend.text=element_text(size=15,colour="black"),
legend.title=element_text(size=15,colour="black"),
axis.title.y=element_text(size=20,colour="black"),
axis.line = element_line(colour = "black"),
axis.title.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x=element_blank())+
annotate("text",x=0.75, y=20, label= "6",size=6)+
annotate("text",x=1, y=29, label= "6",size=6)+
annotate("text",x=1.25, y=34, label= "6",size=6)+
annotate("text",x=1.75, y=20, label= "21",size=6)+
annotate("text",x=2, y=35, label= "11",size=6)+
annotate("text",x=2.25, y=38, label= "5",size=6)+
annotate("text",x=2.75, y=16, label= "2",size=6)+
annotate("text",x=3, y=10, label= "2",size=6)+
annotate("text",x=3.25, y=33, label= "2",size=6)+
annotate("text",x=0.6, y=49, label= "a",size=10)+
scale_y_continuous(breaks=seq(0,50,10),expand = c(0, 0), limits = c(0, 50))+
theme(plot.margin=margin(0.5, 0.5, 0.5, 0.5, unit = "cm"))
DA_plot
Medium calanoids
sqrt transformation
Cal_mod1 <- aov(sqrt(copepod_density) ~ Location+Sample.Month+Sample.Month:Location, data = Field_Cal)check normality assumptions
qqnorm(Cal_mod1$residuals)
qqline(Cal_mod1$residuals)
shapiro.test(Cal_mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: Cal_mod1$residuals
## W = 0.97514, p-value = 0.2489looks good
summary(Cal_mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## Location 2 764 381.9 2.775 0.071575 .
## Sample.Month 2 2361 1180.7 8.581 0.000602 ***
## Location:Sample.Month 4 347 86.8 0.631 0.642738
## Residuals 52 7155 137.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(Cal_mod1, which = "Sample.Month")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sqrt(copepod_density) ~ Location + Sample.Month + Sample.Month:Location, data = Field_Cal)
##
## $Sample.Month
## diff lwr upr p adj
## Jul-Aug 2.539252 -5.593352 10.671855 0.7330054
## Sep-Aug -19.666952 -32.121934 -7.211970 0.0010607
## Sep-Jul -22.206204 -35.546919 -8.865488 0.0005527get sample sizes
#count sample sizes per category
Field_Cal %>%
group_by(Sample.Month,Location) %>%
summarize(count = n()) ## # A tibble: 9 × 3
## # Groups: Sample.Month [3]
## Sample.Month Location count
## <chr> <chr> <int>
## 1 Aug IN 21
## 2 Aug OB 5
## 3 Aug OS 11
## 4 Jul IN 6
## 5 Jul OB 6
## 6 Jul OS 6
## 7 Sep IN 2
## 8 Sep OB 2
## 9 Sep OS 2make plot
#reorder factors
Field_Cal$Sample.Month <- factor(Field_Cal$Sample.Month,
c("Jul", "Aug", "Sep"))
Field_Cal$Location <- factor(Field_Cal$Location,
c("IN", "OS", "OB"))
Cal_plot<- ggplot(Field_Cal, aes(x=Sample.Month, y=sqrt(copepod_density), fill=Location)) +
geom_boxplot() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+theme_classic()+
scale_fill_manual(values = c("IN" = "#619CFF", "OS" = "#F8766D", "OB" = "#00BA38"))+
theme(axis.text=element_text(size=15,colour="black"),
legend.position = "none",
axis.title=element_text(size=20,colour="black"),
axis.line = element_line(colour = "black"))+
xlab(bquote('Month'))+
ylab(bquote('sqrt ( Med. Calanoid Density ) ( # /'~m^3~')'))+
annotate("text",x=0.75, y=54, label= "6",size=6)+
annotate("text",x=1, y=54, label= "6",size=6)+
annotate("text",x=1.25, y=58, label= "6",size=6)+
annotate("text",x=1.75, y=56, label= "21",size=6)+
annotate("text",x=2, y=65, label= "11",size=6)+
annotate("text",x=2.25, y=53, label= "5",size=6)+
annotate("text",x=2.75, y=14, label= "2",size=6)+
annotate("text",x=3, y=32, label= "2",size=6)+
annotate("text",x=3.25, y=32, label= "2",size=6)+
annotate("text",x=0.6, y=68, label= "b",size=10)+
scale_y_continuous(breaks=seq(0,70,10),expand = c(0, 0), limits = c(0, 70))+
theme(plot.margin=margin(0.5, 0.5, 0.5, 0.5, unit = "cm"))
Cal_plot
Final Plot
ANOVA_plot<-ggarrange(DA_plot,Cal_plot,
ncol = 1, nrow = 2)
ANOVA_plot
save plot
setwd("/Users/hailaschultz/Dropbox/Other studies/Aurelia project/Data Analysis/output")
ggsave(filename = "Field_Copepod_ANOVAs.png", plot = ANOVA_plot, width = 164, height = 328, units="mm", device='png', dpi=600)
ggsave(filename = "Field_Copepod_ANOVAs.tif", plot = ANOVA_plot, width = 164, height = 328, units="mm", device='tiff', dpi=600)