back to overview    

Preparation: Please open your RStudio project and download the new data (‘Phylogeny_ConservativeCon_Checklist.nex’, ‘palms_specsxsites_phylo.csv’ and ‘palmtree_pruned.nex’) for today from Stud-IP and copy them into your .Rproj folder ‘/data/’. You can use getwd() to locate your current working directory, which should be your project folder. Please install the following R-packages using install.packages():

 

If you want to visualize this tutorial in the viewer inside RStudio (to save space on your screen) run the following chunk of code:

#install.packages("rstudioapi") # install an R-package required for this step
dir <- tempfile()
dir.create(dir)
download.file("https://gift.uni-goettingen.de/mcmmb/index.html", destfile = file.path(dir, "index.html"))
download.file("https://gift.uni-goettingen.de/mcmmb/Day7.html", destfile = file.path(dir, "Day7.html"))
htmlFile <- file.path(dir, "Day7.html")
rstudioapi::viewer(htmlFile)

Now you can conveniently copy code from the viewer into your script.
 

Load R packages & island data set

library("dplyr") # data manipulation

library("sf") #  Geospatial library

library("TNRS") # Taxonomic standardisation of plant species names
# library("rWCVP") # Taxonomic standardisation of plant species names
# library("rWCVPdata") # According data for standardization
library("ape") # Analyses of phylogenetics and evolution
library("picante") # Phylogenies and ecology
library("pez") # Phylogenies and ecology

library("psych") # basic stats stuff
library("ggplot2") # for plotting
library("cowplot") # for multi-panel plots
library("GGally") # pairs plot
library("RColorBrewer") # colour palettes

library("bioregion") # bioregionalization package used to switch data formats

 

1. Load data into R

1.1 Load a phylogeny & look at structure of ‘multiPhylo’ object

Load multi-tree object into R
Data from Faurby, S., Eiserhardt, W.L., Baker, W.J. & Svenning, J.-C. (2016).
An all-evidence species-level supertree for the palms (Arecaceae).
Molecular Phylogenetics and Evolution, 100, 57-69.

palmtrees <- read.nexus("data/Phylogeny_ConservativeCon_Checklist.nex")

palmtrees
## 1000 phylogenetic trees
palmtrees[[1]]
## 
## Phylogenetic tree with 2539 tips and 2538 internal nodes.
## 
## Tip labels:
##   Nypa_fruticans, Sabal_mauritiiformis, Sabal_pumos, Sabal_domingensis, Sabal_minor, Sabal_palmetto, ...
## 
## Rooted; includes branch length(s).
palmtree <- palmtrees[[1]]
str(palmtree)
## List of 4
##  $ edge       : int [1:5076, 1:2] 2540 2541 2541 2542 2543 2544 2545 2546 2546 2547 ...
##  $ edge.length: num [1:5076] 8.89 95.72 9.44 9.78 12.25 ...
##  $ Nnode      : int 2538
##  $ tip.label  : chr [1:2539] "Nypa_fruticans" "Sabal_mauritiiformis" "Sabal_pumos" "Sabal_domingensis" ...
##  - attr(*, "class")= chr "phylo"
##  - attr(*, "order")= chr "cladewise"

Explanation of structure of a phylogeny in R

edge: a two-column matrix where each row represents a branch (or edge) of the tree. The nodes and the tips are symbolized with integers. The n tips are numbered from 1 to n, and the m (internal) nodes from n+1 to n+m (the root being n + 1). For each row, the first column contains the ancestral node indices (where the branch starts). The second column contains the descendant node or tip indices (where the branch ends).

edge.length (optional): A numeric vector giving the lengths of the branches given by edge.

tip.label: A vector of mode character giving the labels of the tips. The order of these labels corresponds to the integers 1 to n in edge.

Nnode: An integer value giving the number of nodes in the tree (m).

node.label (optional): A vector of mode character giving the labels of the nodes.

root.edge (optional): A numeric value giving the length of the branch at the root if it exists.

# Edges
head(palmtree$edge)
##      [,1] [,2]
## [1,] 2540 2541
## [2,] 2541    1
## [3,] 2541 2542
## [4,] 2542 2543
## [5,] 2543 2544
## [6,] 2544 2545
#Edge length
head(palmtree$edge.length)
## [1]  8.890634 95.723409  9.438997  9.780648 12.246395 24.072388
# Number of nodes (m) (n is the number of tips)
palmtree$Nnode
## [1] 2538
# Tip labels
head(palmtree$tip.label)
## [1] "Nypa_fruticans"       "Sabal_mauritiiformis" "Sabal_pumos"         
## [4] "Sabal_domingensis"    "Sabal_minor"          "Sabal_palmetto"
length(palmtree$tip.label)
## [1] 2539
is.rooted(palmtree) # does the phylogeny have a root? i.e. can the most basal ancestor be identified?
## [1] TRUE

 

We can also plot the phylogenetic tree.

plot(palmtree, type = "fan", cex = 0.3,
     edge.color = "gray70", tip.color = "#ef8a62",
     main = paste("Phylogenetic tree for", length(palmtree$tip.label), "palm species"))

 

1.2. Loading palm distribution data

Data from Kreft, H., Sommer, J.H. & Barthlott, W. (2006).The significance of geographic range size for spatial diversity
patterns in Neotropical palms. Ecography, 29, 21-30.

