inst/doc/ecospace.R

In ecospace: Simulating Community Assembly and Ecological Diversification Using Ecospace Frameworks

## -----------------------------------------------------------------------------
library(ecospace)
nchar <- 9
ecospace <- create_ecospace(nchar = nchar, char.state = rep(3, nchar),
 char.type = rep(c("factor", "ord.fac", "ord.num"), nchar / 3))

## -----------------------------------------------------------------------------
Smax <- 50
set.seed(3142)
neutral_sample <- neutral(Sseed = 5, Smax = Smax, ecospace = ecospace)
head(neutral_sample, 10)

## -----------------------------------------------------------------------------
set.seed(3142)
Sseed = 5
redund_sample <- redundancy(Sseed = Sseed, Smax = Smax, ecospace = ecospace)

## -----------------------------------------------------------------------------
unique(redund_sample)

## -----------------------------------------------------------------------------
set.seed(3142)
redund_sample2 <- redundancy(Sseed = Sseed, Smax = Smax, 
                             ecospace = ecospace, strength = 0.95)

## ---- fig.width = 5, fig.height = 5-------------------------------------------
library(FD, quietly = TRUE)
pc <- prcomp(FD::gowdis(redund_sample))
plot(pc$x, type = "n", 
     main = paste("Redundancy model,\n", Smax, "species"))
text(pc$x[,1], pc$x[,2], labels = seq(Smax), col = c(rep("red", Sseed),
     rep("black", 5), rep("slategray", (Smax - Sseed - 5))), 
     pch = c(rep(19, Sseed), rep(21, (Smax - Sseed))), cex = .8)

pc.r <- prcomp(FD::gowdis(redund_sample2))
plot(pc.r$x, type = "n", main = 
       paste("Redundancy model (95% identical),\n", Smax, "species"))
text(pc.r$x[,1], pc.r$x[,2], labels = seq(Smax), col = c(rep("red",
     Sseed), rep("black", 5), rep("slategray", (Smax - Sseed - 5))),
     pch = c(rep(19, Sseed), rep(21, (Smax - Sseed))), cex = .8)

## -----------------------------------------------------------------------------
set.seed(3142)
Sseed = 5
partS_sample <- partitioning(Sseed = Sseed, Smax = Smax, 
                             ecospace = ecospace)

## -----------------------------------------------------------------------------
set.seed(3142)
Sseed = 5
partR_sample <- partitioning(Sseed = Sseed, Smax = Smax, 
                             ecospace = ecospace, rule = "relaxed")

## ---- fig.width = 5, fig.height = 5-------------------------------------------
pc.ps <- prcomp(FD::gowdis(partS_sample))
plot(pc.ps$x, type = "n", 
     main = paste("'Strict' partitioning model,\n", Smax, "species"))
text(pc.ps$x[,1], pc$x[,2], labels = seq(Smax), 
     col = c(rep("red", Sseed), rep("black", 5), 
             rep("slategray", (Smax - Sseed - 5))), 
     pch = c(rep(19, Sseed), rep(21, (Smax - Sseed))), cex = .8)

pc.pr <- prcomp(FD::gowdis(partR_sample))
plot(pc.pr$x, type = "n", 
     main = paste("'Relaxed' partitioning model,\n", Smax, "species"))
text(pc.pr$x[,1], pc.pr$x[,2], labels = seq(Smax), 
     col = c(rep("red", Sseed), rep("black", 5), 
             rep("slategray", (Smax - Sseed - 5))),
     pch = c(rep(19, Sseed), rep(21, (Smax - Sseed))), cex = .8)

## -----------------------------------------------------------------------------
set.seed(3142)
Sseed = 5
exp_sample <- expansion(Sseed = Sseed, Smax = Smax, 
                        ecospace = ecospace)

## ---- fig.width = 5, fig.height = 5-------------------------------------------
pc.e <- prcomp(FD::gowdis(exp_sample))
plot(pc.e$x, type = "n", 
     main = paste("Expansion model,\n", Smax, "species"))
text(pc.e$x[,1], pc$x[,2], labels = seq(Smax), 
     col = c(rep("red", Sseed), rep("black", 5), 
             rep("slategray", (Smax - Sseed - 5))), 
     pch = c(rep(19, Sseed), rep(21, (Smax - Sseed))), cex = .8)

