Functional diversity with trait probability densities

Introduction

Species richness counts how many species live in a place. Functional diversity asks what those species are like: how tall they grow, whether they are trees or herbs, how big their seeds are. Two places can have the same number of species but very different traits, and a place can keep its number of species while its traits change.

phyloregion measures functional diversity with trait probability densities (TPD, Carmona et al. 2016). Each species is a small hill of probability placed at its traits in trait space. A grid cell’s trait spectrum is the sum of its species’ hills. From these spectra we get functional alpha diversity for each grid cell, functional beta diversity between grid cells, functional regions, and, for a whole group of species, the trait spectrum of threatened versus non-threatened species and its erosion under extinction scenarios (Carmona et al. 2021).

The functions work with any mix of continuous and categorical traits and use the same sparse-matrix approach as phylobeta, so they are fast for large community matrices.

1. Data

We use the africa data: a grid-cell × species matrix for woody plants of southern Africa, with traits for each species. We use one continuous trait, maximum height (log-transformed, because heights are skewed), and one categorical trait, growth habit.

library(phyloregion)
library(terra)

data(africa)
p <- vect(system.file("ex/sa.json", package = "phyloregion"))

tr <- data.frame(species = africa$trait$species,
                 height = log10(africa$trait$h_max),
                 habit = africa$trait$habit)
head(tr)
#>                            species     height habit
#> 1        Abutilon_angulatum_OM1934  0.6020600  herb
#> 2    Abutilon_sonneratianum_LTM034  0.3010300 shrub
#> 3       Acalypha_chirindica_OM2341  0.6989700 shrub
#> 4 Acalypha_glabrata_glabrata_OM441  0.7853298  tree
#> 5       Acalypha_sonderiana_OM2163 -1.2941363  tree
#> 6  Acokanthera_oblongifolia_OM2240  1.1760913 shrub

The first column must contain the species names. Every other column is used as a trait: numeric columns are treated as continuous and character or factor columns as categorical. Species with missing values are left out.

2. How the trait space is built

Continuous traits are standardized and become the axes of trait space (a ruler for one trait, a sheet for two, a box for three).
Each combination of categories gets its own trait space, so a tree’s density never spills into the trait space of herbs.
With more than three continuous traits or more than eight combinations of categories, all traits are summarized by two mixed principal component axes.
Each species’ hill is a normal kernel with a width chosen by bw.nrd0, cut so that it keeps 99% of its density (quantile = 0.99).

3. Functional alpha diversity

fa <- tpd_alpha(africa$comm, tr)
head(fa)
#>   grids  SR      FRic      FEve       FRed CWM_height share_herb share_liana
#> 1 v3635  80 12.210164 0.4877628 23.5452355  1.1423663 0.01250000  0.06250000
#> 2 v3636 114 12.698571 0.4697282 30.9120155  1.1453329 0.00877193  0.05263158
#> 3 v3637  74 11.884560 0.4741508 19.3920281  1.1239405 0.01351351  0.02702703
#> 4 v3638  56 11.721758 0.4831587 13.8707198  1.0943056 0.01785714  0.03571429
#> 5 v3639  44 12.047362 0.4968894 11.7807683  1.1115378 0.02272727  0.06818182
#> 6 v3640   8  6.023681 0.6876728  0.9854933  0.9423217 0.00000000  0.12500000
#>   share_shrub share_succulent share_tree
#> 1   0.2625000               0  0.6625000
#> 2   0.3157895               0  0.6228070
#> 3   0.3108108               0  0.6486486
#> 4   0.3214286               0  0.6250000
#> 5   0.2727273               0  0.6363636
#> 6   0.5000000               0  0.3750000

tpd_alpha returns one row per grid cell:

SR: species richness
FRic: functional richness, the volume of trait space occupied
FEve: functional evenness, how evenly the density is spread
FRed: functional redundancy, the average number of extra species sharing the same traits
CWM_height: community-weighted mean of each continuous trait
share_: share of each category

Map them with the grid-cell polygons:

m <- merge(p, fa, by = "grids")
par(mfrow = c(1, 3), mar = c(1, 1, 2, 3))
plot(m, "FRic", main = "FRic", axes = FALSE, border = NA,
     col = hcl.colors(20, "YlOrRd", rev = TRUE))
plot(m, "FEve", main = "FEve", axes = FALSE, border = NA,
     col = hcl.colors(20, "YlOrRd", rev = TRUE))
plot(m, "FRed", main = "FRed", axes = FALSE, border = NA,
     col = hcl.colors(20, "YlOrRd", rev = TRUE))

4. Functional beta diversity

tpd_beta compares the trait spectra of every pair of grid cells. The total dissimilarity is one minus their overlap. It is split into two parts (Carmona et al. 2016):

beta.unique: density in trait space used by only one of the two grid cells (like turnover)
beta.nested: different amounts of density in trait space used by both (like nestedness)

fb <- tpd_beta(africa$comm, tr)
sapply(fb, function(d) round(mean(d), 3))
#> beta.unique beta.nested  beta.total 
#>       0.025       0.272       0.297

tpd_beta_core returns the shared and unshared parts of each pair, in the same way as phylobeta_core.