species <- read.csv("data/palms_species_per_gridcell.csv",
                    sep = ",", stringsAsFactors = FALSE)
head(species)
##   grid_id spp_Id    GENUS EPITHET
## 1   13735    550 Ammandra natalia
## 2   13736    550 Ammandra natalia
## 3   13737    550 Ammandra natalia
## 4   13738    550 Ammandra natalia
## 5   13910    550 Ammandra natalia
## 6   13911    550 Ammandra natalia
length(unique(species$grid_id)) # number of grid cells
## [1] 6638

 

Add a column for the full species name with “” as separator: ”species”. The ”” separator is necessary because we want the species names to be coded the same as the names in the tip.labels of the phylogenetic tree.

species$species_name  <- paste(species$GENUS, species$EPITHET, sep = "_")
unique(species$species_name)[1:5]
## [1] "Ammandra_natalia"        "Ammandra_decasperma"    
## [3] "Ammandra_dasyneura"      "Phytelephas_tumacana"   
## [5] "Phytelephas_tenuicaulis"

 

How many species in the palm distribution dataset are missing from the phylogeny?
 

#Species in the palm dataset
length(unique(species$species_name)) # 550 species
## [1] 550
table(unique(species$species_name) %in% palmtree$tip.label) # need for unique species
## 
## FALSE  TRUE 
##    80   470
# Alternative ways to count the number of species NOT in the phylogenetic tree
# length(unique(species$species[which(!species$species %in% palmtree$tip.label)]))
# sum(!unique(species$species) %in% palmtree$tip.label)

Which species are missing?

# Species that are NOT in the phylogeny
unique(species$species_name)[!(unique(species$species_name) %in% palmtree$tip.label)]
##  [1] "Ammandra_natalia"             "Ammandra_dasyneura"          
##  [3] "Geonoma_weberbaueri"          "Geonoma_rubescens"           
##  [5] "Geonoma_polyandra"            "Geonoma_poeppiginana"        
##  [7] "Geonoma_paraguensis"          "Geonoma_myriantha"           
##  [9] "Geonoma_longevaginata"        "Geonoma_longepedunculata"    
## [11] "Geonoma_linearis"             "Geonoma_jussieuana"          
## [13] "Geonoma_gastoniana"           "Geonoma_gamiova"             
## [15] "Geonoma_densa"                "Geonoma_brevispatha"         
## [17] "Geonoma_arundinacea"          "Geonoma_appuniana"           
## [19] "Calyptrogyne_kunaria"         "Calyptrogyne_condesata"      
## [21] "Astrocaryum_gynacanthum"      "Desmoncus_phoenicocarpus"    
## [23] "Desmoncus_cirrhiferus"        "Bactris_trailiana"           
## [25] "Bactris_macana"               "Bactris_horrispatha"         
## [27] "Bactris_hondurensis"          "Bactris_grayumi"             
## [29] "Bactris_corrossilla"          "Aiphanes_leidostachys"       
## [31] "Aiphanes_aculeata"            "Gastrococcus_crispa"         
## [33] "Attalea_septuagenta"          "Attalea_luetzelbergii"       
## [35] "Polyandrococcus_caudescens"   "Syagrus_shmithii"            
## [37] "Syagrus_leptospatha"          "Hypospathe_altissima"        
## [39] "Hypospathe_macrorachis"       "Hypospathe_elegans"          
## [41] "Prestoea_roseospadix"         "Euterpe_longebracteata"      
## [43] "Leopoldinia_major"            "Wettinia_angusta"            
## [45] "Wettinia_aquatorialis"        "Socratea_hecantonandra"      
## [47] "Iriartella_stenocarpua"       "Dictyocarpum_ptarianum"      
## [49] "Dictyocarpum_lamarckianum"    "Dictyocarpum_fuscum"         
## [51] "Chamaedorea_vocanensis"       "Chamaedorea_tuerkheimii"     
## [53] "Chamaedorea_sullivaniorum"    "Chamaedorea_selvae"          
## [55] "Chamaedorea_orephila"         "Chamaedorea_keeleriorum"     
## [57] "Chamaedorea_hopperiana"       "Chamaedorea_ernesti-angustii"
## [59] "Ceroxylon_werberbaueri"       "Sabal_miamiensis"            
## [61] "Sabal_guatemalensis"          "Sabal_gretheriae"            
## [63] "Brahea_nitida"                "Acoelorraphe_wrightii"       
## [65] "Coplothrinax_wrightii"        "Coplothrinax_cookii"         
## [67] "Coccothrinax_hiorami"         "Thrinax_rivularis"           
## [69] "Thrinax_morrisii"             "Thrinax_ekmaniana"           
## [71] "Thrinax_compacta"             "Crysophila_williamsii"       
## [73] "Crysophila_warscewiczii"      "Crysophila_stauracantha"     
## [75] "Crysophila_nana"              "Crysophila_macrocarpa"       
## [77] "Crysophila_kalbreyeri"        "Crysophila_guagara"          
## [79] "Crysophila_grayumii"          "Crysophila_cookii"
#  Alternative way with dplyr (finds all rows in x that aren't in y)
setdiff(species$species_name, palmtree$tip.label)
##  [1] "Ammandra_natalia"             "Ammandra_dasyneura"          
##  [3] "Geonoma_weberbaueri"          "Geonoma_rubescens"           
##  [5] "Geonoma_polyandra"            "Geonoma_poeppiginana"        
##  [7] "Geonoma_paraguensis"          "Geonoma_myriantha"           
##  [9] "Geonoma_longevaginata"        "Geonoma_longepedunculata"    
## [11] "Geonoma_linearis"             "Geonoma_jussieuana"          
## [13] "Geonoma_gastoniana"           "Geonoma_gamiova"             
## [15] "Geonoma_densa"                "Geonoma_brevispatha"         
## [17] "Geonoma_arundinacea"          "Geonoma_appuniana"           
## [19] "Calyptrogyne_kunaria"         "Calyptrogyne_condesata"      
## [21] "Astrocaryum_gynacanthum"      "Desmoncus_phoenicocarpus"    
## [23] "Desmoncus_cirrhiferus"        "Bactris_trailiana"           
## [25] "Bactris_macana"               "Bactris_horrispatha"         
## [27] "Bactris_hondurensis"          "Bactris_grayumi"             
## [29] "Bactris_corrossilla"          "Aiphanes_leidostachys"       
## [31] "Aiphanes_aculeata"            "Gastrococcus_crispa"         
## [33] "Attalea_septuagenta"          "Attalea_luetzelbergii"       
## [35] "Polyandrococcus_caudescens"   "Syagrus_shmithii"            
## [37] "Syagrus_leptospatha"          "Hypospathe_altissima"        
## [39] "Hypospathe_macrorachis"       "Hypospathe_elegans"          
## [41] "Prestoea_roseospadix"         "Euterpe_longebracteata"      
## [43] "Leopoldinia_major"            "Wettinia_angusta"            
## [45] "Wettinia_aquatorialis"        "Socratea_hecantonandra"      
## [47] "Iriartella_stenocarpua"       "Dictyocarpum_ptarianum"      
## [49] "Dictyocarpum_lamarckianum"    "Dictyocarpum_fuscum"         
## [51] "Chamaedorea_vocanensis"       "Chamaedorea_tuerkheimii"     
## [53] "Chamaedorea_sullivaniorum"    "Chamaedorea_selvae"          
## [55] "Chamaedorea_orephila"         "Chamaedorea_keeleriorum"     
## [57] "Chamaedorea_hopperiana"       "Chamaedorea_ernesti-angustii"
## [59] "Ceroxylon_werberbaueri"       "Sabal_miamiensis"            
## [61] "Sabal_guatemalensis"          "Sabal_gretheriae"            
## [63] "Brahea_nitida"                "Acoelorraphe_wrightii"       
## [65] "Coplothrinax_wrightii"        "Coplothrinax_cookii"         
## [67] "Coccothrinax_hiorami"         "Thrinax_rivularis"           
## [69] "Thrinax_morrisii"             "Thrinax_ekmaniana"           
## [71] "Thrinax_compacta"             "Crysophila_williamsii"       
## [73] "Crysophila_warscewiczii"      "Crysophila_stauracantha"     
## [75] "Crysophila_nana"              "Crysophila_macrocarpa"       
## [77] "Crysophila_kalbreyeri"        "Crysophila_guagara"          
## [79] "Crysophila_grayumii"          "Crysophila_cookii"

 

