ProgInRBCIfunc: species_richness_1

#' ---
#' title: 02 Species Richness and Simpson Diversity
#' author: Edie Bishop
#' ---

#' # Introduction
#' In this package we will use the ProgInRBCIfunc package to call on two functions, species_richness2 and
#' simpson_diversity2, which will calculate the species richness and the simpson diversity of a population.
#' species_richness2 and simpson_diversity2 work slightly differently to species_richness_1 and
#' simpson_diversity_1 in that they themselves call on a third function, universal_diversity2, which
#' calculates a D (value of diversity) based on the population and the value of q which is supplied to it.
#' Species richness has a q value of 0, and simpson diversity has a q value of 2, which are defined within
#' the species_richness2 and simpson_diversity2 functions. The general diversity function can be represented
#' by the equation below:
#'
#' $$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.*
#'
#' The universal diversity function calculates this equation, for a given value of q. The species_richness
#' and simpson_diversity functions the modify D appropriately, whereby species richness = $D_{0}$, and
#' Simpson Diversity = $1/D_{2}$.
#'
#' We will also call on the script "create_datasets.R" 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 the species richness and simpson doversity, and then calculate the richness and simpson
#' diversity values. We will then check that the fucntions created are calculating diversity measures
#' correctly by comparing answers using the rdiversity package.
#'
#' 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

#' First set the working directory to the Run_scripts folder within the package so we we van access the
#' create_datasets.R file, which generatres the populations.
devtools::wd(".", "Run_Scripts")
# The ProgInRBCIfunc function contains functions which calculate various indices of diversity
library(ProgInRBCIfunc)
# Also load the rdiversity package so that we can check that our functions are calculating correct values
# of Dq
library(rdiversity)
# 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_richness2(tree.pop) # As we know that there are 301 different species in the community, the species richness = 301 as expected

simpson_diversity2(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.

#'
#'------------------------------------------------------------------------------------
#'
#' #### 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_richness2(quadrat.pop) # There are less species within the randomly selected quadrat than the whole of the community. This valuewill change each time quadrat.pop is created.

simpson_diversity2(quadrat.pop) # This value will change each time quadrat.pop is created


#'
#'--------------------------------------------------------------------------------------
#'
#' #### 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_richness2(quadrat10.pop) # More species than quadrat.pop, which makes sense. This value will change each time quadrat50.pop is created.

simpson_diversity2(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.


#'
#'-------------------------------------------------------------------------------------
#'
#' #### 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_richness2(one.pop) # Vlaue is 1, as there is only one species

simpson_diversity2(one.pop) # Value is 1, since two randomly selected individuals are always going to be from the same species


#'
#'-------------------------------------------------------------------------------------
#'
#' #### 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_richness2(uneven.pop) # Species richness is again 301, as there are again 301 species present in the population

simpson_diversity2(uneven.pop) # There is around a 25% chance of selecting two individuals from the same species


#'
#'-------------------------------------------------------------------------------------
#'
#' #### 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_richness2(mixed.pop) # Again there are 301 species present in the population

simpson_diversity2(mixed.pop)
#' And now we see that there is a very small chance of random selecting two individuals from the same species

#'
#'----------------------------------------------------------------------------------
#'
#' #### 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_richness2(even.pop) # Same species richness

simpson_diversity2(even.pop) # And an even smaller chance of selecting two individuals from the same species


#'
#'----------------------------------------------------------------------------------
#'
#' #### 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_richness2(rand.pop) # also 301

simpson_diversity2(rand.pop) # Also small number


#'
#'----------------------------------------------------------------------------------
#'
#' #### 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_richness2(rand50.pop)# also 301
simpson_diversity2(rand50.pop)# also a very small number
universal_diversity2(rand50.pop, 2)
# = 281

#' I will now check that the same values are calculated using the rdiversity package
meta <- metacommunity(rand50.pop) # create an object which stores information which the rdiversity package uses to calculate diversity
res <- meta_gamma(meta, qs = c(0,2)) # Create a dataframe of the biodiversity of the whole metacommunitiy at values of q 1-5
res
#' The raw diversity values for q = 0 and q = 2 are the same when they are calculated by the rdiversity
#' package. So that means that my functions are working! YAY!!
#'
#' ----------------------------------------------------------------------------------
#'
#' # 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!