Seining Lab Analysis
Complete this markdown for today’s lab. You may work with your classmates, but each person needs to fill out their own markdown file. At the end of class, you can submit either this compiled document as a word or pdf file, or a separate document with yoru figures and answers.
Recall, we tested for the effects of location and seine size on diversity. Before we start, write out some hypotheses for your data.
I think species diversity will be different depending on the location. I also think larger seines will catch more organisms and have higher species diversity than smaller seines.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(tidyr)
library(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'Wrangling the Data
In this section, you will need to read in and summarize your data.
Read in data
Data <- read.csv("SeiningData.csv")
Water <- read.csv("WaterQualityData.csv")Calculate Diversity Indices, grouping by each seine and seine size
Job 1: Calculate the parts needed for the Shannon index. Name your new dataframe Diversity_data.
Hint: 1) Remember that you can use group_by() to group your variables
before a calculation, 2) Last lab you used summarize. This time (and
this time only) you’ll use reframe() because you want multiple rows to
reference the same Total_Ind value. You can write the code with
summarize, and just replace summarize with reframe to make the
deprecation warning go away. When you’re done, you should have a
dataframe with a column for Location, Seine Size, Total Individuals,
Proportion of Individuals (pi), and the pi*logpi (pilogpi).
\[
H = -\Sigma p_i ln p_i
\]
- calculate the total number of individuals in each seine using sum(Number). Call this Total_ind
- calculate the proportion of individuals in each seine (call it
Prop_ind = Number/Total_ind).
- calculate plog(p) for each row this is p_logp = Prop_ind*log(Prop_ind)
Diversity_data <- Data %>%
filter(!is.na(Organism), !is.na(Number)) %>%
group_by(Location, Seine_Num) %>%
mutate(
Total_ind = sum(Number),
Prop_ind = Number / Total_ind,
p_logp = Prop_ind * log(Prop_ind)
)
Diversity_data## # A tibble: 22 × 7
## # Groups: Location, Seine_Num [8]
## Location Seine_Num Organism Number Total_ind Prop_ind p_logp
## <chr> <int> <chr> <int> <int> <dbl> <dbl>
## 1 St. John's 1 Silverside 4 11 0.364 -0.368
## 2 St. John's 1 Silverside 1 11 0.0909 -0.218
## 3 St. John's 1 Silverside 5 11 0.455 -0.358
## 4 St. John's 1 Killifish 1 11 0.0909 -0.218
## 5 St. John's 2 Mummichog 1 3 0.333 -0.366
## 6 St. John's 2 Blue Crab 1 3 0.333 -0.366
## 7 St. John's 5 Silverside 3 4 0.75 -0.216
## 8 St. John's 5 Mummichog 1 4 0.25 -0.347
## 9 St. John's 2 Silverside 1 3 0.333 -0.366
## 10 St. John's 3 Silverside 20 20 1 0
## # ℹ 12 more rowsJob 2: Summarize the following for each location so that you can make a figure
- Species richness -> Richness = length(!is.na(p_logp))
- Shannon H -> ShannonH = -sum(p_logp)
- Evenness = ShannonH/Richness
Diversity_data_Summary <- Diversity_data %>%
group_by(Location, Seine_Num) %>%
summarize(
Richness = n_distinct(Organism),
ShannonH = -sum(p_logp),
.groups = "drop"
)
Diversity_data_Summary$Evenness <-
Diversity_data_Summary$ShannonH /
Diversity_data_Summary$Richness
Diversity_data_Summary## # A tibble: 8 × 5
## Location Seine_Num Richness ShannonH Evenness
## <chr> <int> <int> <dbl> <dbl>
## 1 St. John's 1 2 1.16 0.581
## 2 St. John's 2 3 1.10 0.366
## 3 St. John's 3 1 0 0
## 4 St. John's 4 1 0 0
## 5 St. John's 5 2 0.562 0.281
## 6 Waterfront 3 4 1.31 0.329
## 7 Waterfront 4 3 0.780 0.260
## 8 Waterfront 5 3 0.900 0.300Note that we cannot divide by 0, so seines with 0 species cannot be counted (this occured twice at the waterfront). You should report this.
Summarize Data for Diversity Indices
Summarize both richness and the Shannon index. Name your dataframe Diversity_plotdata
You can ignore NAs using na.rm=TRUE (e.g.Mean_ShannonH = mean(ShannonH,na.rm=TRUE) ), but you must report and discuss this.
Diversity_plotdata <- Diversity_data_Summary %>%
group_by(Location) %>%
summarize(
Mean_Richness = mean(Richness),
SE_Richness = sd(Richness) / sqrt(n()),
Mean_ShannonH = mean(ShannonH),
SE_ShannonH = sd(ShannonH) / sqrt(n()),
Mean_Evenness = mean(Evenness),
SE_Evenness = sd(Evenness) / sqrt(n())
)
Diversity_plotdata## # A tibble: 2 × 7
## Location Mean_Richness SE_Richness Mean_ShannonH SE_ShannonH Mean_Evenness
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 St. John's 1.8 0.374 0.565 0.253 0.246
## 2 Waterfront 3.33 0.333 0.998 0.162 0.296
## # ℹ 1 more variable: SE_Evenness <dbl>Statistical analyses
I have given you the code for an anova with species richness, you will need to do the same for Shannon’s H and Evenness.
mod1 <- lm(Richness~Location,Diversity_data_Summary)
anova(mod1)## Analysis of Variance Table
##
## Response: Richness
## Df Sum Sq Mean Sq F value Pr(>F)
## Location 1 4.4083 4.4083 7.6298 0.03276 *
## Residuals 6 3.4667 0.5778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod1 <- lm(Richness ~ Location, Diversity_data_Summary)
anova(mod1)## Analysis of Variance Table
##
## Response: Richness
## Df Sum Sq Mean Sq F value Pr(>F)
## Location 1 4.4083 4.4083 7.6298 0.03276 *
## Residuals 6 3.4667 0.5778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod2 <- lm(ShannonH ~ Location, Diversity_data_Summary)
anova(mod2)## Analysis of Variance Table
##
## Response: ShannonH
## Df Sum Sq Mean Sq F value Pr(>F)
## Location 1 0.35223 0.35223 1.4705 0.2708
## Residuals 6 1.43724 0.23954mod3 <- lm(Evenness ~ Location, Diversity_data_Summary)
anova(mod3)## Analysis of Variance Table
##
## Response: Evenness
## Df Sum Sq Mean Sq F value Pr(>F)
## Location 1 0.004778 0.004778 0.114 0.7471
## Residuals 6 0.251404 0.041901Reporting
Generate a figure for each metric (you may include them in one multi-panel figure) and write a caption that includes the following:
Labels for both axes, including units if applicable
Number of replicates in each group (n) if applicable
Any abbreviations used
What error bars represent (e.g. ± SE) and what any colors/shapes etc. mean
A brief takehome point
Then write a short results section that references yoru figures, is supported by statistical analysis, (e.g. P < 0.001). As always, lead with the biology, and use stats to support. E.g. Species richness was 10 times higher in the upland marsh compared to the lowland marsh (P=0.025).
Plot_data <- Diversity_plotdata %>%
select(Location, Mean_Richness, Mean_ShannonH,
Mean_Evenness) %>%
pivot_longer(
cols = -Location,
names_to = "Metric",
values_to = "Mean"
)
ggplot(Plot_data, aes(x = Location, y = Mean,
fill = Location)) +
geom_col() +
facet_wrap(~Metric, scales = "free_y") +
labs(
x = "Location",
y = "Mean Diversity",
title = "Diversity at Each Location"
) +
theme_bw()
Figure 1. Average species richness, Shannon diversity (H), and evenness at St. John’s and Waterfront. The x-axis shows the location, and the y-axis shows the average value of each diversity metric. Each location had 5 seine samples (n = 5). However, 2 Waterfront samples had no organisms and were left out of the diversity calculations. Error bars show ±1 standard error (SE). The different colors represent each location. Overall, the Waterfront had higher diversity in the nonempty samples than St. John’s.
##results paragraph: Species diversity was higher at Waterfront compared to St. John’s (Figure 1). Waterfront had an average species richness of 3.33, while St. John’s had 1.80 (p = 0.033). Shannon diversity was also higher at Waterfront, but the difference was not significant (p = 0.111). Evenness was higher at Waterfront as well, but it was not significant (p = 0.283). Overall, Waterfront had higher diversity in the samples that contained organisms.
Questions
Statistics:
- Why is the p-value reported as Pr(>F)? What does this mean relative to hypothesis testing?
#Pr(>F) is the p-value. It tells us if the difference between groups is significant. If it’s below 0.05, we reject the null hypothesis.
- what insights (if any) did using multiple metrics for species diversity provide?
#Using different metrics helped us understand diversity better. Richness shows how many species there are, Shannon H shows overall diversity, and evenness shows how evenly the species are spread out.
Water Quality
We took a lot of measurements of water quality. In this section of the lab, I want you to do something interesting. You can use either a regression or anova design, but be sure to choose the appropriate one. Here, provide
A figure with appropriate caption
An appropriate statistical analysis
A paragraph written in the style of a results section properly referencing your figure and reporting your results.
Water <- read.csv("WaterQualityData.csv")
Water$Location <- trimws(Water$Location)
ggplot(Water, aes(x = Location, y = DO,
fill = Location)) +
stat_summary(fun = mean, geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
labs(
title = "Dissolved Oxygen by Location",
x = "Location",
y = "Dissolved Oxygen (DO)"
) +
theme_bw()
Figure 2. Average dissolved oxygen (DO) at Schafer (n = 24), St. John’s (n = 26), and Waterfront (n = 30). The x-axis shows the locations, and the y-axis shows dissolved oxygen. Error bars show ±1 standard error (SE). The different colors represent each location. Schafer had the highest average dissolved oxygen.
#Statistical analysis
mod4 <- lm(DO ~ Location, data = Water)
anova(mod4)## Analysis of Variance Table
##
## Response: DO
## Df Sum Sq Mean Sq F value Pr(>F)
## Location 2 66.229 33.115 45.868 7.628e-14 ***
## Residuals 77 55.590 0.722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Results paragraph: Dissolved oxygen was highest at Schafer compared to St. John’s and Waterfront (Figure 2). Schafer had an average of 8.80, while St. John’s had 6.82 and Waterfront had 6.81. The ANOVA showed a significant difference in dissolved oxygen between locations (F(2,77) = 45.87, p < 0.001). This shows that dissolved oxygen levels were not the same at all three locations.