2. Taxonomic standardization

Missing species from the phylogeny may be synonyms of (current) accepted species names, and are therefore undetected in the phylogenetic tree.
Let’s try standardizing the missing species names.

Write species missing from phylogeny into a vector

specmissing <- unique(species[which(!species$species_name %in% palmtree$tip.label), c(2:5)])
head(specmissing)
##      spp_Id    GENUS      EPITHET         species_name
## 1       550 Ammandra      natalia     Ammandra_natalia
## 173     548 Ammandra    dasyneura   Ammandra_dasyneura
## 808     541  Geonoma  weberbaueri  Geonoma_weberbaueri
## 3948    529  Geonoma    rubescens    Geonoma_rubescens
## 4035    528  Geonoma    polyandra    Geonoma_polyandra
## 4249    526  Geonoma poeppiginana Geonoma_poeppiginana
nrow(specmissing) #check the number
## [1] 80

 

Match palm species names to those in World Checklist of Vascular Plants, WCVP.
It is always important to check whether you have access to an updated taxonomy. We rely on the TNRS R package.

tnrs_names <- TNRS(taxonomic_names = specmissing$species_name)

# Alternative with package rWCVP (nb. the new version of the package is currently in active development so it's not on CRAN, you have to manually download it from GitHub, https://matildabrown.github.io/rWCVP/)
# taxstand <- wcvp_match_names(specmissing,
#                              name_col = "species",
#                              id_col = "spp_Id")
# head(taxstand)

 

Replace old names by new names

tnrs_names$Accepted_species <- gsub(" ", "_", tnrs_names$Accepted_species)

# update the combined species name
for(i in 1:nrow(tnrs_names)){
  species[which(species$species_name == tnrs_names[i, "Name_submitted"]), "species_name"] <-
    tnrs_names[i, "Accepted_species"]
}

 

Check names that are not in the phylogeny

