LargeScaleClass.Rmd

---
title: 'Large scale biodiversity patterns'
author: Jesús N. Pinto-Ledezma
output:
  pdf_document: default
  html_document: 
    theme: readable
    toc: yes
---

Today you will learn basic tools in R for visualizing species distributions, build geographical ranges, testing drivers of gradients of biodiversity under different approaches. 
You will need four datasets, that will be provided for you:
  1. Species occurrence data points - live.oaks.txt
  2. Species geograhical ranges - Furnarii_ranges_geo.shp
  3. Phylogenetic tree - furnariiMCC.nex
  4. Environmental predictors - bio1.bil and bio12.bil

# Set up your data and your working directory

Set up a working directory and put the data files in that directory. Tell R that this is the directory you will be using, and read in your data:

```{r, warnings = FALSE, message = FALSE, eval = FALSE}
setwd("path/for/your/directory")
```

Install and load the following packages

```{r, eval = FALSE}
require(maptools) # Tools for Handling Spatial Objects
require(rgdal) # Bindings for the 'Geospatial' Data Abstraction Library
require(raster) # Geographic Data Analysis and Modeling
require(sp) # Classes and Methods for Spatial Data
require(rangeBuilder) # not installed on this machine
library(spdep) # Spatial Dependence: Weighting Schemes, Statistics
library(ncf) # Spatial Covariance Functions 
require(geiger) # Analysis of Evolutionary Diversification
require(dismo) # Species Distribution Modeling
library(letsR) # Data Handling and Analysis in Macroecology
library(rworldmap) # Mapping Global Data
require(spatialreg) # Spatial Regression Analysis 
require(picante) # Integrating Phylogenies and Ecology
require(ape) # Analyses of Phylogenetics and Evolution
```
# From point occurrences to range maps

Load species occurrences data points. We will use occurrences from Live oaks, that were obtained from iDigBio between 20 and 24 July 2018 by Jeannine Cavender-Bares. Notice that these occurrence data points were visually examined and any localities that were outside the known range of the species, or in unrealistic locations (e.g., water bodies, crop fields) were discarded. 

```{r, warnings = FALSE, message = FALSE, eval = FALSE}
oaks <- read.table("Data/OCC/live.oaks.txt", header = TRUE)
```

```{r, warnings = FALSE, message = FALSE, eval = FALSE}
head(oaks)

tail(oaks)
```

Plot the points (x = Logitude, y = Latitude) and a world map, for reference. We need to load a data object from the {rworldmap} package

```{r, eval = FALSE, fig.keep='all'}
{plot(oaks[c(2:3)], col = "blue", pch = 19)
plot(countriesCoarse, add = TRUE)}
```

Lets try with a single species, in this case, Quercus virginiana
```{r, eval = FALSE}
unique(oaks$Species)

que_vir <- subset(oaks, oaks$Species == "Quercus_virginiana")
```

```{r, eval = FALSE}
{plot(que_vir$Longitude, que_vir$Latitude, pch = 15)
plot(countriesCoarse, add = TRUE)}
```

Cool, right? 

# Data checking and cleaning

Check if there are any duplicated points. We will check if exist any duplicate occurrence data point, if so, then remove all duplicates.

```{r, eval = FALSE}
oaks_dups <- duplicated(oaks[, c(2:3)])
### NOTE: the function "duplicated" returns the results of a logical test (e.g. TRUE or FALSE)
# How many are duplicates?
length(which(oaks_dups == TRUE))
# How many are NOT duplicates?
length(which(oaks_dups == FALSE))
# Keep only those lines that are not duplicates
oaks_dups_row <- which(oaks_dups == TRUE)
# What's the size? That is, how many points are duplicates
length(oaks_dups_row)
# Create another object withoyt the duplicate records
oaks_nodups <- oaks[-oaks_dups_row,]
# What are the dimensions of the new object?
dim(oaks_nodups)
# Take a look at the first rows of data
head(oaks_nodups)
```

What was the result? Are there duplicated occurrences?

Lets plot the results!

```{r, eval = FALSE}
{plot(oaks$Longitude, oaks$Latitude, pch = 19, col = "red", cex = 2)
points(oaks_nodups$Longitude, oaks_nodups$Latitude, pch = 16, col = "black")}
```

# Range maps from point data

