Unit 3: final
Katerina Georgiopoulos
Unit 3: Final, for this exam you will be using the file “fouling.csv” found on canvas:
This data come from work I did with a few undergraduate students in 2022. Fouling communities are marine communities that foul or cover hard surfaces (like dock pilings, boats, rocks, etc.). They are comprised of sessile (non-moving) invertebrates. They also happen to contain a lot of invasive species. This project collected data on fouling communities that adhered themselves to plastic plates hanging in 2 locations in MA - Beverly Harbor and Gloucester Harbor. Plates were left in the water for 1 month, 2 months, 6 months, or 12 months. Plates were assessed for species present at each location and time point (there are multiple plates for each location and time point). Data was organized for us to make statistical comparisons across time and space.
The data is set up as follows: “UniqueID” - this was an ID created to represent individual plates, they are given a “G” or a “B” indicating Gloucester or Beverly, a number for months in water (1,2,6, or 12), and a final letter (A,B,C,D,E, or F) indicating a single plate.
“Location” - another indicator of where the plate was deployed (Beverly or Gloucester)
All of the other columns are species observed on the plates and the values below them indicate abundance/count
PART 1: Visualizing community similarities and differences
library(factoextra)## Warning: package 'factoextra' was built under R version 4.4.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.4.3## Welcome to factoextra!## Want to learn more? See two factoextra-related books at https://www.datanovia.com/en/product/practical-guide-to-principal-component-methods-in-r/library(FactoMineR)## Warning: package 'FactoMineR' was built under R version 4.4.3library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.4.3library(RColorBrewer)
fouling <- read.csv("~/Downloads/fouling.csv")
pca.data <- PCA(fouling[, 3:20],
scale.unit = TRUE,
graph = FALSE)
fviz_eig(pca.data,
addlabels = TRUE,
ylim = c(0, 70))
fviz_pca_var(pca.data,
col.var = "cos2",
gradient.cols = c("gold", "orange", "red", "darkred"),
repel = TRUE)
pca.data <- PCA(fouling[, 3:20],
scale.unit = TRUE,
ncp = 2,
graph = FALSE)
fouling$Location <- as.factor(fouling$Location)
nb.cols <- 2
mycolors <- colorRampPalette(brewer.pal(3, "Set1"))(nb.cols)
a <- fviz_pca_ind(pca.data,
col.ind = fouling$Location,
palette = mycolors,
addEllipses = TRUE)
ggpar(a,
title = "Principal Component Analysis",
xlab = "PC1",
ylab = "PC2",
legend.title = "Location",
legend.position = "top",
ggtheme = theme_minimal())## Ignoring unknown labels:
## • linetype : "Location"
##Question: how would you describe the scree plot? is most of the
influence from a handful of species or is it more spreadout? The scree
plot shows that the variation is fairly spread out across the different
points. There is not a huge drop after the first few components. This
can mean that the patterns are influenced by multiple species instead of
just a few.
- Question: With the final PCA plot, tell me about how similar/different the two locations are, do you spot any outliers? The two locations overlap a bit, this could mean that the communities are relatively similar. There are also a few differenes between the locaations. Gloucester has some samples that are far from the main group (33,35,and 36), as well as Beverly (29,31,and 32).
PART 2: create a ggplot bar plot for the data
- Code: using ggplot, create a barplot representing the mean Ciona value (this is a dominant organism that has shown an ability to supress diversity) and add error bars representing standard deviation. NOTE: you will need to follow similar steps to R Lab 6 to calculate, store, and merge values…
library(ggplot2)
Ciona.mean <- aggregate(Ciona ~ Location, data = fouling, mean)
Ciona.sd <- aggregate(Ciona ~ Location, data = fouling, sd)
Ciona.summary <- merge(Ciona.mean, Ciona.sd, by = "Location")
names(Ciona.summary) <- c("Location", "Mean", "SD")
Ciona.summary## Location Mean SD
## 1 Beverly 32.25 47.38518
## 2 Gloucester 37.95 34.27823- Question: How does average Ciona abundance differ between Beverly and Gloucester? The average Ciona abundance was higher in Gloucester than in Beverly. Glouster had an average abdundance of 37.95, while Beverly had an average fo 32.25.
PART 3: creating a different ggplot
- Code: using ggplot, create a scatter plot with a linear line showing Ciona counts vs. Ascidiella counts across all sites (Ascidiella is a closely related species to Ciona and the two likely compete for space) and add the CI visual to the line.
ggplot(fouling, aes(x = Ciona, y = Ascidiella)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
xlab("Ciona abundance") +
ylab("Ascidiella abundance") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
- Question: how would you describe the relationship between Ciona and Ascidiella abundance? There is a weak negative relationship between Ciona and Ascidiella abundance. As the Ciona abundance increase the Ascidiella starts decreasing slightly. The points are also fairly spread out around the line suggesting that the realtionship might not be very strong.
BONUS: Ciona and Ascidiella are both scientific names and therefore should be written in italics! You will need to look up how to do this, but if you can create the plot above and have the axes labels show the italicized names you can earn 3 bonus points!
ggplot(fouling, aes(x = Ciona, y = Ascidiella)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
labs(x = expression(italic("Ciona") ~ "abundance"),
y = expression(italic("Ascidiella") ~ "abundance")) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
PART 4: altering the data and running a statistical test
- Add a column to the fouling.csv file, this column should represent the number of months a plate has been in the water. Remember, the UniqueID column contains information about number of months - use this to determine the value needed for the new column (should be a 1, 2, 6, or 12). You can add this column in R or using sheets/excel (whatever is easiest for you)
- With the column added, make sure the updated dataset is back in R.
- Code: convert the new column to a factor
- Code: run a two-way ANOVA that looks at how the interaction of the new column AND Location influence the abundance of Ciona.
fouling$Unique.ID## [1] "G1A" "G1B" "G1C" "G1D" "B1A" "B1B" "B1C" "B1D" "G2A" "G2B"
## [11] "G2C" "G2D" "G2E" "B2A" "B2B" "B2C" "B2D" "B2E" "G6A" "G6B"
## [21] "G6C" "G6D" "G6E" "B6A" "B6B" "B6C" "B6D" "B6E" "G12A" "G12B"
## [31] "G12C" "G12D" "G12E" "G12F" "B12A" "B12B" "B12C" "B12D" "B12E" "B12F"fouling$Months <- ifelse(grepl("1", fouling$Unique.ID), 1,
ifelse(grepl("2", fouling$Unique.ID), 2,
ifelse(grepl("6", fouling$Unique.ID), 6, 12)))
fouling[, c("Unique.ID", "Months")]## Unique.ID Months
## 1 G1A 1
## 2 G1B 1
## 3 G1C 1
## 4 G1D 1
## 5 B1A 1
## 6 B1B 1
## 7 B1C 1
## 8 B1D 1
## 9 G2A 2
## 10 G2B 2
## 11 G2C 2
## 12 G2D 2
## 13 G2E 2
## 14 B2A 2
## 15 B2B 2
## 16 B2C 2
## 17 B2D 2
## 18 B2E 2
## 19 G6A 6
## 20 G6B 6
## 21 G6C 6
## 22 G6D 6
## 23 G6E 6
## 24 B6A 6
## 25 B6B 6
## 26 B6C 6
## 27 B6D 6
## 28 B6E 6
## 29 G12A 1
## 30 G12B 1
## 31 G12C 1
## 32 G12D 1
## 33 G12E 1
## 34 G12F 1
## 35 B12A 1
## 36 B12B 1
## 37 B12C 1
## 38 B12D 1
## 39 B12E 1
## 40 B12F 1fouling$Months <- as.factor(fouling$Months)
str(fouling)## 'data.frame': 40 obs. of 21 variables:
## $ Unique.ID : chr "G1A" "G1B" "G1C" "G1D" ...
## $ Location : Factor w/ 2 levels "Beverly","Gloucester": 2 2 2 2 1 1 1 1 2 2 ...
## $ D..vexillum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Anemone : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Semibalanus : int 0 0 1 0 0 0 0 0 0 0 ...
## $ Mytilus : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Styela : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Juv..Mussel : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Diplosoma : int 2 10 18 7 0 2 0 0 9 16 ...
## $ Encrusting.bryozoan: int 0 0 0 0 1 0 0 0 0 0 ...
## $ Oyster : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ascidiella : int 19 25 64 13 25 32 11 20 4 13 ...
## $ Botryllus : int 0 0 2 0 0 0 0 0 0 0 ...
## $ Eudendrium : int 29 58 0 58 145 145 319 261 58 29 ...
## $ Limpet : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Corella : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Bugula : int 20 20 7 4 60 40 0 0 20 20 ...
## $ Ribbed.Mussel : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Ciona : int 40 53 53 38 10 12 0 0 69 70 ...
## $ Botrylloides : int 0 6 0 1 2 4 0 1 2 3 ...
## $ Months : Factor w/ 3 levels "1","2","6": 1 1 1 1 1 1 1 1 2 2 ...anova.model <- aov(Ciona ~ Months * Location, data = fouling)
summary(anova.model)## Df Sum Sq Mean Sq F value Pr(>F)
## Months 2 10459 5230 4.283 0.02194 *
## Location 1 3 3 0.002 0.96344
## Months:Location 2 13332 6666 5.459 0.00879 **
## Residuals 34 41518 1221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Question: what does your output tell you about Ciona abundance? Ciona abundance was not significantly affected by number of months in the water (p = 0.02194). Location did not have a significant effect on Ciona abundance (p = 0.96344). There was however a significant relationship between months and location (p = 0.00879), this meant that time had a different effect between Beverly and Gloucester.
- Code: create a visual that helps represent how the new column and Location relate to abundance.
library(ggplot2)
ggplot(fouling, aes(x = Months, y = Ciona, fill = Location)) +
geom_boxplot() +
scale_fill_manual(values = c("Beverly" = "purple", "Gloucester" = "blue")) +
xlab("Months in water") +
ylab("Ciona abundance") +
theme_minimal()