specmissing <- unique(species$species[which(!species$species %in% palmtree$tip.label)])
specmissing
## [1] "Geonoma_longepedunculata" ""                        
## [3] "Leopoldinia_major"        "Acoelorraphe_wrightii"
length(specmissing)
## [1] 4

 

Instead of 80 missing species, we only have 4 after taxonomic harmonization.

Remaining missing species would have to be looked up and added manually to the phylogeny by an expert.
If missing species cannot be added reliably, we need to remove them from the dataset for analysis.

# Removal of the species not in the phylogeny
length(unique(species$species_name))
## [1] 527
species <- species[which(!species$species_name %in% specmissing), ]
length(unique(species$species_name))
## [1] 523

 

3. Spatial species composition at coarse resolution

3.1. Loading data

Read in polygon shapefile

# gridcells corresponding to palm distribution data and americas coastline
grid <- st_read("data/30min_grid_select50%.shp")
## Reading layer `30min_grid_select50%' from data source 
##   `P:\POSTDOC\Teaching\MCMM - Macroecology SE\MCMMB\data\30min_grid_select50%.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6038 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -116.0007 ymin: -35 xmax: -35.00066 ymax: 36.5
## CRS:           NA
americas <- st_read("data/americas.shp")
## Reading layer `americas' from data source 
##   `P:\POSTDOC\Teaching\MCMM - Macroecology SE\MCMMB\data\americas.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.7158 ymin: -55.91972 xmax: -29.84 ymax: 49.37666
## CRS:           NA
# remember to define the Coordinate Reference System (CRS)
st_crs(grid) <- "EPSG:4326" # "+proj=longlat +ellps=WGS84 +no_defs"
st_crs(americas) <- "EPSG:4326" # "+proj=longlat +ellps=WGS84 +no_defs"

# Rename ID column for matching
names(grid)[1] <- "grid_id"

 

3.2. Coarsening (i.e. reducing) the resolution


We will now convert the grid shapefile to a coarser resolution to have fewer grid cells.

To do so, we first create an empty grid of a 2-degree resolution.

# st_make_grid creates a grid covering the bounding box (i.e. the extent) of the geometry of an sf object (in this case, our "grid" object) 
grid_2degrees <- st_make_grid(x = grid, cellsize = 2)%>%
  st_sf() %>% 
  mutate(ID = row_number()) # add a unique ID to each grid cell

ggplot() +
  geom_sf(data = grid_2degrees, fill = "white", color = "black") +
  geom_sf(data = grid, fill = NA, color = "orange")

As we can see on the previous plot, the aggregated grid has been made after the extent of the grid. We therefore mask it to remove empty cells:

# Let's first "unite" the grid cells in our grid shapefile
grid_union <- st_union(grid)

sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
grid_2degrees <- st_intersection(grid_2degrees, grid_union)
## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
par(mfrow = c(1,1))
plot(st_geometry(grid_2degrees))

ggplot() +
  geom_sf(data = grid_2degrees, fill = "white", color = "black") +
  geom_sf(data = grid, fill = NA, color = "orange")+
  geom_sf(data = grid_2degrees, fill = NA, color = "black", linewidth = 1) +
  scale_x_continuous(limits = c(-90, -60)) +
  scale_y_continuous(limits = c(20, -5))

We now make a spatial join between the original grid and the coarser one. We have to specify largest = TRUE in the function st_join() to only keep the coarser cell with the bigger overlap.

grid <- st_join(grid, grid_2degrees, largest = TRUE) # grid now as a new column ID corresponding to the coarser cell size
## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

 

3.3. Add new ID also to species data.frame

species <- inner_join(species, st_drop_geometry(grid), by = "grid_id")
head(species)
##   grid_id spp_Id    GENUS EPITHET     species_name AREA PORTION  ID
## 1   13735    550 Ammandra natalia Aphandra_natalia 0.25     100 675
## 2   13736    550 Ammandra natalia Aphandra_natalia 0.25     100 676
## 3   13737    550 Ammandra natalia Aphandra_natalia 0.25     100 676
## 4   13738    550 Ammandra natalia Aphandra_natalia 0.25     100 676
## 5   13910    550 Ammandra natalia Aphandra_natalia 0.25     100 675
## 6   13911    550 Ammandra natalia Aphandra_natalia 0.25     100 676

We can now remove species duplicates per new ID.

species_coarse <- species %>% 
  group_by(ID) %>% 
  distinct(species_name, .keep_all = TRUE) %>% 
  ungroup() %>% 
  dplyr::select(-grid_id, -spp_Id, -GENUS, -EPITHET, -AREA, -PORTION)
Alternative way: manual aggregation with a for-loop and the grid centroids
Alternative

Extract the grid centroids

# convert into points shapefile
grid_centroids <- st_centroid(grid)

ggplot() +
  geom_sf(data = grid, fill = NA) +
  geom_sf(data = grid_centroids, shape = 20) +
  geom_sf(data = americas, fill = NA, color = "black") +
  theme_void() +
  scale_x_continuous(limits = c(-80, -60)) +
  scale_y_continuous(limits = c(20, 5))

Add coordinates to SpatialPolygonsDataframe

