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] 4

When 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.

library(ggplot2)
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, union

library(tidyr)

theme_set(theme_classic()) #this is me being picky - I like the classic or theme_bw plots best

Reading 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.csv")
Assignment_dat <- read.csv("Assignment_final.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 like

Now 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 every 8 hours.
## 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.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.

# Create a dataframe with clustering measures at the quadrat level

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 grouped output by 'Site'. You can override using the
## `.groups` argument.

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.

Clustering_plotdat <- Summary_df %>%
  group_by(Site) %>%
  summarise(Mean_ID = mean(ID_count,na.rm=TRUE),
            SD_ID = sd(ID_count,na.rm=TRUE),
            N_ID = length(!is.na(ID_count)),
            SE_ID = SD_ID/sqrt(N_ID))

Clustering_plotdat

ggplot(Clustering_plotdat,aes(x=Site,y=Mean_ID,fill=Site))+
  geom_bar(stat="identity")+
  geom_errorbar(aes(ymin=Mean_ID-SE_ID,
                    ymax=Mean_ID+SE_ID),
                width=0.3,size=1)+
  theme_classic()+
  theme(legend.position="none")+
  xlab("Site")+
  ylab("Index of Dispersion")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_bar()`).

Figure 1. Mean snail clustering at each sampling site, measured using the index of dispersion. Higher values represent greater snail clustering. Error bars represent ± SE.

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

Density_plotdat <- Snail_Survey %>%
  group_by(Site) %>%
  summarise(Mean_Snails = mean(Snail_Num,na.rm=TRUE),
            SD_Snails = sd(Snail_Num,na.rm=TRUE),
            N_Snails = length(!is.na(Snail_Num)),
            SE_Snails = SD_Snails/sqrt(N_Snails))

Density_plotdat

ggplot(Density_plotdat,aes(x=Site,y=Mean_Snails,fill=Site))+
  geom_bar(stat="identity")+
  geom_errorbar(aes(ymin=Mean_Snails-SE_Snails,
                    ymax=Mean_Snails+SE_Snails),
                width=0.3,size=1)+
  theme_classic()+
  theme(legend.position="none")+
  xlab("Site")+
  ylab("Snail Density (snails/m²)")

Figure 2. Mean snail density at each sampling site in snails per square meter. Each quadrat measured 1 m × 1 m. Error bars represent ± 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?

**CP had the most snails and also the nost clustering. This is what I expected because a site with more snails has more oppurtunities for the snails to be grouped together.

Q2: What else might affect the snail counts that should be included in a more rigorous statistical analysis?

**Other things that could affect snail counts include the amount of spartina, food availability, predators, water level and other differences in habitat between the sites.

What is the size-structure snails? How do we expect it to affect Spartina primary production?

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.

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_data

Now 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).

Size_plotdat <- Size_data %>%
  group_by(Site) %>%
  summarise(Mean_Length = mean(Length,na.rm=TRUE),
            SD_Length = sd(Length,na.rm=TRUE),
            N_Length = length(!is.na(Length)),
            SE_Length = SD_Length/sqrt(N_Length))

Size_plotdat

ggplot(Size_plotdat,aes(x=Site,y=Mean_Length,fill=Site))+
  geom_bar(stat="identity")+
  geom_errorbar(aes(ymin=Mean_Length-SE_Length,
                    ymax=Mean_Length+SE_Length),
                width=0.3,size=1)+
  theme_classic()+
  theme(legend.position="none")+
  xlab("Site")+
  ylab("Mean Snail Length (cm)")

Figure 3. Mean snail length in centimeters at each sampling site. Error bars represent ± SE.

Q3: Does snail size appear to differ between sites? If so, what do you think might explain the differences?

**Yes, snail size appears to differ somewhat between sites. This could be caused by differences in food availability, habitat conditions, or the ages of the snails at each site.

#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 scaling

You 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

# Find average snail size at each site
Interaction_df <- Size_data %>%
  group_by(Site) %>%
  summarise(Length = mean(Length, na.rm=TRUE))

# Add snail density
Interaction_df <- merge(Interaction_df, Density_plotdat, by="Site")

# Calculate biomass
Interaction_df <- Interaction_df %>%
  mutate(IB = 6E-05*Length^3.0155)

# Calculate metabolic biomass
Interaction_df <- Interaction_df %>%
  mutate(MBM = IB^.75)

# Calculate total metabolic biomass
Interaction_df <- Interaction_df %>%
  mutate(total_MBM = MBM*Mean_Snails)

# Calculate interaction strength
Interaction_df <- Interaction_df %>%
  mutate(Interaction_Strength = .45*-0.047*total_MBM)

# Show final results
Interaction_df

Q4: Do snails benefit or harm S. alterniflora at SMCM? How can you tell?

**The snails seem to negatively affect S. alterniflora at SMCM because the interaction strength values are negative. This means that Littoraria reduces the primary production of Spartina. Areas with more snails and higher metabolic biomass would likely experience a greater negative effect on Spartina growth.

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!

Snail Analysis

Halimah Hajhassan