This part can be used to create “simple” range maps based on geometry (e.g. minimum convex polygons, etc.), without considering environmental variables (no ENMs or SDMs).

## Convex hull (minimum convex polygon)

This model draws a convex hull around all ’presence’ points.

```{r, eval = FALSE}
# create a polygon around the species' records
oaks_hull <- convHull(oaks_nodups[, c(2:3)])
```

```{r, eval = FALSE}
# Plot the created polygon
{plot(oaks_hull)
points(oaks_nodups$Longitude, oaks_nodups$Latitude, pch = 16, col = "black")}
```

Ok, now let's try with a single species,

```{r, eval = FALSE}
que_vir <- subset(oaks_nodups, oaks_nodups$Species == "Quercus_virginiana")
```


```{r, eval = FALSE}
que_vir_hull <- convHull(que_vir[, c(2:3)])
```

```{r, eval = FALSE}
{plot(que_vir_hull)
points(que_vir$Longitude, que_vir$Latitude, pch = 16, col = "black")}
```

Now lets plot all ranges.

```{r, eval = FALSE}
# Plot all live oaks
{plot(oaks_hull)
points(oaks_nodups$Longitude, oaks_nodups$Latitude, pch = 16, col = "black")

# Plot only Quecus virginiana
plot(que_vir_hull, add = TRUE)
points(que_vir$Longitude, que_vir$Latitude, pch = 16, col = "red")

# Add world maps
plot(countriesCoarse, add = TRUE)}
```

Using simple convex hull seems not to be a good tool for live oaks, now lets try another approach. Now, we will use the dynamic alpha hull from the package {rangeBuilder}.

## Dynamic Alpha hull

```{r, eval = FALSE}
que_vir_alphahull <- getDynamicAlphaHull(que_vir, fraction = 0.95, 
                                         coordHeaders = c("Longitude", "Latitude"),
                                         clipToCoast = 'no')[[1]]
```

Now plot and see the differences.

```{r, eval = FALSE}
{plot(que_vir_alphahull, lwd = 2, col = "red") 
plot(que_vir_hull, add = TRUE, lwd = 2, lty = 2)
points(que_vir$Longitude, que_vir$Latitude, pch = 16, col = "green")
plot(countriesCoarse, add = TRUE, lwd = 2)}
```

Is there any difference? Please, explain the difference!

Until here we explored how to plot, clean and build species geographical ranges using occurrences. Now will use species geographical ranges of the largest continental endemic radiation (Furnariides) to explore the geographical gradients of species diversity.

# Diversity gradients

## Prepare data and mapping