# Add latitude and longitude of grid centroids as extra columns
head(grid)
grid <- cbind(grid, st_coordinates(grid_centroids))
names(grid)[4:5] <- c("longitude", "latitude")
head(grid)

 

Make new IDs for neighboring cells
Two steps: i) we build an empty matrix with the new IDs, ii) with a for-loop, we add the new IDs into the grid object.

range(grid$longitude)
range(grid$latitude)

longitude <- seq(-116, -34, by = 2)
latitude  <- seq(-35, 35, by = 2)

new_ID_matrix <- matrix(c(1:(length(longitude)*length(latitude))),
                        nrow = length(longitude), ncol = length(latitude),
                        dimnames = list(longitude, latitude), byrow = TRUE)
dim(new_ID_matrix)
new_ID_matrix[1:5, 1:5]

# Filling
grid$new_ID <- NA

for(i in 1:nrow(grid)){
  grid$new_ID[i] <- new_ID_matrix[which(longitude < grid$longitude[i] & 
                                          (longitude + 2) > grid$longitude[i]),
                                  which(latitude < grid$latitude[i] &
                                          (latitude + 2) > grid$latitude[i])]
}

head(grid)

 

Join neighboring cells based on new IDs & aggregate species distribution data using the new grid

grid_2degrees <- group_by(grid, new_ID)
grid_2degrees <- summarise(grid_2degrees, do_union = TRUE, is_coverage = TRUE)

head(grid_2degrees)

ggplot(grid_2degrees) +
  geom_sf(fill = "gray90") +
  geom_sf(data = americas, fill = NA, color = "black") +
  theme_void() + ylim(-40, 40)

# Add new ID also to species data.frame
species <- inner_join(species, st_drop_geometry(grid), by = "grid_id")
head(species)

 

 

4. Calculating phylogenetic community metrics

4.1. Faith’s PD

Prune phylogeny to exclude species not in palm distribution data
Make sure that the same number of species in your data set are in your phylogeny

palmtree_pruned <- drop.tip(palmtree,
                            palmtree$tip.label[which(!palmtree$tip.label %in% unique(species$species_name))])
length(palmtree$tip.label)
## [1] 2539
length(palmtree_pruned$tip.label)
## [1] 518
length(unique(species$species_name))
## [1] 518
# You can also save the pruned tree 
#write.nexus(palmtree_pruned, file = "data/palmtree_pruned.nex")

 

Convert long table to species by grid-cell table (a community matrix)

NB. We need to do this operation because the functions we will use to estimate phylogenetic diversity requires a community matrix (i.e. species as columns and plots/grid-cells as rows)

library("bioregion")
species_grid <- net_to_mat(species_coarse[, c("ID", "species_name")])

# Alternative without the package bioregion
# species <- unique(species[, c("new_ID", "species")])
# species_grid <- as.data.frame.matrix(table(species_coarse))
# species_grid <- as.matrix(species_grid)

dim(species_grid)
## [1] 479 518
species_grid[1:5, 1:5]
##     Aphandra_natalia Ammandra_decasperma Phytelephas_tumacana
## 675                1                   0                    0
## 676                1                   0                    0
## 677                1                   0                    0
## 634                1                   0                    0
## 635                1                   0                    0
##     Phytelephas_tenuicaulis Phytelephas_seemannii
## 675                       1                     0
## 676                       1                     0
## 677                       1                     0
## 634                       1                     0
## 635                       1                     0
# Save the community matrix
# write.csv(species_grid, file = "data/palms_specsxsites_phylo.csv")

 

# if you have a hard time reading in the large palm phylogeny or
# in case you got stuck above
palmtree_pruned <- read.nexus("data/palmtree_pruned.nex")
species_grid <- as.matrix(read.csv("data/palms_specsxsites_phylo.csv", row.names = 1))

Calculate phylogenetic diversity

We here calculate the Faith’s PD (phylogenetic diversity index), which can be written like this:
\[ PD = \sum_{branches}l(b) \] with \(l\) the length of a branch \(b\).Basically, this index sums the length of branches connecting all species in the assemblage.

# ?pd
pd_palms <- pd(species_grid, palmtree_pruned)
dim(pd_palms) # one value of PD per grid cell (plus a value for Species Richness, SR)
## [1] 479   2
head(pd_palms)
##           PD SR
## 675 1575.824 69
## 676 1465.944 60
## 677 1524.880 62
## 634 1467.002 61
## 635 1525.790 62
## 636 1642.958 71
# Adding the ID of grid cells
pd_palms <- data.frame(ID = as.integer(row.names(species_grid)),
                       pd_palms)

head(pd_palms)
##      ID       PD SR
## 675 675 1575.824 69
## 676 676 1465.944 60
## 677 677 1524.880 62
## 634 634 1467.002 61
## 635 635 1525.790 62
## 636 636 1642.958 71
# Alternative way using tydiverse:
# pd_palms <- pd_palms %>% 
#   tibble::rownames_to_column(var = "ID")

 

pd() calculates Faith’s PD, and a corrected version of it being the residuals of a regression with total abundance or species richness per plot. A positive residual means the sample has more phylogenetic diversity than expected given its species count (i.e., its species are more distantly related than average). A negative residual means less diversity than expected. Let’s look into this:

