[This article was first published on ** R Code – Geekcologist**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

*This post has been written by **André L. Luza **& **Vinicius A. G. Bastazini**

In the face of the current biodiversity crisis, understanding how ecological communities respond to species loss is more critical than ever. In a recent study, *“ Coping with Collapse: Functional Robustness of Coral-Reef Fish Network to Simulated Cascade Extinction“* (see also a previous post on this study) , we shed light on this issue by examining the cascading effects of species extinctions within coral-fish networks in the Southwestern Atlantic. In this post, we will show how to calculate functional network robustness in R, using the methods we developed in our study. The example below is just illustrative, and it is not intend to produce the same results as those reported in our manuscript (see also the GitHub tutorial).

The basis of the framework is the species-habitat network approach (Marini et al.2019; Bastazini et al. in press), however it can be directly applied to any type of bipartite networks formed by interacting species. The illustrated framework is depicted in the figure below. In step (1), we modeled the occupancy probability of fish as a function of coral and turf algae cover using site occupancy modeling. Based on the model output, we classified species as either coral-associated or co-occurring fish. In step (2), coral and fish species were connected based on the predicted site occupancy probability of each coral-associated fish (fish with different colors in the center of the network) relative to the cover of each coral species (corals with different colors, in the left). Then, coral-associated and co-occurring fish in subnetwork 2 were connected based on Pearson’s correlation values between fish site occupancy probability. Once the networks are constructed, we applied a removal algorithm that eliminated corals and subsequently calculated the direct and indirect effects of coral species removal on network robustness at each elimination step (t=0,t=1,…,|A|). Lost links are shown in red. In step (3), we related fish species composition and species traits at each elimination step. In step (4), we computed the loss in trait space area following corals and fish removal. The area delimited by the black polygon depicts the trait space area att=0, and the area delimited by the red polygon depicts the trait space area att=1.Finally, in step (5), we applied a hyperbolic function (non-linear model) to the simulated species elimination data, analyzing both the remaining taxonomic diversity (TD, represented on the first y-axis with a solid curve) and functional diversity (FD, represented on the second y-axis with a dashed curve) along the gradient of coral elimination (x-axis). Note that the approach is based on the analysis of two bipartite networks. Thus, the approach described here can be directly applied to assess the robustness of bipartite networks.

Analytical framework. Please see Luza et al. (2024) for further details

We will simulate some data to illustrate our approach. Consider a set of eight species in the Partite A, 42 species in the Partite B, and 21 species in the Partite C. The subnetwork 1 will be composed by the Partites A and B, and the subnetwork 2 will be composed by the Partites B and C. To fill matrix cells with values, we will gather values constrained between[0,1], resembling fish occupancy probability (ψ) (subnetwork 1) and co-occurrence (Pearson’s correlationρ) (subnetwork 2) from aBetadistribution with shape parametersa=0.5andb=1. This will produce a distribution of values high density close to zero.

set.seed(2456)nspA <-8 # sp in Partite AnspB<-42 # sp in Partite BnspC<-21 # sp in Partite Cnspp <- nspB+nspC # all species# create subnetwork 1subnetwork1 <- matrix(rbeta(nspA*nspB, shape1=.5, shape2 = 1),byrow=F, nrow=nspA, ncol=nspB)subnetwork1 [subnetwork1>1]<-1 # set to 1 values larger than 1# create subnetwork 2subnetwork2 <- matrix(rbeta(nspB*nspC,shape1=.5, shape2 = 1),byrow=F, nrow=nspB, ncol=nspC)subnetwork2 [subnetwork2>1]<-1 # set to 1 values larger than 1# set rownames() & colnames()# subnetwork 1rownames(subnetwork1) <- paste0 ("coral", seq(1,nspA))colnames(subnetwork1) <- paste0 ("fish", seq(1,nspB))# subnetwork 2rownames(subnetwork2) <- paste0 ("fish", seq(1,nspB))colnames(subnetwork2) <- paste0 ("fish", seq(nspB+1,nspB+nspC))# establish a threshold of rho = 0.8 for a coocurring speciessubnetwork2[subnetwork2<0.8]<-0# rm empty colssubnetwork2 <- subnetwork2 [, which(colSums(subnetwork2)>0)]# arrange the matrices# ordering (degree) --------------------require(bipartite); require(igraph)subnetwork1<-sortweb(subnetwork1,sort.order="dec")subnetwork2<-sortweb(subnetwork2,sort.order="dec")# define affinitiesrequire(reshape)subnetwork1_df <- melt(subnetwork1,as.is=T)subnetwork1_df$aff <- ifelse (subnetwork1_df$value > 0.9, 1,0)colnames(subnetwork1_df) <- c("coral", "fish", "value", "aff")

Note that in the last step we melted the object “subnetwork_1” and created the object “subnetwork1_df”. We did so to make this new object the base for filter data and to establish a binary variable depicting the relationship between species of partites A and B. Now that we produced the matrices and sorted them out, we can illustrate the tripartite network.

plotweb2(data.matrix(subnetwork1), data.matrix(subnetwork2), method = "normal", empty=T, col.interaction = "gray", ybig=1, labsize = 0.75, spacing=0.01, lab.space =0.1, method2="normal", spacing2=0.01, empty2=T, col.interaction2 = "orange" , col.pred2 = "orange", col.prey2 = "gray80")

As we already established the links between species in Partites B and C (fish in our empirical data), we will simulate one trait data set to be used in the trait-based analyzes. The trait values (n=6traits) of species in the Partites B and C were gathered from Normal distribution with averageμ=0and standard deviationσ=0.2.

# produce trait datatrait_data <- as.data.frame ( replicate (6, rnorm (n = nspp, mean = 0,sd = 0.2) ) )# set rownames() and colnames()rownames(trait_data) <- paste0 ("fish", seq(1,nspp))colnames(trait_data) <- paste0 ("trait", seq(1,ncol(trait_data)))

After getting the trait values, we can start the trait space analyses based on the Partites B and C. First we will calculate the distances between species based on their traits, and run one ordination analysis (Principal Coordinate Analysis, ‘dudi.pco’ function of the package*ade4*) to summarize these distances into different vectors. We will use the Gower distance (package*cluster*) as it enables us to handle several types of traits (continuous, categorical, binary traits).

# trait distance matrix (gower)require(cluster)gower_matrix <- daisy (trait_data, metric=("gower"))# Building the functional space based on a PCOA require(ade4)pco<-dudi.pco(quasieuclid(gower_matrix), scannf=F, nf=10) # quasieuclid() transformation to make the gower matrix as euclidean. nf= number of axis

Now we will build the complete trait space. The area of this polygon is the Functional Diversity-FD measure for(t=0), when no species elimination has taken place. We will apply the function ‘chull’ (package*grDevices*) to the data of the two first PCoA axes, which will link species located in the boundaries of the trait space. Then the function ‘Polygon’ (package*sp*) will be used to get the polygon, enabling the calculation of the polygon area (FD). These two functions will also be used to group species and obtain FD for other groups (e.g.coral-associated fish, co-occurring fish).

## complete trait spaceall_pool <- cbind (pco$li[,1:2], ext = F, sp = rownames(trait_data))# convex hulla_pool <- all_pool [chull(all_pool[,1:2], y = NULL),] # its convex hull

Then we can start the species removal process. Within the ‘lapply’ loop, there is a sequential selection of rows of the matrix ‘subnetwork1’. As this matrix was already sorted, the selection obeys the degree criterion. The ‘rm_corals’ object define which species in partite A will be selected (from the species highest to the lower degree). Then we selected the species of partite B which have affinity with the selected species (‘aff==1 & coral %in% rm_corals’), define the group of species (using ‘chull’), and calculate FD (using ‘Polygon’ and extracting the ‘@area’). After we create, the ‘RFS’ object, which contains the result of the division of the difference between the complete and simplified trait space area by the complete trait space area. If no trait space was lost thenRFS=1. As we were interested in how much remains of the trip space, we calculated1−RFSto be stored in our output ‘res’.

# --------------------------------------# removal based on degree of Partite A# calculate functional trait space loss# based on RFS (Luza et al. 2022)require(sp)RFS_corals <- lapply (seq (1, nspA), function (ncoral) { # choose the spp of Partite A to remove rm_corals <- rownames (subnetwork1)[1:ncoral] # coral-associated fish coral_associated <- subnetwork1_df[which(subnetwork1_df$aff == 1 & subnetwork1_df$coral %in% rm_corals ),] # reduced space setB<-cbind(all_pool, ext1=ifelse(all_pool$sp %in% unique(coral_associated$fish),T,F)) pk <-setB[which(setB$ext1==F),] f <- pk [chull(pk, y = NULL),] # hull # quantifying reduction in functional space # https://chitchatr.wordpress.com/2015/01/23/calculating-the-area-of-a-convex-hull/ chull.poly.complete <- Polygon(a_pool[,1:2], hole=F) # complete space chull.area.complete <- chull.poly.complete@area # polygon area # coral chull.poly.coral <- if (nrow (f)==0) {0} else {Polygon(f[,1:2], hole=F)} # reduced space chull.area.coral <- if (nrow (f)==0) {0} else {chull.poly.coral@area} # area # calculating reductions across all corals RFS<-data.frame (corals=(chull.area.complete-chull.area.coral)/chull.area.complete) # reduction in functional space area # results into a dataframe res <- list (remaining_associated_SR = 1-length(unique(coral_associated$fish))/nspB, # proportion of species remaining in partite B remaining_SR_all = 1-length(unique(coral_associated$fish))/nspp, # proportion of species remaining relative to the total RFS = 1-(RFS), # 1-RFS coral.associated = coral_associated$fish[order(coral_associated$fish)], # removed fish corals.removed = rm_corals, # removed corals space.coral = f # trait space to plot ) ; # return res})

By running this we will end up with a list with length equal to the number of species in partite A. This output comprise the simulated losses in partite B (coral-associated) produced by the direct losses of species in partite A (corals). Then we melt the list to have the relevant results.

# total functional lossFD_loss_total <- lapply (RFS_corals, function (i) i$RFS)FD_loss_total<- unlist(FD_loss_total)# loss species richness of partite B (coral associated fish)CA_loss_total <- lapply (RFS_corals, function (i) i$remaining_associated_SR)CA_loss_total<- do.call (rbind, CA_loss_total)# total SR loss ALL_loss_total <- lapply (RFS_corals, function (i) i$remaining_SR_all)ALL_loss_total<- do.call (rbind, ALL_loss_total)# removal of partite A species (Corals)loss_corals <- lapply (RFS_corals, function (i) 1-length(i$corals.removed)/nspA)loss_corals<- do.call (rbind, loss_corals)

We now proceed by creating a data frame for the robustness analyzes. First, we create data fort=0,where no loss has taken place yet. Subsequently, we bind to this data frame, the results of the direct species loss.

# analysis dataset# first step of no loss (t=0)analysis_dataset<- (data.frame (remain_all=1, remain_fish_associated=1, remain_coral=1, remain_FD = 1))# bind removalsanalysis_dataset<-rbind (analysis_dataset, data.frame (remain_all=ALL_loss_total, remain_fish_associated=CA_loss_total, remain_coral=loss_corals, remain_FD = FD_loss_total))# show the headinghead(analysis_dataset)## remain_all remain_fish_associated remain_coral remain_FD## 1 1.0000000 1.0000000 1.000 1.0000000## 2 0.9523810 0.9285714 0.875 1.0000000## 3 0.9047619 0.8571429 0.750 1.0000000## 4 0.8571429 0.7857143 0.625 0.8387945## 5 0.8412698 0.7619048 0.500 0.8387945## 6 0.8095238 0.7142857 0.375 0.8387945

Now we will simulate the influence of the indirect loss of corals. We start by selecting the species in partite A (corals) and B (coral-associated fish) to be removed. Then we select associated fish plus those associated with them by checking which species has the sum ofρ>0. These species will be removed and the remaining proportion of species richness (SR) and FD will be calculated.

# simulating secondary (indirect) extinctions -----------------RFS_corals_secondary_extinctions <- lapply (seq (1, nspA), function (ncoral) { # choose the spp of Partite A to remove rm_corals <- rownames (subnetwork1)[1:ncoral] # coral-associated fish coral_associated <- subnetwork1_df[which(subnetwork1_df$aff == 1 & subnetwork1_df$coral %in% rm_corals ),] # select associated fish plus those associated with them ass_cs <- subnetwork2[unique(coral_associated$fish),] ass_cs <- ass_cs[,colSums(ass_cs)>0] all_to_remove <- unlist(dimnames(ass_cs)) # reduced space setB<-cbind(all_pool, ext1=ifelse(all_pool$sp %in% all_to_remove,T,F)) pk <-setB[which(setB$ext1==F),] f <- pk [chull(pk, y = NULL),] # quantifying reduction in functional space # https://chitchatr.wordpress.com/2015/01/23/calculating-the-area-of-a-convex-hull/ chull.poly.complete <- Polygon(a_pool[,1:2], hole=F) # complete space chull.area.complete <- chull.poly.complete@area # area # coral chull.poly.coral <- if (nrow (f)==0) {0} else {Polygon(f[,1:2], hole=F)} # reduced space chull.area.coral <- if (nrow (f)==0) {0} else {chull.poly.coral@area} # area # calculating reductions across all corals RFS<-data.frame (corals=(chull.area.complete-chull.area.coral)/chull.area.complete) # reduction in functional space area # results res <- list (all_to_remove = all_to_remove, richness_remaining = length(all_to_remove)/nspp, richness_remaining_associated = length(unique(coral_associated$fish))/nspp, RFS = 1-RFS, coral.associated = coral_associated$fish[order(coral_associated$fish)], corals.removed = rm_corals, space.coral = f ) ; # return res})

As before, we selected the results we want, melt the lists, and bind to the data frame for robustness analyses.

# total lossloss_total_secondary <- lapply (RFS_corals_secondary_extinctions, function (i) i$RFS)loss_total_secondary<- do.call (rbind, loss_total_secondary)# richnessloss_total_secondary_SR <- lapply (RFS_corals_secondary_extinctions, function (i) i$richness_remaining)loss_total_secondary_SR<- do.call (rbind,loss_total_secondary_SR)# richness associated fishloss_total_secondary_SR_CA <- lapply (RFS_corals_secondary_extinctions, function (i) i$richness_remaining_associated)loss_total_secondary_SR_CA<- do.call (rbind,loss_total_secondary_SR_CA)# bind indirect extinctions in the analysis datasetanalysis_dataset<-cbind(analysis_dataset, data.frame (remain_RFS_secondary = c(1,unlist(loss_total_secondary)), remain_SR_secondary = c(1,1-unlist (loss_total_secondary_SR)), remain_SR_secondary_associated = c(1,1-unlist (loss_total_secondary_SR_CA)), ext.lower=1, no = seq(0,nspA)/nspA ) )

Next we present the function and fit the non-linear model to the data. The object that the function requires is a data frame with the proportion of species remaining in partite A and the remaining SR or FD in the other partites. The function also have a function to plot.

# run function of hyperbolic curve -------------------------------------------------# the function# Fit hyperbolic functionfit.hyperbolica<-function (object, plot.it = T, ...) { #if (class(object) != "bipartite") # stop("This function cannot be meaningfully applied to objects of this class.") N <- colSums(object) #if (all(object[-nrow(object), 2] == 1)) { # y <- -object[, 3]} else {y <- -object[, 2]} y <- object[,3]#(sum(y) - cumsum(y))/sum(y) x <- (object[, "no"]/max(object[, "no"])) fit <- try(nls(y ~ 1 - x^a, start = list(a = 0.5))) # originally set to 1, but used here 0.5 to make the example work if (class(fit) == "try-error") fit <- nls((y + rnorm(length(y), s = 0.001)) ~ 1 - x^a, start = list(a = 1)) # option to plot if (plot.it) { par(mar = c(5, 5, 1, 1)) plot(x, y, xlab = "Fraction of eliminated corals", ylab = "Fraction of surviving fish", axes = TRUE, type = "n", cex.lab = 1) points(x, y)#, ...) lines(seq(0, 1, length=9), predict(fit, newdata = data.frame(x = seq(0, 1, length=9))), col = "red", lwd = 2) } return(list(exponent = as.numeric(coef(fit)[1]), x=x, y=y, model=fit, preds=predict(fit, newdata = data.frame(x = seq(0, 1, length=9)))))}par(mfrow=c(2,2))# run functionhyper_curve_secondary <- fit.hyperbolica (analysis_dataset[,c("no", "ext.lower", "remain_all")], plot.it = F)# you can plot if you wanthyper_curve_associated <- fit.hyperbolica (analysis_dataset[,c("no", "ext.lower", "remain_SR_secondary")], plot.it = F)hyper_curve_associated_RFS <- fit.hyperbolica (analysis_dataset[,c("no", "ext.lower", "remain_FD")], plot.it = F)hyper_curve_secondary_RFS <- fit.hyperbolica (analysis_dataset[,c("no", "ext.lower", "remain_RFS_secondary")], plot.it = F)

After fitting the model to the data, we bind the predictions to the analysis dataset, and estimated network robustness by integrating (summing up) infinitesimally small values of the spline interpolated using the fitted hyperbolic function applied to the minimum and maximum of the proportion of removed partite A species/corals.

# predictionsanalysis_dataset$pred_remain_SR_all <- hyper_curve_secondary$predsanalysis_dataset$pred_remain_SR_secondary <- hyper_curve_associated$predsanalysis_dataset$pred_remain_RFS_associated <- hyper_curve_associated_RFS$predsanalysis_dataset$pred_remain_RFS_secondary <- hyper_curve_secondary_RFS$preds# robustness# srSR1 <- (integrate(splinefun(hyper_curve_secondary_RFS$x, analysis_dataset$pred_remain_SR_all), min(hyper_curve_secondary_RFS$x), max(hyper_curve_secondary_RFS$x), subdivisions = max(100L, length(hyper_curve_secondary_RFS$x))))# SR secondarySR2 <- (integrate(splinefun(hyper_curve_secondary_RFS$x, analysis_dataset$pred_remain_SR_secondary), min(hyper_curve_secondary_RFS$x), max(hyper_curve_secondary_RFS$x), subdivisions = max(100L, length(hyper_curve_secondary_RFS$x))))# FD cora lassociatedFD1 <- (integrate(splinefun(hyper_curve_secondary_RFS$x, analysis_dataset$pred_remain_RFS_associated), min(hyper_curve_secondary_RFS$x), max(hyper_curve_secondary_RFS$x), subdivisions = max(100L, length(hyper_curve_secondary_RFS$x))))# fd otehrFD2 <- (integrate(splinefun(hyper_curve_secondary_RFS$x, analysis_dataset$pred_remain_RFS_secondary), min(hyper_curve_secondary_RFS$x), max(hyper_curve_secondary_RFS$x), subdivisions = max(100L, length(hyper_curve_secondary_RFS$x))))# create a table to show robustness resultsrequire(dplyr); require(knitr)data.frame(rbind(SR_direct = SR1$value, SR_indirect = SR2$value, FD_direct = FD1$value, FD_indirect = FD2$value )) %>% kable(format = "pipe",col.names = "Robustness (R)")

Finally, we can now plot the robustness curves for direct and indirect effects of the removal of partite A species on the other partities.

# plotrequire(ggplot2)# plotggplot(analysis_dataset[order(analysis_dataset$remain_coral,decreasing=F),], aes (x= (1-remain_coral), y = remain_fish_associated)) + # functional diversity # direct loss geom_line( aes(y=pred_remain_RFS_associated),linetype=1,size=1.2,col="orange") + geom_ribbon( aes(ymax=pred_remain_RFS_associated), ymin=0, linetype=1, size=1.2, col = NA, alpha=0.3, fill="orange") + geom_point( aes(y=remain_FD), size=3, stroke=2, col="orange", fill="white", shape=23, alpha=1) + # indirect loss geom_line( aes(y=pred_remain_RFS_secondary),linetype=1,size=1.2, col = "orange3") + geom_ribbon( aes(ymax=pred_remain_RFS_secondary), ymin=0, linetype=1, size=1.2, col = NA, alpha=0.4, fill="orange3") + geom_point( aes(y=remain_RFS_secondary),size=3, stroke=2, col="orange3", fill="white", shape=21, alpha=1) + # taxonomic diversity geom_ribbon( aes(ymax=pred_remain_SR_all), ymin=0, linetype=1, size=1.2, col = NA, alpha=0.3, fill="gray50") + geom_line( aes(y=pred_remain_SR_all), linetype=1,size=1.2,col="gray50") + geom_point( aes(y=remain_all),size=3, stroke=2, col="gray50", fill="white", shape=23, alpha=1) + #geom_line( aes(y=remain_fish_associated*100),linetype=1,size=1.2,col="orange") + geom_line( aes(y=pred_remain_SR_secondary),linetype=1,size=1.2,col="black") + geom_ribbon( aes(ymax=pred_remain_SR_secondary), ymin=0, linetype=1, size=1.2, col = NA, alpha=0.8, fill="black") + geom_point( aes(y=remain_SR_secondary), size=3, stroke=2, col="black", fill="white", shape=21, alpha=1) + scale_y_continuous( # Features of the first axis name = "Remaining Fish Taxonomic Diversity", # Add a second axis and specify its features sec.axis = sec_axis(trans=~.*1, name="Remaining Fish Functional Diversity") ) + xlab ("Proportion of corals eliminated") + theme_bw() + theme (axis.title=element_text(size=14), axis.text = element_text(size=12), plot.title = element_text(size=18))

Attack Tolerance Curves (ATC). The shaded area below each curve depicts the hyperbolic function curve fitted to the data shown in the two Y-axes. Yellow tones represent losses in functional diversity, while grey and black represent losses in taxonomic diversity. Diamonds and circles denote the direct and indirect effects of coral species extinction on such biodiversity dimensions.

**References**

Bastazini VAG, Gianuca AT, Vizentin-Bugoni J, Gonçalves MSS & Dias RA . *In press*. Identificando o potencial de conservação de áreas úmidas usando ferramentas da teoria de redes complexas. In: Gonçalves, MSS,Bastazini VAG, Andretti C, Lanés LE, Volcan M (orgs).*Biodiversidade e conservação de ecossistemas aquáticos do sul do Brasil*. Editora USEB.

Luza AL, Bender MG, Ferreira CEL, Floeter SR, Francini-Filho RB, Longo GO, Pinheiro HT, Quimbayo JP,Bastazini VAG**.**2024. Coping with collapse: *Functional robustness of coral-reef fish network to simulated cascade extinction*. Global Change Biology 30: e17513. https://doi.org/10.1111/gcb.17513

Marini L, Bartomeus I, Rader R, Lami F. 2019. *Species–habitat networks: A tool to improve landscape management for conservation*. Journal of Applied Ecology 56: 923–928.https://doi.org/10.1111/1365-2664.13337

*Related*

To **leave a comment** for the author, please follow the link and comment on their blog: ** R Code – Geekcologist**.

R-bloggers.com offers **daily e-mail updates** about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.