## ----setup, 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").') ## ----install, eval=FALSE------------------------------------------------------ # install.packages("ecoregime") # devtools::install_github(repo = "MSPinillos/ecoregime", dependencies = T, build_vignettes = T) ## ----ecoregime---------------------------------------------------------------- library(ecoregime) ## ----citation----------------------------------------------------------------- citation("ecoregime") ## ----data--------------------------------------------------------------------- # ID of the sampling units and observations ID_sampling <- LETTERS[1:10] ID_obs <- 1:3 # Define initial species abundances set.seed(789) initial <- data.frame(sampling_unit = ID_sampling, sp1 = round(runif(10), 2), sp2 = round(runif(10), 2)) # Simulate artificial dynamics simulated_abun <- lapply(1:nrow(initial), function(isampling){ sp1 <- initial$sp1[isampling] sp2 <- initial$sp2[isampling] while (length(sp1) < 3) { set.seed(isampling) sp1 <- c(sp1, sp1[length(sp1)]*sample(seq(1, 2, by = 0.01), 1)) sp2<- c(sp2, sp2[length(sp2)]*sample(seq(0, 1, by = 0.01), 1)) } data.frame(sampling_unit = initial$sampling_unit[isampling], time = 1:3, sp1 = sp1/(sp1+sp2), sp2 = sp2/(sp1+sp2)) }) # Compile species abundances of all sampling units in the same data frame abundance <- do.call(rbind, simulated_abun) ## ----abundance---------------------------------------------------------------- head(abundance) ## ----state_dissim------------------------------------------------------------- # Generate a matrix containing dissimilarities between states state_dissim <- vegan::vegdist(abundance[, c("sp1", "sp2")], method = "bray") as.matrix(state_dissim)[1:6, 1:6] ## ----state_mds, fig = TRUE, fig.height = 5, fig.width = 5, fig.align = "center"---- # Multidimensional scaling state_mds <- smacof::smacofSym(state_dissim, ndim = 2) state_mds <- data.frame(state_mds$conf) # Define different colors for each trajectory traj.colors <- grDevices::palette.colors(10, palette = "Classic Tableau") # Plot the distribution of the states in the state space plot(state_mds$D1, state_mds$D2, col = rep(traj.colors, each = 3), # use different colors for each sampling unit pch = rep(ID_sampling, each = 3), # use different symbols for each sampling unit xlab = "MDS D1", ylab = "MDS D2", main = "State space") ## ----traj_state_mds, fig = TRUE, fig.height = 5, fig.width = 5, fig.align = "center"---- # Plot the EDR trajectories in the state space plot_edr(state_mds, trajectories = abundance$sampling_unit, states = as.integer(abundance$time), traj.colors = traj.colors, xlab = "MDS D1", ylab = "MDS D2", main = "State space") legend("bottomright", unique(abundance$sampling_unit), lty = 1, col = traj.colors, ncol = 5, bty = "n", cex = 0.8) ## ----traj_dissim-------------------------------------------------------------- # Generate a matrix containing dissimilarities between trajectories traj_dissim <- ecotraj::trajectoryDistances( ecotraj::defineTrajectories(state_dissim, sites = rep(ID_sampling, each = 3), surveys = rep(ID_obs, 10)), distance.type = "DSPD" ) as.matrix(traj_dissim)[1:6, 1:6] ## ----traj_mds, fig = TRUE, fig.height = 5, fig.width = 5, fig.align = "center"---- # Multidimensional scaling traj_mds <- smacof::smacofSym(traj_dissim, ndim = 2) traj_mds <- data.frame(traj_mds$conf) # Plot the distribution of the trajectories in the trajectory space plot(traj_mds$D1, traj_mds$D2, col = traj.colors, # use different colors for each sampling unit pch = ID_sampling, # use different symbols for each sampling unit xlab = "MDS D1", ylab = "MDS D2", main = "Trajectory space") ## ----plot_edr, fig = TRUE, fig.height = 5, fig.width = 10, fig.align = "center"---- par(mfrow = c(1, 2)) plot_edr(state_mds, trajectories = abundance$sampling_unit, states = as.integer(abundance$time), type = "gradient", variable = abundance$sp1, state.colors = hcl.colors(5, "Viridis", rev = T), initial = TRUE, xlab = "MDS D1", ylab = "MDS D2", main = "sp1 relative abundance") legend("bottomright", legend = c(paste0("sp1 = ", round(max(abundance$sp1), 2)), rep(NA, 28), paste0("sp1 = ", round(min(abundance$sp1), 2))), fill = hcl.colors(30, "Viridis"), border = NA, y.intersp = 0.2, cex = 0.8, bty = "n") plot_edr(state_mds, trajectories = abundance$sampling_unit, states = as.integer(abundance$time), type = "gradient", variable = abundance$sp2, state.colors = hcl.colors(5, "Viridis", rev = T), initial = TRUE, xlab = "MDS D1", ylab = "MDS D2", main = "sp2 relative abundance") legend("bottomright", legend = c(paste0("sp2 = ", round(max(abundance$sp2), 2)), rep(NA, 28), paste0("sp2 = ", round(min(abundance$sp2), 2))), fill = hcl.colors(30, "Viridis"), border = NA, y.intersp = 0.2, cex = 0.8, bty = "n") ## ----RETRA-EDR---------------------------------------------------------------- # Use set.seed to obtain reproducible results of the segment space in RETRA-EDR set.seed(123) # Apply RETRA-EDR repr_traj <- retra_edr(d = state_dissim, trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10), minSegs = 2) ## ----retra_summary------------------------------------------------------------ summary(repr_traj) ## ----retra_segments----------------------------------------------------------- lapply(repr_traj, "[[", "Segments") ## ----plot_retra, fig = TRUE, fig.height = 5, fig.width = 5, fig.align = "center"---- # Plot the representative trajectories of an EDR plot(repr_traj, # <-- This is a RETRA object returned by retra_edr() # data to generate the state space d = state_mds, trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10), # use the colors previously used for individual trajectories. # (make them more transparent to highlight representative trajectories) traj.colors = adjustcolor(traj.colors, alpha.f = 0.3), # display representative trajectories in blue RT.colors = "royalblue", # select T2 to be displayed with a different color (black) select_RT = "T1", sel.color = "black", # Identify artificial links using dashed lines (default) and a different color (red) link.lty = 2, link.color = "red", # We can use other arguments in plot() xlab = "MDS D1", ylab = "MDS D2", main = "Representative trajectories") # Add a legend legend("bottomright", c("T1", "T2", "Link"), col = c("black", "royalblue", "red"), lty = c(1, 1, 2)) ## ----link_length-------------------------------------------------------------- repr_traj$T1$Link_distance ## ----define_retra_ls---------------------------------------------------------- # List including the sequence of segments for each new trajectory new_traj_ls <- list( # Exclude the first segment in T1 repr_traj$T1$Segments[2:length(repr_traj$T1$Segments)], # Exclude the first segment in T2 repr_traj$T2$Segments[2:length(repr_traj$T2$Segments)], # First segment of T1 and T2: segment composed of states 2 and 3 of the sampling unit C ("C[2-3]") "C[2-3]" ) new_traj_ls ## ----define_retra_df---------------------------------------------------------- # Generate a data frame indicating the states forming the new trajectories new_traj_df <- data.frame( # name of the new trajectories (as many times as the number of states) RT = c(rep("T1.1", 6), rep("T2.1", 6), rep("T1.2", 2)), # name of the trajectories (sampling units) RT_traj = c(rep("B", 2), rep("I", 2), rep("D", 2), # for the first trajectory (T1.1) rep("B", 2), rep("I", 2), rep("G", 2), # for the second trajectory (T2.1) rep("C", 2)), # for the third trajectory (T1.2) # states in each sampling unit RT_states = c(2:3, 2:3, 2:3, # for the first trajectory (T1.1) 2:3, 2:3, 1:2, # for the second trajectory (T2.1) 2:3), # for the third trajectory (T1.2) # representative trajectories obtained in retra_edr() RT_retra = c(rep("T1", 6), rep("T2", 6), rep("T1", 2)) # The last segment belong to both (T1, T2), choose one ) new_traj_df ## ----define_retra------------------------------------------------------------- # Define representative trajectories using a data frame new_repr_traj <- define_retra(data = new_traj_df, # Information of the state space d = state_dissim, trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10), # RETRA object returned by retra_edr() retra = repr_traj) # Define representative trajectories using a list with sequences of segments new_repr_traj_ls <- define_retra(data = new_traj_ls, # Information of the state space d = state_dissim, trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10), # RETRA object returned by retra_edr() retra = repr_traj) if (all.equal(new_repr_traj, new_repr_traj_ls)) { print("Yes, both are equal!") } ## ----plot_newretra, fig = TRUE, fig.height = 5, fig.width = 5, fig.align = "center"---- plot(new_repr_traj, # <-- This is the RETRA object returned by define_retra() # data to generate the state space d = state_mds, trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10), # display individual trajectories in light blue traj.colors = "lightblue", # display representative trajectories in dark blue RT.colors = "darkblue", # select T1.2 to be displayed in a different color (red) select_RT = "T1.2", sel.color = "coral", # Identify artificial links using dashed lines (default), but use the same # color than the representative trajectories (default) link.lty = 2, link.color = NULL, # We can use other arguments in plot() xlab = "MDS D1", ylab = "MDS D2", main = "Defined representative trajectories") # Add a legend legend("bottomright", c("T1.1, T2.1", "T1.2", "Link"), col = c("darkblue", "coral", "darkblue"), lty = c(1, 1, 2)) ## ----sp_composition----------------------------------------------------------- # Set an ID in the abundance matrix abundance$ID <- paste0(abundance$sampling_unit, abundance$time) # Identify the states forming both representative trajectories traj_states <- lapply(repr_traj, function(iRT){ segments <- iRT$Segments seg_components <- strsplit(gsub("\\]", "", gsub("\\[", "-", segments)), "-") traj_states <- vapply(seg_components, function(iseg){ c(paste0(iseg[1], iseg[2]), paste0(iseg[1], iseg[3])) }, character(2)) traj_states <- unique(as.character(traj_states)) traj_states <- data.frame(ID = traj_states, RT_states = 1:length(traj_states)) }) sp_comp <- lapply(traj_states, function(iRT){ data <- merge(iRT, abundance, by = "ID", all.x = T) data <- data[order(data$RT_states), ] data$Shannon <- vegan::diversity(data[, c("sp1", "sp2")], index = "shannon") return(data) }) sp_comp$T1 ## ----plot_sp_diversity, fig = TRUE, fig.height = 5, fig.width = 10, fig.align = "center"---- par(mfrow = c(1, 2)) # Plot the variation of sp1 in T2 plot(x = sp_comp$T2$RT_states, y = sp_comp$T2$sp1, type = "l", col = "royalblue", ylim = c(0, 1), xlab = "RT state", ylab = "Relative abundance", main = "Variation of species composition") # Add the variation of sp1 in T1 lines(x = sp_comp$T1$RT_states, y = sp_comp$T1$sp1, col = "black") # Add the variation of sp2 in T2 lines(x = sp_comp$T2$RT_states, y = sp_comp$T2$sp2, col = "royalblue", lty = 2) # Add the variation of sp2 in T1 lines(x = sp_comp$T1$RT_states, y = sp_comp$T1$sp2, col = "black", lty = 2) # Add a legend legend("bottom", c("T1", "T2", "sp1", "sp2"), col = c("black", "royalblue", "darkgrey", "darkgrey"), pch = c(20, 20, NA, NA), lty = c(NA, NA, 1, 2), ncol = 2, cex = 0.9, bty = "n") # Plot the variation of species diversity in T1 plot(x = sp_comp$T2$RT_states, y = sp_comp$T2$Shannon, type = "l", col = "royalblue", ylim = c(0, 1), xlab = "RT state", ylab = "Shannon index", main = "Variation of species diversity") # Add the variation of species diversity in T2 lines(x = sp_comp$T1$RT_states, y = sp_comp$T1$Shannon, col = "black") # Add a legend legend("topright", c("T1", "T2"), col = c("black", "royalblue"), lty = 1, cex = 0.9, bty = "n") ## ----dDis_data---------------------------------------------------------------- # Change the trajectory identifier and the number of the state abundance_T1 <- sp_comp$T1 abundance_T1$sampling_unit <- "T1" abundance_T1$time <- abundance_T1$RT_states # Add representative trajectories' data to the abundance matrix abundance_T1 <- rbind(abundance, abundance_T1[, names(abundance)]) # Calculate state dissimilarities including the representative trajectory state_dissim_T1 <- vegan::vegdist(abundance_T1[, c("sp1", "sp2")], method = "bray") ## ----dDis_T2------------------------------------------------------------------ # Compute dDis taking T2 as reference dDis_T1 <- dDis(d = state_dissim_T1, d.type = "dStates", trajectories = abundance_T1$sampling_unit, states = abundance_T1$time, reference = "T1") dDis_T1 ## ----dDis_D_C----------------------------------------------------------------- # dDis: reference D dDis_D <- dDis(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling, reference = "D") # dDis: reference C dDis_C <- dDis(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling, reference = "C") dDis_D dDis_C ## ----dDis_w------------------------------------------------------------------- # Define w values initial_sp2 <- abundance[which(abundance$time == 1), ]$sp2 # Identify the index of the reference trajectories ind_D <- which(ID_sampling == "D") ind_C <- which(ID_sampling == "C") # Compute dDis with weighted trajectories: # Considering I as reference dDis_D_w <- dDis(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling, reference = "D", w.type = "precomputed", w.values = initial_sp2[-ind_D]) # Considering F as reference dDis_C_w <- dDis(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling, reference = "C", w.type = "precomputed", w.values = initial_sp2[-ind_C]) dDis_D_w dDis_C_w ## ----dBD---------------------------------------------------------------------- # Calculate dBD dBD(d = state_dissim, d.type = "dStates", trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10)) ## ----dEve--------------------------------------------------------------------- # Calculate dEve dEve(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling) ## ----dEvew-------------------------------------------------------------------- # Calculate dEve weighting trajectories by the initial abundance of sp2 dEve(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling, w.type = "precomputed", w.values = initial_sp2) ## ----data2-------------------------------------------------------------------- # ID of the sampling units for EDR2 initial$sampling_unit <- paste0("inv_", LETTERS[1:10]) # Define initial species abundances for sp3 set.seed(654) initial$sp3 <- round(runif(10, 0, 0.1), 2) # Simulate artificial dynamics simulated_abun2 <- lapply(1:nrow(initial), function(isampling){ sp1 <- initial$sp1[isampling] sp2 <- initial$sp2[isampling] sp3 <- initial$sp3[isampling] while (length(sp1) < 3) { set.seed(isampling) sp1 <- c(sp1, sp1[length(sp1)]*sample(seq(0.5, 1, by = 0.01), 1)) sp2<- c(sp2, sp2[length(sp2)]*sample(seq(0.2, 0.5, by = 0.01), 1)) sp3<- c(sp3, sp3[length(sp3)]*sample(seq(1.7, 2, by = 0.01), 1)) } data.frame(sampling_unit = initial$sampling_unit[isampling], time = 1:3, sp1 = sp1/(sp1+sp2+sp3), sp2 = sp2/(sp1+sp2+sp3), sp3 = sp3/(sp1+sp2+sp3)) }) # Compile species abundances of all sampling units in the same data frame abundance2 <- do.call(rbind, simulated_abun2) ## ----data3-------------------------------------------------------------------- # ID of the sampling units and observations ID_sampling3 <- LETTERS[11:20] # Define initial species abundances set.seed(987) initial3 <- data.frame(sampling_units = ID_sampling3, # sp1 = round(runif(10, 0, 0.05), 2), sp4 = round(runif(10, 0, 0.3), 2)) # Simulate artificial dynamics simulated_abun3 <- lapply(1:nrow(initial), function(isampling){ # sp1 <- initial3$sp1[isampling] sp2 <- initial$sp2[isampling] sp4 <- initial3$sp4[isampling] while (length(sp2) < 3) { set.seed(isampling) sp2<- c(sp2, sp2[length(sp2)]*sample(seq(0.8, 1, by = 0.01), 1)) sp4 <- c(sp4, sp4[length(sp4)]*sample(seq(1, 2, by = 0.01), 1)) } data.frame(sampling_unit = initial3$sampling_unit[isampling], time = 1:3, sp2 = sp2/(sp2+sp4), sp4 = sp4/(sp2+sp4)) }) # Compile species abundances of all sampling units in the same data frame abundance3 <- do.call(rbind, simulated_abun3) ## ----traj_dissim_allEDR------------------------------------------------------- # Bind all abundance matrices abundance_allEDR <- data.table::rbindlist(list(abundance, abundance2, abundance3), fill = T) abundance_allEDR[is.na(abundance_allEDR)] <- 0 # Calculate state dissimilarities including states in the three EDRs state_dissim_allEDR <- vegan::vegdist(abundance_allEDR[, paste0("sp", 1:4)], method = "bray") # Calculate trajectory dissimilarities including trajectories in the three EDRs traj_dissim_allEDR <- ecotraj::trajectoryDistances( ecotraj::defineTrajectories(state_dissim_allEDR, sites = abundance_allEDR$sampling_unit, surveys = abundance_allEDR$time)) ## ----state_mds_allEDR, fig = TRUE, fig.height = 5, fig.width = 5, fig.align = "center"---- # Multidimensional scaling st_mds_all <- smacof::smacofSym(state_dissim_allEDR, ndim = nrow(as.matrix(state_dissim_allEDR))-1) st_mds_all <- data.frame(st_mds_all$conf) # Plot ecological trajectories in the state space state.colorsEDRs <- rep(grDevices::palette.colors(9, palette = "Set 1"), each = 30) # Set an empty plot plot(st_mds_all$D1, st_mds_all$D2, type = "n", xlab = "MDS D1", ylab = "MDS D2", main = "EDRs in the state space") # Add arrows for (isampling in seq(1, 90, 3)) { # From observation 1 to observation 2 shape::Arrows(st_mds_all[isampling, 1], st_mds_all[isampling, 2], st_mds_all[isampling + 1, 1], st_mds_all[isampling + 1, 2], col = state.colorsEDRs[isampling], arr.adj = 1) # From observation 2 to observation 3 shape::Arrows(st_mds_all[isampling + 1, 1], st_mds_all[isampling + 1, 2], st_mds_all[isampling + 2, 1], st_mds_all[isampling + 2, 2], col = state.colorsEDRs[isampling], arr.adj = 1) } # Add a legend legend("topleft", paste0("EDR", 1:3), col = unique(state.colorsEDRs), lty = 1, cex = 0.9, bty = "n") ## ----EDR_dissim--------------------------------------------------------------- # Compute the dissimilarities between EDRs EDR_dissim <- dist_edr(d = traj_dissim_allEDR, d.type = "dTraj", edr = rep(c("EDR1", "EDR2", "EDR3"), each = 10), metric = "dDR") round(EDR_dissim, 2)