ProgInRBCIfunc: species_richness_1

#' ---
#' title: 03 Species Richness, Simpson Diversity, Shannon entropy and Berger-Parker index
#' author: Edie Bishop
#' ---

#' # Introduction
#' In this package we will use the ProgInRBCIfunc package to call on four different functions, species_richness3,
#' simpson_diversity3, shannon_entropy3 and berger-parker3, each call on universal_diversity3 to calculate
#' a diversity measure. universal_diversity3 calculates a diversity measure, represented by the following
#' equation:
#'
#' $$D_{q} = \left(\sum_{i\in\{1...N_{s}, p_{i}\neq0 } p_{i}^q\right)^{1/1-q}$$
#'
#' *Where $N_{s}$ is the number of species we have observed at least one of, and $p_{i}$ is the proportion
#' of individuals that are in the $i^{th}$ species. $D_{q}$ is then the $q^{th}$ diversity measure.*
#'
#' Each individual diversity function will pass a specific value for q, and then modify the calculated
#' diversity measure appropriately such that: \
#' **Species Richness** = $D_{0}$ *(0 - $N_{s}$)* = total number of species within a community. \
#' **Simpson DIversity** = $1/D_{2}$ *(1-0)* = probability that two randomly selected individuals will be from the same
#' species. \
#' **Shannon Entropy** = $\log(D_{1})$ *(0-10)* = how surprising the species of the next individual observed
#' in a sequence is likely to be.\
#' **Berger-Parker Index** = 1/$D_{\infty}$ *(0-1)* = tells us how dominant the most dominant species is. \
#'
#' We will also call on the script "create_datasets.R" in this project, which creates a number of vector
#' files, which contain populations generated from the BCI_2010 dataset (which we obtained from the BCI
#' package).
#'
#' For each population, we will first plot the population, so that we know what to expect when we
#' calculate each diversity measure.
#'
#' Species richness is a count of the total number of species within a community. Simpson diversity is
#' a measure of the probability that two individuals randomly selected from a sample will be from the
#' same species, and therefore will provide a value between zero and one.
#'
#' # Load required packages and scripts

# The ProgInRBCIfunc function contains functions which calculate various indices of diversity
library(ProgInRBCIfunc)
# the "create_datasets" script produces a list of populations
source("create_datasets.R")

#' #### Tree.pop
#' tree.pop is a vector which contains the count for each species in the entire BCI_2010 dataset
plot(sort(tree.pop))
#' There are a lot of species in the community which only have a few individuals, and there are only a few
#' species which have a lot of individuals within the community
species_richness3(tree.pop)
# As we know that there are 301 different species in the community, the species richness = 301 as expected
simpson_diversity3(tree.pop)
# 0.449 indicates that there is an almost 50% chance that two randomly selected individuals from the
# community will be from the same species.
shannon_entropy3(tree.pop) # = 4.019
berger_parker3(tree.pop) # = 0.135
#'
#'------------------------------------------------------------------------------------
#'
#' #### quadrat.pop
#' quadrat.pop was generated by selected one quadrat at random, and counting the number of individuals
#' from each species within that quadrat.
plot(sort(quadrat.pop))
#' As we might expect, the graph looks similar to that of tree-pop. There were a few species with many
#' individuals within this quadrat, and many species with only a few individuals.
species_richness3(quadrat.pop)
# There are less species within the randomly selected quadrat than the whole of the community. This value
# will change each time quadrat.pop is created.
simpson_diversity3(quadrat.pop)
# This value will change each time quadrat.pop is created
shannon_entropy3(quadrat.pop) # = 3.29
berger_parker3(quadrat.pop) # = 0.109
#'
#'--------------------------------------------------------------------------------------
#'
#' #### quadrat10.pop
#' quadrat10.pop is similar to quadrat.pop, although this time it has drawn from 10 random quadrats and
#' summed the number of each species within those 10 quadrats
plot(sort(quadrat10.pop))
#' The plot still looks fairly similar to the previous two
species_richness3(quadrat10.pop)
# More species than quadrat.pop, which makes sense. This value will change each time quadrat50.pop is
# created.
simpson_diversity3(quadrat10.pop)
# And were back to around 50% chance of randomly selecting two individuals from the same species. This
# value will change each time quadrat50.pop is created.
shannon_entropy3(quadrat10.pop) # = 3.904
berger_parker3(quadrat10.pop) # = 0.1708
#'
#'-------------------------------------------------------------------------------------
#'
#' #### one.pop
#' one.pop is a population with are large number of individuals of only one sepcies
plot(sort(one.pop))
#' Values for all species are 0, except for the one species which contains 222718 individuals
species_richness3(one.pop)
# Vlaue is 1, as there is only one species
simpson_diversity3(one.pop)
# Value is 1, since two randomly selected individuals are always going to be from the same species
shannon_entropy3(one.pop) # = 0
berger_parker3(one.pop) # = 1
#'
#'-------------------------------------------------------------------------------------
#'
#' #### uneven.pop
#' uneven.pop contains one species with 111728 individuals, and 300 species with 369 individuals each
plot(sort(uneven.pop))
# The plot looks similar to the last one, since 369 is so much smaller a number than 111728
species_richness3(uneven.pop)
# Species richness is again 301, as there are again 301 species present in the population
simpson_diversity3(uneven.pop)
# There is around a 25% chance of selecting two individuals from the same species
shannon_entropy3(uneven.pop) # = 3.535
berger_parker3(uneven.pop) # = 0.5016
#'
#'-------------------------------------------------------------------------------------
#'
#' #### mixed.pop
#' mixed.pop contains a population with a variable number of individuals for each species, no two species
#' have the same number of individuals
plot(sort(mixed.pop))
#' The plot shows a linear increase in the number of individuals of each species, when they are sorted
#' from lowest to huighest
species_richness3(mixed.pop)
# Again there are 301 species present in the population
simpson_diversity3(mixed.pop)
#' And now we see that there is a very small chance of random selecting two individuals from the same species
shannon_entropy3(mixed.pop) # = 5.516
berger_parker3(mixed.pop) # = 0.0066
#'
#'----------------------------------------------------------------------------------
#'
#' #### even.pop
#' even.pop contains a population which has 739 individuals of every species
plot(sort(even.pop))
#' As we would expect, the graph shows a straight horizontal line, since all species have the same number
#' of individuals.
species_richness3(even.pop)
# Same species richness
simpson_diversity3(even.pop)
# And an even smaller chance of selecting two individuals from the same species
shannon_entropy3(even.pop) # = 5.71
berger_parker3(even.pop) # = 0.003322
#'
#'----------------------------------------------------------------------------------
#'
#' #### rand.pop
#' rand.pop contains a random number of individuals of each species, drawn once from the BCI_2010 data
plot(sort(rand.pop))
#' This will change each time the population is generated
species_richness3(rand.pop)
# also 301
simpson_diversity3(rand.pop)
# Also small number
shannon_entropy3(rand.pop) # = 5.706
berger_parker3(rand.pop) # = 0.00361
#'
#'----------------------------------------------------------------------------------
#'
#' #### rand50.pop
#' rand50.pop contains a population with a random number of individuals of each species, drawn from the
#' BCI_2010 data 50 times.
plot(sort(rand50.pop))
#' Another random curve which will change each time the population is generated
species_richness3(rand50.pop)
# also 301
simpson_diversity3(rand50.pop)
# also a very small number
shannon_entropy3(rand50.pop) # = 4.019
berger_parker3(rand50.pop) # = 0.135
#'
#' ----------------------------------------------------------------------------------
#'
#' # Conclsuion
#' All of the values of species richness and simpson diversity (apart from the randomly generated
#' populations) are the same as in the first script using the first functions. So assuming the initial
#' functions worked correctly, we can see that our new functions are also working correctly!
#'
#'
universal_diversity3(tree.pop, 0) # = 301
universal_diversity3(tree.pop, 1) # = 55.66
universal_diversity3(tree.pop, 2) # = 22.25
universal_diversity3(tree.pop, Inf) # = 7.39