## ---- fig.width = 5, fig.height = 5-------------------------------------------
library(vegan, quietly = TRUE)
start <- neutral_sample[1:Sseed,]
neu <- neutral_sample[(Sseed + 1):Smax,]
red <- redund_sample2[(Sseed + 1):Smax,]
par <- partR_sample[(Sseed + 1):Smax,]
exp <- exp_sample[(Sseed + 1):Smax,]
nmds.data <- rbind(start, neu, red, par, exp)
all <- metaMDS(gowdis(nmds.data), zerodist = "add", k = 2, trymax = 10)

plot(all$points[,1], all$points[,2], 
     col = c(rep("red", Sseed), rep("orange", nrow(neu)), 
             rep("red", nrow(red)), rep("blue", nrow(par)),
             rep("purple", nrow(exp))), 
     pch = c(rep(19, Sseed), rep(21, nrow(neu)), 
             rep(22, nrow(red)), rep(23, nrow(par)), 
             rep(24, nrow(exp))), 
     main = paste("Combined models,\n", Smax, "species per model"),
     xlab = "Axis 1", ylab = "Axis 2", cex = 2, cex.lab = 1.5, lwd = 1)

leg.txt <- c("seed", "neutral", "redundancy", "partitioning",
             "expansion")
leg.col <- c("red", "orange", "red", "blue", "purple")
leg.pch <- c(19, 21, 22, 23, 24)
legend("topright", inset = .02, legend = leg.txt, pch = leg.pch, 
       col = leg.col, cex = .75)

## -----------------------------------------------------------------------------
# Using Smax = 10 here to illustrate calculation for first 25 species in neutral assemblage
options(warn = -1)
metrics <- calc_metrics(samples = neutral_sample, Smax = 10, 
                        Model = "Neutral")
metrics

## -----------------------------------------------------------------------------
# Calculate statistics for just the entire sample
options(warn = -1)
metrics <- calc_metrics(samples = neutral_sample, increm = FALSE)
metrics

## -----------------------------------------------------------------------------
# Using Smax = 10 here to illustrate calculation for first 10 species in neutral assemblage
options(warn = -1)
metrics <- calc_metrics(samples = neutral_sample, Smax = 10, 
                        Model = "Neutral", increm = TRUE)
metrics

## ---- fig.width = 5, fig.height = 5-------------------------------------------
nreps <- 1:25 # A sequence of the samples to be simulated
n.samples <- lapply(X = nreps, FUN = neutral, Sseed = 3, Smax = 20,
                    ecospace)

# Calculate functional diversity metrics for simulated samples
n.metrics <- lapply(X = nreps, FUN = calc_metrics, samples = n.samples,
                    Model = "neutral", Param = "NA")

# Combine lists together into a single dataframe (the function is new to this package, 
# but the newer 'rbindlist' function in 'data.table' package is even faster)
all <- rbind_listdf(n.metrics)

# Calculate mean dynamics across simulations
means <- n.metrics[[1]]
for(n in 1:20) {
  means[n,4:11] <- apply(all[which(all$S == means$S[n]),4:11], 2, mean,                          na.rm = TRUE)
}

# Plot statistics as function of species richness, overlaying mean dynamics
par(mfrow = c(2,4), mar = c(4, 4, 1, .3))
attach(all)

plot(S, H, type = "p", cex = .75, col = "gray")
lines(means$S, means$H, type = "l", lwd = 2)
plot(S, D, type = "p", cex = .75, col = "gray")
lines(means$S, means$D, type = "l", lwd = 2)
plot(S, M, type = "p", cex = .75, col = "gray")
lines(means$S, means$M, type = "l", lwd = 2)
plot(S, V, type = "p", cex = .75, col = "gray")
lines(means$S, means$V, type = "l", lwd = 2)
plot(S, FRic, type = "p", cex = .75, col = "gray")
lines(means$S, means$FRic, type = "l", lwd = 2)
plot(S, FEve, type = "p", cex = .75, col = "gray")
lines(means$S, means$FEve, type = "l", lwd = 2)
plot(S, FDiv, type = "p", cex = .75, col = "gray")
lines(means$S, means$FDiv, type = "l", lwd = 2)
plot(S, FDis, type = "p", cex = .75, col = "gray")
lines(means$S, means$FDis, type = "l", lwd = 2)