pd_model <- lm(PD ~ SR, data = pd_palms)
pd_palms$residuals <- residuals(pd_model)
head(pd_palms)
##      ID       PD SR  residuals
## 675 675 1575.824 69 -29.973347
## 676 676 1465.944 60  44.204392
## 677 677 1524.880 62  62.238817
## 634 634 1467.002 61  24.812397
## 635 635 1525.790 62  63.148715
## 636 636 1642.958 71  -3.741623
plot_grid(
  nrow = 1, ncol = 2,
  ggplot(pd_palms, aes(SR, PD)) +
    geom_point(color = "black", alpha = 0.3) +
    stat_smooth(method = "lm", formula = y ~ x, color = "red") +
    labs(x = "Species_number", y = "Faith's PD") +
    theme_bw(),
  ggplot(pd_palms, aes(SR, residuals)) +
    geom_point(color = "black", alpha = 0.3) +
    labs(x = "Species_number", y = "Faith's PD residuals") +
    stat_smooth(method = "lm", formula = y ~ x, color = "red") +
    theme_bw()
)

4.2. Null models and other PD indices

Calculate generic phylogenetic metrics (absolute and standardized)

Phylogenetic diversity metrics are standardized by first using a null model to generate expected values.
The principle of null models is illustrated in the following figure:

Then, the observed value is compared to the null, usually using standardized effect sizes (SES):
\[SES = (observed - mean.null)/sd.null\] where mean.null is the mean of the null values and sd.null is the standard deviation of the null values.
A positive SES value indicates phylogenetic overdispersion (more diverse than expected by chance). A negative SES indicates phylogenetic clustering.

Null models can also apply to other phylogenetic diversity indices like the Mean Pairwise Distance (MPD) and the Mean Nearest Taxon Distance (MNTD).

MPD consists in calculating all the pairwise distances between species in a grid cell and then takes the mean.
MNTD calculates the mean distance separating each species in the grid cell from its closest relative. The primary difference is that MPD measures the average relatedness of all species to one another, while MNTD measures how close species are to their closest relative (evaluating recent evolutionary divergence at the tips of the tree). These indices are described in Webb el al. (2002).

We will use some functions from the pez package to calculate the observed and standardized effect size (SES) of these metrics:

c.data <- comparative.comm(palmtree_pruned , species_grid)
pd.null <- generic.null(c.data, c(.pd, .mpd, .mntd),
                        null.model = "richness",
                        comp.fun = .ses, permute = 100)

colnames(pd.null) <- c("Faith_PD", "Corrected_FaithPD", "MPD", "MNTD")
rownames(pd.null) <- sites(c.data)
pd.null <- data.frame(pd.null)
pd.null$ID <- as.numeric(rownames(pd.null))

dim(pd.null); head(pd.null)
## [1] 479  21
##      Faith_PD.observed Corrected_FaithPD.observed MPD.observed MNTD.observed
## 1000         1161.9594                  210.59122     126.3598      44.62270
## 1001         1020.6110                  253.30102     123.3111      51.51361
## 1034          475.7734                   97.03081     141.8437      64.31507
## 1035          543.9522                   42.50417     122.3554      44.78403
## 1036          862.3636                   13.25003     109.3301      34.82029
## 1037         1101.8782                  -53.99906     105.3146      28.25024
##      Faith_PD.null.mean Corrected_FaithPD.null.mean MPD.null.mean
## 1000          1215.3054                   64.096509      136.6966
## 1001          1032.2442                   98.695874      137.1126
## 1034           495.3685                   21.325890      134.7540
## 1035           699.2178                   80.068061      136.6099
## 1036          1105.0137                   74.727306      135.6249
## 1037          1399.4184                    6.364399      135.4209
##      MNTD.null.mean Faith_PD.SE Corrected_FaithPD.SE    MPD.SE   MNTD.SE
## 1000       39.80326    9.245731             9.195371 0.5962760 0.4555122
## 1001       44.95946    8.594374             8.582740 0.6584741 0.7079973
## 1034       74.98508    5.564996             5.541187 1.1250735 1.7148519
## 1035       59.27905    7.146856             7.155512 0.9184345 1.1266196
## 1036       42.09754    8.285158             8.255822 0.5656292 0.5896826
## 1037       35.59100    9.426104             9.417544 0.4814656 0.3781952
##      Faith_PD.SES Corrected_FaithPD.SES    MPD.SES   MNTD.SES Faith_PD.rank
## 1000    -5.769804             15.931355 -17.335686  10.580261          0.31
## 1001    -1.353588             18.013496 -20.959817   9.257300          0.48
## 1034    -3.521137             13.662222   6.301544  -6.222118          0.35
## 1035   -21.725009             -5.249644 -15.520396 -12.865938          0.04
## 1036   -29.287320             -7.446536 -46.487730 -12.340950          0.01
## 1037   -31.565556             -6.409681 -62.530439 -19.409973          0.02
##      Corrected_FaithPD.rank MPD.rank MNTD.rank   ID
## 1000                   0.95     0.07      0.83 1000
## 1001                   0.97     0.03      0.83 1001
## 1034                   0.92     0.73      0.25 1034
## 1035                   0.31     0.06      0.13 1035
## 1036                   0.29     0.01      0.12 1036
## 1037                   0.22     0.01      0.04 1037
# select columns of interest for final data set
pd.out <- pd.null[, c("ID", "Faith_PD.observed", "Faith_PD.SES",
                      "MPD.observed", "MPD.SES", "MNTD.observed", "MNTD.SES")] 
