Unit 3: final

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

NOTE: for each prompt, I need to see code in order to give you credit! NOTE: you can perform all of these tasks in one code chunk or many code chunks, this is up to you

Code: bring in the file “fouling.csv”, call it whatever you want

fouling <- read.csv("fouling.csv", header = TRUE)
str(fouling)

## 'data.frame':    40 obs. of  20 variables:
##  $ Unique.ID          : chr  "G1A" "G1B" "G1C" "G1D" ...
##  $ Location           : chr  "Gloucester" "Gloucester" "Gloucester" "Gloucester" ...
##  $ 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 ...

Code: Create a series of PCA plots that help visualize the patterns by location (IMPORTANT - you should be able to utilize the exact same code from the PCA tutorial - just be aware that your dataset will likely have a different name…)

library(FactoMineR)
library(factoextra)

## Loading required package: ggplot2

## Welcome to factoextra!

## Want to learn more? See two factoextra-related books at https://www.datanovia.com/library/principal-component-methods

library(ggplot2)

species <- fouling[,-c(1,2)]

pca.data <- PCA(
  species,
  scale.unit = TRUE,
  ncp = 18,
  graph = TRUE
)

fviz_eig(
  pca.data,
  addlabels = TRUE,
  ylim = c(0, 70)
)

fviz_pca_var(
  pca.data,
  col.var = "cos2",
  gradient.cols = c(
    "pink",
    "hotpink",
    "orchid",
    "darkorchid"
  ),
  repel = TRUE
)

fviz_pca_ind(
  pca.data,
  habillage = fouling$Location,
  addEllipses = TRUE,
  ellipse.type = "confidence",
  repel = TRUE
)

Question: how would you describe the scree plot? is most of the influence from a handful of species or is it more spreadout?

Based on the scree plot, the variation of this data is fairly spread out. There is not one of two variables that completely dominate the variation.

Question: With the final PCA plot, tell me about how similar/different the two locations are, do you spot any outliers?

The two locations seem to be pretty similar. There is a lot of overlap and a few outliers, the most outlying from Gloucester but some from Beverly as well.

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…

ciona.mean <- aggregate(
  Ciona ~ Location,
  data = fouling,
  FUN = mean
)
ciona.mean

##     Location Ciona
## 1    Beverly 32.25
## 2 Gloucester 37.95

ciona.sd <- aggregate(
  Ciona ~ Location,
  data = fouling,
  FUN = sd
)
ciona.sd

##     Location    Ciona
## 1    Beverly 47.38518
## 2 Gloucester 34.27823

names(ciona.mean)[2] <- "Mean"
names(ciona.sd)[2] <- "SD"

ciona.summary <- merge(
  ciona.mean,
  ciona.sd,
  by = "Location"
)
ciona.summary

##     Location  Mean       SD
## 1    Beverly 32.25 47.38518
## 2 Gloucester 37.95 34.27823

ggplot(
  ciona.summary,aes(x = Location, y = Mean)) +
  geom_col() +
  geom_errorbar(
    aes(
      ymin = Mean - SD, ymax = Mean + SD),
    width = 0.2
  ) +
  labs(
    x = "Location",
    y = "Mean Ciona abundance",
    title = "Mean Ciona Abundance by Location"
  ) +
  theme_minimal()

Question: How does average Ciona abundance differ between Beverly and Gloucester?

Average Ciona abundance is slightly higher in Gloucester (Mean = 37.95) than in Beverly (Mean = 32.25). However, there is more variability in Beverly (SD = 47.385) than in Gloucester (SD = 34.278). Therefore, there seems to be little visual difference between the bars themselves, with more difference in the error bars.

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
  ) +
  labs(
  x = expression(italic(Ciona) ~ "abundance"),
  y = expression(italic(Ascidiella) ~ "abundance"),
  title = "Relationship Between Ciona and Ascidiella"
  ) +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

Question: how would you describe the relationship between Ciona and Ascidiella abundance?

The fitted line has a slightly negative slope, but plot points are widely scattered and the relationship is weak. Therefore, there appears to be little to no relationship between the abundance of Ciona and Ascidiella.

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!

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)

fouling$Months <- ifelse(
  grepl("12", fouling$`Unique.ID`), 12,
  ifelse(
    grepl("6", fouling$`Unique.ID`), 6,
    ifelse(
      grepl("2", fouling$`Unique.ID`), 2,
      1
    )
  )
)

With the column added, make sure the updated dataset is back in R.

str(fouling)

## 'data.frame':    40 obs. of  21 variables:
##  $ Unique.ID          : chr  "G1A" "G1B" "G1C" "G1D" ...
##  $ Location           : chr  "Gloucester" "Gloucester" "Gloucester" "Gloucester" ...
##  $ 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             : num  1 1 1 1 1 1 1 1 2 2 ...

Code: convert the new column to a factor

fouling$Months <- as.factor(fouling$Months)
str(fouling$Months)

##  Factor w/ 4 levels "1","2","6","12": 1 1 1 1 1 1 1 1 2 2 ...

Code: run a two-way ANOVA that looks at how the interaction of the new column AND Location influence the abundance of Ciona.

ciona.anova <- aov(
  Ciona ~ Months * Location,
  data = fouling
)
summary(ciona.anova)

##                 Df Sum Sq Mean Sq F value Pr(>F)  
## Months           3  10637    3546   2.845 0.0531 .
## Location         1      0       0   0.000 0.9990  
## Months:Location  3  14792    4931   3.956 0.0166 *
## Residuals       32  39882    1246                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Question: what does your output tell you about Ciona abundance?

The two-way ANOVA indicates a link between Months and Location on Ciona abundance (p = 0.0166). Therefore, Ciona abundance changes over time depending on location.

Code: create a visual that helps represent how the new column and Location relate to abundance.

ciona.mean.time <- aggregate(
  Ciona ~ Months + Location,
  data = fouling,
  FUN = mean
)
ciona.mean.time

##   Months   Location     Ciona
## 1      1    Beverly   5.50000
## 2      2    Beverly  18.50000
## 3      6    Beverly 100.00000
## 4     12    Beverly  18.62500
## 5      1 Gloucester  46.00000
## 6      2 Gloucester  45.16667
## 7      6 Gloucester  36.16667
## 8     12 Gloucester  21.75000