5. Functional regions

Functional beta diversity can be used like any other beta diversity in phyloregion. Grid cells with very few species look completely different from everything else, so we keep grid cells with at least 10 species before clustering.

keep <- fa$grids[fa$SR >= 10]
fbk <- as.dist(as.matrix(fb$beta.total)[keep, keep])
fr <- phyloregion(fbk, k = 5, pol = p)
plot(fr)

6. Trait spectra of groups: threatened versus non-threatened

Instead of grid cells, we can compare groups of species, as Carmona et al. (2021) did for the global spectra of plants and vertebrates. tpd_groups turns a vector of group labels into a groups × species matrix that can be used in tpd_alpha and tpd_beta.

status <- ifelse(africa$IUCN$IUCN %in% c("CR", "EN", "VU"),
                 "threatened", "not_threatened")
names(status) <- africa$IUCN$Species
g <- tpd_groups(status)
rowSums(g)
#> not_threatened     threatened            all 
#>           1021             91           1112

tpd_alpha(g, tr)[, c("grids", "SR", "FRic", "FRed", "CWM_height")]
#>            grids  SR     FRic      FRed CWM_height
#> 1 not_threatened 737 23.82952 125.80932  0.8904009
#> 2     threatened  62 12.46710  11.21385  0.5593309
#> 3            all 799 24.46077 133.15548  0.8647109
round(as.matrix(tpd_beta(g, tr)$beta.total), 3)
#>                not_threatened threatened   all
#> not_threatened          0.000      0.399 0.031
#> threatened              0.399      0.000 0.368
#> all                     0.031      0.368 0.000

In this example, threatened species are on average shorter (a smaller CWM_height) and occupy a part of the trait space that differs from the spectrum of non-threatened species.

7. Extinction scenarios

tpd_extinction removes species at random according to their extinction probabilities, many times, and compares the trait spectrum before and after. As a null model, the same number of species is removed completely at random. The probabilities below are hypothetical values for each IUCN category; use probabilities suited to your question.

pr <- c(LC = 0.001, NT = 0.01, VU = 0.1, EN = 0.5, CR = 0.9)
prob <- setNames(pr[africa$IUCN$IUCN], africa$IUCN$Species)
set.seed(1)
ex <- tpd_extinction(tr, prob, nsim = 50)
ex$summary[ex$summary$quantity %in% c("lost", "FRic", "overlap"), ]
#>      scenario quantity       mean      lower      upper
#> 1  extinction     lost 26.2200000 22.0000000 31.0000000
#> 2  extinction     FRic  0.9994839  1.0000000  1.0000000
#> 5  extinction  overlap  0.9846655  0.9811280  0.9877163
#> 6      random     lost 26.2200000 22.0000000 31.0000000
#> 7      random     FRic  0.9935484  0.9483871  1.0000000
#> 10     random  overlap  0.9904504  0.9879078  0.9934542

If the overlap with the current spectrum is lower for extinctions than for random losses of the same number of species, extinctions are reorganizing the trait spectrum more than expected by chance.

Session info

sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: aarch64-apple-darwin23
#> Running under: macOS Tahoe 26.6.2
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/Los_Angeles
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] terra_1.9-46      Matrix_1.7-5      ape_5.8-1         phyloregion_2.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] magic_1.6-1        sass_0.4.10        generics_0.1.4     maptpx_1.9-7      
#>  [5] slam_0.1-55        rcdd_1.6-1         lattice_0.22-9     digest_0.6.39     
#>  [9] magrittr_2.0.5     evaluate_1.0.5     grid_4.6.0         RColorBrewer_1.1-3
#> [13] iterators_1.0.14   picante_1.8.2      betapart_1.6.1     fastmap_1.2.0     
#> [17] foreach_1.5.2      jsonlite_2.0.0     doSNOW_1.0.20      combinat_0.0-8    
#> [21] mgcv_1.9-4         itertools_0.1-3    permute_0.9-10     codetools_0.2-20  
#> [25] jquerylib_0.1.4    abind_1.4-8        cli_3.6.6          rlang_1.3.0       
#> [29] splines_4.6.0      cachem_1.1.0       yaml_2.3.12        vegan_2.7-3       
#> [33] otel_0.2.0         geometry_0.5.2     tools_4.6.0        predicts_0.2-2    
#> [37] parallel_4.6.0     minpack.lm_1.2-4   dplyr_1.2.1        colorspace_2.1-2  
#> [41] fastmatch_1.1-8    vctrs_0.7.3        R6_2.6.1           lifecycle_1.0.5   
#> [45] smoothr_1.3.0      MASS_7.3-65        cluster_2.1.8.2    pkgconfig_2.0.3   
#> [49] pillar_1.11.1      bslib_0.12.0       glue_1.8.1         phangorn_2.12.1   
#> [53] Rcpp_1.1.2         xfun_0.62          tibble_3.3.1       tidyselect_1.2.1  
#> [57] knitr_1.51         htmltools_0.5.9    nlme_3.1-169       snow_0.4-4        
#> [61] igraph_2.3.3       rmarkdown_2.32     clustMixType_0.4-2 compiler_4.6.0    
#> [65] quadprog_1.5-8