# First we load the necessary libraries
library(behaviouR)
library(ggpubr)
# Then we read in our data
data("CensusData")
# Part 1. Population density estimation. First let's focus on your data NOTE you
# need to change 'A' to the name you used
MyCensusData <- subset(CensusData, Partner == "A")
# Now we will subset the data to focus on your focal species for density
# estimation NOTE in this example we subsetting the data so that we only have the
# squirrels; you will need to change the code to subset based on your focal
# species
MyCensusDataFocal <- subset(MyCensusData, Species == "squirrel")
# Here we will determine the width of our sampling area by creating a histogram
# of perpendicular detection distances.
hist(MyCensusDataFocal$PerpendicularDistance, xlab = "Perpendicular distance (D)",
ylab = "Number of observations", main = "")
# In this example there is a clear break between 15 and 20 meters, so we will
# only use observations that were within 15 meters. NOTE: Your data will look
# different than this!
# Change the value to the cutoff point indicated in your data
CutOffPoint <- 15
MyCensusDataFocalAdjusted <- subset(MyCensusDataFocal, PerpendicularDistance < CutOffPoint)
# Now we need to subset by each site NOTE: You will need to change to the site
# names you used
MyCensusDataFocalSiteA <- subset(MyCensusDataFocalAdjusted, Site == "SiteA")
MyCensusDataFocalSiteB <- subset(MyCensusDataFocalAdjusted, Site == "SiteB")
# Now we will calculate the population density based on our two surveys Change
# the following to indicate the distance (in meters) of your survey for site A
SiteACensusDistance <- 500
# Now we will calculate the area of our census. The sample area (a) is equal to
# the length of the transect multiplied by twice the width or a= 2wl. We divide
# by 1000 to convert our answer to square kilometers.
SiteACensusArea <- 2 * SiteACensusDistance * CutOffPoint/1000
# Now we need to calculate the number of animals using the following code
NumberFocalAnimalsSiteA <- nrow(MyCensusDataFocalSiteA)
PopulationDensitySiteA <- NumberFocalAnimalsSiteA/SiteACensusArea
PopulationDensitySiteA
SiteBCensusDistance <- 500
# Now we will calculate the area of our census. The sample area (a) is equal to
# the length of the transect multiplied by twice the width or a= 2wl
SiteBCensusArea <- 2 * SiteBCensusDistance * CutOffPoint/1000
# Now we need to calculate the number of animals using the following code
NumberFocalAnimalsSiteB <- nrow(MyCensusDataFocalSiteB)
PopulationDensitySiteB <- NumberFocalAnimalsSiteB/SiteBCensusArea
PopulationDensitySiteB
# Now we can compare population density using the following code. Was the
# population density of site A higher than site B?
PopulationDensitySiteA > PopulationDensitySiteB
# Was the population density of site B higher than site A?
PopulationDensitySiteB > PopulationDensitySiteA
# **Question 1**. What were the population density estimates (reported as number
# of individuals per square kilometer) for your two sites? Do your results of
# the population density estimates match your predictions? Why or why not?
# Part 2. Comparing biodiversity
# First, we will estimate the alpha diversity, or the diversity within a
# particular area or ecosystem. The alpha diversity is simply the number of
# different species present at each site Here we subset by partner and site A
# (you may need to change these!)
MyCensusDataSiteA <- subset(MyCensusData, Partner == "A" & Site == "SiteA")
# What were the unique species present?
unique(MyCensusDataSiteA$Species)
# How many unique species were there?
SiteANumberSpecies <- length(unique(MyCensusDataSiteA$Species))
SiteANumberSpecies
# Here we subset by partner and site B (you may need to change these!)
MyCensusDataSiteB <- subset(MyCensusData, Partner == "A" & Site == "SiteB")
# What were the unique species present?
unique(MyCensusDataSiteB$Species)
# How many unique species were there?
SiteBNumberSpecies <- length(unique(MyCensusDataSiteB$Species))
SiteBNumberSpecies
# What was the species richness for both sites sampled?
unique(c(MyCensusDataSiteA$Species, MyCensusDataSiteB$Species))
# **Question 2**. Which of your sites had higher species richness (i.e. number of
# species)?
# Now we will estimate beta diversity, which estimates changes in species
# diversity between ecosystems or along environmental gradients
# Beta diversity = 2c / S1 + S2
# This code tells us which species both sites have in common
intersect(unique(MyCensusDataSiteA$Species), unique(MyCensusDataSiteB$Species))
# Now we calculate the number of species in common
SpeciesInCommonBothSites <- length(intersect(unique(MyCensusDataSiteA$Species), unique(MyCensusDataSiteB$Species)))
SpeciesInCommonBothSites
# Scientists are often interested in the community similarity between sites. For
# this we will calculate Sørenson’s index; a value of 1 means exactly the same
# number of species a value of 0 means no overlap
2 * SpeciesInCommonBothSites/(SiteANumberSpecies + SiteBNumberSpecies)
# Now we will use your partner's data
MyPartnersCensusDataSiteA <- subset(CensusData, Partner == "B" & Site == "SiteC")
unique(MyPartnersCensusDataSiteA$Species)
length(unique(MyPartnersCensusDataSiteA$Species))
SiteANumberSpecies <- length(unique(MyPartnersCensusDataSiteA$Species))
SiteANumberSpecies
MyPartnersCensusDataSiteB <- subset(CensusData, Partner == "B" & Site == "SiteD")
unique(MyPartnersCensusDataSiteB$Species)
length(unique(MyPartnersCensusDataSiteB$Species))
SiteBNumberSpecies <- length(unique(MyPartnersCensusDataSiteB$Species))
SiteBNumberSpecies
# What was the species richness for both sites sampled?
unique(c(MyPartnersCensusDataSiteA$Species, MyPartnersCensusDataSiteB$Species))
# **Question 3**. How did the alpha diversity of each of your sites compare with
# that of your partner?
# Gamma diversity is total species over a large area or region; there are many
# different ways that this can be measured. The way we will do it is a bit of an
# oversimplification by simply comparing the number of species seen during the
# census at both locations.
# First we subset by the first partner
MyCensusData <- subset(CensusData, Partner == "A")
unique(MyCensusData$Species)
# Then we plot the results
gghistogram(data = MyCensusData, x = "Species", stat = "count")
# Now we subset by the second partner
MyPartnersCensusData <- subset(CensusData, Partner == "B")
unique(MyPartnersCensusData$Species)
# Then we plot the results
gghistogram(data = MyPartnersCensusData, x = "Species", stat = "count")
# Now we can plot all the data together, separated by partner
gghistogram(data = CensusData, x = "Species", stat = "count", facet.by = "Partner",
fill = "Partner", x.text.angle = 90) + xlab("Species") + ylab("Number of individuals")
# **Question 4**. How did the gamma diversity of your site compare with that of
# your partners?
# There are special packages in R that can measure different diversity indices,
# as this is a tool many ecologists use. We will use the 'vegan' package.
library(vegan)
# First we need to convert our data into a table thate can be used to calculate
# the indices.
BiodiversityTable <- table(CensusData$Site, CensusData$Species)
# Simpson's Index (D) measures the probability that two individuals randomly
# selected from a sample will belong to the same species (or some category other
# than species). With this index, 1 represents infinite diversity and 0 means no
# diversity.
H <- diversity(BiodiversityTable, index = "simpson")
H
# Species evenness refers to how close in numbers each species in an environment
# is. We can calculate that using the following code. The value is constrained
# between 0 and 1, with communities that that have a more even representation of
# species having values closer to 1.
J <- H/log(specnumber(BiodiversityTable))
J
# **Question 5**. Which of the four sites was most diverse? Which was the most
# even? Why is it important to consider diversity and evenness when studying
# biodiversity?