head(pd.out)
##        ID Faith_PD.observed Faith_PD.SES MPD.observed    MPD.SES MNTD.observed
## 1000 1000         1161.9594    -5.769804     126.3598 -17.335686      44.62270
## 1001 1001         1020.6110    -1.353588     123.3111 -20.959817      51.51361
## 1034 1034          475.7734    -3.521137     141.8437   6.301544      64.31507
## 1035 1035          543.9522   -21.725009     122.3554 -15.520396      44.78403
## 1036 1036          862.3636   -29.287320     109.3301 -46.487730      34.82029
## 1037 1037         1101.8782   -31.565556     105.3146 -62.530439      28.25024
##        MNTD.SES
## 1000  10.580261
## 1001   9.257300
## 1034  -6.222118
## 1035 -12.865938
## 1036 -12.340950
## 1037 -19.409973
# change NAs to zeroes; zero phylogenetic diversity if grid has only 1 spp!
pd.out[is.na(pd.out)] <- 0

head(pd.out)
##        ID Faith_PD.observed Faith_PD.SES MPD.observed    MPD.SES MNTD.observed
## 1000 1000         1161.9594    -5.769804     126.3598 -17.335686      44.62270
## 1001 1001         1020.6110    -1.353588     123.3111 -20.959817      51.51361
## 1034 1034          475.7734    -3.521137     141.8437   6.301544      64.31507
## 1035 1035          543.9522   -21.725009     122.3554 -15.520396      44.78403
## 1036 1036          862.3636   -29.287320     109.3301 -46.487730      34.82029
## 1037 1037         1101.8782   -31.565556     105.3146 -62.530439      28.25024
##        MNTD.SES
## 1000  10.580261
## 1001   9.257300
## 1034  -6.222118
## 1035 -12.865938
## 1036 -12.340950
## 1037 -19.409973
head(pd_palms)
##      ID       PD SR  residuals
## 675 675 1575.824 69 -29.973347
## 676 676 1465.944 60  44.204392
## 677 677 1524.880 62  62.238817
## 634 634 1467.002 61  24.812397
## 635 635 1525.790 62  63.148715
## 636 636 1642.958 71  -3.741623
palmdiversity <- left_join(pd_palms, pd.out, by = "ID")
head(palmdiversity)
##    ID       PD SR  residuals Faith_PD.observed Faith_PD.SES MPD.observed
## 1 675 1575.824 69 -29.973347          1575.824    -22.64146     125.8718
## 2 676 1465.944 60  44.204392          1465.944    -20.16233     126.7775
## 3 677 1524.880 62  62.238817          1524.880    -14.64466     127.4341
## 4 634 1467.002 61  24.812397          1467.002    -17.51154     130.6003
## 5 635 1525.790 62  63.148715          1525.790    -14.27851     129.7774
## 6 636 1642.958 71  -3.741623          1642.958    -16.91324     128.3756
##     MPD.SES MNTD.observed   MNTD.SES
## 1 -34.21616      26.52567 -13.819731
## 2 -30.35257      27.22376 -15.956642
## 3 -29.61332      26.70619 -14.736498
## 4 -20.95976      28.74301  -8.141586
## 5 -21.56188      28.67301 -10.356458
## 6 -24.68063      25.54811 -15.702618
# Plot
plot_grid(
  nrow = 1, ncol = 2,
  
  ggplot(palmdiversity, aes(SR, Faith_PD.SES)) +
    geom_point(color = "black", alpha = 0.3) +
    stat_smooth(method = "lm", formula = y ~ x, color = "red") +
    labs(x = "Species_number", y = "Faith's PD standardized") +
    theme_bw(),
  
  ggplot(palmdiversity, aes(Faith_PD.observed, Faith_PD.SES)) +
    geom_point(color = "black", alpha = 0.3) +
    stat_smooth(method = "lm", formula = y ~ x, color = "red") +
    labs(x = "Faith's PD observed", y = "Faith's PD standardized") +
    theme_bw()
)

 

Plot correlations among the diversity metrics in the palm diversity data.frame (use ggpairs()).

ggpairs(palmdiversity[, c(3:10)],
        upper = list(continuous = wrap(ggally_cor, digits = 1))) +
  theme_bw()

# Alternative:
# psych::pairs.panels(palmdiversity[, c(3:10)], method = "pearson",
#             density = FALSE, ellipses = FALSE, hist.col = "white")

 

Join the palmdiversity data.frame to the grid_2degrees shapefile to be able to link the diversity metrics and the poylgons for plotting.

