Modern concepts and methods in Macroecology and Biogeography
Day 8: Functional biogeography
Pierre Denelle, Nicola Pavanetto & Holger Kreft
June 10th, 2026
 back to overview  Â
Preparation: Please open your RStudio project and
download the new data (PalmTraits_1.0.txt,
palm_species_per_gridcell_2degrees.rds,
and the shapefile from the 2degrees_grid
folder) 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():
dplyrRColorBrewerpsychade4fundiversitysfggplot2GGallyfactoextracowplotbioregionfunspaceleaflet
Â
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/Day8.html", destfile = file.path(dir, "Day8.html"))
htmlFile <- file.path(dir, "Day8.html")
rstudioapi::viewer(htmlFile)Now you can conveniently copy code from the viewer into your script.
Load R packages & island data set
library("ade4") # Dissimilarity matrix and PCoA
library("dplyr") # data.frame operations
library("factoextra") # PCA
library("fundiversity") # Functional diversity
library("ggplot2") # plot
library("GGally") # pair plot
library("psych") # pair plot
library("RColorBrewer") # color gradients
library("sf") # spatial analyses
library("cowplot") # combine several plots
library("bioregion") # conversion between long and wide formats
library("funspace") # create and plot functional trait spaces
library("leaflet") # Create interactive mapsÂ
1. Trait dataset
1.1. Trait data
For today’s practical we will use a database of palm traits from Kissling, W. D., Balslev, H., Baker, W. J., Dransfield, J., Göldel, B., Lim, J. Y., Onstein, R. E., & Svenning, J.-C. (2019). PalmTraits 1.0, a species-level functional trait database of palms worldwide. Scientific Data, 6(1), 1–13.. You can download the dataset as a .txt format from StudIP by clicking here.
# Importing the dataset
tra <- read.table("data/PalmTraits_1.0.txt",
stringsAsFactors = FALSE, sep = "\t", header = TRUE)
dim(tra)## [1] 2557 29head(tra)## SpecName accGenus accSpecies PalmTribe
## 1 Acanthophoenix crinita Acanthophoenix crinita Areceae
## 2 Acanthophoenix rousselii Acanthophoenix rousselii Areceae
## 3 Acanthophoenix rubra Acanthophoenix rubra Areceae
## 4 Acoelorrhaphe wrightii Acoelorrhaphe wrightii Trachycarpeae
## 5 Acrocomia aculeata Acrocomia aculeata Cocoseae
## 6 Acrocomia crispa Acrocomia crispa Cocoseae
## PalmSubfamily Climbing Acaulescent Erect StemSolitary StemArmed LeavesArmed
## 1 Arecoideae 0 0 1 1 1 1
## 2 Arecoideae 0 0 1 1 1 1
## 3 Arecoideae 0 0 1 1 1 1
## 4 Coryphoideae 0 0 1 0 0 1
## 5 Arecoideae 0 0 1 1 1 1
## 6 Arecoideae 0 0 1 1 1 1
## MaxStemHeight_m MaxStemDia_cm UnderstoreyCanopy MaxLeafNumber
## 1 10.0 20 canopy 15
## 2 25.0 30 canopy NA
## 3 15.0 18 canopy 20
## 4 9.1 15 canopy 25
## 5 12.0 50 canopy 30
## 6 18.0 35 canopy 15
## Max_Blade_Length_m Max_Rachis_Length_m Max_Petiole_length_m
## 1 2.3 NA NA
## 2 3.0 NA NA
## 3 3.1 3.0 NA
## 4 1.3 0.7 0.65
## 5 3.5 2.5 NA
## 6 3.0 NA 0.65
## AverageFruitLength_cm MinFruitLength_cm MaxFruitLength_cm
## 1 0.65 0.6 0.7
## 2 2.00 NA NA
## 3 1.00 NA NA
## 4 0.70 NA NA
## 5 4.25 3.5 5.0
## 6 2.50 NA NA
## AverageFruitWidth_cm MinFruitWidth_cm MaxFruitWidth_cm FruitSizeCategorical
## 1 0.5 NA NA small
## 2 0.8 NA NA small
## 3 0.7 NA NA small
## 4 0.7 0.5 0.9 small
## 5 4.6 3.8 5.4 large
## 6 1.8 NA NA small
## FruitShape FruitColorDescription MainFruitColors Conspicuousness
## 1 black black cryptic
## 2 ovoid black black cryptic
## 3 ovoid black black cryptic
## 4 ovoid orange-brown; becomming black brown; black cryptic
## 5 ovoid yellowish green green cryptic
## 6 globose yellow to orange yellow; orange conspicuousWhen using a functional trait table, we should look at:
- trait type (continuous, ordinal, categorical, binary)
- trait coverage (percentage of species having a trait value)
- trait distribution (if too asymmetrical, the trait can be transformed)
- trait correlation (are some traits carrying the same information?)
- ecological significance (is the trait varying against environment? A biotic variable? A physiological process?)
Let’s go through these points.
1.1.1. Trait type
The traits available relate to growth form and habit, armature (e.g. spinescence), stem size, leaves and fruits. They can be either binary traits (0 or 1, sometimes 2 when one species varies in that trait), categorical or continuous.
str(tra)## 'data.frame': 2557 obs. of 29 variables:
## $ SpecName : chr "Acanthophoenix crinita" "Acanthophoenix rousselii" "Acanthophoenix rubra" "Acoelorrhaphe wrightii" ...
## $ accGenus : chr "Acanthophoenix" "Acanthophoenix" "Acanthophoenix" "Acoelorrhaphe" ...
## $ accSpecies : chr "crinita" "rousselii" "rubra" "wrightii" ...
## $ PalmTribe : chr "Areceae" "Areceae" "Areceae" "Trachycarpeae" ...
## $ PalmSubfamily : chr "Arecoideae" "Arecoideae" "Arecoideae" "Coryphoideae" ...
## $ Climbing : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Acaulescent : int 0 0 0 0 0 0 1 0 1 0 ...
## $ Erect : int 1 1 1 1 1 1 0 1 0 1 ...
## $ StemSolitary : int 1 1 1 0 1 1 1 1 1 1 ...
## $ StemArmed : int 1 1 1 0 1 1 1 1 1 1 ...
## $ LeavesArmed : int 1 1 1 1 1 1 1 1 1 1 ...
## $ MaxStemHeight_m : num 10 25 15 9.1 12 18 0 NA 0 NA ...
## $ MaxStemDia_cm : num 20 30 18 15 50 35 NA NA NA NA ...
## $ UnderstoreyCanopy : chr "canopy" "canopy" "canopy" "canopy" ...
## $ MaxLeafNumber : int 15 NA 20 25 30 15 NA NA 6 NA ...
## $ Max_Blade_Length_m : num 2.3 3 3.1 1.3 3.5 3 NA NA 0.9 NA ...
## $ Max_Rachis_Length_m : num NA NA 3 0.7 2.5 NA NA NA 0.54 NA ...
## $ Max_Petiole_length_m : num NA NA NA 0.65 NA 0.65 NA NA 0.51 NA ...
## $ AverageFruitLength_cm: num 0.65 2 1 0.7 4.25 2.5 2 NA 2.25 4.6 ...
## $ MinFruitLength_cm : num 0.6 NA NA NA 3.5 NA NA NA 1.5 3.8 ...
## $ MaxFruitLength_cm : num 0.7 NA NA NA 5 NA NA NA 3 5.4 ...
## $ AverageFruitWidth_cm : num 0.5 0.8 0.7 0.7 4.6 1.8 2 NA 2.25 4.6 ...
## $ MinFruitWidth_cm : num NA NA NA 0.5 3.8 NA NA NA 1.5 3.8 ...
## $ MaxFruitWidth_cm : num NA NA NA 0.9 5.4 NA NA NA 3 5.4 ...
## $ FruitSizeCategorical : chr "small" "small" "small" "small" ...
## $ FruitShape : chr "" "ovoid" "ovoid" "ovoid" ...
## $ FruitColorDescription: chr "black" "black" "black" "orange-brown; becomming black" ...
## $ MainFruitColors : chr "black" "black" "black" "brown; black" ...
## $ Conspicuousness : chr "cryptic" "cryptic" "cryptic" "cryptic" ...1.1.2. Trait coverage
Percentage of species having a trait value informed per trait.
nrow(tra) # total number of species## [1] 2557summary(tra)## SpecName accGenus accSpecies PalmTribe
## Length:2557 Length:2557 Length:2557 Length:2557
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## PalmSubfamily Climbing Acaulescent Erect
## Length:2557 Min. :0.0000 Min. :0.0000 Min. :0.0000
## Class :character 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Mode :character Median :0.0000 Median :0.0000 Median :1.0000
## Mean :0.2198 Mean :0.1091 Mean :0.7302
## 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :2.0000 Max. :2.0000 Max. :2.0000
##
## StemSolitary StemArmed LeavesArmed MaxStemHeight_m
## Min. :0.0000 Min. :0.00000 Min. :0.0000 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.: 2.50
## Median :1.0000 Median :0.00000 Median :0.0000 Median : 6.00
## Mean :0.8167 Mean :0.08193 Mean :0.4165 Mean : 10.86
## 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.: 15.00
## Max. :2.0000 Max. :1.00000 Max. :1.0000 Max. :170.00
## NA's :375 NA's :55 NA's :120 NA's :446
## MaxStemDia_cm UnderstoreyCanopy MaxLeafNumber Max_Blade_Length_m
## Min. : 0.00 Length:2557 Min. : 4.00 Min. : 0.150
## 1st Qu.: 2.00 Class :character 1st Qu.: 8.00 1st Qu.: 1.000
## Median : 5.00 Mode :character Median :11.00 Median : 1.695
## Mean : 12.38 Mean :14.37 Mean : 2.374
## 3rd Qu.: 17.00 3rd Qu.:18.00 3rd Qu.: 3.000
## Max. :175.00 Max. :75.00 Max. :25.000
## NA's :602 NA's :1251 NA's :659
## Max_Rachis_Length_m Max_Petiole_length_m AverageFruitLength_cm
## Min. : 0.050 Min. :0.0000 Min. : 0.300
## 1st Qu.: 0.750 1st Qu.:0.2500 1st Qu.: 1.050
## Median : 1.500 Median :0.5500 Median : 1.500
## Mean : 1.972 Mean :0.8517 Mean : 2.196
## 3rd Qu.: 2.700 3rd Qu.:1.2500 3rd Qu.: 2.500
## Max. :18.500 Max. :6.7500 Max. :45.000
## NA's :1026 NA's :1347 NA's :505
## MinFruitLength_cm MaxFruitLength_cm AverageFruitWidth_cm MinFruitWidth_cm
## Min. : 0.300 Min. : 0.500 Min. : 0.200 Min. : 0.200
## 1st Qu.: 1.000 1st Qu.: 1.400 1st Qu.: 0.750 1st Qu.: 0.700
## Median : 1.500 Median : 2.000 Median : 1.050 Median : 1.000
## Mean : 2.181 Mean : 3.102 Mean : 1.594 Mean : 1.479
## 3rd Qu.: 2.500 3rd Qu.: 3.500 3rd Qu.: 1.800 3rd Qu.: 1.800
## Max. :40.000 Max. :50.000 Max. :20.000 Max. :13.000
## NA's :1651 NA's :1641 NA's :563 NA's :1563
## MaxFruitWidth_cm FruitSizeCategorical FruitShape FruitColorDescription
## Min. : 0.220 Length:2557 Length:2557 Length:2557
## 1st Qu.: 1.000 Class :character Class :character Class :character
## Median : 1.500 Mode :character Mode :character Mode :character
## Mean : 2.125
## 3rd Qu.: 2.500
## Max. :20.000
## NA's :1555
## MainFruitColors Conspicuousness
## Length:2557 Length:2557
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## # Summary example for continuous traits
summary(tra$AverageFruitWidth_cm)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.200 0.750 1.050 1.594 1.800 20.000 563summary(tra$MaxStemHeight_m)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 2.50 6.00 10.86 15.00 170.00 446# Binary traits
summary(tra$Climbing)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2198 0.0000 2.0000# Categorical traits
summary(tra$FruitSizeCategorical)## Length Class Mode
## 2557 character character# Not very informative, but we can use a continency table to obtain a more useful summary:
table(tra$FruitSizeCategorical) # 505 observation has a "empty" value ( ≠from NAs!)##
## large small
## 505 251 18011.1.3. Trait distribution
To visualize the distribution of a specific trait value, we can plot
a histogram, a boxplot or a violin plot.
Side note
RE: The R Graph Gallery provides
many nice examples of plot types along with R code.
# Example with Maximum stem height
plot_grid(
ggplot(tra) +
geom_boxplot(aes(x = 1, y = MaxStemHeight_m)),
ggplot(tra) +
geom_violin(aes(x = 1, y = MaxStemHeight_m)),
ggplot(tra) +
geom_histogram(aes(MaxStemHeight_m)),
nrow = 1)## Warning: Removed 446 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 446 rows containing non-finite outside the scale range
## (`stat_ydensity()`).## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.## Warning: Removed 446 rows containing non-finite outside the scale range
## (`stat_bin()`).
For categorical and binary traits, we can use the function
table() to retrieve the proportion of each value.
The next chunk illustrates for example how to plot the joined proportion
of palm species with armed leaves and stems.
Example of non-armed stem (Cocos nucifera) & armed-stem (Aiphanes minima)
table(tra$StemArmed); table(tra$LeavesArmed)##
## 0 1
## 2297 205##
## 0 1
## 1422 1015# Combined contingency table
table(tra$StemArmed, tra$LeavesArmed)##
## 0 1
## 0 1406 775
## 1 16 189# Plot
plot(table(tra$StemArmed, tra$LeavesArmed),
main = "Frequency of Armature classes",
xlab = "Stem Armed", ylab = "Leaves Armed")
prpX <- prop.table(table(tra$StemArmed))
prpY <- prop.table(table(tra$StemArmed, tra$LeavesArmed), margin = 1)
text(prpX/2 * c(1,-1) + 0:1,
prpY/2 * c(-1,-1,1,1) + c(1,1,0,0),
paste("n=", table(tra$StemArmed, tra$LeavesArmed), sep=""))
1.1.4. Trait correlation
Some traits can carry a very similar information. To see this, you
can look at the correlation between two traits.
One way to look at
this is to plot one trait directly against another. For example, in our
dataset we have two traits related to stem size: maximum height and
diameter of stem, and we can see how much they are correlated.
ggplot(tra, aes(MaxStemHeight_m, MaxStemDia_cm)) +
geom_point() +
annotate(x = 100, y = 100, geom = "label",
label = paste0("Pearson correlation coefficient = ",
round(cor(tra$MaxStemHeight_m,
tra$MaxStemDia_cm,
use = "complete.obs"), 2))) +
labs(x = "Maximum Stem Height (m)",
y = "Maximum Stem Diameter (cm)") +
theme_bw()## Warning: Removed 647 rows containing missing values or values outside the scale range
## (`geom_point()`).
The correlation value is not very high, but considering other categorical trait can add further insight into such plots. In the example below, using climbing feature as an additional grouping variable provides more details.
ggplot(tra, aes(MaxStemHeight_m, MaxStemDia_cm)) +
geom_point(alpha = 0.5,
aes(color = as.factor(Climbing))) +
scale_color_viridis_d("Climbing") +
labs(x = "Maximum Stem Height (m)",
y = "Maximum Stem Diameter (cm)") +
theme_bw() +
theme(legend.position = "bottom")## Warning: Removed 647 rows containing missing values or values outside the scale range
## (`geom_point()`).
When you have more than two traits to look at, the best way to visualize the correlation between multiple traits is to use a correlation pair plot. Let’s make one for all the traits related to leaves and fruits.
# Pair plot for leaves traits
ggpairs(tra[, c("MaxLeafNumber", "Max_Blade_Length_m",
"Max_Rachis_Length_m", "Max_Petiole_length_m")])
#Alternative:
# pairs.panels(tra[, c("MaxLeafNumber", "Max_Blade_Length_m",
# "Max_Rachis_Length_m", "Max_Petiole_length_m")],
# density = FALSE, ellipses = FALSE, hist.col = "grey")
# Pair plot for fruits traits
ggpairs(tra[, c("AverageFruitLength_cm", "MinFruitLength_cm",
"MaxFruitLength_cm", "AverageFruitWidth_cm",
"MinFruitWidth_cm", "MaxFruitWidth_cm")])
For the leaves traits, the last three traits are highly correlated. We can just keep the blade length and the number of leaves. Regarding the fruit traits, they are all highly correlated, we’ll take only AverageFruitLength_cm. Rule of thumb: usually keep traits that have a correlation r <0.65
1.2. Trait selection
Conclusion:
- to keep it easy: we just keep four continuous
traits
- stem: height and diameter are highly correlated, we just take the
height
- leaf: many leaf traits are correlated, we’ll take only Blade length
and the maximum number of leaves.
- fruit: many fruit traits are correlated, we’ll take only AverageFruitLength_cm
We here create a vector with our selected traits.
# Vector of traits we want to keep
tra_select <- c("MaxStemHeight_m", "MaxLeafNumber", "Max_Blade_Length_m",
"AverageFruitLength_cm")
And then we make a new species x trait table (species in rows
and traits in columns).
# Only species column and traits of interest
sp_tra <- tra[, c("SpecName", tra_select)]Let’s make a pair plot once again to ensure our selected traits are not highly correlated.
# Distribution of continuous traits
ggpairs(sp_tra[, tra_select])
The correlation coefficients are all below 65%. However, some traits
are highly asymmetrical.
This can heavily affect the construction of functional spaces and the
calculation of functional metrics.
We therefore
log-transform all the traits.
# Log-transforming all the continuous traits
sp_tra$log_height <- log(sp_tra$MaxStemHeight_m + 1) # adding 1 here to avoid -Inf values
sp_tra$log_leaf_nb <- log(sp_tra$MaxLeafNumber)
sp_tra$log_blade <- log(sp_tra$Max_Blade_Length_m)
sp_tra$log_fruit_length <- log(sp_tra$AverageFruitLength_cm)
# Example of log-transformed fruit length values
ggplot(sp_tra, aes(AverageFruitLength_cm, log_fruit_length)) +
geom_point() +
labs(x = "Fruit length (cm)", y = "Log-transformed fruit length") +
theme_bw()## Warning: Removed 505 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Change vector of traits of interest
tra_select <- c("log_height", "log_leaf_nb", "log_blade", "log_fruit_length")1.3. Final trait table
We here merge the trait data with the distribution 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.
We here import a transformed version of the data to speed up the
calculations. Now each cell is at a 2 degrees resolution, and species
occurrences have been aggregated into these bigger cells. See practical
7 for more details.
species <- readRDS("data/palm_species_per_gridcell_2degrees.rds")
length(unique(species$new_ID)) # number of grid cells## [1] 479n_distinct(species$species_name) # number of species## [1] 545We can remove from our trait table all the species that are not
present in our species distribution table.
# Look at what species names in the trait dataset are present in the grid distribution dataset
dim(sp_tra[which(sp_tra$SpecName %in% unique(species$species_name)), ])## [1] 0 9head(species$species_name)## [1] "Ammandra_natalia" "Ammandra_natalia" "Ammandra_natalia" "Ammandra_natalia"
## [5] "Ammandra_natalia" "Ammandra_natalia"head(sp_tra$SpecName)## [1] "Acanthophoenix crinita" "Acanthophoenix rousselii"
## [3] "Acanthophoenix rubra" "Acoelorrhaphe wrightii"
## [5] "Acrocomia aculeata" "Acrocomia crispa"# Replace space with underscore in trait table before merging
sp_tra$species_name <- gsub(" ", "_", sp_tra$SpecName)
table(sp_tra$species_name %in% species$species_name)##
## FALSE TRUE
## 2095 462sp_tra <- sp_tra[which(sp_tra$species_name %in% species$species_name), ]
dim(sp_tra)## [1] 462 10head(sp_tra)## SpecName MaxStemHeight_m MaxLeafNumber Max_Blade_Length_m
## 5 Acrocomia aculeata 12 30 3.50
## 9 Acrocomia hassleri 0 6 0.90
## 19 Aiphanes acaulis 0 10 1.83
## 22 Aiphanes chiribogensis 3 9 1.45
## 23 Aiphanes deltoidea 2 12 2.95
## 24 Aiphanes duquei 5 9 1.35
## AverageFruitLength_cm log_height log_leaf_nb log_blade log_fruit_length
## 5 4.25 2.564949 3.401197 1.2527630 1.44691898
## 9 2.25 0.000000 1.791759 -0.1053605 0.81093022
## 19 1.00 0.000000 2.302585 0.6043160 0.00000000
## 22 1.10 1.386294 2.197225 0.3715636 0.09531018
## 23 1.30 1.098612 2.484907 1.0818052 0.26236426
## 24 1.05 1.791759 2.197225 0.3001046 0.04879016
## species_name
## 5 Acrocomia_aculeata
## 9 Acrocomia_hassleri
## 19 Aiphanes_acaulis
## 22 Aiphanes_chiribogensis
## 23 Aiphanes_deltoidea
## 24 Aiphanes_duquei2. Trait maps
Since our distribution data is linked to a georeferenced grid cell, we can map the distribution of our traits.
We first import the georeferenced grid (at 2 degrees resolution, you can find the shapefile here) and the outline of Americas. The CRS is missing, and has to be filled manually. Latitude and longitude centroids have already been calculated.
# Spatial data
grid.2 <- st_read("data/2degrees_grid/2degrees_grid.shp")
americas <- st_read("data/americas.shp")## Reading layer `2degrees_grid' from data source
## `P:\POSTDOC\Teaching\MCMM - Macroecology SE\MCMMB\data\2degrees_grid\2degrees_grid.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 494 features and 1 field (with 3 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -116.0007 ymin: -35 xmax: -35.00066 ymax: 36.5
## Geodetic CRS: WGS 84## 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# coordinate system
st_crs(grid.2) <- "EPSG:4326" # "+proj=longlat +ellps=WGS84 +no_defs"
st_crs(americas) <- "EPSG:4326" # "+proj=longlat +ellps=WGS84 +no_defs"
plot(st_geometry(grid.2))
Â
2.1. Trait geographical distribution
We here aim at mapping the distribution of one trait. We will do it
using maximal stem height.
First, we merge the species x grid cell table with our trait table.
# Merging palm height with spatial data frame
species_height <- left_join(
species[, c("new_ID", "species_name")],
sp_tra[, c("species_name", "MaxStemHeight_m")],
by = "species_name")
head(species_height)## # A tibble: 6 × 3
## new_ID species_name MaxStemHeight_m
## <int> <chr> <dbl>
## 1 675 Ammandra_natalia NA
## 2 676 Ammandra_natalia NA
## 3 677 Ammandra_natalia NA
## 4 634 Ammandra_natalia NA
## 5 635 Ammandra_natalia NA
## 6 636 Ammandra_natalia NAWe can then calculate the average species height per grid cell.
# Calculating the average height per cell With dplyr
grid_mean_height <- species_height %>%
group_by(new_ID) %>%
summarise(mean_height = mean(MaxStemHeight_m, na.rm = TRUE))
dim(grid_mean_height); head(grid_mean_height)## [1] 479 2## # A tibble: 6 × 2
## new_ID mean_height
## <int> <dbl>
## 1 23 25
## 2 26 6
## 3 27 6
## 4 28 6
## 5 29 6
## 6 30 17We then merge these values with our spatial grid object. Once this merge has been done, we can map the average height per cell!
# Plot
grid.2 <- left_join(grid.2, grid_mean_height, by = "new_ID")
ggplot(grid.2) +
geom_sf(aes(fill = mean_height), color = NA) +
geom_sf(data = americas, color = "black", fill = NA) +
scale_fill_viridis_c("Average height (m)") +
theme_bw()
2.2. Trait coverage
We saw that not all species have a trait value. This lack of
information can be spatially distributed in an uneven way. We therefore
need to make an assessment of the spatial trait coverage.
# Coverage calculation with dplyr:
grid_coverage_height <- species_height %>%
group_by(new_ID) %>%
summarise(coverage_height = sum(!is.na(MaxStemHeight_m))/
length(MaxStemHeight_m))
dim(grid_coverage_height); tail(grid_coverage_height)## [1] 479 2## # A tibble: 6 × 2
## new_ID coverage_height
## <int> <dbl>
## 1 1412 1
## 2 1413 1
## 3 1414 1
## 4 1436 1
## 5 1437 1
## 6 1455 1# Joining the coverage column and generating the plot
grid.2 <- left_join(grid.2, grid_coverage_height, by = "new_ID")
ggplot(grid.2) +
geom_sf(aes(fill = coverage_height), color = NA) +
geom_sf(data = americas, color = "black", fill = NA) +
scale_fill_viridis_c("Coverage (%)", option = "E") +
theme_bw()
To quickly visualize or access the informations in maps (e.g. height values in each grid cell/polygons) you could also generate interactive maps, e.g. with the leaflet package..
library(leaflet)
# Interactive plot of average palm height.
# Creating a color palette
pal <- colorNumeric(palette = "viridis", domain = grid.2$mean_height)
leaflet(data = grid.2)%>%
addPolygons(color = pal(grid.2$mean_height ),
popup = paste("Mean Height:",
grid.2$mean_height, "<br>",
"Trait coverage:", grid.2$coverage_height, "<br>",
"Cell ID:", grid.2$new_ID))# Investigating "suspicious" cells
species_height %>%
filter(new_ID == "1322")## # A tibble: 1 × 3
## new_ID species_name MaxStemHeight_m
## <int> <chr> <dbl>
## 1 1322 Sabal_minor 03. Functional diversity
When we have several functional traits, we can compute a set of
functional diversity indices that can inform us on the assembly
processes within a given assemblage of species.
We here focus on the following indices: functional richness, functional evenness, dispersion and divergence (see the lecture for definition of these indices).
3.1. FD per grid cell
We here at calculating and then mapping the functional richness per
grid cell using the fundiversity package.
# Before using the functions from fundiversity, we need to build a
# species-by-sites matrix
library("bioregion")
cell_sp <- net_to_mat(species[, c("new_ID", "species_name")])
# cell_sp <- as.matrix(table(species$new_ID, species$species))
cell_sp[1:5, 1:5] # sites as rows and species as columns## Ammandra_natalia Ammandra_decasperma Ammandra_dasyneura
## 675 1 0 0
## 676 1 0 0
## 677 1 0 0
## 634 1 0 0
## 635 1 0 0
## Phytelephas_tumacana Phytelephas_tenuicaulis
## 675 0 1
## 676 0 1
## 677 0 1
## 634 0 1
## 635 0 1# And we use the selected traits from our species x trait table
rownames(sp_tra) <- sp_tra$species_name
sp_tra[1:5, tra_select]## log_height log_leaf_nb log_blade log_fruit_length
## Acrocomia_aculeata 2.564949 3.401197 1.2527630 1.44691898
## Acrocomia_hassleri 0.000000 1.791759 -0.1053605 0.81093022
## Aiphanes_acaulis 0.000000 2.302585 0.6043160 0.00000000
## Aiphanes_chiribogensis 1.386294 2.197225 0.3715636 0.09531018
## Aiphanes_deltoidea 1.098612 2.484907 1.0818052 0.26236426# Functional richness
fd_palm <- fd_fric(traits = sp_tra[, tra_select],
sp_com = cell_sp)## Removed 45 species with missing trait(s)## Differing number of species between trait dataset and site-species matrix
## Taking subset of species## Warning: package 'future' was built under R version 4.5.3## Warning in fd_fric(traits = sp_tra[, tra_select], sp_com = cell_sp): Some sites
## had less species than traits so returned FRic is 'NA'Many cells (the ones with species number less than the number of traits) are not having a Frich value.
head(fd_palm)## site FRic
## 1 675 7.943896
## 2 676 10.795302
## 3 677 16.218305
## 4 634 6.537796
## 5 635 11.512031
## 6 636 16.681220summary(fd_palm$FRic)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000207 1.457130 4.794811 5.479266 8.866735 17.403105 108We can now add to this dataframe also other indices of functional diversity.
# Functional dispersion
fd_palm <- left_join(fd_palm,
fd_fdis(traits = sp_tra[, tra_select],
sp_com = cell_sp),
by = "site")## Removed 45 species with missing trait(s)## Differing number of species between trait dataset and site-species matrix
## Taking subset of species# Functional divergence
fd_palm <- left_join(fd_palm,
fd_fdiv(traits = sp_tra[, tra_select],
sp_com = cell_sp),
by = "site")## Removed 45 species with missing trait(s)
## Differing number of species between trait dataset and site-species matrix
## Taking subset of species# Functional evenness
fd_palm <- left_join(fd_palm,
fd_feve(traits = sp_tra[, tra_select],
sp_com = cell_sp),
by = "site")%>%
mutate(across(site, as.integer))## Removed 45 species with missing trait(s)
## Differing number of species between trait dataset and site-species matrix
## Taking subset of specieshead(fd_palm)## site FRic FDis FDiv FEve
## 1 675 7.943896 1.313072 0.7576665 0.8657444
## 2 676 10.795302 1.424242 0.7674720 0.8607406
## 3 677 16.218305 1.528642 0.7683450 0.8703439
## 4 634 6.537796 1.284308 0.7331159 0.8687230
## 5 635 11.512031 1.466312 0.7735406 0.8769224
## 6 636 16.681220 1.555033 0.7703221 0.8795797By doing a pair plot, we can have a look at the correlations between these different indices, and also with species richness.
# Adding species richness to grid
palm_SR <- species %>%
group_by(new_ID) %>%
summarise(SR = n())
fd_palm <- left_join(fd_palm, palm_SR, by = c("site" = "new_ID"))ggpairs(fd_palm[, c("SR", "FRic", "FDis", "FDiv", "FEve")])
Like for phylogenetic diversity, there is a positive correlation
between functional richness and species richness. To take this bias into
account, we should perform a null model. We won’t do it here, but the
principle is to reshuffle the trait values per cell a certain number of
times and then to compute the Standardized Effect Size (see practical
7).
3.2. Map of FD
Like for individual traits, we can also map the indeces of FD per cell.
# Merge results with grid
grid.2 <- left_join(grid.2, fd_palm, by = c("new_ID" = "site"))
plot_grid(
ggplot(grid.2) +
geom_sf(aes(fill = SR), color = NA) +
geom_sf(data = americas, color = "black", fill = NA) +
scale_fill_viridis_c("Species Richness") +
labs(title = "Species richness") +
theme_bw() +
theme(legend.position = "bottom"),
ggplot(grid.2) +
geom_sf(aes(fill = FRic), color = NA) +
geom_sf(data = americas, color = "black", fill = NA) +
scale_fill_viridis_c("Functional Richness") +
labs(title = "Functional richness") +
theme_bw() +
theme(legend.position = "bottom"),
ggplot(grid.2) +
geom_sf(aes(fill = FDis), color = NA) +
geom_sf(data = americas, color = "black", fill = NA) +
scale_fill_viridis_c("Functional dispersion") +
labs(title = "Functional dispersion") +
theme_bw() +
theme(legend.position = "bottom"),
nrow = 1)
4. Trait space
Another classical analysis in functional ecology is to build trait spaces that summarize the main axes of variations between species.
The main objective is to reduce dimensionality, in order to get the main axes of variation between species, through the use of multivariate analyses.
Since we only have continuous traits, we can build this space using a PCA.
Side-node
If there are some categorical traits, you can first calculate a distance
matrix and then compute a Principal Coordinates Analysis (PCoA). You can
find more details here: Applied
Multivariate Statistics in R, Jonathan D. Bakker.
We use the princomp function to to build the PCA.
# Selecting the trait columns (RE: we are using the log-transformed trait values)
sp_tra_pca <- sp_tra[, tra_select]
# Removing NAs
sp_tra_pca <- sp_tra_pca[complete.cases(sp_tra_pca), ]
# Perfoming the Principal Component Analysis
pca_palm <- princomp(sp_tra_pca, cor = FALSE)
# Plotting the results
factoextra::fviz_pca_biplot(pca_palm,
geom = "point",
col.var = "black")
# Look at percentage of variation explained by the PCA axes
eigenvalues <- factoextra::get_eigenvalue(pca_palm)Another option to visually assess (and compare) trait combinations is
by using the funspace package, which relies on kernel
density estimation to calculate the actual boundaries of the trait space
and some indices of functional diversity.
library(funspace)
# The main function of the package, funspace(), accept as a main argument any
# object containing trait information, eg. a PCA like the one we just created.
palm_func_space <- funspace(PCs = c(1, 2),
x = pca_palm,
n_divisions = 300)
# The summary functions returns the"loadings", which are, in a very simplified way, the weights
# (i.e. coefficients) that show how much each original variable contributes to
# a specific Principal Component (PC, also sometimes called Dimensions, Dim.)
summary(palm_func_space)##
## Functional space based on a PCA with 4 dimensions
## Dimensions 1 and 2 are considered in analyses
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4
## log_height 0.743 0.557 0.051 0.114
## log_leaf_nb 0.269 0.083 -0.421 -0.213
## log_blade 0.785 -0.309 0.208 -0.214
## log_fruit_length 0.493 -0.393 -0.179 0.285
##
## Percentage of variance explained for each trait:
## Comp.1 Comp.2 Overall_explained
## log_height 55.27 31.05 86.33
## log_leaf_nb 7.22 0.68 7.91
## log_blade 61.66 9.54 71.20
## log_fruit_length 24.35 15.45 39.81
##
## ------------------------------------------------------------
##
## Functional diversity indicators:
##
## ---> For the global set of species:
##
## Functional richness (99.9% probability threshold) = 29.59
##
## Functional divergence = 0.55# The summary also tells us the percentage of variance of each trait that is explained
# by each principal component (Comp.1 and Comp.2 in this example) and across components (Overall_explained).
# We can then plot the functional space:
plot(palm_func_space,
quant.plot = TRUE,
pnt = TRUE,
arrows = T,
colors = c("white", "yellow", "orange", "red"))
References
Jonathan Bakker (2026), Applied Multivariate Statistics in R Copyright, University of Washington
Pictures of palms: http://idtools.org/id/palms/palmid/gallery.php
The world’s tallest monocotyledon:
tra$SpecName[which(tra$Climbing == 0)][which.max(tra$MaxStemHeight_m[which(tra$Climbing == 0)])]## [1] "Ceroxylon quindiuense"
Ceroxylon quindiuense