## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = requireNamespace("vegan", quietly = TRUE) ) ## ----vegan, eval = !requireNamespace("vegan", quietly = TRUE), echo = FALSE, comment = NA---- # message('This vignette requires the "vegan" package. Please, install it: install.packages("vegan").') ## ----draw_arc, echo=FALSE----------------------------------------------------- draw_arc <- function(x, y, r, ini, fin, ...){ ini_rad <- ini*pi/180 fin_rad <- fin*pi/180 angle <- seq(ini_rad, fin_rad, by = 0.01) arc_x <- x + r * cos(angle) arc_y <- y + r * sin(angle) lines(arc_x, arc_y, ...) } ## ----install, eval=FALSE------------------------------------------------------ # install.packages("ecoregime") # devtools::install_github(repo = "MSPinillos/ecoregime", dependencies = T, build_vignettes = T) ## ----setup-------------------------------------------------------------------- library(ecoregime) ## ----citation----------------------------------------------------------------- citation("ecoregime") ## ----EDR_data----------------------------------------------------------------- abun_EDR1 <- EDR_data$EDR1$abundance abun_EDR2 <- EDR_data$EDR2$abundance abun_EDR3 <- EDR_data$EDR3$abundance ## ----dist_traj---------------------------------------------------------------- head(EDR_data$EDR3_disturbed$abundance) ## ----st_landscape, fig = TRUE, fig.height=5, fig.width=6---------------------- library(data.table) # Include all abundances in one matrix abun_all <- rbindlist(list( EDR_data$EDR1$abundance, EDR_data$EDR2$abundance, EDR_data$EDR3$abundance, EDR_data$EDR3_disturbed$abundance), fill = TRUE) # Compute state dissimilarities d_all <- vegan::vegdist(abun_all[, paste0("sp", 1:12)], method = "bray") # We can apply mMDS to the dissimilarity matrix mds_all <- smacof::mds(d_all, ndim = 3) # Visualize the stability landscape plot_edr(x = mds_all$conf, trajectories = paste0(abun_all$EDR, "_", abun_all$traj), # provide a unique ID to each trajectory states = as.integer(abun_all$state), traj.colors = c(rep(palette.colors(4, "Table")[c(1:2, 4)], each= 30), rep(palette.colors(4, "Table")[3], 3)), xlab = "MDS D1", ylab = "MDS D2", main = "Stability landscape and disturbed trajectories") legend("bottomleft", c(paste0("EDR", 1:3), "Disturbed trajectories"), lty = 1, lwd = 2, col = palette.colors(4, "Table")[c(1,2,4,3)], bty = "n", cex = 0.8) ## ----dDis--------------------------------------------------------------------- # Species abundances in pre-disturbance states of the disturbed trajectories # (Pre-disturbance states are identified by disturbed_states = 0) abun_predist <- EDR_data$EDR3_disturbed$abundance[disturbed_states == 0] selcols <- names(EDR_data$EDR1$abundance) ## EDR1 ------------------------------------------------------------------------ # Species abundances in EDR1 and the predisturbed states of disturbed trajectories abun1_predist <- rbind(EDR_data$EDR1$abundance, abun_predist[, ..selcols]) # State dissimilarities in EDR1 and the predisturbed states of disturbed trajectories d1_predist <- vegan::vegdist(x = abun1_predist[, paste0("sp", 1:12)], method = "bray") # dDis of the disturbed trajectories in relation to EDR1 dDis1 <- sapply(unique(abun_predist$traj), function(ipredist){ dDis(d = d1_predist, d.type = "dStates", trajectories = abun1_predist$traj, states = abun1_predist$state, reference = as.character(ipredist)) }) ## EDR2 ------------------------------------------------------------------------ # Species abundances in EDR2 and the predisturbed states of disturbed trajectories abun2_predist <- rbind(EDR_data$EDR2$abundance, abun_predist[, ..selcols]) # State dissimilarities in EDR2 and the predisturbed states of disturbed trajectories d2_predist <- vegan::vegdist(x = abun2_predist[, paste0("sp", 1:12)], method = "bray") # dDis of the disturbed trajectories in relation to EDR2 dDis2 <- sapply(unique(abun_predist$traj), function(ipredist){ dDis(d = d2_predist, d.type = "dStates", trajectories = abun2_predist$traj, states = abun2_predist$state, reference = as.character(ipredist)) }) ## EDR3 ------------------------------------------------------------------------ # Species abundances in EDR3 and the predisturbed states of disturbed trajectories abun3_predist <- rbind(EDR_data$EDR3$abundance, abun_predist[, ..selcols]) # State dissimilarities in EDR3 and the predisturbed states of disturbed trajectories d3_predist <- vegan::vegdist(x = abun3_predist[, paste0("sp", 1:12)], method = "bray") # dDis of the disturbed trajectories in relation to EDR3 dDis3 <- sapply(unique(abun_predist$traj), function(ipredist){ dDis(d = d3_predist, d.type = "dStates", trajectories = abun3_predist$traj, states = abun3_predist$state, reference = as.character(ipredist)) }) ## Compare dynamic dispersion -------------------------------------------------- # Compare dDis values for the three EDRs dDis_df <- data.frame(EDR1 = dDis1, EDR2 = dDis2, EDR3 = dDis3) ## ----dDis_comparison, echo=FALSE---------------------------------------------- knitr::kable(dDis_df, digits = 3) ## ----EDR_dist, fig = TRUE, fig.height=5, fig.width=6-------------------------- # Index of the states in the reference EDR iref <- which(abun_all$EDR == 3 & is.na(abun_all$disturbed_states)) # Index of the states in disturbed trajectories idist <- which(!is.na(abun_all$disturbed_states)) # Index of the states in the reference EDR and disturbed trajectories iref_idist <- c(iref, idist) # Filter the dissimilarity matrix to select the states of the reference EDR and # the disturbed trajectories d_all <- as.matrix(d_all) d_EDR3_dist <- d_all[iref_idist, iref_idist] # Apply mMDS for the new matrix mds_EDR3_dist <- smacof::mds(d_EDR3_dist, ndim = 3) # Vector of colors state.colors <- abun_all[iref_idist, # color for reference states state.colors := ifelse(is.na(disturbed_states), palette.colors(8, "R4")[8], # color for pre-disturbance states ifelse(disturbed_states == 0, palette.colors(8, "R4")[7], # color for disturbed states ifelse(traj == 31, palette.colors(8, "R4")[2], ifelse(traj == 32, palette.colors(8, "R4")[3], palette.colors(8, "R4")[4]))))]$state.colors state.colors <- state.colors[!is.na(state.colors)] # Plot the reference EDR and disturbed trajectories plot_edr(x = mds_EDR3_dist$conf, trajectories = paste0(abun_all$EDR[iref_idist], "_", abun_all$traj[iref_idist]), states = as.integer(abun_all$state[iref_idist]), type = "states", # it allows assigning different colors to the states state.colors = state.colors, xlab = "MDS D1", ylab = "MDS D2", main = "Reference EDR and disturbed trajectories") legend("topleft", c("Reference EDR", "Pre-disturbed states", "Disturbed states T31", "Disturbed states T32", "Disturbed states T33"), lty = 1, lwd = 2, col = palette.colors(8, "R4")[c(8, 7, 2:4)], cex = 0.9, bty = "n") ## ----retra-------------------------------------------------------------------- # State dissimilarities for EDR3 (considering only the undisturbed trajectories) d_EDR3 <- vegan::vegdist(EDR_data$EDR3$abundance[, paste0("sp", 1:12)]) # Representative trajectories retra <- retra_edr(d = d_EDR3, trajectories = EDR_data$EDR3$abundance$traj, states = EDR_data$EDR3$abundance$state, minSegs = 5) ## ----summary_retra------------------------------------------------------------ # Summarize retra summary(retra) # Define T4 as the unique representative trajectory and generate an object of class 'RETRA' retra_ref <- define_retra(data = retra$T4$Segments, d = d_EDR3, trajectories = EDR_data$EDR3$abundance$traj, states = EDR_data$EDR3$abundance$state, retra = retra) ## ----plot_retra, fig.dim = c(6,5)--------------------------------------------- # Plot EDR3 and its representative trajectories plot(retra, d = d_EDR3, trajectories = EDR_data$EDR3$abundance$traj, states = EDR_data$EDR3$abundance$state, select_RT = "T4", xlab = "MDS D1", ylab = "MDS D2", main = "Representative trajectories in EDR3") legend("topleft", c("Representative trajectory 'T4'", "Other representative trajectories", "Individual trajectories in EDR3"), lty = 1, col = c("red", "black", "grey"), bty = "n") ## ----graphical_indices, echo=FALSE, fig.dim=c(6,4.8)-------------------------- par(mar = c(0,0,0,0)) # Representative trajectory plot(x = 0:10, y = c(0,0,0:8), type = "n", axes = F, xlab = "", ylab = "") for (i in 0:9) { shape::Arrows(x0 = i, y0 = 1, x1 = i + 1, y1 = 1, lwd = 2, arr.adj = 1) } text(x = 0.2, y = 1.3, "RT") # Reference lines(x = c(0, 10), y = c(2.5, 2.5), col = "#E41A1C", lty = 3, lwd = 2) # Disturbed trajectory for (i in 0.5:1.5) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 2.5, y0 = 2.5, x1 = 4, y1 = 7, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 4, y0 = 7, x1 = 5, y1 = 6, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 5, y0 = 6, x1 = 5.5, y1 = 4, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 5.5, y0 = 4, x1 = 6.5, y1 = 3.5, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 6.5, y0 = 3.5, x1 = 7.5, y1 = 3.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 0.7, y = 2.8, "DT", col = "#E41A1C") # Resistance lines(x = c(2.5, 4)-0.2, y = c(2.5, 7), col = "#377EB8") text(x = 2.8, y = 5, "Rt", col = "#377EB8") # Amplitude draw_arc(x = 2.5, y = 2.5, r = 0.7, ini = 0, fin = 70, lwd = 2, col = "#4DAF4A") text(x = 3.3, y = 3.1, "A", col = "#4DAF4A") # Recovery lines(x = c(4, 5), y = c(7, 7), lty = 3, lwd = 2, col = "#984EA3") draw_arc(x = 4, y = 7, r = 0.6, ini = -45, fin = 0, lwd = 2, col = "#984EA3") text(x = 4.9, y = 6.7, "Rc", col = "#984EA3") # Net change lines(x = c(2.5, 6.5), y = c(2.5, 3.5), lty = 3, lwd = 2, col = "#999999") draw_arc(x = 2.5, y = 2.5, r = 1.5, ini = 0, fin = 15, lwd = 2, col = "#999999") text(x = 4.3, y = 2.7, "NC", col = "#999999") # Legend legend("topright", c("RT: Representative trajectory", "DT: Disturbed trajectory", "Rt: Resistance", "A: Amplitude", "Rc: Recovery", "NC: Net change"), col = c("black", "#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#999999"), lty = 1, text.col = c("black", "#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#999999"), bty = "n") # States text(x = c(2.5, 4, 5.2, 5.3, 6.5, 7.5), y = c(2.2, 7.4, 6, 4, 3.8, 3.8), 0:5, cex = 0.9, col = "#E41A1C") ## ----graphical_resistance, echo=FALSE, fig.dim=c(6,4.8)----------------------- par(mar = c(0,0,0,0)) # Representative trajectory plot(x = c(0, 10), y = c(0, 4.5), type = "n", axes = F, xlab = "", ylab = "") for (i in 0:9) { shape::Arrows(x0 = i, y0 = 0, x1 = i + 1, y1 = 0, lwd = 2, arr.adj = 1) } text(x = 0.2, y = 0.2, "RT") # Reference lines(x = c(0, 10), y = c(2.5, 2.5), col = "#E41A1C", lty = 3, lwd = 2) # Disturbed trajectories for (i in 5.5:6.5) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 7.5, y0 = 2.5, x1 = 8.5, y1 = 4, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 5.7, y = 2.7, "DT2", col = "#E41A1C") for (i in 0.5:1.5) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 2.5, y0 = 2.5, x1 = 8, y1 = 1, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 0.7, y = 2.7, "DT1", col = "#E41A1C") # Resistance lines(x = c(7.5, 8.5)-0.2, y = c(2.5, 4)+0.05, lwd = 2, col = "#377EB8") text(x = 7.2, y = 3.4, "High Rt", col = "#377EB8") lines(x = c(2.5, 8)-0.1, y = c(2.5, 1)-0.2, lwd = 2, col = "#377EB8") text(x = 5, y = 1.3, "Low Rt", col = "#377EB8") # Legend legend("topleft", c("RT: Representative trajectory", "DT: Disturbed trajectories", "Rt: Resistance"), col = c("black", "#E41A1C", "#377EB8"), lty = 1, text.col = c("black", "#E41A1C", "#377EB8"), bty = "n") # States text(x = c(2.5, 7.5, 8.7, 8.3), y = c(2.7, 2.3, 4.2, 1), c(0, 0, 1, 1), cex = 0.9, col = "#E41A1C") ## ----resistance--------------------------------------------------------------- # To calculate resistance, we need a state dissimilarity matrix for the disturbed trajectories d_disturbed <- vegan::vegdist(EDR_data$EDR3_disturbed$abundance[, paste0("sp", 1:12)], method = "bray") # Compute resistance # Note that the disturbed states are identified by disturbed_states = 1 Rt <- resistance(d = d_disturbed, trajectories = EDR_data$EDR3_disturbed$abundance$traj, states = EDR_data$EDR3_disturbed$abundance$state, disturbed_trajectories = unique(EDR_data$EDR3_disturbed$abundance$traj), disturbed_states = EDR_data$EDR3_disturbed$abundance[disturbed_states == 1]$state) ## ----Rt_results, echo=FALSE--------------------------------------------------- knitr::kable(Rt, row.names = F, col.names = c("Disturbed trajectories", "Rt"), digits = 3) ## ----graphical_amplitude, echo=FALSE, fig.dim=c(6,3.8)------------------------ par(mar = c(0,0,0,0)) # Representative trajectory plot(x = c(0, 16), y = c(0, 8), type = "n", axes = F, xlab = "", ylab = "") for (i in 0:14) { shape::Arrows(x0 = i, y0 = 0, x1 = i + 1, y1 = 0, lwd = 2, arr.adj = 1) } text(x = 0.2, y = 0.3, "RT") # Reference lines(x = c(0, 15), y = c(2.5, 2.5), col = "#E41A1C", lty = 3, lwd = 2) # Disturbed trajectory 1 for (i in 0.5:1.5) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 2.5, y0 = 2.5, x1 = 4, y1 = 0.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 0.7, y = 2.8, "DT1", col = "#E41A1C") # Disturbed trajectory 2 for (i in 5:6) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 7, y0 = 2.5, x1 = 14, y1 = 4.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 5, y = 2.8, "DT2", col = "#E41A1C") # # Amplitude lines(x = c(4, 4), y = c(0.5, 2.5), lwd = 2, col = "#4DAF4A") text(x = 5.05, y = 1.5, expression(A[abs1] ~ "< 0"), col = "#4DAF4A") # plotrix::draw.arc(x = 2.5, y = 2.5, radius = 0.5, deg2 = -70, lwd = 2, col = "#4DAF4A") draw_arc(x = 2.6, y = 2.5, r = 0.4, ini = -60, fin = 0, lwd = 2, col = "#4DAF4A") text(x = 3.85, y = 2.2, expression(A[rel1] ~ "< 0"), col = "#4DAF4A") lines(x = c(14, 14), y = c(2.5, 4.5), lwd = 2, col = "#4DAF4A") text(x = 15.05, y = 3.5, expression(A[abs2] ~ "> 0"), col = "#4DAF4A") # plotrix::draw.arc(x = 7, y = 2.5, radius = 1.5, deg2 = 32, lwd = 2, col = "#4DAF4A") draw_arc(x = 7, y = 2.5, r = 1.5, ini = 0, fin = 17, lwd = 2, col = "#4DAF4A") text(x = 9.5, y = 2.75, expression(A[rel2] ~ "> 0"), col = "#4DAF4A") text(x = c(0), y = c(6), adj = 0, expression("|"~A[abs1]~"| = |"~A[abs2]~"|"), cex = 0.9, col = "#4DAF4A") text(x = c(0), y = c(5.5), adj = 0, expression("|"~A[rel1]~"| > |"~A[rel2]~"|"), cex = 0.9, col = "#4DAF4A") # Legend legend("topleft", c("RT: Representative trajectory", "DT: Disturbed trajectories", "A: Amplitude"), col = c("black", "#E41A1C", "#4DAF4A"), lty = 1, text.col = c("black", "#E41A1C", "#4DAF4A"), bty = "n") # States text(x = c(2.5, 7, 4, 14), y = c(2.7, 2.3, 0.3, 4.7), c(0, 0, 1, 1), cex = 0.9, col = "#E41A1C") ## ----amplitude---------------------------------------------------------------- # We need a state dissimilarity matrix containing the states of the disturbed # trajectories and the representative trajectory taken as the reference abun <- rbind(EDR_data$EDR3$abundance, EDR_data$EDR3_disturbed$abundance, fill = T) d <- vegan::vegdist(abun[, paste0("sp", 1:12)], method = "bray") # Compute amplitude A <- amplitude(d = d, trajectories = abun$traj, states = abun$state, disturbed_trajectories = abun[disturbed_states == 1]$traj, disturbed_states = abun[disturbed_states == 1]$state, reference = retra_ref, method = "nearest_state") ## ----A_results, echo=FALSE---------------------------------------------------- knitr::kable(A, col.names = c("Disturbed trajectories", "Reference", "A~abs~", "A~rel~"), digits = 3) ## ----graphical_recovery, echo=FALSE, fig.dim=c(6,3.8)------------------------- par(mar = c(0,0,0,0)) # Representative trajectory plot(x = c(0, 16), y = c(0, 8), type = "n", axes = F, xlab = "", ylab = "") for (i in 0:14) { shape::Arrows(x0 = i, y0 = 1, x1 = i + 1, y1 = 1, lwd = 2, arr.adj = 1) } text(x = 0.2, y = 1.3, "RT") # Reference lines(x = c(0, 15), y = c(2.5, 2.5), col = "#E41A1C", lty = 3, lwd = 2) # Disturbed trajectory 1 for (i in 0.5:1.5) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 2.5, y0 = 2.5, x1 = 4, y1 = 4.5, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 4, y0 = 4.5, x1 = 5, y1 = 3.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 0.7, y = 2.8, "DT1", col = "#E41A1C") # Disturbed trajectory 2 for (i in 6:7) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 8, y0 = 2.5, x1 = 9.5, y1 = 4.5, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 9.5, y0 = 4.5, x1 = 14, y1 = 3.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 6.2, y = 2.8, "DT2", col = "#E41A1C") # States text(x = c(2.5, 3.5, 5, 8, 9.5, 14), y = c(2.3, 4.5, 3.3, 2.3, 4.7, 3.3), c(0:2, 0:2), cex = 0.9, col = "#E41A1C") # Recovery DT1 lines(x = c(4, 5), y = c(4.5, 4.5), lty = 3, col = "#984EA3") lines(x = c(5, 5), y = c(3.5, 4.5), lwd = 2, col = "#984EA3") text(x = 5.2, y = 4, adj = 0, expression(Rc[abs1] ~ "> 0"), col = "#984EA3") draw_arc(x = 4, y = 4.5, r = 0.5, ini = -45, fin = 0, lwd = 2, col = "#984EA3") text(x = 4, y = 4.7, adj = 0, expression(Rc[rel1] ~ "> 0"), col = "#984EA3") # Recovery DT2 lines(x = c(9.5, 14), y = c(4.5, 4.5), lty = 3, col = "#984EA3") lines(x = c(14, 14), y = c(3.5, 4.5), lwd = 2, col = "#984EA3") text(x = 15.3, y = 4, expression(Rc[abs2] ~ "> 0"), col = "#984EA3") draw_arc(x = 9.5, y = 4.5, r = 1, ini = -13, fin = 0, lwd = 2, col = "#984EA3") text(x = 11, y = 4.3, adj = 0, expression(Rc[rel2] ~ "> 0"), col = "#984EA3") text(x = c(0), y = c(6), adj = 0, expression(Rc[abs1]~" = "~Rc[abs2]), cex = 0.9, col = "#984EA3") text(x = c(0), y = c(5.5), adj = 0, expression(Rc[rel1]~" > "~Rc[rel2]), cex = 0.9, col = "#984EA3") # Legend legend("topleft", c("RT: Representative trajectory", "DT: Disturbed trajectories", "Rc: Recovery"), col = c("black", "#E41A1C", "#984EA3"), lty = 1, text.col = c("black", "#E41A1C", "#984EA3"), bty = "n") ## ----graphical_recovery2, echo=FALSE, fig.dim=c(6,3.8)------------------------ par(mar = c(0,0,0,0)) # Representative trajectory plot(x = c(0, 16), y = c(0, 8), type = "n", axes = F, xlab = "", ylab = "") for (i in 0:14) { shape::Arrows(x0 = i, y0 = 0, x1 = i + 1, y1 = 0, lwd = 2, arr.adj = 1) } text(x = 0.2, y = 0.3, "RT") # Reference lines(x = c(0, 15), y = c(2.5, 2.5), col = "#E41A1C", lty = 3, lwd = 2) # Disturbed trajectory 3 for (i in 0.5:1.5) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 2.5, y0 = 2.5, x1 = 4, y1 = 0.5, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 4, y0 = 0.5, x1 = 5, y1 = 1.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 0.7, y = 2.8, "DT3", col = "#E41A1C") # Disturbed trajectory 4 for (i in 6:7) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 8, y0 = 2.5, x1 = 9.5, y1 = 4.5, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 9.5, y0 = 4.5, x1 = 14, y1 = 5.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 6.2, y = 2.8, "DT4", col = "#E41A1C") # States text(x = c(2.5, 3.4, 5.1, 8, 9.2, 14.2), y = c(2.7, 0.6, 1.7, 2.3, 4.5, 5.6), c(0:2, 0:2), cex = 0.9, col = "#E41A1C") # Recovery DT3 lines(x = c(4, 5), y = c(0.5, 0.5), lty = 3, col = "#984EA3") lines(x = c(5, 5), y = c(0.5, 1.5), lwd = 2, col = "#984EA3") text(x = 5.2, y = 1, adj = 0, expression(Rc[abs3] ~ "< 0"), col = "#984EA3") draw_arc(x = 4, y = 0.5, r = 0.5, ini = 0, fin = 45, lwd = 2, col = "#984EA3") text(x = 3.8, y = 0.3, adj = 0, expression(Rc[rel3] ~ "< 0"), col = "#984EA3") # Recovery DT4 lines(x = c(9.5, 14), y = c(4.5, 4.5), lty = 3, col = "#984EA3") lines(x = c(14, 14), y = c(5.5, 4.5), lwd = 2, col = "#984EA3") text(x = 15.3, y = 5, expression(Rc[abs4] ~ "< 0"), col = "#984EA3") draw_arc(x = 9.5, y = 4.5, r = 1.3, ini = 0, fin = 13, lwd = 2, col = "#984EA3") text(x = 11.1, y = 4.7, adj = 0, expression(Rc[rel4] ~ "< 0"), col = "#984EA3") text(x = c(0), y = c(6.5), adj = 0, expression(Rc[abs3]~" = "~Rc[abs4]), cex = 0.9, col = "#984EA3") text(x = c(0), y = c(6), adj = 0, expression(Rc[rel3]~" < "~Rc[rel4]~" < "~Rc[rel2]~" < "~Rc[rel1]), cex = 0.9, col = "#984EA3") ## ----recovery----------------------------------------------------------------- # Compute recovery using the same data used for amplitude Rc <- recovery(d = d, trajectories = abun$traj, states = abun$state, disturbed_trajectories = abun[disturbed_states == 1]$traj, disturbed_states = abun[disturbed_states == 1]$state, reference = retra_ref, method = "nearest_state") ## ----plot_recovery, fig.height=5, fig.width=10-------------------------------- # Number of states after the disturbed state Rc <- data.table::data.table(Rc) Rc[, ID_post := 1:(.N), by = disturbed_trajectories] par(mfrow = c(1, 2)) # Plot absolute recovery over time plot(x = range(Rc$ID_post), y = range(Rc$Rc_abs), type = "n", xlab = "Nb. states after disturbance", ylab = "Absolute recovery", main = "Variation of absolute recovery") for (i in unique(Rc$disturbed_trajectories)) { lines(Rc[disturbed_trajectories == i, c("ID_post", "Rc_abs")], col = which(unique(Rc$disturbed_trajectories) %in% i) + 1) } legend("bottomleft", legend = unique(Rc$disturbed_trajectories), lty = 1, col = seq_along(unique(Rc$disturbed_trajectories)) + 1, bty = "n") # Plot relative recovery over time plot(x = range(Rc$ID_post), y = range(Rc$Rc_rel), type = "n", xlab = "Nb. states after disturbance", ylab = "Relative recovery", main = "Variation of relative recovery") for (i in unique(Rc$disturbed_trajectories)) { lines(Rc[disturbed_trajectories == i, c("ID_post", "Rc_rel")], col = which(unique(Rc$disturbed_trajectories) %in% i) + 1) } legend("topright", legend = unique(Rc$disturbed_trajectories), lty = 1, col = seq_along(unique(Rc$disturbed_trajectories)) + 1, bty = "n") ## ----graphical_NetChange, echo=FALSE, fig.dim=c(6,3.8)------------------------ par(mar = c(0,0,0,0)) # Representative trajectory plot(x = c(0, 16), y = c(0, 8), type = "n", axes = F, xlab = "", ylab = "") for (i in 0:14) { shape::Arrows(x0 = i, y0 = 1, x1 = i + 1, y1 = 1, lwd = 2, arr.adj = 1) } text(x = 0.2, y = 1.3, "RT") # Reference lines(x = c(0, 15), y = c(2.5, 2.5), col = "#E41A1C", lty = 3, lwd = 2) # Disturbed trajectory 1 for (i in 0.5:1.5) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 2.5, y0 = 2.5, x1 = 4, y1 = 4.5, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 4, y0 = 4.5, x1 = 5, y1 = 3.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 0.7, y = 2.8, "DT1", col = "#E41A1C") # Disturbed trajectory 2 for (i in 6:7) { shape::Arrows(x0 = i, y0 = 2.5, x1 = i + 1, y1 = 2.5, lwd = 2, arr.adj = 1, col = "#E41A1C") } shape::Arrows(x0 = 8, y0 = 2.5, x1 = 9.5, y1 = 3.5, lwd = 2, arr.adj = 1, col = "#E41A1C") shape::Arrows(x0 = 9.5, y0 = 3.5, x1 = 14, y1 = 1.5, lwd = 2, arr.adj = 1, col = "#E41A1C") text(x = 6.2, y = 2.8, "DT2", col = "#E41A1C") # States text(x = c(2.5, 4, 5.3, 8, 9.5, 14), y = c(2.3, 4.7, 3.7, 2.3, 3.7, 1.3), c(0:2, 0:2), cex = 0.9, col = "#E41A1C") # Net change DT1 lines(x = c(2.5, 5), y = c(2.5, 3.5), lty = 3, col = "#999999") lines(x = c(5, 5), y = c(2.5, 3.5), lwd = 2, col = "#999999") text(x = 5.2, y = 3.2, adj = 0, expression(NC[abs1] ~ "> 0"), col = "#999999") draw_arc(x = 2.5, y = 2.5, r = 0.6, ini = 0, fin = 25, lwd = 2, col = "#999999") text(x = 3.2, y = 2.7, adj = 0, expression(NC[rel1] ~ "> 0"), col = "#999999") # Net change DT2 lines(x = c(8, 14), y = c(2.5, 1.5), lty = 3, col = "#999999") lines(x = c(14, 14), y = c(2.5, 1.5), lwd = 2, col = "#999999") text(x = 15.2, y = 2, expression(NC[abs2] ~ "< 0"), col = "#999999") draw_arc(x = 8, y = 2.5, r = 1.7, ini = -10, fin = 0, lwd = 2, col = "#999999") text(x = 10, y = 2.3, adj = 0, expression(NC[rel2] ~ "< 0"), col = "#999999") text(x = c(0), y = c(6), adj = 0, expression("|"~NC[abs1]~"| = |"~NC[abs2]~"|"), cex = 0.9, col = "#999999") text(x = c(0), y = c(5.5), adj = 0, expression("|"~NC[rel1]~"| > |"~NC[rel2]~"|"), cex = 0.9, col = "#999999") # Legend legend("topleft", c("RT: Representative trajectory", "DT: Disturbed trajectories", "NC: Net change"), col = c("black", "#E41A1C", "#999999"), lty = 1, text.col = c("black", "#E41A1C", "#999999"), bty = "n") ## ----net_change--------------------------------------------------------------- # Compute net change using the same data used for amplitude NC <- net_change(d = d, trajectories = abun$traj, states = abun$state, disturbed_trajectories = abun[disturbed_states == 1]$traj, disturbed_states = abun[disturbed_states == 1]$state, reference = retra_ref, method = "nearest_state") ## ----plot_NetChange, fig.height=5, fig.width=10------------------------------- # ID post-disturbance states NC <- data.table::data.table(NC) NC[, ID_post := 1:(.N), by = disturbed_trajectories] par(mfrow = c(1, 2)) # Plot absolute net change over time plot(x = range(NC$ID_post), y = range(NC$NC_abs), type = "n", xlab = "Nb. states after disturbance", ylab = "Absolute net change", main = "Variation of absolute net change") for (i in unique(NC$disturbed_trajectories)) { lines(NC[disturbed_trajectories == i, c("ID_post", "NC_abs")], col = which(unique(NC$disturbed_trajectories) %in% i) + 1) } legend("topleft", legend = unique(NC$disturbed_trajectories), lty = 1, col = seq_along(unique(NC$disturbed_trajectories)) + 1, bty = "n") # Plot relative net change over time plot(x = range(NC$ID_post), y = range(NC$NC_rel), type = "n", xlab = "Nb. states after disturbance", ylab = "Relative net change", main = "Variation of relative net change") for (i in unique(NC$disturbed_trajectories)) { lines(NC[disturbed_trajectories == i, c("ID_post", "NC_rel")], col = which(unique(NC$disturbed_trajectories) %in% i) + 1) } legend("bottomleft", legend = unique(NC$disturbed_trajectories), lty = 1, col = seq_along(unique(NC$disturbed_trajectories)) + 1, bty = "n") ## ----A_Rc_NC------------------------------------------------------------------ # Merge the results for resistance, amplitude, recovery, and net change results <- Reduce(function(x, y) merge(x, y, all = T), list(Rt, A, Rc, NC)) results <- results[which(results$ID_post == 14), c("disturbed_trajectories", "Rt", "A_abs", "Rc_abs", "NC_abs")] ## ----all_results, echo=FALSE-------------------------------------------------- knitr::kable(results, row.names = F, digits = 3, col.names = c("Disturbed trajectories", "Rt", "A~abs~", "Rc~abs~", "NC~abs~")) ## ----dDis_post---------------------------------------------------------------- # Species abundances in post-disturbance states of the disturbed trajectories # The states after the release of the disturbance are identified by disturbed_states > 1 abun_post <- EDR_data$EDR3_disturbed$abundance[disturbed_states > 1] selcols <- names(EDR_data$EDR3$abundance) # Species abundances in EDR3 and the post-disturbance states of disturbed trajectories abun3_post <- rbind(EDR_data$EDR3$abundance, abun_post[, ..selcols]) # State dissimilarities in EDR3 and the post-disturbance states of disturbed trajectories d3_post <- as.matrix(vegan::vegdist(x = abun3_post[, paste0("sp", 1:12)], method = "bray")) # dDis of each disturbed trajectory dDis_dist <- sapply(unique(abun_post$traj), function(idist){ rm_disturbed <- unique(abun_post$traj[which(abun_post$traj != idist)]) irm <- which(abun3_post$traj %in% rm_disturbed) dDis(d = d3_post[-irm, -irm], d.type = "dStates", trajectories = abun3_post$traj[-irm], states = abun3_post$state[-irm], reference = as.character(idist)) }) ## ----dDis_post_tb, echo = F--------------------------------------------------- knitr::kable(data.frame(disturbed_trajectory = 31:33, dDis = dDis_dist), digits = 3, row.names = F, col.names = c("Disturbed trajectories", "dDis")) ## ----plot_disturbed, fig.dim=c(7,6)------------------------------------------- # Plot EDR3 and its representative trajectories plot(retra_ref, d = d, trajectories = abun$traj, states = abun$state, traj.colors = c(rep("grey", 30), 2:4), xlab = "MDS D1", ylab = "MDS D2", main = "Comparison of disturbed trajectories and EDR3") legend("topleft", c("Representative trajectory", "Individual trajectories in EDR3", "Trajectory 31", "Trajectory 32", "Trajectory 33"), lty = 1, col = c("black", "grey", 2:4), bty = "n")