universal_diversity3(quadrat.pop, 0) # = 52
universal_diversity3(quadrat.pop, 1) # = 26
universal_diversity3(quadrat.pop, 2) # = 15
universal_diversity3(quadrat.pop, Inf) # = 5.8

universal_diversity3(quadrat10.pop, 0) # = 131
universal_diversity3(quadrat10.pop, 1) # = 39.84
universal_diversity3(quadrat10.pop, 2) # = 17.76
universal_diversity3(quadrat10.pop, Inf) # = 6.05

universal_diversity3(one.pop, 0) # = 1
universal_diversity3(one.pop, 1) # = 1
universal_diversity3(one.pop, 2) # = 1
universal_diversity3(one.pop, Inf) # = 1

universal_diversity3(uneven.pop, 0) # = 301
universal_diversity3(uneven.pop, 1) # = 34.3
universal_diversity3(uneven.pop, 1.01) # = 32.9

universal_diversity2(uneven.pop, 1.01) # = 32.9

universal_diversity2(uneven.pop, 1) # = 1 - NOT CORRECT

universal_diversity3(uneven.pop, 2) # = 3.96
universal_diversity3(uneven.pop, Inf) # = 1.99

universal_diversity3(mixed.pop, 0) # = 301
universal_diversity3(mixed.pop, 1) # = 248.5
universal_diversity3(mixed.pop, 0.99) # = 248.85
universal_diversity3(mixed.pop, 1.01) # = 248.23
universal_diversity3(mixed.pop, 2) # = 226.1
universal_diversity3(mixed.pop, Inf) # = 151
universal_diversity3(mixed.pop, 20) # = 170.6
universal_diversity3(mixed.pop, 30) # = 165.69

universal_diversity3(even.pop, 0) # = 301
universal_diversity3(even.pop, 1) # = 301
universal_diversity3(even.pop, 0.99) # = 301
universal_diversity3(even.pop, 2) # = 301
universal_diversity3(even.pop, Inf) # = 301

universal_diversity3(rand.pop, 0) # = 301
universal_diversity3(rand.pop, 1) # = 300.8
universal_diversity3(rand.pop, 2) # = 300.6
universal_diversity3(rand.pop, Inf) # = 272.3

universal_diversity3(rand50.pop, 0) # = 301
universal_diversity3(rand50.pop, 1) # = 291.6
universal_diversity3(rand50.pop, 2) # = 283.8
universal_diversity3(rand50.pop, Inf) # = 178.16