Pre Requirements: We have to create 3 separate CSV files with following data structure.

Note: 

# ------------------- 1. INSTALL & LOAD PACKAGES -------------------

install.packages("mFD")    # Install mFD (only run once)

library(mFD)               # Main package for Functional Diversity calculations

library(knitr)             # For displaying clean tables

library(dplyr)             # For data manipulation like select(), filter(), mutate()

library(tidyverse)         # Collection of data science packages (includes ggplot2, readr, etc.)

library(ggplot2)           # For plotting

library(patchwork)         # To combine multiple ggplot2 plots

# ------------------- 2. LOAD DATA -------------------

# Load species x traits table

raptor <- read.csv("C:/.../raptor.csv", row.names = 1)

# Ensures that species names are row names, required by mFD

if (is.null(rownames(raptor))) {

  stop("Species names are missing in the traits data.")

}

str(raptor)  # Shows the structure of the traits dataset

knitr::kable(head(raptor), caption = "Species x traits matrix")  # Shows first few rows

# Load species x assemblages (sites) matrix

raptorocc <- read.csv("C:/.../occh.csv", row.names = 1)

raptorocc <- as.matrix(raptorocc)  # Convert to matrix

storage.mode(raptorocc) <- "numeric"  # Force numerical mode (if it's 0s and 1s)

str(raptorocc)

knitr::kable(as.data.frame(raptorocc[1:3, 1:5]), caption = "Species x assemblages matrix")

# Load trait metadata (trait types: nominal, ordinal, etc.)

raptorcat <- read.csv("C:/.../category.csv", stringsAsFactors = FALSE)

str(raptorcat)

# Check expected columns

if (!all(c("trait_name", "trait_type") %in% colnames(raptorcat))) {

  stop("Trait category file must have 'trait_name' and 'trait_type' columns.")

}

# Check that all traits in raptor exist in raptorcat

if (!all(colnames(raptor) %in% raptorcat$trait_name)) {

  stop("Mismatch between traits in raptor and raptorcat.")

}

# ------------------- 3. CONVERT TRAIT COLUMNS -------------------

# Convert categorical variables into factor format

raptor$Habitat   <- as.factor(raptor$Habitat)

raptor$Niche <- as.factor(raptor$Niche)

raptor$Lifestyle <- as.factor(raptor$Lifestyle)

raptor$Active <- as.factor(raptor$Active)

str(raptor)

# ------------------- 4. TRAIT SUMMARY -------------------

raptorsumm <- mFD::sp.tr.summary(tr_cat = raptorcat, sp_tr = raptor, stop_if_NA = TRUE)

raptorsumm$tr_types  # Shows how each trait was classified (quantitative, ordinal, nominal)

raptorsumm$mod_list  # Details on trait transformation and coding

# ------------------- 5. OCCURRENCE MATRIX SUMMARY -------------------

raptor_summ <- mFD::asb.sp.summary(asb_sp_w = raptorocc)

raptor_occ <- raptor_summ$asb_sp_occ  # Binary presence/absence version

# ------------------- 6. FUNCTIONAL DISTANCE MATRIX -------------------

# This calculates pairwise distances between species based on their traits.

sp_dist_raptor <- mFD::funct.dist(

  sp_tr = raptor,        # Trait data

  tr_cat = raptorcat,    # Trait types

  metric = "gower",      # Distance metric used

  #  What is "Gower"?:

  # Gower distance can handle mixed data types:

  # - Quantitative (continuous)

  # - Ordinal (ordered categorical)

  # - Nominal (unordered categorical)

  # It standardizes variables so all contribute equally.

  scale_euclid = "scale_center",  # Center and scale numeric traits before Euclidean distances

  ordinal_var = "classic",    # Uses classic method for ordinal traits (rank then Euclidean)

  weight_type = "equal",      # All traits contribute equally to distance

  stop_if_NA = TRUE           # Stop if any missing data found

)

round(sp_dist_raptor, 3)  # View the resulting distance matrix (rounded)

# Save the distance matrix

write.csv(as.data.frame(as.matrix(sp_dist_raptor)),

          "C:/.../raptorsdistance.csv")

# ------------------- 7. FUNCTIONAL SPACE QUALITY -------------------

# Creates multi-dimensional trait space using Principal Coordinates Analysis (PCoA)

fspaces_raptor <- mFD::quality.fspaces(

  sp_dist = sp_dist_raptor,  # Functional distance matrix

  maxdim_pcoa = 10,          # Compute up to 10 dimensions

  deviation_weighting = "absolute",  # Use absolute deviation for space quality

  fdist_scaling = FALSE,     # Do not scale distances again

  fdendro = "average"        # Use average linkage for dendrogram

)

round(fspaces_raptor$quality_fspaces, 3)  # Quality for 2D to 10D spaces