The geographical ranges correspond to the Infraorder Furnariides (Aves). This data is available thorough BirdLife International (http://datazone.birdlife.org/species/requestdis) and you can use any other group available on IUCN or BIEN (for plants in the Americas). In any case, you first need to download the polygons in shapefile format.

To load the Furnariides geographical ranges we will use the function **{readOGR}** from the package **{rgdal}**. There are other functions

```{r, eval = FALSE}
franges <- readOGR(dsn = "Data/Franges", layer = "Furnarii_ranges_geo")
```

Now explore the data inside the ranges. Notice that to access to the information, we will use **@** instead of **$**. 

```{r, eval = FALSE}
head(franges@data)
```

## Raster of species richness

First create an empty raster for the Neotropics using the extent of the furnariides ranges under a spatial resolution of 1º long-lat.

```{r, eval = FALSE}
neo_ras <- raster()
# Set the raster "extent" 
extent(neo_ras) <- extent(franges)
res(neo_ras) <- 1

neo_ras

values(neo_ras) <- 0
```

Now using the empty raster we will **rasterize** the species identities in each cell or pixel. The resulting raster will be the Furnariides' species richness across the Neotropics.

```{r, eval = FALSE}
f_sr_raster <- rasterize(franges, neo_ras, field = "SCINAME", 
                         fun = function(x,...){length(unique(na.omit(x)))})
```

Plot the raster.

```{r, eval = FALSE}
{plot(f_sr_raster)
plot(countriesCoarse, add = T)}
```

Let's try changing the colors.

```{r, eval = FALSE}
#change the color scale
colfuncYellows <- colorRampPalette(c("#d7191c", "#fdae61", "#ffffbf", "#abd9e9", 
                                     "#2c7bb6"))
````

```{r, eval = FALSE}
{plot(f_sr_raster, col = rev(colfuncYellows(100)), axes = FALSE, box = FALSE, 
     zlim = c(minValue(f_sr_raster), maxValue(f_sr_raster)), 
     xlab = "Furnariides richness", legend.width = 2)
plot(countriesCoarse, add = T)}
```

Awesome, right?. Now, please describe the observed pattern!

Let's try to rasterize other information from the polygon data set. We will use the information in the column **RD**, this data correspond to the numbers of nodes from the tips to the root of a phylogenetic tree or just **root distance**, thus, will use the RD to calculate the MRD metric **(mean root distance)** that measures the evolutionary derivedness of species within an assemblage (Kerr & Currie, 1999) and can be used to determine whether a local fauna is constituted primarily by early-diverged or by recently originated species (Hawkins et al., 2012). 

```{r, eval = FALSE}
head(franges@data)

f_MRD_raster <- rasterize(franges, neo_ras, field = "RD", fun = mean)
```

```{r, eval = FALSE}
{plot(f_MRD_raster)
plot(countriesCoarse, add = T)}
```

Let's try changing the colors.

```{r, eval = FALSE}
{plot(f_MRD_raster, col = rev(colfuncYellows(100)), axes = FALSE, box = FALSE, 
     zlim = c(minValue(f_MRD_raster), maxValue(f_MRD_raster)), 
     xlab = "Furnariides mean root distance", legend.width = 2)
plot(countriesCoarse, add = T)}
```

Let’s plot the two raster.

```{r, eval = FALSE}
{par(mfrow = c(1, 2))
plot(f_sr_raster, col = rev(colfuncYellows(100)), axes = FALSE, box = FALSE, 
     zlim = c(minValue(f_sr_raster), maxValue(f_sr_raster)), 
     xlab = "Furnariides richness", legend.width = 2)

plot(f_MRD_raster, col = rev(colfuncYellows(100)), axes = FALSE, box = FALSE, 
     zlim = c(minValue(f_MRD_raster), maxValue(f_MRD_raster)), 
     xlab = "Furnariides mean root distance", legend.width = 2)}
```

See if there is a relationship...

```{r, eval = FALSE}
obj <- lm(values(f_MRD_raster) ~ values(f_sr_raster))
cor.test(values(f_sr_raster), values(f_MRD_raster))
```

```{r, eval = FALSE}
{plot(values(f_MRD_raster) ~ values(f_sr_raster), xlab = "SR", ylab = "MRD")
abline(obj, col = "red", lwd = 2)}
```

Hmmm. What happened in here? So, please answer the next questions.

*From the mean root distance map, it is possible to explain the Furnariides diversity gradient? If so, please explain from an evolutionary perspective* 

Let's be more quantitative in describing that relationship with a linear model.

## Scale dependency

Now we will explore one of the oldest problems in ecology and evolution, the **scale dependency** in the data. So, to explore this scale dependence we will rasterize the Furnariides ranges, but using different spatial resolutions from 2º to 4º degrees of long-lat.

```{r, eval = FALSE}
# 2º degrees
neo_ras_2dg <- raster()
# Set the raster "extent" 
extent(neo_ras_2dg) <- extent(franges)
res(neo_ras_2dg) <- 2
neo_ras_2dg
values(neo_ras_2dg) <- 0

# 4º degrees
neo_ras_4dg <- raster()
# Set the raster "extent" 
extent(neo_ras_4dg) <- extent(franges)
res(neo_ras_4dg) <- 4
neo_ras_4dg
values(neo_ras_4dg) <- 0
```

```{r, eval = FALSE}
# Furnariides at 2º of long-lat
f_sr_2dg_raster <- rasterize(franges, neo_ras_2dg, field = "SCINAME", fun = function(x,...){length(unique(na.omit(x)))})

# Furnariides at 4º of long-lat
f_sr_4dg_raster <- rasterize(franges, neo_ras_4dg, field = "SCINAME", fun = function(x,...){length(unique(na.omit(x)))})
```

Plot the three maps.

```{r, eval = FALSE}
{par(mfrow = c (1, 3))
plot(f_sr_raster, main = "Furnariides richness 1dg")
plot(countriesCoarse, add = T)

plot(f_sr_2dg_raster, main = "Furnariides richness 2dg")
plot(countriesCoarse, add = T)

plot(f_sr_4dg_raster, main = "Furnariides richness 4dg")
plot(countriesCoarse, add = T)}
```

So, is there a scale effect?

*Explain the differences between the three maps*

## Correlative relationships

Load the environmental variables that correspond to bio1 (Annual Mean Temperature) and bio12 (Annual Precipitation). These data correspond to two variables out of 19 from WorldClim (http://www.worldclim.org/current). We will use these two variables just for educational purposes, rather to make a complete evaluation of the species-environmental relationships.

```{r, eval = FALSE}
bio1 <- raster("Data/Envi/bio1.bil")
bio1
bio12 <- raster("Data/Envi/bio12.bil")
bio12
```

Plot the environmental variables

```{r, eval = FALSE}
{par(mfrow = c(2, 1))
plot(bio1)
plot(bio12)}
```

Ok, the bio1 and bio12 layers are at global scale, so now will need to crop to the extent of the Neotropics.

```{r, eval = FALSE}
bio1_neo <- crop(bio1, extent(franges))
bio12_neo <- crop(bio12, extent(franges))
```

```{r, eval = FALSE}
{par(mfrow = c(1, 2))
plot(bio1_neo, main = "Annual Mean Temperature")
plot(bio12_neo, main = "Annual Precipitation")}
```

Much better!

Now we will obtain the coordinates from the Furnariides diversity raster. These coordinates then will be used to extract the information from the bio1 and bio12 climatic layers.

```{r, eval = FALSE}
f_ras_coords <- xyFromCell(f_sr_raster, 1:length(values(f_sr_raster)))
```

Obtain the values from bio1 and bio12 for each cell or pixel using the coordinates.

```{r, eval = FALSE}
f_ras_bio1 <- extract(bio1_neo, f_ras_coords)
f_ras_bio12 <- extract(bio12_neo, f_ras_coords)

f_ras_rich <- values(f_sr_raster)
f_ras_mrd <- values(f_MRD_raster)
fdata <- na.omit(data.frame(f_ras_coords, f_ras_rich, f_ras_mrd, f_ras_bio1, f_ras_bio12))

#f_ras_rich_noNA <- ifelse(is.na(f_ras_rich), 0, f_ras_rich)
#f_ras_bio1_noNA <- ifelse(is.na(f_ras_bio1), 0, f_ras_bio1)
#f_ras_bio12_noNA <- ifelse(is.na(f_ras_bio12), 0, f_ras_bio12)

#fdata_noNA <- data.frame(f_ras_coords, f_ras_rich_noNA, f_ras_bio1_noNA,
#                         f_ras_bio12_noNA)
```

Now make a simple correlation between the Furnariides richness and bio1 and bio12.

```{r, eval = FALSE}
cor.test(fdata$f_ras_rich, fdata$f_ras_bio1)
```

```{r, eval = FALSE}
cor.test(fdata$f_ras_rich, fdata$f_ras_bio12)
```

*Which environmental variable is more related with Furnariides richness?*

*Please explain the relationship from an ecological perspective*

```{r, eval = FALSE}
{par(mfrow = c(1, 2))
plot(fdata$f_ras_bio1, fdata$f_ras_rich, xlab = "Bio 1", ylab = "Richness")

plot(fdata$f_ras_bio12, fdata$f_ras_rich, xlab = "Bio 12", ylab = "Richness")}
```

## Considering spatial autocorrelation

We now that a correlation is not causation, so to explore the true relationships between a variable and N covariables we need to build a model or fit a model. To explore this relationships we will first explore a simple Ordinary Least Square regression or OLS.

```{r, eval = FALSE}
fols <- lm(f_ras_rich ~ f_ras_bio1 + f_ras_bio12, data = fdata)
summary(fols)
```

Let's complicate our model a little bit...

```{r, eval = FALSE}
fols2 <- lm(f_ras_rich ~ f_ras_bio1 + f_ras_bio12 + f_ras_mrd, data = fdata)
summary(fols2)
```

*What is telling us this OLS?*

Now, explore the spatial autocorrelation in the Furnariides richness gradient.

```{r, eval = FALSE}
autocor_SR <- ncf::correlog(fdata$x, fdata$y, z = fdata$f_ras_rich, na.rm = T, 
                         increment = 1, resamp = 1)
```

Let's use an correlogram to explore the spatial autocorrelation.

```{r, eval = FALSE}
{plot(autocor_SR$correlation[1:50], type = "b", pch = 1, cex = 1.2, lwd = 1.5,
     ylim = c(-1, 1), xlab = "Distance class", ylab = "Moran's I", cex.lab = 1.2, 
     cex.axis = 1.2)
abline(h = 0)}
```

*Is there a spatial autocorrelation in the data?*

What about the residuals? Now, we will explore the spatial autocorrelation in the residuals.

```{r, eval = FALSE}
coords <- fdata[1:2]
coords <- as.matrix(coords)
```

Build a neighbourhood contiguity by distance

```{r, eval = FALSE}
nb1.5 <- dnearneigh(coords, 0, 1.5)
```

Using the neighbourhood contiguity build a spatial weights for neighbours lists.

```{r, eval = FALSE}
nb1.5.w <- nb2listw(nb1.5, glist = NULL, style = "W", zero.policy = TRUE)
```

Extract the residuals from the OLS model

```{r, eval = FALSE}
residuals_ols <- residuals(fols2)
plot(residuals_ols)
```

Calculate a univariate spatial correlogram.

```{r, eval = FALSE}
autocor_ols_res <- correlog(fdata$x, fdata$y, z = residuals(fols), 
                        increment = 1, resamp = 1)
```

plot the autocorrelagram for the residuals

```{r, eval = FALSE}
{plot(autocor_ols_res$correlation[1:50], type = "b", pch = 1, cex = 1.2, lwd = 1.5,
     ylim = c(-0.5, 1), xlab = "distance", ylab = "Moran's I", cex.lab = 1.5, 
     cex.axis = 1.2)
abline(h = 0)
title(main = "OLS residuals", cex = 1.5)}
```

Ohhh, seems that the residuals have a strong spatial autocorrelation, that is a problem because if we found autocorrelation in the residuals much of the explanation that we obtained can be biased.

Let's inspect two two autocorrelograms.

```{r, eval = FALSE}
{par(mfrow = c(2, 1))

plot(autocor_SR$correlation[1:50], type = "b", pch = 1, cex = 1.2, lwd = 1.5,
     ylim = c(-1, 1), xlab = "Distance class", ylab = "Moran's I", cex.lab = 1.2, 
     cex.axis = 1.2)
abline(h = 0)
title(main = "OLS model", cex = 1.5)

plot(autocor_ols_res$correlation[1:50], type = "b", pch = 1, cex = 1.2, lwd = 1.5,
     ylim = c(-0.5, 1), xlab = "Distance class", ylab = "Moran's I", cex.lab = 1.5, 
     cex.axis = 1.2)
abline(h = 0)
title(main = "OLS residuals", cex = 1.5)}
```

Hmmm, seems that there is a strong spatial autocorrelation, thus any conclusion based on the use of this OLS model can be biased.

To try to solve this important issue, we will use **spatial simultaneous autoregressive error model estimation (Aka SAR model)**, this model account for spatial autocorrelation by adding an extra term (autoregressive) in the form of a spatial-weight matrix that specifies the neighborhood of each cell or pixel and the relative weight of each neighbor.

Let's fit the SAR model.

```{r, eval = FALSE}
sar_nb1.5.w <- errorsarlm(fols2, listw = nb1.5.w, data = fdata, quiet = FALSE, 
                          zero.policy = TRUE, na.action = na.exclude)
summary(sar_nb1.5.w)
residuals_sar_nb1.5.w <- residuals(sar_nb1.5.w)
```

Now estimate the spatial autocorrelation of the SAR model.

```{r, eval = FALSE}
autocor_sar_nb1.5.w <- correlog(fdata$x, fdata$y, z = residuals(sar_nb1.5.w), 
                                na.rm = T, increment = 1, resamp = 1)
```

Plot the autocorrelogram under the SAR model.

```{r, eval = FALSE}
{plot(autocor_sar_nb1.5.w$correlation[1:50], type = "b", pch = 4, cex = 1.2, lwd = 1.5,
     ylim = c(-0.5, 1), xlab = "distance", ylab = "Moran's I", cex.lab = 1.5, 
     cex.axis = 1.2)
abline(h = 0)
title(main = "Correlogram SARerr", cex = 1.5)}
```

Ohhh, where is the autocorrelation in the residuals? Now compare the two autocorrelograms.

```{r, eval = FALSE}
{par(mfrow = c(2, 1))
plot(autocor_ols_res$correlation[1:50], type = "b", pch = 1, cex = 1.2, lwd = 1.5,
     ylim = c(-0.5, 1), xlab = "distance", ylab = "Moran's I", cex.lab = 1.5, 
     cex.axis = 1.2)
abline(h = 0)
title(main = "OLS residuals", cex = 1.5)

plot(autocor_sar_nb1.5.w$correlation[1:50], type = "b", pch = 4, cex = 1.2, lwd = 1.5,
     ylim = c(-0.5, 1), xlab = "distance", ylab = "Moran's I", cex.lab = 1.5, 
     cex.axis = 1.2)
abline(h = 0)
title(main = "Correlogram SARerr", cex = 1.5)}
```

Ok, now we know that the SAR model can solve the problem in the spatial autocorrelation in the residuals, let's try to make some inferences.

```{r, eval = FALSE}
summary(sar_nb1.5.w)
```

```{r, eval = FALSE}
summary(fols2)
```

Looking to the summary of the SAR and OLS models, explain the differences in the coefficients between both models.

Now let's compare the prediction of both models. To calculate a R2 to the SAR model, we will use the function **SARr2()** from GitHub.

```{r, eval = FALSE}
source("https://raw.githubusercontent.com/jesusNPL/BetaDivNA/master/SARr2.R")
```

```{r, eval = FALSE}
SARr2(Lfull = sar_nb1.5.w$LL, Lnull = sar_nb1.5.w$logLik_lm.model, N = nrow(fdata))
```

Comparing the two models (OLS and SAR), please answer the following questions:

*1. Which model have the best explanation?*

*2. What can we conclude from these results?*

# Part 2
The main goal of the second part is to present basic understanding about measuring phylogenetic diversity with in communities or best known as the analysis of community phylogenetics. The community phylogenetics integrates ecological and evolutionary concepts and explores the mecanisms (e.g., biotic interactions or environemntal filters) governing the assembly of ecological communities. 

There are different sources of information and web pages with a lot of information about this field. The most common and useful are the web pages of the books: “Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology” (Garamszegi, 2014) (http://www.mpcm-evolution.org/) and “Phylogenies in Ecology” (Cadotte and Davies, 2016) (https://www.utsc.utoronto.ca/~mcadotte/page-3/).

# Prepare data
## Load phylogentic and community data

```{r, results = "hide", warnings = FALSE, message = FALSE}
ftree <- ape::read.nexus("Data/Phylo/furnariiMCC.nex")
```

Lets inspect the tree.

```{r, results = "hide", eval = FALSE}
length(ftree$tip.label)
head(ftree$tip.label)
```

Now load the community data. Well we don't have a community data, but we can build one ;p...To do so, using the furnariides ranges from the part 1, we will construct a presence-absence matrix (PAM) in which species are columns and rows correspond to grid cells or biological assemblages.

```{r, results = "hide", eval = FALSE}
fpam <- lets.presab(franges, xmn = -110.0817, xmx = -34.78897, ymn = -55.98222, ymx = 29.0965, remove.cells = TRUE, remove.sp = TRUE, resol = 1, count = FALSE)
```

Inspect the data

```{r, results = "hide", eval = FALSE}
class(fpam)
fpam$Presence_and_Absence_Matrix[1:10, 1:10]
dim(fpam$Presence_and_Absence_Matrix)
```

## Prepare data - match phylogenetic and assemblage data
Now we will match both information that we need for further analyses.

```{r, results = "hide", eval = FALSE}
fcom <- fpam$Presence_and_Absence_Matrix
dim(fcom)
sppNames <- fpam$Species_name
sppNames <- gsub(" ", "_", sppNames)

colnames(fcom) <- c("long", "lat", sppNames)
fcom[1:10, 1:10]
dim(fcom)
```

Now match both phylogentic and assemblages information

```{r, results = "hide", eval = FALSE}
matched <- match.phylo.comm(phy = ftree, comm = fcom[, 3:654])
```

Explore the resulting match dataset
```{r, results = "hide", eval = FALSE}
matched$phy

dim(matched$comm)
```

Ok, we now have all the information we need for exploring some phylogenetic diversity metrics.

# Phylogenetic diversity metrics
## Phylogenetic diversity - Faith's PD
```{r, results = "hide", warnings = FALSE, message = FALSE, eval = FALSE}
Fpd <- pd(matched$comm, matched$phy, include.root = FALSE) # Faith's PD
head(Fpd)
```
Explore the relationship between phylogenetic diversity and species richness within each assemblage.
```{r, results = "hide", warnings = FALSE, message = FALSE, eval = FALSE}
cor.test(Fpd$SR, Fpd$PD)

plot(Fpd$SR, Fpd$PD)
```

## Mean pairwise distance

```{r, results = "hide", warnings = FALSE, message = FALSE, eval = FALSE}
Fmpd <- mpd(matched$comm, cophenetic(matched$phy)) # Faith's PD
head(Fmpd, 10)
```

## Mean nearest-pairwise distance

```{r, results = "hide", warnings = FALSE, message = FALSE, eval = FALSE}
Fmntd <- mntd(matched$comm, cophenetic(matched$phy))
head(Fmntd, 10)
```

# Community diversity metrics

The analyses of community phylogenetic started making inferences about the mechanisms structuring the local communities through the evaluation of phylogenetic arrangements in local communities (see Cavender-Bares et al. 2009 for an initial criticism). However, new methods are now available, such that more complex balance between ecological and historical processes at local and regional scales can be incorporated into the analyses (Pigot and Etienne 2015, Pinto-Ledezma et al. 2019).

Now, lets calulate some of the most common metrics. But for practical porpuses we will just run the null models 99 times. In real life you should run the null model at least with 999 repetitions...

MPD - Mean pairwise distance separating taxa in a community

```{r, results = "hide", warnings = FALSE, message = FALSE, eval = FALSE}
# SESMPD
Fsesmpd <- ses.mpd(matched$comm, cophenetic(matched$phy), runs = 99)
head(Fsesmpd, 15)
```

MNTD - Mean nearest taxon distance for taxa in a community
 
```{r, results = "hide", warnings = FALSE, message = FALSE, eval = FALSE}
# SESMNTD
Fsesmntd <- ses.mntd(matched$comm, cophenetic(matched$phy), runs = 99)
head(Fsesmntd, 15)
```

The end!!! for now...

# References

Cadotte, M. W. and Davies, T. J. (2016). Phylogenies in Ecology: A Guide to Concepts and Methods. Princeton: Princeton University Press.

Cavender-Bares, J., Kozak, K. H., Fine, P. V. A. and Kembel, S. W. (2009). The merging of community ecology and phylogenetic biology. Ecology letters 12, 693–715.

Garamszegi, L. Z. (2014). Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology. (ed. Garamszegi, L. Z.) Berlin: Springer-Verlag.
1. Hawkins, B.A., McCain, Ch.M., Davies, T.J., Buckley, L.B., Anacker, B.L., Cornell, H.V., Damschen, E.I., Grytness, J.A., Harrison, S., Holt, R.D., Kraft, N.J.B. & Stephens, P.R. (2012) Different evolutionary histories underlie congruent species richness gradients on birds and mammals. Journal of Biogeography, 39, 825–841.

Kerr, J.T. & Currie, D.J. (1999) The relative importance of evolutionary and environmental controls on broad-scale patterns of species richness in North America. Ecoscience, 6, 329–337.

Pigot, A. L., & Etienne, R. S. (2015). A new dynamic null model for phylogenetic community structure. Ecology Letters, 18(2), 153–163.

Pinto-Ledezma, J. N., Simon, L. M., Diniz-Filho, J. A. F., & Villalobos, F. (2017). The geographical diversification of Furnariides: the role of forest versus open habitats in driving species richness gradients. Journal of Biogeography, 44(8), 1683–1693.

Pinto-Ledezma, J. N., Jahn, A. E., Cueto, V. R., Diniz-Filho, J. A. F., & Villalobos, F. (2019). Drivers of Phylogenetic Assemblage Structure of the Furnariides, a Widespread Clade of Lowland Neotropical Birds. The American Naturalist, 193(2), E41–E56.