Snail Analysis
Hope Parrish
Lab overview
This document will walk you through the creation of two figures for the snail distribution lab. You will need to install the packages “tidyr”, “dplyr” and “ggplot”.
Before the start of lab next week, you will need to upload either:
This compiled document (html) with figures, captions, and answers to the post-lab questions.
A word document with figures, captions, and answers to the questions (3) in the document.
As a reminder, a figure caption should include (from the SMCM Biology Style Manual):
Figure captions should largely stand alone, and the reader should be able to understand what a figure is showing without reading the results section. Captions are placed below the figure, and represent a concise description of the data being displayed, including the following components:
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)
If a statement is made that is supported by statistical analysis, include the test and P value here in addition to in the text of the Results section (e.g. t-test: P < 0.001) Note: for this lab we will be ignoring a statistical test - but keep this in mind for future labs.
If more than one dataset is presented, provide legend details (e.g. gray circles = elevated temperature treatment, black circles = control)
In today’s lab we are going to use the data collected to answer the following questions:
How does the overall number of snails, snail clustering, and snail size vary across sites?
How do we expect Littoraria to affect Spartina on St. Mary’s campus?
Introduction and R Review
For this lab, we will be using ggplot and dplyr to look at the water quality data we collected.
This is an R Markdown file, which “chunks” of code included. To run a chunk, click the green arrow in the upper right hand corner. You will see the result below the chunk. Notice the # that indicates a “comment”, or a note about the code that isn’t part of the code itself.
2+2 # add 2 and 2 together## [1] 4When you are completley finished with the lab, click “Knit” at the top, and submit the html file to Blackboard. Don’t reinvent the wheel, use the examples I give you to make your own figures.
Example figure and analysis
Loading the packages
We will use a couple of different packages - ggplot and dplyr. Let’s make sure they are installed and loaded.
To install a package go to the “Tools -> Install Packages” and type in “ggplot2”, “dplyr”, “tidyr”. They might already be installed on your machine.
install.packages("ggplot2")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)library(ggplot2)
install.packages("dplyr")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)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, unioninstall.packages("tidyr")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)library(tidyr)
theme_set(theme_classic()) #this is me being picky - I like the classic or theme_bw plots bestReading in the datasets
To read in the data, make sure that your data file and this R markdown file are saved in the same directory. Then set your working directory by clicking “Session -> Set Working Directory -> To Source File Location” with this file up. Now your working directory is wherever this file is saved, and R knows where to look for the datafiles we ask it to get for us. We are going to read in two datasets - one comparing the scores of three color groups, and the other looking at a score on an assignment and the final score in a course.
Color_dat <- read.csv("Color_groups (1).csv")
Assignment_dat <- read.csv("Assignment_final (1).csv")We can look at the structure of each dataset using the “str” command.
str(Color_dat)## 'data.frame': 30 obs. of 2 variables:
## $ Color: chr "Blue" "Blue" "Blue" "Blue" ...
## $ Score: num 0.798 0.402 0.717 1.483 1.107 ...str(Assignment_dat)## 'data.frame': 30 obs. of 2 variables:
## $ Assignment: num 3.218 0.422 5.892 4.84 5.775 ...
## $ Final : num 9.66 1.16 36.91 22.64 30.71 ...Example barplot with errorbars
Making our barplot for the color data will have two components. First we will summarize the data. Then, we will create the plot
To summarize the data, we’ll use what’s called a “pipe” from dplyr. This will let us group the data by color, and get the mean and standard error of each group. This is a fantastic way to organize your data and have it be traceable. You can see the tricks here: https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf
Color_plotdat <- Color_dat %>% #this line tells R we're making a new dataframe and to start with the Color_dat file. The "%>%" is the pipe
group_by(Color) %>% # group the data by color
summarise(Mean_Score = mean(Score,na.rm=TRUE), # mean of each group
SD_Score = sd(Score,na.rm=TRUE), # standard deviation of each group
N_Score = length(!is.na(Score)), # number of replicates in each group
SE_Score = SD_Score/sqrt(N_Score)) # standard error of each group
Color_plotdat #display the data so you can see what the result looks likeNow that we have the data formatted, we can make our figure.
Make the plot
Now that we have our means, we can make a plot using ggplot, including errorbars.
ggplot(Color_plotdat,aes(x=Color,y=Mean_Score,fill=Color))+ # gives the data for the x and y axes, and the fill as the Color color
geom_bar(stat="identity")+ # adds the bars to the plot
geom_errorbar(aes(ymin=Mean_Score-SE_Score,ymax=Mean_Score+SE_Score),width=0.3,linewidth=1) # adds the error bars to the plot
Since these colors are a bit confusing, we can change the color palette.
ggplot(Color_plotdat,aes(x=Color,y=Mean_Score,fill=Color))+ # gives the data for the x and y axes, and the fill as the Color color
geom_bar(stat="identity")+ # adds the bars to the plot
geom_errorbar(aes(ymin=Mean_Score-SE_Score,ymax=Mean_Score+SE_Score),width=0.3,size=1) + # adds the error bars to the plot
scale_fill_manual(values = c("darkblue", "darkgreen", "darkred")) #sets a custom color scheme## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
The legend is redundant, and the axes labels are straight from the dataframe, so let’s clean it up. The gray background is also ugly, so let’s use theme_classic() to clean it up.
ggplot(Color_plotdat,aes(x=Color,y=Mean_Score,fill=Color))+ # gives the data for the x and y axes, and the fill as the Color color
geom_bar(stat="identity")+ # adds the bars to the plot
geom_errorbar(aes(ymin=Mean_Score-SE_Score,ymax=Mean_Score+SE_Score),width=0.3,size=1) + # adds the error bars to the plot
scale_fill_manual(values = c("darkblue", "darkgreen", "darkred"))+ #sets custom colors
theme_classic()+theme(legend.position="none")+
xlab("Color")+
ylab("Score")
> Add Caption Here
Reading in the marsh survey data and loading packages (don’t forget to set your working directory! You should save your csv and this markdown in the same folder)
library(dplyr)
library(ggplot2)
library(tidyr)
Snail_Survey <- read.csv("Snail_survey_data (1).csv")At what site is there the most spatial variation in snail density?
Today we’re using dplyr again to summarize the data, with a new commands: summarize. Summarize will allow you to do calculations based on the values in other rows.
Below, I will give you summary code for the quadrats nested within sites. You will get a value of “NaN” if there are no snails at a site, since you can’t divide by 0.
Now, your turn is to plot snail clustering at each site. You’ll need to use summary and plotting code from the R review above.
# Create a dataframe with clustering measures at the quadrat level
str(Snail_Survey)## 'data.frame': 81 obs. of 16 variables:
## $ Group : chr "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" ...
## $ Site : chr "CP" "CP" "CP" "CP" ...
## $ Quadrat : int 1 1 1 2 2 2 3 3 3 1 ...
## $ Sub_Quadrat: int 1 2 3 1 2 3 1 2 3 1 ...
## $ Stalk_Num : int 12 13 25 20 12 15 4 5 6 5 ...
## $ Snail_Num : int 200 200 124 140 109 90 42 10 27 23 ...
## $ Length_cm : num 1.3 1.4 1.7 1 1.7 2 1.8 2.1 2.3 1.6 ...
## $ Length2_cm : num 1.8 1.9 1.8 1.9 1.9 2.2 1.9 1.7 2.4 2 ...
## $ Length3_cm : num 1.5 1.6 0.9 2 1.8 1.9 2.2 1.9 2.2 2.3 ...
## $ Length4_cm : num 1.5 2.2 1.4 1 1.8 1.9 2.1 1.8 2.2 2.5 ...
## $ Length5_cm : num 2.5 1.2 1.6 0.8 0.7 1.7 2 2 2.1 2 ...
## $ Length6_cm : num 2 1.3 1.9 2.3 1.2 1.8 1.9 1.8 2.2 2 ...
## $ Length7_cm : num 0.6 1 0.9 1.7 1.7 0.9 1 0.9 1.9 2.2 ...
## $ Length8_cm : num 2 2.2 1.8 1.9 1.6 1.3 1.8 1.7 2.1 1.9 ...
## $ Length9_cm : num 1.5 1.7 2 0.9 1 2.1 2 2.2 2.2 1.7 ...
## $ Length10_cm: num 1.4 1.3 2.1 1.5 2 1.3 1.9 1.9 2.2 2 ...Summary_df <- Snail_Survey %>% # start with the snail survey dataframe and store it as "Quadrat_df"
group_by(Site,Quadrat) %>% # group by the quadrat (I already made each group's quadrats unique)
summarize(ID_count = var(Snail_Num)/mean(Snail_Num))# compute the index of dispersion for each quadrat within a site## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by Site and Quadrat.
## ℹ Output is grouped by Site.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(Site, Quadrat))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead. Site_plotdat <- Summary_df %>% #this line tells R we're making a new dataframe and to start with the Color_dat file. The "%>%" is the pipe
group_by(Site) %>% # group the data by color
summarise(Mean_ID_count = mean(ID_count,na.rm=TRUE), # mean of each group
SD_ID_count = sd(ID_count,na.rm=TRUE), # standard deviation of each group
N_ID_count = length(!is.na(ID_count)), # number of replicates in each group
SE_ID_count = SD_ID_count/sqrt(N_ID_count)) # standard error of each group
Site_plotdatggplot(Site_plotdat,aes(x=Site,y=Mean_ID_count))+ # gives the data for the x and y axes, and the fill as the Color color
geom_bar(stat="identity")+ # adds the bars to the plot
geom_errorbar(aes(ymin=Mean_ID_count-SE_ID_count,ymax=Mean_ID_count+SE_ID_count),width=0.3,size=1) + # adds the error bars to the plot
scale_fill_manual(values = c("darkblue", "darkgreen", "darkred"))+ #sets custom colors
theme_classic()+theme(legend.position="none")+
xlab("Site")+
ylab("Clustering of Snails")## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_bar()`).
**Figure 1. Mean clustering of snails among quadrats at each sampling site. Clustering was quantified using the index of dispersion calculated for each quadrat. Bars represent the mean index of dispersion ± standard error (SE) across quadrats within each site.**
Now repeat to create a figure for the overall number of snails (hint: remember mean() like you used in the index of dispersion?) Since you're grouping at the quadrat (which was 1m by 1m), be sure to display these numbers as density
``` r
str(Snail_Survey)## 'data.frame': 81 obs. of 16 variables:
## $ Group : chr "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" ...
## $ Site : chr "CP" "CP" "CP" "CP" ...
## $ Quadrat : int 1 1 1 2 2 2 3 3 3 1 ...
## $ Sub_Quadrat: int 1 2 3 1 2 3 1 2 3 1 ...
## $ Stalk_Num : int 12 13 25 20 12 15 4 5 6 5 ...
## $ Snail_Num : int 200 200 124 140 109 90 42 10 27 23 ...
## $ Length_cm : num 1.3 1.4 1.7 1 1.7 2 1.8 2.1 2.3 1.6 ...
## $ Length2_cm : num 1.8 1.9 1.8 1.9 1.9 2.2 1.9 1.7 2.4 2 ...
## $ Length3_cm : num 1.5 1.6 0.9 2 1.8 1.9 2.2 1.9 2.2 2.3 ...
## $ Length4_cm : num 1.5 2.2 1.4 1 1.8 1.9 2.1 1.8 2.2 2.5 ...
## $ Length5_cm : num 2.5 1.2 1.6 0.8 0.7 1.7 2 2 2.1 2 ...
## $ Length6_cm : num 2 1.3 1.9 2.3 1.2 1.8 1.9 1.8 2.2 2 ...
## $ Length7_cm : num 0.6 1 0.9 1.7 1.7 0.9 1 0.9 1.9 2.2 ...
## $ Length8_cm : num 2 2.2 1.8 1.9 1.6 1.3 1.8 1.7 2.1 1.9 ...
## $ Length9_cm : num 1.5 1.7 2 0.9 1 2.1 2 2.2 2.2 1.7 ...
## $ Length10_cm: num 1.4 1.3 2.1 1.5 2 1.3 1.9 1.9 2.2 2 ...Summary_df <- Snail_Survey %>% # start with the snail survey dataframe and store it as "Quadrat_df"
group_by(Site,Quadrat) %>% # group by the quadrat (I already made each group's quadrats unique)
summarize(Snail_density = mean(Snail_Num))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by Site and Quadrat.
## ℹ Output is grouped by Site.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(Site, Quadrat))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead. Site_plotdat <- Summary_df %>% #this line tells R we're making a new dataframe and to start with the Color_dat file. The "%>%" is the pipe
group_by(Site) %>% # group the data by color
summarise(Mean_density = mean(Snail_density, na.rm=TRUE), # mean density at each site
SD_density = sd(Snail_density, na.rm=TRUE), # standard deviation
N_density = length(!is.na(Snail_density)), # number of replicates
SE_density = SD_density/sqrt(N_density)) # standard error
Site_plotdatggplot(Site_plotdat,aes(x=Site,y=Mean_density)) +
geom_bar(stat="identity")+ # adds the bars to the plot
geom_errorbar(aes(ymin=Mean_density-SE_density,ymax=Mean_density+SE_density),width=0.3,size=1) + # adds the error bars to the plot
scale_fill_manual(values = c("darkblue", "darkgreen", "darkred"))+ #sets custom colors
theme_classic()+theme(legend.position="none")+
xlab("Site")+
ylab("Snail density")
Figure 2. This graph shows mean snail density across sampling sites. Snail density was calculated as the mean number of snails per quadrat, and bars represent the mean density across quadrats within each site. Error bars represent ± standard error (SE).
Answer the following questions:
Q1: Which site has the most snails? Which site appears has the most clustering? Is that what you expected? Why or why not? Church point has the most snails and the most clustering. This is what I expected because there was not many snails found at the other locations and they were more spread out if there was any. The results made sense based off of observations made during the lab.
Q2: What else might affect the snail counts that should be included in a more rigorous statistical analysis?
Factors could include environmental conditions such as temperature, moisture, vegetation, food availability, habitat type, and substrate. Other factors such as time of day, weather conditions, season, and sampling effort could also influence the number of snails observed.
What is the size-structure snails? How do we expect it to affect Spartina primary production?
The size structure desribes the proportion of large, medium, and small snails in the population. The size structure could influence Spartina primary production becuase a population with more large snails could have a stronger effect on plant growth and productivity than a population with smaller snails.
First, we need to reformat our data so that we have only one column of lengths. We can do that with a command called “pivot longer,” which takes wide format data and makes it long.
str(Snail_Survey)## 'data.frame': 81 obs. of 16 variables:
## $ Group : chr "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" ...
## $ Site : chr "CP" "CP" "CP" "CP" ...
## $ Quadrat : int 1 1 1 2 2 2 3 3 3 1 ...
## $ Sub_Quadrat: int 1 2 3 1 2 3 1 2 3 1 ...
## $ Stalk_Num : int 12 13 25 20 12 15 4 5 6 5 ...
## $ Snail_Num : int 200 200 124 140 109 90 42 10 27 23 ...
## $ Length_cm : num 1.3 1.4 1.7 1 1.7 2 1.8 2.1 2.3 1.6 ...
## $ Length2_cm : num 1.8 1.9 1.8 1.9 1.9 2.2 1.9 1.7 2.4 2 ...
## $ Length3_cm : num 1.5 1.6 0.9 2 1.8 1.9 2.2 1.9 2.2 2.3 ...
## $ Length4_cm : num 1.5 2.2 1.4 1 1.8 1.9 2.1 1.8 2.2 2.5 ...
## $ Length5_cm : num 2.5 1.2 1.6 0.8 0.7 1.7 2 2 2.1 2 ...
## $ Length6_cm : num 2 1.3 1.9 2.3 1.2 1.8 1.9 1.8 2.2 2 ...
## $ Length7_cm : num 0.6 1 0.9 1.7 1.7 0.9 1 0.9 1.9 2.2 ...
## $ Length8_cm : num 2 2.2 1.8 1.9 1.6 1.3 1.8 1.7 2.1 1.9 ...
## $ Length9_cm : num 1.5 1.7 2 0.9 1 2.1 2 2.2 2.2 1.7 ...
## $ Length10_cm: num 1.4 1.3 2.1 1.5 2 1.3 1.9 1.9 2.2 2 ...Size_data <- Snail_Survey %>%
select(-Group) %>% # remove group members names
pivot_longer(cols = starts_with("Length"),names_to = "SnailID",values_to = "Length", values_drop_na = TRUE) # this line takes everything that starts with length and puts it in a column, removing NAs (where there was no snail)
Size_dataNow Your job is to summarize the data and make a figure of the average snail size at each site. You should include error bars and a caption. (Figure 2).
str(Snail_Survey)## 'data.frame': 81 obs. of 16 variables:
## $ Group : chr "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" "MaraHopeLauren" ...
## $ Site : chr "CP" "CP" "CP" "CP" ...
## $ Quadrat : int 1 1 1 2 2 2 3 3 3 1 ...
## $ Sub_Quadrat: int 1 2 3 1 2 3 1 2 3 1 ...
## $ Stalk_Num : int 12 13 25 20 12 15 4 5 6 5 ...
## $ Snail_Num : int 200 200 124 140 109 90 42 10 27 23 ...
## $ Length_cm : num 1.3 1.4 1.7 1 1.7 2 1.8 2.1 2.3 1.6 ...
## $ Length2_cm : num 1.8 1.9 1.8 1.9 1.9 2.2 1.9 1.7 2.4 2 ...
## $ Length3_cm : num 1.5 1.6 0.9 2 1.8 1.9 2.2 1.9 2.2 2.3 ...
## $ Length4_cm : num 1.5 2.2 1.4 1 1.8 1.9 2.1 1.8 2.2 2.5 ...
## $ Length5_cm : num 2.5 1.2 1.6 0.8 0.7 1.7 2 2 2.1 2 ...
## $ Length6_cm : num 2 1.3 1.9 2.3 1.2 1.8 1.9 1.8 2.2 2 ...
## $ Length7_cm : num 0.6 1 0.9 1.7 1.7 0.9 1 0.9 1.9 2.2 ...
## $ Length8_cm : num 2 2.2 1.8 1.9 1.6 1.3 1.8 1.7 2.1 1.9 ...
## $ Length9_cm : num 1.5 1.7 2 0.9 1 2.1 2 2.2 2.2 1.7 ...
## $ Length10_cm: num 1.4 1.3 2.1 1.5 2 1.3 1.9 1.9 2.2 2 ...Summary_df <- Size_data %>%
group_by(Site,Quadrat) %>%
summarize(Mean_Length = mean(Length, na.rm=TRUE))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by Site and Quadrat.
## ℹ Output is grouped by Site.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(Site, Quadrat))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.Site_plotdat <- Summary_df %>%
group_by(Site) %>% # group the data by site
summarise(Avg_Length = mean(Mean_Length, na.rm=TRUE), # mean snail size at each site
SD_Size = sd(Mean_Length, na.rm=TRUE), # standard deviation
N_Size = sum(!is.na(Mean_Length)), # number of replicates
SE_Size = SD_Size/sqrt(N_Size))
Site_plotdatggplot(Site_plotdat, aes(x=Site, y=Avg_Length))+
geom_bar(stat="identity")+
geom_errorbar(aes(ymin=Avg_Length-SE_Size,
ymax=Avg_Length+SE_Size),
width=0.3, size=1)+
theme_classic()+
theme(legend.position="none")+
xlab("Site")+
ylab("Average Snail Size")
**Figure 3. Mean snail size across sampling sites. Snail size was
measured as mean snail length within each quadrat and then averaged
across quadrats at each site. Bars represent mean snail length, and
error bars represent ± standard error (SE).
Q3: Does snail size appear to differ between sites? If so, what do you think might explain the differences?
Yes, snail size appears different between the sitesof water front/church point and Wherrits pond. The most was found at teh water front with church point being a close second. Factors like food availability, habitat conditons, predation, competition, and the environment could all play a role in snail size.
#Effect of Littoraria on Spartina at SMCM
You can calculate the metabolic biomass of Littoraria with this code
mutate(IB = 6E-05*Length^3.0155, #biomass conversion
MBM = IBM^{.75}) #with Metabolic scalingYou can calculate the Interaction Strength using this equation:
mutate(Interaction_Strength = .45*-0.047*total_MBM)Calculate the average interaction strength at each site. You may assume that each snail is the same size within a site
library(dplyr)
Mean_Size <- Size_data %>%
group_by(Site) %>%
summarise(Mean_Length = mean(Length, na.rm = TRUE))
Snail_Count <- Snail_Survey %>%
group_by(Site) %>%
summarise(Total_Snails = sum(Snail_Num, na.rm = TRUE))
Interaction <- Mean_Size %>%
left_join(Snail_Count, by = "Site") %>%
mutate(
IB = 6E-05 * Mean_Length^3.0155,
MBM = IB^0.75,
total_MBM = MBM * Total_Snails,
Interaction_Strength = 0.45 * -0.047 * total_MBM
)
InteractionQ4: Do snails benefit or harm S. alterniflora at SMCM? How can you tell?
Snails are predicted to harm Spartina alterniflora at SMCM. We can tell because the calculated interaction strength is negative. The interaction-strength equation includes a negative effect of Littoraria on Spartina primary production, so more snail biomass results in a more negative interaction strength.
All done! Make sure you have
- 3 Figures and captions (snail density, clustering, and size by site) (please keep these bold so they’re easy for me to find) (4 points each)
- 4 Questions and their answers throughout the markdown document (please keep these bold so they’re easy for me to find) (1-3 are 2 points each, 4 is 4 points)
- Knitted and submitted to Blackboard!