# Plot quality for different spaces (tree, 2D, 3D, etc.)

mFD::quality.fspaces.plot(

  fspaces_quality = fspaces_raptor,

  quality_metric = "mad",  # Mean absolute deviation (lower = better)

  fspaces_plot = c("tree_average", "pcoa_2d", "pcoa_3d", "pcoa_4d")

)

# ------------------- 8. TRAIT–AXIS CORRELATION -------------------

sp_faxes_coord_raptor <- fspaces_raptor$details_fspaces$sp_pc_coord

# Optional: Reduce number of traits

raptor_subset <- raptor[, 1:10]

# Check correlation between original traits and PCoA axes (functional space)

raptor_faxes <- mFD::traits.faxes.cor(

  sp_tr = raptor_subset,

  sp_faxes_coord = sp_faxes_coord_raptor[, 1:4],  # PC1–PC4

  plot = TRUE

)

# Traits significantly related to any axes

raptor_faxes$tr_faxes_stat[which(raptor_faxes$tr_faxes_stat$p.value < 0.05), ]

# ------------------- 9. FUNCTIONAL SPACE PLOT -------------------

big_plot <- mFD::funct.space.plot(

  sp_faxes_coord = sp_faxes_coord_raptor,

  plot_ch = TRUE,          # Plot convex hulls (boundaries)

  plot_vertices = TRUE     # Show vertices (species at the edges)

)

big_plot$patchwork  # Displays the functional space plot

# ------------------- 10. FUNCTIONAL ALPHA DIVERSITY -------------------

# Alpha diversity = diversity within each assemblage

alpha_fd_indices_raptor <- mFD::alpha.fd.multidim(

  sp_faxes_coord = sp_faxes_coord_raptor[, 1:4],

  asb_sp_w = raptorocc,     # Assemblage × species matrix

  ind_vect = c("fdis", "fmpd", "fnnd", "feve", "fric",

               "fdiv", "fori", "fspe", "fide"),  # List of FD indices

  #  What do these indices mean?

  # fdis  = Functional dispersion (spread of species)

  # fmpd  = Mean pairwise distance

  # fnnd  = Mean nearest neighbor distance

  # feve  = Functional evenness

  # fric  = Functional richness (volume of convex hull)

  # fdiv  = Functional divergence

  # fori  = Functional originality

  # fspe  = Specialization

  # fide  = Identity

  scaling = TRUE,           # Standardize metrics

  details_returned = TRUE   # Return intermediate details

)

# Save FD index values

write.csv(alpha_fd_indices_raptor$functional_diversity_indices,

          "C:/.../Test functional diversity indices.csv")

# Save the details list (e.g., centroid coordinates)

write.csv(as.data.frame(alpha_fd_indices_raptor$details),

          "C:/.../2 detail list.csv")

# ------------------- 11. PLOT FUNCTIONAL ALPHA DIVERSITY -------------------

# Plot comparison between 2 sites (assemblages), e.g., "A" and "D"

plots_alpha <- mFD::alpha.multidim.plot(

  output_alpha_fd_multidim = alpha_fd_indices_raptor,

  plot_asb_nm = c("A", "D"),  # Site names

  ind_nm = c("fdis", "fide", "fnnd", "feve", "fric",

             "fdiv", "fori", "fspe")

)

# Save selected plots as high-resolution images

ggsave("C:/.../fric_plot.png", plot = plots_alpha$fric$patchwork, width = 10, height = 8, dpi = 600)

# Repeat for other indices (fdis, fdiv, etc.)

# ------------------- 12. FUNCTIONAL BETA DIVERSITY -------------------

# Beta diversity = dissimilarity between communities

beta_fd_indices_raptors <- mFD::beta.fd.multidim(

  sp_faxes_coord = sp_faxes_coord_raptor[, 1:4],

  asb_sp_occ = raptor_occ,   # Presence/absence

  beta_family = "Jaccard",   # Dissimilarity index

  #  What is "Jaccard"?

  # Jaccard index = a/b+c+a

  # Measures dissimilarity between two sets (communities)

  # Based on shared vs. unique species

  details_returned = TRUE

)

str(beta_fd_indices_raptors)  # Explore results

# ------------------- 13. PLOT FUNCTIONAL BETA DIVERSITY -------------------

# Visualize difference between two assemblages (e.g., "Mix" and "Open")

beta_plot_raptor <- mFD::beta.multidim.plot(

  output_beta_fd_multidim = beta_fd_indices_raptors,

  plot_asb_nm = c("Mix", "Open"),

  plot_sp_nm = c("Aegypius_monachus", "Clanga_clanga", "Gyps_himalayensis"),  # Optional focus species

  beta_family = "Jaccard",

  faxes = c("PC1", "PC2")

)

print(beta_plot_raptor)  # Show the beta diversity plot