head(grid_2degrees); head(palmdiversity)
## Simple feature collection with 6 features and 1 field
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -72.00066 ymin: -35 xmax: -58.00066 ymax: -33
## Geodetic CRS:  WGS 84
##    ID                       geometry
## 22 22 MULTILINESTRING ((-72.00066...
## 23 23 POLYGON ((-71.50066 -35, -7...
## 26 26 GEOMETRYCOLLECTION (POLYGON...
## 27 27 POLYGON ((-64.00066 -33, -6...
## 28 28 POLYGON ((-62.00066 -33, -6...
## 29 29 MULTIPOLYGON (((-58.00066 -...
##    ID       PD SR  residuals Faith_PD.observed Faith_PD.SES MPD.observed
## 1 675 1575.824 69 -29.973347          1575.824    -22.64146     125.8718
## 2 676 1465.944 60  44.204392          1465.944    -20.16233     126.7775
## 3 677 1524.880 62  62.238817          1524.880    -14.64466     127.4341
## 4 634 1467.002 61  24.812397          1467.002    -17.51154     130.6003
## 5 635 1525.790 62  63.148715          1525.790    -14.27851     129.7774
## 6 636 1642.958 71  -3.741623          1642.958    -16.91324     128.3756
##     MPD.SES MNTD.observed   MNTD.SES
## 1 -34.21616      26.52567 -13.819731
## 2 -30.35257      27.22376 -15.956642
## 3 -29.61332      26.70619 -14.736498
## 4 -20.95976      28.74301  -8.141586
## 5 -21.56188      28.67301 -10.356458
## 6 -24.68063      25.54811 -15.702618
grid_2degrees <- left_join(grid_2degrees, palmdiversity, by = "ID")
head(grid_2degrees)
## Simple feature collection with 6 features and 10 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -72.00066 ymin: -35 xmax: -58.00066 ymax: -33
## Geodetic CRS:  WGS 84
##   ID      PD SR residuals Faith_PD.observed Faith_PD.SES MPD.observed MPD.SES
## 1 22      NA NA        NA                NA           NA           NA      NA
## 2 23 104.614  1 -110.5213           104.614    0.8066271            0       0
## 3 26 104.614  1 -110.5213           104.614    6.1761888            0       0
## 4 27 104.614  1 -110.5213           104.614    5.3427762            0       0
## 5 28 104.614  1 -110.5213           104.614    6.2064961            0       0
## 6 29 104.614  1 -110.5213           104.614    7.9920893            0       0
##   MNTD.observed MNTD.SES                       geometry
## 1            NA       NA MULTILINESTRING ((-72.00066...
## 2             0        0 POLYGON ((-71.50066 -35, -7...
## 3             0        0 GEOMETRYCOLLECTION (POLYGON...
## 4             0        0 POLYGON ((-64.00066 -33, -6...
## 5             0        0 POLYGON ((-62.00066 -33, -6...
## 6             0        0 MULTIPOLYGON (((-58.00066 -...

 

Plot maps with colours according to the diversity metrics.

Species richness, PD, MPD & MNTD (absolute)

plot_grid(
  nrow = 2, ncol = 2,
  ggplot(grid_2degrees[complete.cases(grid_2degrees$SR), ]) +
    geom_sf(color = NA, aes(fill = SR)) +
    geom_sf(data = americas, fill = NA, color = "black") +
    labs(title = "Palm species richness") +
    scale_fill_viridis_c("Species\nnumber") +
    theme_void() +
    theme(plot.title = element_text(margin = margin(0, 0, 10, 0))) + 
    ylim(-40, 40),
  
  ggplot(grid_2degrees[complete.cases(grid_2degrees$SR), ]) +
    geom_sf(color = NA, aes(fill = PD)) +
    geom_sf(data = americas, fill = NA, color = "black") +
    labs(title = "Faith's Phylogenetic diversity") +
    scale_fill_viridis_c("PD (my)") +
    theme_void() +
    theme(plot.title = element_text(margin = margin(0, 0, 10, 0))) + 
    ylim(-40, 40),
  
  ggplot(grid_2degrees[complete.cases(grid_2degrees$SR), ]) +
    geom_sf(color = NA, aes(fill = MPD.observed)) +
    geom_sf(data = americas, fill = NA, color = "black") +
    labs(title = "Mean pairwise distance") +
    scale_fill_viridis_c("MPD (my)") +
    theme_void() +
    theme(plot.title = element_text(margin = margin(0, 0, 10, 0))) + 
    ylim(-40, 40),
  
  ggplot(grid_2degrees[complete.cases(grid_2degrees$SR), ]) +
    geom_sf(color = NA, aes(fill = MNTD.observed)) +
    geom_sf(data = americas, fill = NA, color = "black") +
    labs(title = "Mean nearest taxon distance") +
    scale_fill_viridis_c("MNTD (my)") +
    theme_void() +
    theme(plot.title = element_text(margin = margin(0, 0, 10, 0))) + 
    ylim(-40, 40)
)

 

PD, MPD & MNTD (standardized)
Create maps.

 

References

Faurby S, Eiserhardt WL, Baker WJ, Svenning JC. An all-evidence species-level supertree for the palms (Arecaceae). Mol Phylogenet Evol. 2016 Jul;100:57-69.

Faith D.P. (1992) Conservation evaluation and phylogenetic diversity. Biological Conservation, 61, 1-10

Kreft, H., Sommer, J.H. and Barthlott, W. (2006), The significance of geographic range size for spatial diversity patterns in Neotropical palms. Ecography, 29: 21-30.

World Checklist of Vascular Plants

Phylogenies and Community Ecology Campbell O. Webb, David D. Ackerly, Mark A. McPeek, Michael J. Donoghue Annual Review of Ecology and Systematics 2002 33:1, 475-505