# Title: Functional and phylogenetic structure of forest bird assemblages # along an Afrotropical elevational gradient # Corresponding authors: Riccardo Pernice & David Hořák # This script includes codes used for estimating functional and # phylogenetic diversity metrics, as well as for statistical analyses. # This script can be run in Rstudio. # Note: (i) Standardized effect sizes of functional and phylogenetic diversity # metrics, may very slightly when running the script, because of randomization # procedures; (ii) the computation of the functional richness, after # removing iteratively species from the species-richer assemblage, can take # several minutes (Lines 1858-4685). # For the analysis, the following files need to be imported: # the Morpho_traits.xlsx (Line 55), # the Dietary_traits.xlsx (Line 275), the Strata_traits.xlsx (Line 480), # the Tactics_traits.xlsx (Line 665), the All_traits.xlsx (Line 848), # the Presence_absence_table.xlsx (Line 71-72), # the Phylo file (Line 1264-1265), Presence_absence_table_2.xlsx (Line 1279-1282), # Morpho_traits_2.xlsx file (Line 1290-1293). # These files have also been uploaded to this repository. # Sections: ## Lines 38-50: Required packages ## Lines 53-1076: Functional diversity metrics and relative null models ## Lines 1078-1258: Comparison trait space made in mFD package and trait space made in gawdis package ## Lines 1262-1359: Phylogenetic signal, Phylogenetic diversity metrics and relative null models ## Lines 1363-1581: Segmented regression of observed FD metrics ## Lines 1583-1611: Segmented regression of observed PD metrics ## Lines 1614-1828: Segmented regression of SES functional diversity metrics ## Lines 1831-1858: Segmented regression of SES phylogenetic diversity metrics ## Lines 1861-4162: Niche packing/expansion contribution to increments in species richness # Packages: library(mFD) library(ape) library(picante) library(RRphylo) library(phytools) library(gawdis) library(ggplot2) library(gridExtra) library(cowplot) library(segmented) library(scales) # Functional diversity based on morphological traits ------- ## Set the proper working directory and import the Morpho_traits.xlsx file. ## Then, convert the table to a matrix. morphotrait.matrix <- as.matrix(Morpho_traits[,-1]) rownames(morphotrait.matrix) <- Morpho_traits$Species ## Modify colnames. colnames(morphotrait.matrix) <- c("Wing", "Tail", "Culmen","Tarsus", "Weight") ## Log-transform all morphotraits and create a new dataframe. traits.scaled <- apply(morphotrait.matrix,2,log) rownames(traits.scaled) <- Morpho_traits$Species morphotraits <- as.data.frame(traits.scaled) ## Set the proper working directory and import the ## Presence_absence_table.xlsx file. Then, convert the table to a matrix. occurrence.matrix <- as.matrix(Presence_absence_table[,-1]) rownames(occurrence.matrix) <- Checklist_CMR_$Altitude ## Generate the functional space where functional diversity metrics are ## computed. fspace <- tr.cont.fspace(sp_tr = morphotraits, pca = T, nb_dim = 4, scaling = "scale_center", compute_corr = "pearson") ## Check the quality of the space. fspace$quality_metrics ## Check the percentage of variance explained and cumulative ## percentage for each axis. fspace$eigenvalues_percentage_var head(fspace$"sp_faxes_coord") ## Create a dist object containing species euclidean distances in the Fspace. dist_mat <- as.matrix(fspace$sp_dist_multidim$`4D`) dist_mat[1:4, 1:4] ## Check the correlation between traits. fspace$"tr_correl" ## Test correlation between functional axes and traits sp_faxes_coord_morphology <- fspace$sp_faxes_coord sp_tr_faxes <- traits.faxes.cor( sp_tr = morphotraits, sp_faxes_coord = sp_faxes_coord_morphology, plot = T) print(sp_tr_faxes) ## Plot the convex-hull of the species pool for the PCA axes. big_plot <- funct.space.plot(sp_faxes_coord = sp_faxes_coord_morphology, faxes = NULL, name_file = NULL, faxes_nm = NULL, range_faxes = c(NA, NA), color_bg = "grey95", color_pool = "darkgreen", fill_pool = "white", shape_pool = 21, size_pool = 1, plot_ch = TRUE, color_ch = "black", fill_ch = "white", alpha_ch = 0.5, plot_vertices = TRUE, color_vert = "blueviolet", fill_vert = "blueviolet", shape_vert = 23, size_vert = 1, plot_sp_nm = NULL , nm_size = 3, nm_color = "black", nm_fontface = "plain", check_input = TRUE) print(big_plot) ## Compute Functional Diversity metrics. alpha_fd_indices <- alpha.fd.multidim(sp_faxes_coord = sp_faxes_coord_morphology, asb_sp_w = occurrence.matrix, ind_vect = c("fric", "feve", "fnnd"), scaling = F) fd_ind_values_morphology <- alpha_fd_indices$functional_diversity_indices fd_ind_values_morphology ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional richness is then calculated ## across bird assemblages. null.FRic <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_morphology),length(rownames(sp_faxes_coord_morphology)), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fric", scaling = F) } null.fric.rep <- replicate(500,null.FRic(sp_faxes_coord_morphology)) null.fric.rep ## Create a new matrix where to store functional richness ## values of the null distribution. indices <-seq (1,999,2) newset <- NULL for (i in indices) { newset <- cbind (newset, null.fric.rep [[c(i,2)]]) } View(newset) ## Calculate the standardized effect size of the functional richness. average.null.fric. <- apply(newset,1,mean) average.null.fric. sd.null.fric. <- apply(newset,1,sd) sd.null.fric. ses.value.fric. <- (fd_ind_values_morphology$fric - average.null.fric.)/sd.null.fric. ses.value.fric. fric.rank <- apply(cbind(fd_ind_values_morphology$fric,newset),1,rank) [1,] fric.rank p.val.fric <- fric.rank/501 p.val.fric ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional nearest neighbor ## distance is then calculated across bird assemblages. null.FNnd <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_morphology),length(rownames(sp_faxes_coord_morphology)), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = order.czech.community, ind_vect = "fnnd", scaling = F) } null.fnnd.rep <- replicate(500,null.FNnd(sp_faxes_coord_morphology)) null.fnnd.rep ## Create a new matrix where to store functional nearest neighbor ## distance values of the null distribution. newset.fnnd <- NULL for (i in indices) { newset.fnnd <- cbind (newset.fnnd, null.fnnd.rep [[c(i,2)]]) } View(newset.fnnd) ## Calculate the standardized effect size of the functional nearest neighbor ## distance. average.null.fnnd. <- apply(newset.fnnd,1,mean) sd.null.fnnd. <- apply(newset.fnnd,1,sd) ses.value.fnnd. <- (fd_ind_values_morphology$fnnd - average.null.fnnd.)/sd.null.fnnd. fnnd.rank <- apply(cbind(fd_ind_values_morphology$fnnd,newset.fnnd),1,rank) [1,] p.val.fnnd <- fnnd.rank/501 p.val.fnnd ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional evenness is then ## calculated across bird assemblages. null.FEve <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_morphology),length(rownames(sp_faxes_coord_morphology)), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = order.czech.community, ind_vect = "feve", scaling = F) } null.feve.rep <- replicate(500,null.FEve(sp_faxes_coord_morphology)) null.feve.rep ## Create a new matrix where to store functional evenness values ## of the null distribution. newset.feve <- NULL for (i in indices) { newset.feve <- cbind (newset.feve, null.feve.rep [[c(i,2)]]) } View(newset.feve) ## Calculate the standardized effect size of the functional evenness. average.null.feve. <- apply(newset.feve,1,mean) sd.null.feve. <- apply(newset.feve,1,sd) ses.value.feve. <- (fd_ind_values_morphology$feve - average.null.feve.)/sd.null.feve. feve.rank <- apply(cbind(fd_ind_values_morphology$feve,newset.feve),1,rank) [1,] p.val.feve <- feve.rank/501 # Functional diversity based on dietary traits ------- ## Set the proper working directory and import the Dietary_traits.xlsx file. diet.matrix <- as.matrix(Dietary_traits[,-1]) rownames(diet.matrix) <- Dietary_traits$Species diet.traits <- as.data.frame(diet.matrix) ## Create a trait category dataset needed to compute distances between ## species based on functional traits in order to build the functional space ## in which indices will be computed. Trait_category_diet <- data.frame(trait_name = c("Invert", "Fruits", "Leaves", "Nectar", "Seeds", "Vert"), trait_type = c("F","F","F","F","F", "F"), fuzzy_name = c("Consume","Consume", "Consume","Consume", "Consume","Consume")) sp_dist_diet <- funct.dist(sp_tr = diet.traits, tr_cat = Trait_category_diet, metric = "gower", weight_type = "equal") ## In order to create the functional space in which functional diversity ## metrics are computed, we will perform a PCoA using the trait-based ## distances. fspaces_quality_diet <- quality.fspaces(sp_dist = sp_dist_diet, maxdim_pcoa = 10, fdist_scaling = F, fdendro = "average", deviation_weighting = "squared") # The space with the best quality has the lowest quality metric. round(fspaces_quality_diet$quality_fspaces, 3) # Illustrate the quality of the selected functional space. figure.fquality <- quality.fspaces.plot(fspaces_quality = fspaces_quality_diet, quality_metric = "rmsd", fspaces_plot = c("tree_average", "pcoa_1d", "pcoa_2d", "pcoa_3d", "pcoa_4d", "pcoa_5d"), gradient_deviation = c(neg = "darkblue", nul = "grey80", pos = "darkred"), gradient_deviation_quality = c(low = "yellow", high = "red"), x_lab = "Trait-based distance") figure.fquality ## Test correlation between functional axes and traits. ## Print only traits with significant effect. sp_faxes_coord_diet <- fspaces_quality_diet$details_fspaces$sp_pc_coord diet_tr_faxes <- traits.faxes.cor(sp_tr = diet.traits, sp_faxes_coord = sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")], plot = T) diet_tr_faxes$tr_faxes_stat[which(diet_tr_faxes$tr_faxes_stat$p.value < 0.05),] diet_tr_faxes$tr_faxes_plot diet_tr_faxes$tr_faxes_stat ## Calculate Functional diversity metrics. ## Check, first, if rownames of trait matrix and colnames of the community ## matrix are matching. setdiff(rownames(sp_faxes_coord_diet),colnames(occurrence.matrix)) rownames(sp_faxes_coord_diet)[134] <- "Mandingoa nitidula" alpha_fd_indices <- alpha.fd.multidim(sp_faxes_coord = sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")], asb_sp_w = occurrence.matrix, ind_vect = c("fric", "feve", "fnnd"), scaling = F) fd_ind_values_diet <- alpha_fd_indices$functional_diversity_indices fd_ind_values_diet ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional richness is then calculated ## across bird assemblages. null.FRic.diet <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")]) ,length(rownames(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fric", scaling = F) } null.fric.rep.diet <- replicate(500,null.FRic(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")])) null.fric.rep.diet ## Create a new matrix where to store functional richness ## values of the null distribution. newset.fric.diet <- NULL for (i in indices) { newset.fric.diet <- cbind (newset.fric.diet, null.fric.rep.diet [[c(i,2)]]) } View(newset.fric.diet) ## Calculate the standardized effect size of the functional richness. average.null.fric.diet <- apply(newset.fric.diet,1,mean) sd.null.fric.diet <- apply(newset.fric.diet,1,sd) ses.value.fric.diet <- (fd_ind_values_diet$fric - average.null.fric.diet)/sd.null.fric.diet fric.rank.diet <- apply(cbind(fd_ind_values_diet$fric,newset.fric.diet),1,rank) [1,] p.val.fric.diet <- fric.rank.diet/501 ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional nearest neighbor ## distance is then calculated across bird assemblages. null.FNnd.diet <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")]) ,length(rownames(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fnnd", scaling = F) } null.fnnd.rep.diet <- replicate(500,null.FNnd.diet(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")])) ## Create a new matrix where to store unctional nearest neighbor distance ## values of the null distribution. newset.fnnd.diet <- NULL for (i in indices) { newset.fnnd.diet <- cbind (newset.fnnd.diet, null.fnnd.rep.diet [[c(i,2)]]) } ## Calculate the standardized effect size of the functional nearest neighbor ## distance. average.null.fnnd.diet <- apply(newset.fnnd.diet,1,mean) sd.null.fnnd.diet <- apply(newset.fnnd.diet,1,sd) ses.value.fnnd.diet <- (fd_ind_values_diet$fnnd - average.null.fnnd.diet)/sd.null.fnnd.diet ses.value.fnnd.diet fnnd.rank.diet <- apply(cbind(fd_ind_values_diet$fnnd,newset.fnnd.diet),1,rank) [1,] p.val.fnnd.diet <- fnnd.rank.diet/501 p.val.fnnd.diet ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional evenness is then ## calculated across bird assemblages. null.FEve.diet <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")]) ,length(rownames(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "feve", scaling = F) } null.feve.rep.diet <- replicate(500,null.FEve.diet(sp_faxes_coord_diet[, c("PC1", "PC2", "PC3", "PC4", "PC5")])) null.feve.rep.diet ## Create a new matrix where to store functional evenness values of ## the null distribution. newset.feve.diet <- NULL for (i in indices) { newset.feve.diet <- cbind (newset.feve.diet, null.feve.rep.diet [[c(i,2)]]) } ## Calculate the standardized effect size of the functional evenness. average.null.feve.diet <- apply(newset.feve.diet,1,mean) sd.null.feve.diet <- apply(newset.feve.diet,1,sd) ses.value.feve.diet <- (fd_ind_values_diet$feve - average.null.feve.diet)/sd.null.feve.diet ses.value.feve.diet feve.rank.diet <- apply(cbind(fd_ind_values_diet$feve,newset.feve.diet),1,rank) [1,] p.val.feve.diet <- feve.rank.diet/501 p.val.feve.diet # Functional diversity based on traits related to vegetation strata usage ---- ## Set the proper working directory and import the Strata_traits.xlsx file. strata.matrix <- as.matrix(Strata_traits[,-1]) rownames(strata.matrix) <- Strata_traits$Species strata.traits <- as.data.frame(strata.matrix) ## Create a trait category dataset needed to compute distances between ## species based on functional traits in order to build the functional space ## in which indices will be computed. Trait_category_strata <- data.frame(trait_name = c("Ground", "Low", "Middle", "Canopy"), trait_type = c("F","F","F","F"), fuzzy_name = c("Consume","Consume", "Consume","Consume")) sp_dist_strata <- funct.dist(sp_tr = strata.traits, tr_cat = Trait_category_strata, metric = "gower", weight_type = "equal") round(sp_dist_strata, 3) ## In order to create the functional space in which functional diversity ## metrics are computed, we will perform a PCoA using the trait-based ## distances. fspaces_quality_strata <- quality.fspaces(sp_dist = sp_dist_strata, maxdim_pcoa = 10, fdist_scaling = F, fdendro = "average", deviation_weighting = "squared") round(fspaces_quality_strata$quality_fspaces, 3) # Illustrate the quality of the selected functional space. figure.fquality <- quality.fspaces.plot(fspaces_quality = fspaces_quality_strata, quality_metric = "rmsd", fspaces_plot = c("tree_average", "pcoa_1d", "pcoa_2d", "pcoa_3d"), gradient_deviation = c(neg = "darkblue", nul = "grey80", pos = "darkred"), gradient_deviation_quality = c(low = "yellow", high = "red"), x_lab = "Trait-based distance") figure.fquality ## Test correlation between functional axes and traits. ## Print only traits with significant effect. sp_faxes_coord_strata <- fspaces_quality_strata$details_fspaces$sp_pc_coord strata_tr_faxes <- traits.faxes.cor(sp_tr = strata.traits, sp_faxes_coord = sp_faxes_coord_strata[, c("PC1", "PC2","PC3")], plot = T) strata_tr_faxes$tr_faxes_stat[which(strata_tr_faxes$tr_faxes_stat$p.value < 0.05),] # Compute Functional diversity metrics. alpha_fd_indices <- alpha.fd.multidim(sp_faxes_coord = sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")], asb_sp_w = occurrence.matrix, ind_vect = c("fdis", "fric", "feve", "fnnd"), scaling = F) fd_ind_values_strata <- alpha_fd_indices$functional_diversity_indices fd_ind_values_strata ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional richness is then calculated ## across bird assemblages. null.FRic.strata <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")]) ,length(rownames(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fric", scaling = F) } null.fric.rep.strata <- replicate(500,null.FRic.strata(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")])) ## Create a new matrix where to store functional richness ## values of the null distribution. newset.fric.strata <- NULL for (i in indices) { newset.fric.strata <- cbind (newset.fric.strata, null.fric.rep.strata [[c(i,2)]]) } ## Calculate the standardized effect size of the functional richness. average.null.fric.strata <- apply(newset.fric.strata,1,mean) sd.null.fric.strata <- apply(newset.fric.strata,1,sd) ses.value.fric.strata <- (fd_ind_values_strata$fric - average.null.fric.strata)/sd.null.fric.strata ses.value.fric.strata fric.rank.strata <- apply(cbind(fd_ind_values_strata$fric, newset.fric.strata),1,rank) [1,] p.val.fric.strata <- fric.rank.strata/501 ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional nearest neighbor ## distance is then calculated across bird assemblages. null.FNnd.strata <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")]) ,length(rownames(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fnnd", scaling = F) } null.fnnd.rep.strata <- replicate(500,null.FNnd.strata(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")])) ## Create a new matrix where to store functional nearest neighbor distance ## values of the null distribution. newset.fnnd.strata <- NULL for (i in indices) { newset.fnnd.strata <- cbind (newset.fnnd.strata, null.fnnd.rep.strata [[c(i,2)]]) } ## Calculate the standardized effect size of the functional nearest neighbor ## distance. average.null.fnnd.strata <- apply(newset.fnnd.strata,1,mean) sd.null.fnnd.strata <- apply(newset.fnnd.strata,1,sd) ses.value.fnnd.strata <- (fd_ind_values_strata$fnnd - average.null.fnnd.strata)/sd.null.fnnd.strata ses.value.fnnd.strata fnnd.rank.strata <- apply(cbind(fd_ind_values_strata$fnnd,newset.fnnd.strata),1,rank) [1,] p.val.fnnd.strata <- fnnd.rank.strata/501 p.val.fnnd.strata ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional evenness ## is then calculated across bird assemblages. null.FEve.strata <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")]) ,length(rownames(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "feve", scaling = F) } null.feve.rep.strata <- replicate(500,null.FEve.strata(sp_faxes_coord_strata[, c("PC1", "PC2", "PC3")])) ## Create a new matrix where to store functional evenness ## values of the null distribution. newset.feve.strata <- NULL for (i in indices) { newset.feve.strata <- cbind (newset.feve.strata, null.feve.rep.strata [[c(i,2)]]) } ## Calculate the standardized effect size of the functional evenness. average.null.feve.strata <- apply(newset.feve.strata,1,mean) sd.null.feve.strata <- apply(newset.feve.strata,1,sd) ses.value.feve.strata <- (fd_ind_values_strata$feve - average.null.feve.strata)/sd.null.feve.strata feve.rank.strata <- apply(cbind(fd_ind_values_strata$feve,newset.feve.strata),1,rank) [1,] p.val.feve.strata <- feve.rank.strata/501 # Functional diversity based on traits related to foraging tactics ----------- ## Set the proper working directory and import the Tactics_traits.xlsx file. tactics.matrix <- as.matrix(Tactics_traits[,-1]) rownames(tactics.matrix) <- Tactics_traits$Species tactics.traits <- as.data.frame(tactics.matrix) ## Create a trait category dataset needed to compute distances between ## species based on functional traits in order to build the functional space ## in which indices will be computed. Trait_category_tactics <- data.frame(trait_name = c("Foliage gleaning", "Bark probing", "Flycatching", "Sit-and-wait", "Ground search", "Fruit taking", "Nectar taking"), trait_type = c("F","F","F","F","F", "F","F"), fuzzy_name = c("Consume","Consume", "Consume","Consume", "Consume","Consume", "Consume")) sp_dist_tactics <- funct.dist(sp_tr = tactics.traits, tr_cat = Trait_category_tactics, metric = "gower", weight_type = "equal") round(sp_dist_tactics, 3) ## In order to create the functional space in which functional diversity ## metrics are computed, we will perform a PCoA using the trait-based ## distances. fspaces_quality_tactics <- quality.fspaces(sp_dist = sp_dist_tactics, maxdim_pcoa = 10, fdist_scaling = F, fdendro = "average", deviation_weighting = "squared") round(fspaces_quality_tactics$quality_fspaces, 3) ## Test correlation between functional axes and traits ## Print only traits with significant effect. sp_faxes_coord_tactics <- fspaces_quality_tactics$details_fspaces$sp_pc_coord tactics_tr_faxes <- traits.faxes.cor(sp_tr = tactics.traits, sp_faxes_coord = sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")], plot = T) tactics_tr_faxes$tr_faxes_stat[which(tactics_tr_faxes$tr_faxes_stat$p.value < 0.05),] tactics_tr_faxes$tr_faxes_plot ## Compute functional diversity metrics. alpha_fd_indices <- alpha.fd.multidim(sp_faxes_coord = sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")], asb_sp_w = occurrence.matrix, ind_vect = c("fdis", "fric", "feve", "fnnd"), scaling = F) fd_ind_values_tactics <- alpha_fd_indices$functional_diversity_indices fd_ind_values_tactics ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional richness is then calculated ## across bird assemblages. null.FRic.tactics <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")]) ,length(rownames(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fric", scaling = F) } null.fric.rep.tactics <- replicate(500,null.FRic.tactics(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")])) ## Create a new matrix where to store functional richness ## values of the null distribution. newset.fric.tactics <- NULL for (i in indices) { newset.fric.tactics <- cbind (newset.fric.tactics, null.fric.rep.tactics [[c(i,2)]]) } ## Calculate the standardized effect size of the functional richness. average.null.fric.tactics <- apply(newset.fric.tactics,1,mean) sd.null.fric.tactics <- apply(newset.fric.tactics,1,sd) ses.value.fric.tactics <- (fd_ind_values_tactics$fric - average.null.fric.tactics)/sd.null.fric.tactics fric.rank.tactics <- apply(cbind(fd_ind_values_tactics$fric,newset.fric.tactics),1,rank) [1,] p.val.fric.tactics <- fric.rank.tactics/501 ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional nearest neighbor ## distance is then calculated across bird assemblages. null.FNnd.tactics <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")]) ,length(rownames(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fnnd", scaling = F) } null.fnnd.rep.tactics <- replicate(500,null.FNnd.tactics(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")])) null.fnnd.rep ## Create a new matrix where to store functional nearest neighbor distance ## values of the null distribution. newset.fnnd.tactics <- NULL for (i in indices) { newset.fnnd.tactics <- cbind (newset.fnnd.tactics, null.fnnd.rep.tactics [[c(i,2)]]) } ## Calculate the standardized effect size of the functional nearest neighbor ## distance. average.null.fnnd.tactics <- apply(newset.fnnd.tactics,1,mean) sd.null.fnnd.tactics <- apply(newset.fnnd.tactics,1,sd) ses.value.fnnd.tactics <- (fd_ind_values_tactics$fnnd - average.null.fnnd.tactics)/sd.null.fnnd.tactics ses.value.fnnd.tactics fnnd.rank.tactics <- apply(cbind(fd_ind_values_tactics$fnnd,newset.fnnd.tactics),1,rank) [1,] p.val.fnnd.tactics <- fnnd.rank.tactics/501 p.val.fnnd.tactics ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional evenness ## is then calculated across bird assemblages. null.FEve.tactics <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")]) ,length(rownames(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "feve", scaling = F) } null.feve.rep.tactics <- replicate(500,null.FEve.tactics(sp_faxes_coord_tactics[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")])) ## Create a new matrix where to store functional evenness ## values of the null distribution. newset.feve.tactics <- NULL for (i in indices) { newset.feve.tactics <- cbind (newset.feve.tactics, null.feve.rep.tactics [[c(i,2)]]) } ## Calculate the standardized effect size of the functional evenness. average.null.feve.tactics <- apply(newset.feve.tactics,1,mean) sd.null.feve.tactics <- apply(newset.feve.tactics,1,sd) ses.value.feve.tactics <- (fd_ind_values_tactics$feve - average.null.feve.tactics)/sd.null.feve.tactics print(ses.value.feve.tactics) feve.rank.tactics <- apply(cbind(fd_ind_values_tactics$feve,newset.feve.tactics),1,rank) [1,] p.val.feve.tactics <- feve.rank.tactics/501 # Functional diversity based on all traits -------------------------------- ## Set the proper working directory and import the All_traits.xlsx file. all_trait_matrix <- as.matrix(All_traits[,-1]) rownames(all_trait_matrix) <- All_traits$Species all_traits <- as.data.frame(all_trait_matrix) ## Modify columns' name of the morphological traits. new_names <- c("Wing", "Tail", "Culmen","Tarsus", "Weight") colnames(all_traits)[(ncol(all_traits)-4):ncol(all_traits)] <- new_names ## Log-transform all morphological traits. all_traits[,18:22] <- log10(all_traits[,18:22]) ## Create a trait category dataset needed to compute distances between ## species based on functional traits in order to build the functional space ## in which indices will be computed. Trait_category_all_traits <- data.frame(trait_name = c("Invert", "Fruits", "Leaves", "Nectar", "Seeds", "Vert", "Ground", "Low", "Middle", "Canopy", "Foliage gleaning", "Bark probing", "Flycatching", "Sit-and-wait", "Ground search", "Fruit taking", "Nectar taking", "Wing","Tail","Culmen","Tarsus", "Weight"), trait_type = c("F","F","F","F","F", "F","F", "F","F", "F","F","F","F", "F","F","F","F","Q", "Q","Q","Q","Q"), fuzzy_name = c("Diet","Diet", "Diet","Diet", "Diet","Diet", "Strata","Strata", "Strata","Strata", "Tactic","Tactic", "Tactic","Tactic", "Tactic","Tactic", "Tactic","Morph", "Morph","Morph", "Morph","Morph")) sp_dist_all_traits <- funct.dist(sp_tr = all_traits, tr_cat = Trait_category_all_traits, metric = "gower", scale_euclid = "scale_center", weight_type = "equal") round(sp_dist_all_traits, 3) ## In order to create the functional space in which functional diversity ## metrics are computed, we will perform a PCoA using the trait-based ## distances. fspaces_quality_all_traits <- quality.fspaces(sp_dist = sp_dist_all_traits, maxdim_pcoa = 21, fdist_scaling = F, fdendro = "average", deviation_weighting = "squared") round(fspaces_quality_all_traits$quality_fspaces, 3) ## Test correlation between functional axes and traits ## Print only traits with significant effect. sp_faxes_coord_all_traits <- fspaces_quality_all_traits$details_fspaces$sp_pc_coord all_traits_tr_faxes <- traits.faxes.cor(sp_tr = all_traits[,1:9], sp_faxes_coord = sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")], plot = T) all_traits_tr_faxes$tr_faxes_stat[which(all_traits_tr_faxes$tr_faxes_stat$p.value < 0.05),] ## Calculate functional diversity metrics. alpha_fd_indices <- alpha.fd.multidim(sp_faxes_coord = sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")], asb_sp_w = occurrence.matrix, ind_vect = c("fdis", "fric", "feve", "fnnd"), scaling = F) fd_ind_values_all_traits <- alpha_fd_indices$functional_diversity_indices fd_ind_values_all_traits ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional richness is then calculated ## across bird assemblages. null.FRic.all.traits <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")]) ,length(rownames(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fric", scaling = F) } null.fric.rep.all.traits <- replicate(500,null.FRic.all.traits(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")])) null.fric.rep ## Create a new matrix where to store functional richness ## values of the null distribution. newset.fric.all.traits <- NULL for (i in indices) { newset.fric.all.traits <- cbind (newset.fric.all.traits, null.fric.rep.all.traits [[c(i,2)]]) } ## Calculate the standardized effect size of the functional richness. average.null.fric.all.traits <- apply(newset.fric.all.traits,1,mean) sd.null.fric.all.traits <- apply(newset.fric.all.traits,1,sd) ses.value.fric.all.traits <- (fd_ind_values_all_traits$fric - average.null.fric.all.traits)/sd.null.fric.all.traits ses.value.fric.all.traits fric.rank.all.traits <- apply(cbind(fd_ind_values_all_traits$fric, newset.fric.all.traits),1,rank) [1,] p.val.fric.all.traits <- fric.rank.all.traits/501 p.val.fric.all.traits ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional nearest neighbor ## distance is then calculated across bird assemblages. null.FNnd.all.traits <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")]) ,length(rownames(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "fnnd", scaling = F) } null.fnnd.rep.all.traits <- replicate(500,null.FNnd.all.traits(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")])) ## Create a new matrix where to store functional nearest neighbor distance ## values of the null distribution. newset.fnnd.all.traits <- NULL for (i in indices) { newset.fnnd.all.traits <- cbind (newset.fnnd.all.traits, null.fnnd.rep.all.traits [[c(i,2)]]) } ## Calculate the standardized effect size of the functional ## nearest neighbor distance. average.null.fnnd.all.traits <- apply(newset.fnnd.all.traits,1,mean) sd.null.fnnd.all.traits <- apply(newset.fnnd.all.traits,1,sd) ses.value.fnnd.all.traits <- (fd_ind_values_all_traits$fnnd - average.null.fnnd.all.traits)/sd.null.fnnd.all.traits ses.value.fnnd.all.traits fnnd.rank.all.traits <- apply(cbind(fd_ind_values_all_traits$fnnd, newset.fnnd.all.traits),1,rank) [1,] p.val.fnnd.all.traits <- fnnd.rank.all.traits/501 p.val.fnnd.all.traits ## Generate a null model where (i) species names are shuffled across ## PCA scores matrix 500 times, and (ii) functional evenness ## is then calculated across bird assemblages. null.FEve.all.traits <- function(x){ y<- x rownames(y) <- sample(rownames(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")]) ,length(rownames(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")])), replace=F) z <- y[order(rownames(y)),] alpha.fd.multidim(sp_faxes_coord = z, asb_sp_w = occurrence.matrix, ind_vect = "feve", scaling = F) } null.feve.rep.all.traits <- replicate(500,null.FEve.all.traits(sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")])) ## Create a new matrix where to store functional evenness ## values of the null distribution. newset.feve.all.traits <- NULL for (i in indices) { newset.feve.all.traits <- cbind (newset.feve.all.traits, null.feve.rep.all.traits [[c(i,2)]]) } ## Calculate the standardized effect size of the functional evenness. average.null.feve.all.traits <- apply(newset.feve.all.traits,1,mean) sd.null.feve.all.traits <- apply(newset.feve.all.traits,1,sd) ses.value.feve.all.traits <- (fd_ind_values_all_traits$feve - average.null.feve.all.traits)/sd.null.feve.all.traits ses.value.feve.all.traits feve.rank.all.traits <- apply(cbind(fd_ind_values_all_traits$feve,newset.feve.all.traits),1,rank) [1,] p.val.feve.all.traits <- feve.rank.all.traits/501 p.val.feve.all.traits # Assessing whether the mFD package provides an equal contribution of each group of traits into the trait space with its func.dist function ----------- # For this purpose, we test the correlation between functional diversity # metrics, based on all types of traits, using the trait space generated # from the distance matrix with the mFD package and the functional diversity # metrics, based on all types of traits, using the trait space generated # from the distance matrix with the gawdis package. ## Create a distance matrix using the gawdis package. Specify which sets of ## traits are fuzzy. sp_dist_all_traits_2 <- gawdis(all_traits, w.type = "analytic", groups = c(1,1,1,1,1,1, 2,2,2,2,3,3,3,3,3,3,3, 4,4,4,4,4), fuzzy = c(1,2,3)) ## Create the trait space in which functional diversity ## metrics are computed by performing a PCoA using the distnace ## matrix created in gawdis. fspaces_quality_all_traits_2 <- quality.fspaces(sp_dist = sp_dist_all_traits_2, maxdim_pcoa = 21, fdist_scaling = F, fdendro = "average", deviation_weighting = "squared") sp_faxes_coord_all_traits_2 <- fspaces_quality_all_traits_2$details_fspaces$sp_pc_coord round(fspaces_quality_all_traits_2$quality_fspaces, 3) # The space with the best quality is the 7D space ## Computing functional diversity metrics using the trait space generated from ## the distance matrix coming from the gawdis package. alpha_fd_indices_2 <- alpha.fd.multidim(sp_faxes_coord = sp_faxes_coord_all_traits_2[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7")], asb_sp_w = occurrence.matrix, ind_vect = c("fdis", "fric", "feve", "fnnd"), scaling = F) fd_ind_values_all_traits_2 <- alpha_fd_indices_2$functional_diversity_indices fd_ind_values_all_traits_2 ## Functional diversity metrics calculated with mFD package (for the whole ## procedure, see lines 835-933). alpha_fd_indices <- alpha.fd.multidim(sp_faxes_coord = sp_faxes_coord_all_traits[, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")], asb_sp_w = occurrence.matrix, ind_vect = c("fdis", "fric", "feve", "fnnd"), scaling = F) fd_ind_values_all_traits <- alpha_fd_indices$functional_diversity_indices fd_ind_values_all_traits ## Test the correlation between the functional diversity metrics ## obtained from the two packages # dissimilarity methods. ## Functional richness. corr_fric <- cor(fd_ind_values_all_traits$fric, fd_ind_values_all_traits_2$fric) ## Functional nearest neighbor distance. corr_fnnd <- cor(fd_ind_values_all_traits$fnnd, fd_ind_values_all_traits_2$fnnd) ## Functional evenness. corr_feve <- cor(fd_ind_values_all_traits$feve, fd_ind_values_all_traits_2$feve) ## Scatter plots and final image (Fig.S2). df_cor <- data.frame(FRic_mFD = fd_ind_values_all_traits$fric, FRic_gawdis = fd_ind_values_all_traits_2$fric, FNND_mFD = fd_ind_values_all_traits$fnnd, FNND_gawdis = fd_ind_values_all_traits_2$fnnd, FEve_mFD = fd_ind_values_all_traits$feve, FEve_gawdis = fd_ind_values_all_traits_2$feve, Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2350)) p_fric <- ggplot(data = df_cor, aes(x = FRic_mFD, y = FRic_gawdis)) + geom_point(color = "darkred", alpha = 0.6, aes(size = factor(Elevation))) + labs(x = "FRic from mFD package", y = "FRic from gawdis package", size = "Elevation m (a.s.l.)") + theme( axis.text = element_text(color = "black"), axis.title = element_text(face = "bold", color = "black"), panel.background = element_rect(fill = "white"), panel.border = element_blank(), axis.line = element_line(), legend.position = "none") + annotate("text", x = min(df_cor$FRic_mFD), y = max(df_cor$FRic_gawdis), label = paste("r = 0.99"), hjust = -0.5, vjust = 1, color = "darkgreen", size = 5, fontface = "italic")+ annotate("text", x = -Inf, y = Inf, label = "a", hjust = -0.5, vjust = 1, size = 5, fontface = "bold") p_fnnd <- ggplot(data = df_cor, aes(x = FNND_mFD, y = FNND_gawdis)) + geom_point(color = "darkred", alpha = 0.6, aes(size = factor(Elevation))) + labs(x = "FNND from mFD package", y = "FNND from gawdis package", size = "Elevation m (a.s.l.)") + theme( axis.text = element_text(color = "black"), axis.title = element_text(face = "bold", color = "black"), panel.background = element_rect(fill = "white"), panel.border = element_blank(), axis.line = element_line(), legend.position = "none") + annotate("text", x = min(df_cor$FNND_mFD), y = max(df_cor$FNND_gawdis), label = paste("r = 0.98"), hjust = -0.5, vjust = 1, color = "darkgreen", size = 5, fontface = "italic") + annotate("text", x = -Inf, y = Inf, label = "c", hjust = -0.5, vjust = 1, size = 5, fontface = "bold") p_feve <- ggplot(data = df_cor, aes(x = FEve_mFD, y = FEve_gawdis)) + geom_point(color = "darkred", alpha = 0.6, aes(size = factor(Elevation))) + labs(x = "FEve from mFD package", y = "FEve from gawdis package", size = "Elevation m (a.s.l.)") + theme( axis.text = element_text(color = "black"), axis.title = element_text(face = "bold", color = "black"), panel.background = element_rect(fill = "white"), panel.border = element_blank(), axis.line = element_line(), legend.position = "none") + annotate("text", x = min(df_cor$FEve_mFD), y = max(df_cor$FEve_gawdis), label = paste("r = 0.95"), hjust = -0.5, vjust = 1, color = "darkgreen", size = 5, fontface = "italic") + annotate("text", x = -Inf, y = Inf, label = "b", hjust = -0.5, vjust = 1, size = 5, fontface = "bold") ## Create the legend. legend <- get_legend( p_fnnd + guides(size = guide_legend(title = "Elevation m (a.s.l.)", nrow = 1)) + theme(legend.position = "bottom") ) ## Arrange all plots in a grid. plot_grid <- plot_grid( plot_grid(p_fric, p_feve, ncol = 2), plot_grid(NULL, p_fnnd, NULL, ncol = 3, rel_widths = c(1, 2, 1)), ncol = 1, rel_heights = c(1, 1) ) ## Combine the plot grid with the legend. final_plot <- plot_grid( plot_grid, legend, ncol = 1, rel_heights = c(4, 0.5) ) print(final_plot) # Phylogenetic diversity analyses -------------------------------------- ## Import the consensus tree generated in MESQUITE using 10000 trees created ## in Birdtree.org. my.nexus.phylo <- read.nexus("Phylo") my.nexus.phylo ## There is a little polytomy in the consensus tree, but it can be ## resolved randomly. This process can cause a slight mismatch in results ## of phylogenetic diversity metrics. Consensus_tree.1 <- fix.poly(my.nexus.phylo, type = c("resolve")) is.binary(Consensus_tree.1) ## Import the Presence_absence_table_2 for the analyses. ## Note: in this table, species names are changed to match the taxonomy ## in Birdtree.org. community.matrix <- as.matrix(Presence_absence_table_2[,-1]) rownames(community.matrix) <- Presence_absence_table_2$Altitude ## Assessing the power of phylogenetic signal through Blomberg's K for ## each morphological trait, in the following order: wing, tail, culmen, ## tarsus, and body mass. ## For this analysis we need to import the Morpho_traits_2.xlsx file. ## Note: in this table, the taxonomy matches that of the ## Presence_absence_table_2.xlsx file, which is essential ## for carrying out the analysis. traits <- as.matrix(Morpho_traits_2[,-1]) rownames(traits) <- Morpho_traits_2$Species traits ## Create a phylogenetic distance matrix. p.dist.mat <- cophenetic(Consensus_tree.1) ## Traits must be in the same order as the names in ## the phylogenetic distance matrix. traits <- traits[row.names(p.dist.mat),] traits test.wing <- phylosig(Consensus_tree.1, traits[Consensus_tree.1$tip.label, 1], method = "K", test = T, nsim = 1000) test.tail <- phylosig(Consensus_tree.1, traits[Consensus_tree.1$tip.label, 2], method = "K", test = T, nsim = 1000) test.culmen <- phylosig(Consensus_tree.1, traits[Consensus_tree.1$tip.label, 3], method = "K", test = T, nsim = 1000) test.tarsus <- phylosig(Consensus_tree.1,traits[Consensus_tree.1$tip.label, 4], method = "K", test = T, nsim = 1000) test.mass <- phylosig(Consensus_tree.1, traits[Consensus_tree.1$tip.label, 5], method = "K", test = T, nsim = 1000) ## Calculate the phylogenetic mean pairwise distance metric. pairwise <- mpd(community.matrix, cophenetic(Consensus_tree.1), abundance.weighted = F) print(pairwise) ## Randomize tip labels on the phylogenetic tree for the mpd and ## calculate SES and p-value. Because of randomization procedure, ## results may differ slightly. ses.mpd <- ses.mpd(community.matrix,cophenetic(Consensus_tree.1), null.model = "taxa.labels", abundance.weighted = F, runs = 500, iterations = 501) ## Calculate the phylogenetic mean nearest taxonomic distance metric. nearest <- mntd(community.matrix, cophenetic(Consensus_tree.1), abundance.weighted = F) print(nearest) ## Randomize tip labels on the phylogenetic tree for the mntd and ## calculate SES and p-value. Because of randomization procedure, ## results may differ slightly. ses.mntd <- ses.mntd(community.matrix,cophenetic(Consensus_tree.1), null.model = "taxa.labels", abundance.weighted = F, runs = 500, iterations = 501) # Segmented regression analysis of functional diversity based on morphology. ---------- ## Observed functional richness. ## Create a dataframe for the analysis. df.obs.fric <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), Fric = fd_ind_values_morphology$fric) summary(lm.FRic.morph<-lm(Fric~Elevation, data= df.obs.fric)) ## Now, we can test if there is a significant breakpoint in ## the model using a Davies test. davies.test(lm.FRic.morph,seg.Z =~Elevation, k=8) ## A small p-value of the davies test indicates a high probability ## that there is a difference between slopes of a segmented model. ## The function segmented() implements a bootstrap restarting algorithm to ## find the breakpoint in the linear model that with result in the largest # difference in the slopes of the resulting two models. summary(seg.FRic.morph<-segmented(lm.FRic.morph)) slope(seg.FRic.morph) ## Repeate the same procedure for the other metrics. ## Observed functional nearest neighbor distance. ## Create a dataframe for the analysis. df.obs.fnnd <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FNND = fd_ind_values_morphology$fnnd) summary(lm.FNND.morph<-lm(FNND~Elevation, data=df.obs.fnnd)) davies.test(lm.FNND.morph,seg.Z =~Elevation, k=8) summary(seg.FNND.morph<-segmented(lm.FNND.morph)) slope(seg.FNND.morph) ## Observed functional evenness. ## Create a dataframe for the analysis. df.obs.feve <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FEve = fd_ind_values_morphology$feve) summary(lm.FEve.morph<-lm(FEve~Elevation, data= df.obs.feve)) davies.test(lm.FEve.morph,seg.Z =~Elevation, k=8) summary(seg.FEve.morph<-segmented(lm.FEve.morph)) slope(seg.FEve.morph) # Segmented regression analysis of functional diversity based on diet. ---------- ## Observed functional richness. df.obs.fric.diet <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), Fric = fd_ind_values_diet$fric) summary(lm.FRic.diet<-lm(Fric~Elevation, data= df.obs.fric.diet)) ## Davies test. davies.test(lm.FRic.diet, seg.Z =~Elevation, k=8) summary(seg.FRic.diet<-segmented(lm.FRic.diet)) slope(seg.FRic.diet) ## Observed functional nearest neighbor distance. ## Create a dataframe for the analysis. df.obs.fnnd.diet <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FNND = fd_ind_values_diet$fnnd) summary(lm.FNND.diet<-lm(FNND~Elevation, data=df.obs.fnnd.diet)) davies.test(lm.FNND.diet,seg.Z =~Elevation, k=8) summary(seg.FNND.diet<-segmented(lm.FNND.diet)) slope(seg.FNND.diet) ## Observed functional evenness. ## Create a dataframe for the analysis. df.obs.feve.diet <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FEve = fd_ind_values_diet$feve) summary(lm.FEve.diet<-lm(FEve~Elevation, data= df.obs.feve.diet)) davies.test(lm.FEve.diet,seg.Z =~Elevation, k=8) summary(seg.FEve.diet<-segmented(lm.FEve.diet)) slope(seg.FEve.diet) # Segmented regression analysis of functional diversity based on traits related to vegeation strata. ---------- ## Observed functional richness. df.obs.fric.strata <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), Fric = fd_ind_values_strata$fric) summary(lm.FRic.strata<-lm(Fric~Elevation, data= df.obs.fric.strata)) davies.test(lm.FRic.strata, seg.Z =~Elevation, k=8) summary(seg.FRic.strata<-segmented(lm.FRic.strata)) slope(seg.FRic.strata) ## Observed functional nearest neighbor distance. df.obs.fnnd.strata <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FNND = fd_ind_values_strata$fnnd) summary(lm.FNND.strata<-lm(FNND~Elevation, data=df.obs.fnnd.strata)) davies.test(lm.FNND.strata,seg.Z =~Elevation, k=8) summary(seg.FNND.strata<-segmented(lm.FNND.strata)) slope(seg.FNND.strata) ## Observed functional evenness. df.obs.feve.strata <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FEve = fd_ind_values_strata$feve) summary(lm.FEve.strata<-lm(FEve~Elevation, data=df.obs.feve.strata)) davies.test(lm.FEve.strata,seg.Z =~Elevation, k=8) summary(seg.FEve.strata<-segmented(lm.FEve.strata)) slope(seg.FEve.strata) # Segmented regression analysis of functional diversity based on foraging modes ---------- ## Observed functional richness. df.obs.fric.tactics <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), Fric = fd_ind_values_tactics$fric) summary(lm.FRic.tactics<-lm(Fric ~ Elevation, data= df.obs.fric.tactics)) davies.test(lm.FRic.tactics,seg.Z =~Elevation, k=8) summary(seg.FRic.tactics<-segmented(lm.FRic.tactics)) slope(seg.FRic.tactics) ## Observed functional nearest neighbor distance. df.obs.fnnd.tactics <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FNND = fd_ind_values_tactics$fnnd) summary(lm.FNND.tactics<-lm(FNND~Elevation, data=df.obs.fnnd.tactics)) davies.test(lm.FNND.tactics,seg.Z =~Elevation, k=8) summary(seg.FNND.tactics<-segmented(lm.FNND.tactics)) slope(seg.FNND.tactics) ## Observed functional evenness. df.obs.feve.tactics <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FEve = fd_ind_values_tactics$feve) summary(lm.FEve.tactics<-lm(FEve~Elevation, data= df.obs.feve.tactics)) davies.test(lm.FEve.tactics,seg.Z =~Elevation, k=8) summary(seg.FEve.tactics<-segmented(lm.FEve.tactics)) slope(seg.FEve.tactics) # Segmented regression analysis of functional diversity based on all traits ---------- ## Observed functional richness. df.obs.fric.all.traits <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FRic = fd_ind_values_all_traits$fric) summary(lm.FRic.combined<-lm(FRic~Elevation, data= df.obs.fric.all.traits)) davies.test(lm.FRic.combined,seg.Z =~Elevation, k=8) summary(seg.FRic.combined<-segmented(lm.FRic.combined)) slope(seg.FRic.combined) ## Observed functional nearest neighbor distance. df.obs.fnnd.all.traits <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FNND = fd_ind_values_all_traits$fnnd) summary(lm.FNND.combined<-lm(FNND~Elevation, data= df.obs.fnnd.all.traits)) davies.test(lm.FNND.combined,seg.Z =~Elevation, k=8) summary(seg.FNND.combined<-segmented(lm.FNND.combined)) slope(seg.FNND.combined) ## Observed functional evenness. df.obs.feve.all.traits <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), FEve = fd_ind_values_all_traits$feve) summary(lm.FEve.combined<-lm(FEve~Elevation, data= df.obs.feve.all.traits)) davies.test(lm.FEve.combined,seg.Z =~Elevation, k=8) summary(seg.FEve.combined<-segmented(lm.FEve.combined)) slope(seg.FEve.combined) # Segmented regression analysis of phylogenetic diversity metrics ---------- ## Observed mean pairwise distance. df.obs.mpd <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), MPD = ses.mpd$mpd.obs) summary(lm.MPD<-lm(MPD~Elevation, data= df.obs.mpd)) davies.test(lm.MPD,seg.Z =~Elevation, k=8) summary(seg.MPD<-segmented(lm.MPD)) slope(seg.MPD) ## Observed mean nearest taxonomic distance. df.obs.mntd <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), MNTD = ses.mntd$mntd.obs) summary(lm.MNTD<-lm(MNTD~Elevation, data= df.obs.mntd)) davies.test(lm.MNTD,seg.Z =~Elevation, k=8) summary(seg.MNTD<-segmented(lm.MNTD)) slope(seg.MNTD) # Segmented regression analysis of SES functional diversity metrics based on morphology ---------- ## SES.FRic. df.ses.fric <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES.Fric = ses.value.fric.) summary(lm.SES.Fric<-lm(SES.Fric~Elevation, data= df.ses.fric)) davies.test(lm.SES.Fric,seg.Z =~Elevation, k=8) # The small p-value of the davies test indicates a high probability # that there is a difference between slopes of a segmented model. # In other words, it indicates that two seperate linear models might # better describe this ageing pattern than one. summary(seg.SES.Fric<-segmented(lm.SES.Fric)) slope(seg.SES.Fric) ## SES.FNND. df.ses.fnnd <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES.Fnnd = ses.value.fnnd.) summary(lm.SES.Fnnd<-lm(SES.Fnnd~Elevation, data= df.ses.fnnd)) davies.test(lm.SES.Fnnd,seg.Z =~Elevation, k=8) summary(seg.SES.Fnnd<-segmented(lm.SES.Fnnd)) slope(seg.SES.Fnnd) ## SES.FEve. df.ses.feve <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES.FEve = ses.value.feve.) summary(lm.SES.FEve<-lm(SES.FEve~Elevation, data= df.ses.feve)) davies.test(lm.SES.FEve,seg.Z =~Elevation, k=8) summary(seg.SES.FEve<-segmented(lm.SES.FEve)) slope(seg.SES.FEve) # Segmented regression analysis of SES functional diversity metrics based on diet ---------- ## SES.FRic. df.ses.fric.diet <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fric.diet) summary(lm.SES.FRic.diet<-lm(SES~Elevation, data=data.4)) davies.test(lm.SES.FRic.diet,seg.Z =~Elevation, k=8) summary(seg.SES.FRic.diet<-segmented(lm.SES.FRic.diet)) slope(seg.SES.FRic.diet) ## SES.FNND. df.ses.fnnd.diet <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fnnd.diet) summary(lm.SES.FNND.diet<-lm(SES~Elevation, data= df.ses.fnnd.diet)) davies.test(lm.SES.FNND.diet,seg.Z =~Elevation, k=8) summary(seg.SES.FNND.diet<-segmented(lm.SES.FNND.diet)) slope(seg.SES.FNND.diet) ## SES.FEve. df.ses.feve.diet <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.feve.diet) summary(lm.SES.FEve.diet<-lm(SES~Elevation, data=df.ses.feve.diet)) davies.test(lm.SES.FEve.diet,seg.Z =~Elevation, k=8) summary(seg.SES.FEve.diet<-segmented(lm.SES.FEve.diet)) slope(seg.SES.FEve.diet) # Segmented regression analysis of SES functional diversity metrics based on vegetation strata ---------- ## SES.FRic. df.ses.fric.strata <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fric.strata) summary(lm.SES.FRic.strata<-lm(SES~Elevation, data= df.ses.fric.strata)) davies.test(lm.SES.FRic.strata,seg.Z =~Elevation, k=8) summary(seg.SES.FRic.strata<-segmented(lm.SES.FRic.strata)) slope(seg.SES.FRic.strata) ## SES.FNND. df.ses.fnnd.strata <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fnnd.strata) summary(lm.SES.FNND.strata<-lm(SES~Elevation, data= df.ses.fnnd.strata)) davies.test(lm.SES.FNND.strata,seg.Z =~Elevation, k=8) summary(seg.SES.FNND.strata<-segmented(lm.SES.FNND.strata)) slope(seg.SES.FNND.strata) ## SES.FEve. df.ses.feve.strata <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.feve.strata) summary(lm.SES.FEve.strata<-lm(SES~Elevation, data= df.ses.feve.strata)) davies.test(lm.SES.FEve.strata,seg.Z =~Elevation, k=8) summary(seg.SES.FEve.strata<-segmented(lm.SES.FEve.strata)) slope(seg.SES.FEve.strata) # Segmented regression analysis of SES functional diversity metrics based on foraging modes ---------- ## SES.FRic. df.ses.fric.tactics <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fric.tactics) summary(lm.SES.FRic.tactics<-lm(SES~Elevation, data= df.ses.fric.tactics)) davies.test(lm.SES.FRic.tactics,seg.Z =~Elevation, k=8) summary(seg.SES.FRic.tactics<-segmented(lm.SES.FRic.tactics)) slope(seg.SES.FRic.tactics) ## SES.FNND. df.ses.fnnd.tactics <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fnnd.tactics) summary(lm.SES.FNND.tactics<-lm(SES~Elevation, data= df.ses.fnnd.tactics)) davies.test(lm.SES.FNND.tactics,seg.Z =~Elevation, k=8) summary(seg.SES.FNND.tactics<-segmented(lm.SES.FNND.tactics)) slope(seg.SES.FNND.tactics) ## SES.FEve. df.ses.feve.tactics <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.feve.tactics) summary(lm.SES.FEve.tactics<-lm(SES~Elevation, data= df.ses.feve.tactics)) davies.test(lm.SES.FEve.tactics,seg.Z =~Elevation, k=8) summary(seg.SES.FEve.tactics<-segmented(lm.SES.FEve.tactics)) slope(seg.SES.FEve.tactics) # Segmented regression analysis of SES functional diversity metrics based on all traits ---------- ## SES.FRic. df.ses.fric.all.traits <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fric.all.traits) summary(lm.SES.FRic.combined<-lm(SES~Elevation, data= df.ses.fric.all.traits)) davies.test(lm.SES.FRic.combined,seg.Z =~Elevation, k=8) summary(seg.SES.FRic.combined<-segmented(lm.SES.FRic.combined)) slope(seg.SES.FRic.combined) ## SES.FNND. df.ses.fnnd.all.traits <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.fnnd.all.traits) summary(lm.SES.FNND.combined<-lm(SES~ Elevation, data=df.ses.fnnd.all.traits)) davies.test(lm.SES.FNND.combined,seg.Z =~Elevation, k=8) summary(seg.SES.FNND.combined<-segmented(lm.SES.FNND.combined)) slope(seg.SES.FNND.combined) ## SES.FEve. df.ses.feve.all.traits <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.value.feve.all.traits) summary(lm.SES.FEve.combined<-lm(SES~Elevation, data=df.ses.feve.all.traits)) davies.test(lm.SES.FEve.combined,seg.Z =~Elevation, k=8) summary(seg.SES.FEve.combined<-segmented(lm.SES.FEve.combined)) slope(seg.SES.FEve.combined) # Segmented regression analysis of SES phylogenetic diversity metrics ---------- ## SES.MPD. df.ses.mpd <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.mpd$mpd.obs.z) summary(lm.SES.MPD<-lm(SES~Elevation, data= df.ses.mpd)) davies.test(lm.SES.MPD,seg.Z =~Elevation, k=8) summary(seg.SES.MPD<-segmented(lm.SES.MPD)) slope(seg.SES.MPD) ## SES.MNTD. df.ses.mntd <- data.frame(Elevation = c(350, 550, 750, 950, 1150, 1350, 1550, 1750, 1950, 2150), SES = ses.mntd$mntd.obs.z) summary(lm.SES.MNTD<-lm(SES~Elevation, data= df.ses.mntd)) davies.test(lm.SES.MNTD,seg.Z =~Elevation, k=8) summary(seg.SES.MNTD<-segmented(lm.SES.MNTD)) slope(seg.SES.MNTD) # Niche packing/expansion contributions based on morphological traits ------- ## Since the functional volume of the 550-m-assemblage is lower than ## that of the 350-m-assemblage and 750-m-assemblage, it implies that ## all species in the species richer assemblage (550 m) fall inside ## the functional volume of the two adjacent species poorer assemblages ## (100% niche packing). The same holds for the 1350-m-assemblage ## vs the 1550-m-assemblage, and the 1950-m-assemblage vs the ## 2150-m-assemblage. ## We start comparing the 750-m-assemblage against the 950-m-assemblage, ## following the species richness gradient. ## Getting the 750-m-assemblage and the 950-m-assemblage. com.750 <- as.matrix(occurrence.matrix[3,occurrence.matrix[3,]>0]) com.950 <- as.matrix(occurrence.matrix[4,occurrence.matrix[4,]>0]) ## Find unique species. species.750 <- rownames(com.750) species.950 <- rownames(com.950) com_com <- intersect(species.750,species.950) species <- species.750[!species.750 %in% com_com] ## Compute the functional space and extract PCoA scores ## for the 750-m-assemblage. fspace <- tr.cont.fspace(sp_tr = morphotraits, pca = T, nb_dim = 4, scaling = "scale_center", compute_corr = "pearson") sp_faxes_coord_morphology <- fspace$sp_faxes_coord fspace.750m <- sp_faxes_coord_morphology[species.750,] ## Prune the occurrence matrix to include only species present ## in the 750-m-assemblage. community.750m <- occurrence.matrix[,species.750] ## Remove sequentially each species from the list and calculate ## the resulting functional richness (i.e. functional volume). fun <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.750m[-f,] c <- community.750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes <- sapply(1:136,fun) # In total, there are 136 species. ## Extract functional richness values from the list ## and store therm in a matrix. indices.volumes <-seq (1,271,2) # 271 is the number of dataframes in the list. matrix.volumes <- NULL for (i in indices.volumes) { matrix.volumes <- cbind(matrix.volumes, volumes [[c(i,2)]]) } colnames(matrix.volumes)<- rownames(fspace.750m) matrix.volumes <- matrix.volumes[3,] ## Prune the matrix in order to select only unique species composing the ## assemblage. vol.unique.750 <- as.matrix(matrix.volumes[species]) vol.unique.750 ## Function to calculate niche packing and expansion contributions. calculate_niche_contributions <- function(matrix_volumes, species, fd_ind_values, fspace, community, target_value, r) { vol_unique <- as.matrix(matrix_volumes[species]) # Order the matrix and get row indices ord_matrix <- vol_unique[order(vol_unique[, 1]), , drop = FALSE] rows_to_remove <- numeric(0) # Store row indices to remove # Iteratively remove the lowest values for (i in seq_len(nrow(ord_matrix))) { rows_to_remove <- c(rows_to_remove, which(rownames(fspace) == rownames(ord_matrix)[i])) # Calculate the new volume after removing the species new_vol <- alpha.fd.multidim( sp_faxes_coord = fspace[-rows_to_remove, ], asb_sp_w = community[, -rows_to_remove], ind_vect = "fric", scaling = FALSE ) # Check if the new functional richness is less than or equal to the target if (new_vol$functional_diversity_indices$fric[r] <= target_value) { break } } # Calculate the percentage of niche packing and niche expansion niche_expansion <- length(rows_to_remove) / length(rownames(vol_unique)) niche_packing <- 1 - niche_expansion list( niche_expansion = label_percent()(niche_expansion), niche_packing = label_percent()(niche_packing) ) } ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes, species = species, fd_ind_values = fd_ind_values_morphology, fspace = fspace.750m, community = community.750m, target_value = fd_ind_values_morphology[4, 2], # functional volume of the species-poorer assemblage r = 3) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 950-m-assemblage and the 1150-m-assemblage. com.1150 <- as.matrix(occurrence.matrix[5,occurrence.matrix[5,]>0]) ## Find unique species species.1150 <- rownames(com.1150) species.950 <- rownames(com.950) com_com.2 <- intersect(species.1150,species.950) species.2 <- species.950[!species.950 %in% com_com.2] ## Compute the functional space and extract PCoA scores for ## the 950-m-assemblage. fspace <- tr.cont.fspace(sp_tr = morphotraits, pca = T, nb_dim = 4, scaling = "scale_center", compute_corr = "pearson") sp_faxes_coord_morphology <- fspace$sp_faxes_coord fspace.950m <- sp_faxes_coord_morphology[species.950,] ## Prune the community data set to include only species present ## in the 950-m-assemblage. community.950m <- occurrence.matrix[,species.950] ## Function to remove sequentially each species from the list ## and calculate the resulting morphovolume. fun.2 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.950m[-f,] c <- community.950m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.2 <- sapply(1:131,fun.2) # In total, there are 131 species. ## Extract Functional richess values from the list and ## obtain a matrix. indices.2 <-seq (1,261,2) # number of datasets in the list. newset.2 <- NULL for (i in indices.2) { newset.2 <- cbind(newset.2, volumes.2[[c(i,2)]]) } View(newset.2) matrix.volumes.2 <- newset.2 colnames(matrix.volumes.2)<- rownames(fspace.950m) matrix.volumes.2 <- matrix.volumes.2[4,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.2, species = species.2, fd_ind_values = fd_ind_values_morphology, fspace = fspace.950m, community = community.950m, target_value = fd_ind_values_morphology[5, 2], r = 4) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1150-m-assemblage and the 1350-m-assemblage. com.1350 <- as.matrix(occurrence.matrix[6,occurrence.matrix[6,]>0]) ## Find unique species present in the 1150-m-assemblage. species.1150 <- rownames(com.1150) species.1350 <- rownames(com.1350) com_com.3 <- intersect(species.1150,species.1350) species.3 <- species.1150[!species.1150 %in% com_com.3] ## Compute the functional space and extract PCoA scores ## for the 1150-m-assemblage. fspace <- tr.cont.fspace(sp_tr = morphotraits, pca = T, nb_dim = 4, scaling = "scale_center", compute_corr = "pearson") sp_faxes_coord_morphology <- fspace$sp_faxes_coord fspace.1150m <- sp_faxes_coord_morphology[species.1150,] ## Prune the community data set to include only species ## present in the 1150-m-assemblage. community.1150m <- occurrence.matrix[,species.1150] ## Function to remove sequentially each species from the list ## and calculate the resulting functional volume. fun.3 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1150m[-f,] c <- community.1150m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.3 <- sapply(1:105,fun.3) # 105 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.3 <-seq (1,209,2) # 209 data frames in the list. newset.3 <- NULL for (i in indices.3) { newset.3 <- cbind(newset.3, volumes.3[[c(i,2)]]) } View(newset.3) matrix.volumes.3 <- newset.3 colnames(matrix.volumes.3)<- rownames(fspace.1150m) matrix.volumes.3 <- matrix.volumes.3[5,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.3, species = species.3, fd_ind_values = fd_ind_values_morphology, fspace = fspace.1150m, community = community.1150m, target_value = fd_ind_values_morphology[6, 2], r = 5) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1550-m-assemblage and the 1750-m-assemblage. com.1550 <- as.matrix(occurrence.matrix[7,occurrence.matrix[7,]>0]) com.1750 <- as.matrix(occurrence.matrix[8,occurrence.matrix[8,]>0]) ## Find species not shared between the two assemblage. species.1550 <- rownames(com.1550) species.1750 <- rownames(com.1750) com_com.4 <- intersect(species.1550,species.1750) species.4 <- species.1550[!species.1550 %in% com_com.4] ## Generate the functional space and extract PCoA scores ## for the 1550-m-assemblage. fspace <- tr.cont.fspace(sp_tr = morphotraits, pca = T, nb_dim = 4, scaling = "scale_center", compute_corr = "pearson") sp_faxes_coord_morphology <- fspace$sp_faxes_coord fspace.1550m <- sp_faxes_coord_morphology[species.1550,] ## Prune the community data set to ## include only species present in the 1550-m-assemblage. community.1550m <- occurrence.matrix[,species.1550] # Function to remove sequentially each species from the list # and calculate the resulting functional volume. fun.4 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1550m[-f,] c <- community.1550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.4 <- sapply(1:84,fun.4) volumes.4 ## Extract FRic values from the list and obtain a matrix. indices.4 <-seq (1,167,2) # 167 dataframes in the list. newset.4 <- NULL for (i in indices.4) { newset.4 <- cbind(newset.4, volumes.4[[c(i,2)]]) } View(newset.4) matrix.volumes.4 <- newset.4 colnames(matrix.volumes.4)<- rownames(fspace.1550m) matrix.volumes.4 <- matrix.volumes.4[7,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.4, species = species.4, fd_ind_values = fd_ind_values_morphology, fspace = fspace.1550m, community = community.1550m, target_value = fd_ind_values_morphology[8, 2], r = 7) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1750-m-assemblage and the 1950-m-assemblage. com.1950 <- as.matrix(occurrence.matrix[9,occurrence.matrix[9,]>0]) com.1750 <- as.matrix(occurrence.matrix[8,occurrence.matrix[8,]>0]) ## Find species not shared between the two assemblages. species.1950 <- rownames(com.1950) species.1750 <- rownames(com.1750) com_com.5 <- intersect(species.1950,species.1750) species.5 <- species.1750[!species.1750 %in% com_com.5] ## Generate the functional space and extract PCoA scores ## for the 1750-m-assemblage. fspace <- tr.cont.fspace(sp_tr = morphotraits, pca = T, nb_dim = 4, scaling = "scale_center", compute_corr = "pearson") sp_faxes_coord_morphology <- fspace$sp_faxes_coord fspace.1750m <- sp_faxes_coord_morphology[species.1750,] ## Prune the community data set to include only species ## present in the 1750-m-assemblage. community.1750m <- occurrence.matrix[,species.1750] ## Function to remove sequentially each species from the list ## and calculate the resulting functional volume. fun.5 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1750m[-f,] c <- community.1750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.5 <- sapply(1:62,fun.5) # 62 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.5 <-seq (1,123,2) newset.5 <- NULL for (i in indices.5) { newset.5 <- cbind(newset.5, volumes.5[[c(i,2)]]) } View(newset.5) matrix.volumes.5 <- newset.5 colnames(matrix.volumes.5)<- rownames(fspace.1750m) matrix.volumes.5 <- matrix.volumes.5[8,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.5, species = species.5, fd_ind_values = fd_ind_values_morphology, fspace = fspace.1750m, community = community.1750m, target_value = fd_ind_values_morphology[9, 2], r = 8) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) # Niche packing/expansion contributions based on dietary traits ------------- ## The observed functional volume of the 750-m-assemblage is lower than ## that of the 950-m-assemblage, implying that ## all species in the species richer assemblage (750 m) fall inside ## the functional volume of the species poorer assemblage ## (100% niche packing). The same holds for the 1550-m-assemblage ## vs the 1750-m-assemblage, and the 1950-m-assemblage vs the ## 2150-m-assemblage. ## Comparison between the 350-m-assemblage and the 550-m-assemblage. com.350 <- as.matrix(occurrence.matrix[1,occurrence.matrix[1,]>0]) com.550 <- as.matrix(occurrence.matrix[2,occurrence.matrix[2,]>0]) ## Select species not shared between the two assemblages. species.350 <- rownames(com.350) species.550 <- rownames(com.550) com_com <- intersect(species.350,species.550) species <- species.550[!species.550 %in% com_com] ## Generate the functional space and extract PCoA scores ## for the 550-m-assemblage. fspace.550m <- sp_faxes_coord_diet[species.550, c("PC1", "PC2", "PC3", "PC4", "PC5")] ## Prune the community data set to include only species present in ## the 550-m-assemblage. community.550m <- occurrence.matrix[,species.550] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.550m[-f,] c <- community.550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes <- sapply(1:139,fun) # 139 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices <-seq (1,277,2) # 277 dataframes in the list. newset <- NULL for (i in indices) { newset <- cbind(newset, volumes [[c(i,2)]]) } View(newset) matrix.volumes <- newset colnames(matrix.volumes)<- rownames(fspace.550m) matrix.volumes <- matrix.volumes[2,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes, species = species, fd_ind_values = fd_ind_values_diet, fspace = fspace.550m, community = community.550m, target_value = fd_ind_values_diet[1, 2], r = 2) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 550-m-assemblage and the 750-m-assemblage. ## Select species not shared between the two assemblages. species.750 <- rownames(com.750) species.550 <- rownames(com.550) com_com.2 <- intersect(species.750,species.550) species.2 <- species.550[!species.550 %in% com_com.2] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes, species = species.2, fd_ind_values = fd_ind_values_diet, fspace = fspace.550m, community = community.550m, target_value = fd_ind_values_diet[3, 2], r = 2) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 950-m-assemblage and the 1150-m-assemblage. com.950 <- as.matrix(occurrence.matrix[4,occurrence.matrix[4,]>0]) com.1150 <- as.matrix(occurrence.matrix[5,occurrence.matrix[5,]>0]) ## Select species not shared between the two assemblages. species.950 <- rownames(com.950) species.1150 <- rownames(com.1150) com_com.3 <- intersect(species.950,species.1150) species.3 <- species.950[!species.950 %in% com_com.3] ## Generate the functional space and extract PCoA scores ## for the 950-m-assemblage. fspace.950m <- sp_faxes_coord_diet[species.950, c("PC1", "PC2", "PC3", "PC4", "PC5")] ## Prune the community data set to include only species present in ## the 950-m-assemblage. community.950m <- occurrence.matrix[,species.950] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.3 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.950m[-f,] c <- community.950m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.3 <- sapply(1:131,fun.3) # 131 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.3 <-seq (1,261,2) # 261 dataframes in the list. newset.3 <- NULL for (i in indices.3) { newset.3 <- cbind(newset.3, volumes.3 [[c(i,2)]]) } View(newset.3) matrix.volumes.3 <- newset.3 colnames(matrix.volumes.3)<- rownames(fspace.950m) matrix.volumes.3 <- matrix.volumes.3[4,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.3, species = species.3, fd_ind_values = fd_ind_values_diet, fspace = fspace.950m, community = community.950m, target_value = fd_ind_values_diet[5, 2], r = 4) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1150-m-assemblage and the 1350-m-assemblage. com.1350 <- as.matrix(occurrence.matrix[6,occurrence.matrix[6,]>0]) ## Select species not shared between the two assemblages. species.1350 <- rownames(com.1350) species.1150 <- rownames(com.1150) com_com.4 <- intersect(species.1350,species.1150) species.4 <- species.1150[!species.1150 %in% com_com.4] ## Generate the functional space and extract PCoA scores ## for the 1150-m-assemblage. fspace.1150m <- sp_faxes_coord_diet[species.1150, c("PC1", "PC2", "PC3", "PC4", "PC5")] ## Prune the community data set to include only species present in ## the 1150-m-assemblage. community.1150m <- occurrence.matrix[,species.1150] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.4 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1150m[-f,] c <- community.1150m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.4 <- sapply(1:105,fun.4) # 105 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.4 <-seq (1,209,2) # 209 dataframes in the list. newset.4 <- NULL for (i in indices.4) { newset.4 <- cbind(newset.4, volumes.4 [[c(i,2)]]) } View(newset.4) matrix.volumes.4 <- newset.4 colnames(matrix.volumes.4)<- rownames(fspace.1150m) matrix.volumes.4 <- matrix.volumes.4[5,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.4, species = species.4, fd_ind_values = fd_ind_values_diet, fspace = fspace.1150m, community = community.1150m, target_value = fd_ind_values_diet[6, 2], r = 5) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1550-m-assemblage and the 1750-m-assemblage. com.1750 <- as.matrix(occurrence.matrix[8,occurrence.matrix[8,]>0]) com.1550 <- as.matrix(occurrence.matrix[7,occurrence.matrix[7,]>0]) ## Select species not shared between the two assemblages. species.1750 <- rownames(com.1750) species.1550 <- rownames(com.1550) com_com.x <- intersect(species.1750,species.1550) species.x <- species.1550[!species.1550 %in% com_com.x] ## Generate the functional space and extract PCoA scores ## for the 1550-m-assemblage. fspace.1550m <- sp_faxes_coord_diet[species.1550, c("PC1", "PC2", "PC3", "PC4", "PC5")] ## Prune the community data set to include only species present in ## the 1550-m-assemblage. community.1550m <- occurrence.matrix[,species.1550] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.x <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1550m[-f,] c <- community.1550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.x <- sapply(1:84,fun.x) # 84 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.x <-seq (1,167,2) # 167 dataframes in the list. newset.x <- NULL for (i in indices.x) { newset.x <- cbind(newset.x, volumes.x [[c(i,2)]]) } View(newset.x) matrix.volumes.x <- newset.x colnames(matrix.volumes.x)<- rownames(fspace.1550m) matrix.volumes.x <- matrix.volumes.x[7,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.x, species = species.x, fd_ind_values = fd_ind_values_diet, fspace = fspace.1550m, community = community.1550m, target_value = fd_ind_values_diet[8, 2], r = 7) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1750-m-assemblage and the 1950-m-assemblage. com.1750 <- as.matrix(occurrence.matrix[8,occurrence.matrix[8,]>0]) com.1950 <- as.matrix(occurrence.matrix[9,occurrence.matrix[9,]>0]) ## Select species not shared between the two assemblages. species.1750 <- rownames(com.1750) species.1950 <- rownames(com.1950) com_com.5 <- intersect(species.1750,species.1950) species.5 <- species.1750[!species.1750 %in% com_com.5] ## Generate the functional space and extract PCoA scores ## for the 1750-m-assemblage. fspace.1750m <- sp_faxes_coord_diet[species.1750, c("PC1", "PC2", "PC3", "PC4", "PC5")] ## Prune the community data set to include only species present in ## the 1750-m-assemblage. community.1750m <- occurrence.matrix[,species.1750] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.5 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1750m[-f,] c <- community.1750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.5 <- sapply(1:62,fun.5) # 62 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.5 <-seq (1,123,2) # 123 data frames in the list. newset.5 <- NULL for (i in indices.5) { newset.5 <- cbind(newset.5, volumes.5 [[c(i,2)]]) } View(newset.5) matrix.volumes.5 <- newset.5 colnames(matrix.volumes.5)<- rownames(fspace.1750m) matrix.volumes.5 <- matrix.volumes.5[8,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.5, species = species.5, fd_ind_values = fd_ind_values_diet, fspace = fspace.1750m, community = community.1750m, target_value = fd_ind_values_diet[9, 2], r = 8) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) # Niche packing/expansion contributions based on vegetation strata usage ------------- ## The observed functional volume of the 550-m-assemblage is lower than ## that of the 350-m-assemblage, implying that ## all species in the species richer assemblage (550 m) fall inside ## the functional volume of the species poorer assemblage ## (100% niche packing). The same holds for the 550-m-assemblage vs ## the 750-m-assemblgae, the 750-m-assemblage ## vs the 950-m-assemblage, and the 950-m-assemblage vs the ## 1150-m-assemblage. ## Comparison between the 1150-m-assemblage and the 1350-m-assemblage. com.1350 <- as.matrix(occurrence.matrix[6,occurrence.matrix[6,]>0]) ## Select species not shared between the two assemblages. species.1150 <- rownames(com.1150) species.1350 <- rownames(com.1350) com_com.2 <- intersect(species.1150,species.1350) species.2 <- species.1150[!species.1150 %in% com_com.2] ## Generate the functional space and extract PCoA scores ## for the 1150-m-assemblage. fspace.1150m <- sp_faxes_coord_strata[species.1150, c("PC1", "PC2", "PC3")] ## Prune the community data set to include only species present in ## the 1150-m-assemblage. community.1150m <- occurrence.matrix[,species.1150] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.2 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1150m[-f,] c <- community.1150m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.2 <- sapply(1:105,fun.2) # 105 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.2 <-seq (1,209,2) # 209 data frames in the list. newset.2 <- NULL for (i in indices.2) { newset.2 <- cbind(newset.2, volumes.2 [[c(i,2)]]) } View(newset.2) matrix.volumes.2 <- newset.2 colnames(matrix.volumes.2)<- rownames(fspace.1150m) matrix.volumes.2 <- matrix.volumes.2[5,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.2, species = species.2, fd_ind_values = fd_ind_values_strata, fspace = fspace.1150m, community = community.1150m, target_value = fd_ind_values_strata[6, 3], r = 5) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1350-m-assemblage and the 1550-m-assemblage. com.1550 <- as.matrix(occurrence.matrix[7,occurrence.matrix[7,]>0]) com.1350 <- as.matrix(occurrence.matrix[6,occurrence.matrix[6,]>0]) ## Select species not shared between the two assemblages. species.1550 <- rownames(com.1550) species.1350 <- rownames(com.1350) com_com.3 <- intersect(species.1550,species.1350) species.3 <- species.1350[!species.1350 %in% com_com.3] ## Generate the functional space and extract PCoA scores ## for the 1350-m-assemblage. fspace.1350m <- sp_faxes_coord_strata[species.1350, c("PC1", "PC2", "PC3")] ## Prune the community data set to include only species present in ## the 1350-m-assemblage. community.1350m <- occurrence.matrix[,species.1350] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.3 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1350m[-f,] c <- community.1350m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.3 <- sapply(1:85,fun.3) # 85 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.3 <-seq (1,169,2) # 169 data frames in the list. newset.3 <- NULL for (i in indices.3) { newset.3 <- cbind(newset.3, volumes.3 [[c(i,2)]]) } View(newset.3) matrix.volumes.3 <- newset.3 colnames(matrix.volumes.3)<- rownames(fspace.1350m) matrix.volumes.3 <- matrix.volumes.3[6,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.3, species = species.3, fd_ind_values = fd_ind_values_strata, fspace = fspace.1350m, community = community.1350m, target_value = fd_ind_values_strata[7, 3], r = 6) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1550-m-assemblage and the 1750-m-assemblage. ## Select species not shared between the two assemblages. species.1550 <- rownames(com.1550) species.1750 <- rownames(com.1750) com_com.4 <- intersect(species.1550,species.1750) species.4 <- species.1550[!species.1550 %in% com_com.4] ## Generate the functional space and extract PCoA scores ## for the 1550-m-assemblage. fspace.1550m <- sp_faxes_coord_strata[species.1550, c("PC1", "PC2", "PC3")] ## Prune the community data set to include only species present in ## the 1550-m-assemblage. community.1550m <- occurrence.matrix[,species.1550] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.4 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1550m[-f,] c <- community.1550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.4 <- sapply(1:84,fun.4) # 84 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.4 <-seq (1,167,2) # 167 data frames in the list. newset.4 <- NULL for (i in indices.4) { newset.4 <- cbind(newset.4, volumes.4 [[c(i,2)]]) } View(newset.4) matrix.volumes.4 <- newset.4 colnames(matrix.volumes.4)<- rownames(fspace.1550m) matrix.volumes.4 <- matrix.volumes.4[7,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.4, species = species.4, fd_ind_values = fd_ind_values_strata, fspace = fspace.1550m, community = community.1550m, target_value = fd_ind_values_strata[8, 3], r = 7) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1750-m-assemblage and the 1950-m-assemblage. ## Select species not shared between the two assemblages. species.1950 <- rownames(com.1950) species.1750 <- rownames(com.1750) com_com.5 <- intersect(species.1950,species.1750) species.5 <- species.1750[!species.1750 %in% com_com.5] ## Generate the functional space and extract PCoA scores ## for the 1750-m-assemblage. fspace.1750m <- sp_faxes_coord_strata[species.1750, c("PC1", "PC2", "PC3")] ## Prune the community data set to include only species present in ## the 1750-m-assemblage. community.1750m <- occurrence.matrix[,species.1750] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.5 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1750m[-f,] c <- community.1750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.5 <- sapply(1:62,fun.5) # 62 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.5 <-seq (1,123,2) # 123 data frames in the list. newset.5 <- NULL for (i in indices.5) { newset.5 <- cbind(newset.5, volumes.5 [[c(i,2)]]) } View(newset.5) matrix.volumes.5 <- newset.5 colnames(matrix.volumes.5)<- rownames(fspace.1750m) matrix.volumes.5 <- matrix.volumes.5[8,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.5, species = species.5, fd_ind_values = fd_ind_values_strata, fspace = fspace.1750m, community = community.1750m, target_value = fd_ind_values_strata[9, 3], r = 8) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1950-m-assemblage and the 2150-m-assemblage. com.1950 <- as.matrix(occurrence.matrix[9,occurrence.matrix[9,]>0]) com.2150 <- as.matrix(occurrence.matrix[10,occurrence.matrix[10,]>0]) ## Select species not shared between the two assemblages. species.1950 <- rownames(com.1950) species.2150 <- rownames(com.2150) com_com.6 <- intersect(species.1950,species.2150) species.6 <- species.1950[!species.1950 %in% com_com.6] ## Generate the functional space and extract PCoA scores ## for the 1950-m-assemblage. fspace.1950m <- sp_faxes_coord_strata[species.1950, c("PC1", "PC2", "PC3")] ## Prune the community data set to include only species present in ## the 1950-m-assemblage. community.1950m <- occurrence.matrix[,species.1950] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.6 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1950m[-f,] c <- community.1950m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.6 <- sapply(1:57,fun.6) # 57 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.6 <-seq (1,113,2) # 113 data frames in the list. newset.6 <- NULL for (i in indices.6) { newset.6 <- cbind(newset.6, volumes.6 [[c(i,2)]]) } View(newset.6) matrix.volumes.6 <- newset.6 colnames(matrix.volumes.6)<- rownames(fspace.1950m) matrix.volumes.6 <- matrix.volumes.6[9,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.6, species = species.6, fd_ind_values = fd_ind_values_strata, fspace = fspace.1950m, community = community.1950m, target_value = fd_ind_values_strata[10, 3], r = 9) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) # Niche packing/expansion contributions based on foraging tactics ------------- ## The observed functional volume of the 1950-m-assemblage is lower than ## that of the 2150-m-assemblage, implying that ## all species in the species richer assemblage (1950 m) fall inside ## the functional volume of the species poorer assemblage ## (100% niche packing). ## Comparison between the 550-m-assemblage and the 350-m-assemblage. ## Select species not shared between the two assemblages. species.350 <- rownames(com.350) species.550 <- rownames(com.550) com_com <- intersect(species.350,species.550) species <- species.550[!species.550 %in% com_com] ## Generate the functional space and extract PCoA scores ## for the 550-m-assemblage. fspace.550m <- sp_faxes_coord_tactics[species.550, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6")] ## Prune the community data set to include only species present in ## the 550-m-assemblage. community.550m <- occurrence.matrix[,species.550] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.550m[-f,] c <- community.550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes <- sapply(1:139,fun) # 139 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices <-seq (1,277,2) # 277 data frames in the list. newset <- NULL for (i in indices) { newset <- cbind(newset, volumes [[c(i,2)]]) } View(newset) matrix.volumes <- newset colnames(matrix.volumes)<- rownames(fspace.550m) matrix.volumes <- matrix.volumes[2,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes, species = species, fd_ind_values = fd_ind_values_tactics, fspace = fspace.550m, community = community.550m, target_value = fd_ind_values_tactics[1, 3], r = 2) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 550-m-assemblage and the 750-m-assemblage. ## Select species not shared between the two assemblages. species.750 <- rownames(com.750) species.550 <- rownames(com.550) com_com.2 <- intersect(species.750,species.550) species.2 <- species.550[!species.550 %in% com_com.2] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes, species = species.2, fd_ind_values = fd_ind_values_tactics, fspace = fspace.550m, community = community.550m, target_value = fd_ind_values_tactics[3, 3], r = 2) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 750-m-assemblage and the 950-m-assemblage. ## Select species not shared between the two assemblages. species.750 <- rownames(com.750) species.950 <- rownames(com.950) com_com.x <- intersect(species.950,species.750) species.x <- species.750[!species.750 %in% com_com.x] ## Generate the functional space and extract PCoA scores ## for the 750-m-assemblage. fspace.750m <- sp_faxes_coord_tactics[species.750, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6")] ## Prune the community data set to include only species present in ## the 750-m-assemblage. community.750m <- occurrence.matrix[,species.750] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.x <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.750m[-f,] c <- community.750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.x <- sapply(1:136,fun.x) # 136 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.x <-seq (1,271,2) # 271 data frames in the list. newset.x <- NULL for (i in indices.x) { newset.x <- cbind(newset.x, volumes.x [[c(i,2)]]) } View(newset.x) matrix.volumes.x <- newset.x colnames(matrix.volumes.x)<- rownames(fspace.750m) matrix.volumes.x <- matrix.volumes.x[3,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.x, species = species.x, fd_ind_values = fd_ind_values_tactics, fspace = fspace.750m, community = community.750m, target_value = fd_ind_values_tactics[4, 3], r = 3) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 950-m-assemblage and the 1150-m-assemblage. ## Select species not shared between the two assemblages. species.1150 <- rownames(com.1150) species.950 <- rownames(com.950) com_com.3 <- intersect(species.950,species.1150) species.3 <- species.950[!species.950 %in% com_com.3] ## Generate the functional space and extract PCoA scores ## for the 950-m-assemblage. fspace.950m <- sp_faxes_coord_tactics[species.950, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6")] ## Prune the community data set to include only species present in ## the 950-m-assemblage. community.950m <- occurrence.matrix[,species.950] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.3 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.950m[-f,] c <- community.950m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.3 <- sapply(1:131,fun.3) # 131 species present in the assemblage. volumes.3 ## Extract FRic values from the list and obtain a matrix. indices.3 <-seq (1,261,2) # 261 data frames in the list. newset.3 <- NULL for (i in indices.3) { newset.3 <- cbind(newset.3, volumes.3 [[c(i,2)]]) } View(newset.3) matrix.volumes.3 <- newset.3 colnames(matrix.volumes.3)<- rownames(fspace.950m) matrix.volumes.3 <- matrix.volumes.3[4,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.3, species = species.3, fd_ind_values = fd_ind_values_tactics, fspace = fspace.950m, community = community.950m, target_value = fd_ind_values_tactics[5, 3], r = 4) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1150-m-assemblage and the 1350-m-assemblage. ## Select species not shared between the two assemblages. species.1150 <- rownames(com.1150) species.1350 <- rownames(com.1350) com_com.4 <- intersect(species.1350,species.1150) species.4 <- species.1150[!species.1150 %in% com_com.4] ## Generate the functional space and extract PCoA scores ## for the 1150-m-assemblage. fspace.1150m <- sp_faxes_coord_tactics[species.1150, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6")] ## Prune the community data set to include only species present in ## the 1150-m-assemblage. community.1150m <- occurrence.matrix[,species.1150] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.4 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1150m[-f,] c <- community.1150m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.4 <- sapply(1:105,fun.4) # 105 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.4 <-seq (1,209,2) # 209 data frames in the list. newset.4 <- NULL for (i in indices.4) { newset.4 <- cbind(newset.4, volumes.4 [[c(i,2)]]) } View(newset.4) matrix.volumes.4 <- newset.4 colnames(matrix.volumes.4)<- rownames(fspace.1150m) matrix.volumes.4 <- matrix.volumes.4[5,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.4, species = species.4, fd_ind_values = fd_ind_values_tactics, fspace = fspace.1150m, community = community.1150m, target_value = fd_ind_values_tactics[6, 3], r = 5) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1350-m-assemblage and the 1550-m-assemblage. ## Select species not shared between the two assemblages. species.1550 <- rownames(com.1550) species.1350 <- rownames(com.1350) com_com.5 <- intersect(species.1350,species.1550) species.5 <- species.1350[!species.1350 %in% com_com.5] ## Generate the functional space and extract PCoA scores ## for the 1350-m-assemblage. fspace.1350m <- sp_faxes_coord_tactics[species.1350, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6")] ## Prune the community data set to include only species present in ## the 1350-m-assemblage. community.1350m <- occurrence.matrix[,species.1350] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.5 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1350m[-f,] c <- community.1350m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.5 <- sapply(1:85,fun.5) # 85 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.5 <-seq (1,169,2) # 169 data frames in the list. newset.5 <- NULL for (i in indices.5) { newset.5 <- cbind(newset.5, volumes.5 [[c(i,2)]]) } View(newset.5) matrix.volumes.5 <- newset.5 colnames(matrix.volumes.5)<- rownames(fspace.1350m) matrix.volumes.5 <- matrix.volumes.5[6,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.5, species = species.5, fd_ind_values = fd_ind_values_tactics, fspace = fspace.1350m, community = community.1350m, target_value = fd_ind_values_tactics[7, 3], r = 6) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1550-m-assemblage and the 1750-m-assemblage. ## Select species not shared between the two assemblages. species.1550 <- rownames(com.1550) species.1750 <- rownames(com.1750) com_com.6 <- intersect(species.1750,species.1550) species.6 <- species.1550[!species.1550 %in% com_com.6] ## Generate the functional space and extract PCoA scores ## for the 1550-m-assemblage. fspace.1550m <- sp_faxes_coord_tactics[species.1550, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6")] ## Prune the community data set to include only species present in ## the 1550-m-assemblage. community.1550m <- occurrence.matrix[,species.1550] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.6 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1550m[-f,] c <- community.1550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.6 <- sapply(1:84,fun.6) # 84 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.6 <-seq (1,167,2) # 167 data frames in the list. newset.6 <- NULL for (i in indices.6) { newset.6 <- cbind(newset.6, volumes.6 [[c(i,2)]]) } View(newset.6) matrix.volumes.6 <- newset.6 colnames(matrix.volumes.6)<- rownames(fspace.1550m) matrix.volumes.6 <- matrix.volumes.6[7,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.6, species = species.6, fd_ind_values = fd_ind_values_tactics, fspace = fspace.1550m, community = community.1550m, target_value = fd_ind_values_tactics[8, 3], r = 7) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1750-m-assemblage and the 1950-m-assemblage. ## Select species not shared between the two assemblages. species.1950 <- rownames(com.1950) species.1750 <- rownames(com.1750) com_com.y <- intersect(species.1950,species.1750) species.y <- species.1750[!species.1750 %in% com_com.y] ## Generate the functional space and extract PCoA scores ## for the 1750-m-assemblage. fspace.1750m <- sp_faxes_coord_tactics[species.1750, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6")] ## Prune the community data set to include only species present in ## the 1750-m-assemblage. community.1750m <- occurrence.matrix[,species.1750] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.y <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1750m[-f,] c <- community.1750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.y <- sapply(1:62,fun.y) # 62 species present in the assemblage. volumes.y ## Extract FRic values from the list and obtain a matrix. indices.y <-seq (1,123,2) # 123 data frames in the list. newset.y <- NULL for (i in indices.y) { newset.y <- cbind(newset.y, volumes.y [[c(i,2)]]) } View(newset.y) matrix.volumes.y <- newset.y colnames(matrix.volumes.y)<- rownames(fspace.1750m) matrix.volumes.y <- matrix.volumes.y[8,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.y, species = species.y, fd_ind_values = fd_ind_values_tactics, fspace = fspace.1750m, community = community.1750m, target_value = fd_ind_values_tactics[9, 3], r = 8) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) # Niche packing/expansion contributions based on all traits ------------- ## The observed functional volume of the 1350-m-assemblage is lower than ## that of the 1550-m-assemblage, implying that ## all species in the species richer assemblage (1350 m) fall inside ## the functional volume of the species poorer assemblage ## (100% niche packing). The same holds for the comparison between the # 1750-m-assemblage and the 1950-m-assemblage. ## Comparison between the 550-m-assemblage and the 350-m-assemblage. com.350 <- as.matrix(occurrence.matrix[1,occurrence.matrix[1,]>0]) com.550 <- as.matrix(occurrence.matrix[2,occurrence.matrix[2,]>0]) ## Select species not shared between the two assemblages. species.350 <- rownames(com.350) species.550 <- rownames(com.550) com_com <- intersect(species.350,species.550) species <- species.550[!species.550 %in% com_com] ## Generate the functional space and extract PCoA scores ## for the 550-m-assemblage. fspace.550m <- sp_faxes_coord_all_traits[species.550, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8")] ## Prune the community data set to include only species present in ## the 550-m-assemblage. community.550m <- occurrence.matrix[,species.550] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.550m[-f,] c <- community.550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes <- sapply(1:139,fun) # 139 species present in the assemblage. volumes ## Extract FRic values from the list and obtain a matrix. indices <-seq (1,277,2) # 277 data frames in the list. newset <- NULL for (i in indices) { newset <- cbind(newset, volumes [[c(i,2)]]) } View(newset) matrix.volumes <- newset colnames(matrix.volumes)<- rownames(fspace.550m) matrix.volumes <- matrix.volumes[2,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes, species = species, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.550m, community = community.550m, target_value = fd_ind_values_all_traits[1, 3], r = 2) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 550-m-assemblage and the 750-m-assemblage. com.750 <- as.matrix(occurrence.matrix[3,occurrence.matrix[3,]>0]) com.550 <- as.matrix(occurrence.matrix[2,occurrence.matrix[2,]>0]) ## Select species not shared between the two assemblages. species.750 <- rownames(com.750) species.550 <- rownames(com.550) com_com.2 <- intersect(species.750,species.550) species.2 <- species.550[!species.550 %in% com_com.2] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes, species = species.2, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.550m, community = community.550m, target_value = fd_ind_values_all_traits[3, 3], r = 2) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 950-m-assemblage and the 1150-m-assemblage. com.1150 <- as.matrix(occurrence.matrix[5,occurrence.matrix[5,]>0]) com.950 <- as.matrix(occurrence.matrix[4,occurrence.matrix[4,]>0]) ## Select species not shared between the two assemblages. species.1150 <- rownames(com.1150) species.950 <- rownames(com.950) com_com.3 <- intersect(species.950,species.1150) species.3 <- species.950[!species.950 %in% com_com.3] ## Generate the functional space and extract PCoA scores ## for the 950-m-assemblage. fspace.950m <- sp_faxes_coord_all_traits[species.950, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6", "PC7", "PC8")] ## Prune the community data set to include only species present in ## the 950-m-assemblage. community.950m <- occurrence.matrix[,species.950] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.3 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.950m[-f,] c <- community.950m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.3 <- sapply(1:131,fun.3) # 131 species present in the assemblage. volumes.3 ## Extract FRic values from the list and obtain a matrix. indices.3 <-seq (1,261,2) # 261 data frames in the list. newset.3 <- NULL for (i in indices.3) { newset.3 <- cbind(newset.3, volumes.3 [[c(i,2)]]) } View(newset.3) matrix.volumes.3 <- newset.3 colnames(matrix.volumes.3)<- rownames(fspace.950m) matrix.volumes.3 <- matrix.volumes.3[4,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.3, species = species.3, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.950m, community = community.950m, target_value = fd_ind_values_all_traits[5, 3], r = 4) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1150-m-assemblage and the 1350-m-assemblage. com.1150 <- as.matrix(occurrence.matrix[5,occurrence.matrix[5,]>0]) com.1350 <- as.matrix(occurrence.matrix[6,occurrence.matrix[6,]>0]) ## Select species not shared between the two assemblages. species.1150 <- rownames(com.1150) species.1350 <- rownames(com.1350) com_com.4 <- intersect(species.1350,species.1150) species.4 <- species.1150[!species.1150 %in% com_com.4] ## Generate the functional space and extract PCoA scores ## for the 1150-m-assemblage. fspace.1150m <- sp_faxes_coord_all_traits[species.1150, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6", "PC7", "PC8")] ## Prune the community data set to include only species present in ## the 1150-m-assemblage. community.1150m <- occurrence.matrix[,species.1150] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.4 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1150m[-f,] c <- community.1150m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.4 <- sapply(1:105,fun.4) # 105 species present in the assemblage. volumes.4 ## Extract FRic values from the list and obtain a matrix. indices.4 <-seq (1,209,2) # 209 data frames in the list. newset.4 <- NULL for (i in indices.4) { newset.4 <- cbind(newset.4, volumes.4 [[c(i,2)]]) } View(newset.4) matrix.volumes.4 <- newset.4 colnames(matrix.volumes.4)<- rownames(fspace.1150m) matrix.volumes.4 <- matrix.volumes.4[5,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.4, species = species.4, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.1150m, community = community.1150m, target_value = fd_ind_values_all_traits[6, 3], r = 5) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1550-m-assemblage and the 1750-m-assemblage. com.1550 <- as.matrix(occurrence.matrix[7,occurrence.matrix[7,]>0]) com.1750 <- as.matrix(occurrence.matrix[8,occurrence.matrix[8,]>0]) ## Select species not shared between the two assemblages. species.1550 <- rownames(com.1550) species.1750 <- rownames(com.1750) com_com.5 <- intersect(species.1750,species.1550) species.5 <- species.1550[!species.1550 %in% com_com.5] ## Generate the functional space and extract PCoA scores ## for the 1550-m-assemblage. fspace.1550m <- sp_faxes_coord_all_traits[species.1550, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6", "PC7", "PC8")] ## Prune the community data set to include only species present in ## the 1550-m-assemblage. community.1550m <- occurrence.matrix[,species.1550] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.5 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1550m[-f,] c <- community.1550m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.5 <- sapply(1:84,fun.5) # 84 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.5 <-seq (1,167,2) # 167 data frames in the list. newset.5 <- NULL for (i in indices.5) { newset.5 <- cbind(newset.5, volumes.5 [[c(i,2)]]) } View(newset.5) matrix.volumes.5 <- newset.5 colnames(matrix.volumes.5)<- rownames(fspace.1550m) matrix.volumes.5 <- matrix.volumes.5[7,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.5, species = species.5, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.1550m, community = community.1550m, target_value = fd_ind_values_all_traits[8, 3], r = 7) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1750-m-assemblage and the 1950-m-assemblage. ## Select species not shared between the two assemblages. species.1950 <- rownames(com.1950) species.1750 <- rownames(com.1750) com_com.6 <- intersect(species.1950,species.1750) species.6 <- species.1750[!species.1750 %in% com_com.6] ## Generate the functional space and extract PCoA scores ## for the 1750-m-assemblage. fspace.1750m <- sp_faxes_coord_all_traits[species.1750, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6", "PC7", "PC8")] ## Prune the community data set to include only species present in ## the 1750-m-assemblage. community.1750m <- occurrence.matrix[,species.1750] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.6 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1750m[-f,] c <- community.1750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.6 <- sapply(1:62,fun.6) # 62 species in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.6 <-seq (1,123,2) # 123 data frames in the list. newset.6 <- NULL for (i in indices.6) { newset.6 <- cbind(newset.6, volumes.6 [[c(i,2)]]) } View(newset.6) matrix.volumes.6 <- newset.6 colnames(matrix.volumes.6)<- rownames(fspace.1750m) matrix.volumes.6 <- matrix.volumes.6[8,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.6, species = species.6, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.1750m, community = community.1750m, target_value = fd_ind_values_all_traits[9, 3], r = 8) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 1950-m-assemblage and the 2150-m-assemblage. com.1950 <- as.matrix(occurrence.matrix[9,occurrence.matrix[9,]>0]) com.2150 <- as.matrix(occurrence.matrix[10,occurrence.matrix[10,]>0]) ## Select species not shared between the two assemblages. species.1950 <- rownames(com.1950) species.2150 <- rownames(com.2150) com_com.7 <- intersect(species.1950,species.2150) species.7 <- species.1950[!species.1950 %in% com_com.7] ## Generate the functional space and extract PCoA scores ## for the 1950-m-assemblage. fspace.1950m <- sp_faxes_coord_all_traits[species.1950, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6", "PC7", "PC8")] ## Prune the community data set to include only species present in ## the 1950-m-assemblage. community.1950m <- occurrence.matrix[,species.1950] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.7 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.1950m[-f,] c <- community.1950m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.7 <- sapply(1:57,fun.7) # 57 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.7 <-seq (1,113,2) # 113 data frames in the list. newset.7 <- NULL for (i in indices.7) { newset.7 <- cbind(newset.7, volumes.7 [[c(i,2)]]) } View(newset.7) matrix.volumes.7 <- newset.7 colnames(matrix.volumes.7)<- rownames(fspace.1950m) matrix.volumes.7 <- matrix.volumes.7[9,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.7, species = species.7, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.1950m, community = community.1950m, target_value = fd_ind_values_all_traits[10, 3], r = 9) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing) ## Comparison between the 750-m-assemblage and the 950-m-assemblage. com.750 <- as.matrix(occurrence.matrix[3,occurrence.matrix[3,]>0]) com.950 <- as.matrix(occurrence.matrix[4,occurrence.matrix[4,]>0]) ## Select species not shared between the two assemblages. species.750 <- rownames(com.750) species.950 <- rownames(com.950) com_com.8 <- intersect(species.750,species.950) species.8 <- species.750[!species.750 %in% com_com.8] ## Generate the functional space and extract PCoA scores ## for the 1950-m-assemblage. fspace.750m <- sp_faxes_coord_all_traits[species.750, c("PC1", "PC2", "PC3","PC4", "PC5", "PC6", "PC7", "PC8")] ## Prune the community data set to include only species present in ## the 1950-m-assemblage. community.750m <- occurrence.matrix[,species.750] ## Function to remove sequentially each species from the species-rich assemblage ## and calculate the resulting functional volume. fun.8 <- function(x){ f<-0 for(i in 1:x){ f<- i t <- fspace.750m[-f,] c <- community.750m[,-f] v<- alpha.fd.multidim(sp_faxes_coord = t, asb_sp_w = c, ind_vect = "fric", scaling = F) } v } volumes.8 <- sapply(1:136,fun.8) # 136 species present in the assemblage. ## Extract FRic values from the list and obtain a matrix. indices.8 <-seq (1,271,2) # 271 data frames in the list. newset.8 <- NULL for (i in indices.8) { newset.8 <- cbind(newset.8, volumes.8 [[c(i,2)]]) } View(newset.8) matrix.volumes.8 <- newset.8 colnames(matrix.volumes.8)<- rownames(fspace.750m) matrix.volumes.8 <- matrix.volumes.8[3,] ## Estimating niche packing/expansion contribution in percentage. niche_metrics <- calculate_niche_contributions( matrix_volumes = matrix.volumes.8, species = species.8, fd_ind_values = fd_ind_values_all_traits, fspace = fspace.750m, community = community.750m, target_value = fd_ind_values_all_traits[4, 3], r = 3) print(niche_metrics$niche_expansion) print(niche_metrics$niche_packing)