Coalescent-Based Simulation of Ecological Communities

Description

Coalescent-Based Simulation of Ecological Communities as proposed by Munoz et al. (2018) <doi:10.1111/2041-210X.12918>. The package includes a tool for estimating parameters of community assembly by using Approximate Bayesian Computation.

Details

The DESCRIPTION file:

Index of help topics:

abund                   Compute absolute and relative abundances in the
                        local community and the reference pool
coalesc                 Coalescent-based simulation of ecological
                        communities undergoing both neutral and
                        niche-base dynamics
coalesc_abc             Estimation of neutral and non-neutral
                        parameters of community assembly using
                        Approximate Bayesian Computation (ABC)
ecolottery-package      Coalescent-Based Simulation of Ecological
                        Communities
forward                 Simulation of neutral and niche-based community
                        dynamics forward in time
plot_comm               Regional vs. Local trait distributions of
                        abundances
tcor                    Generates Correlated Traits

Further information is available in the following vignettes:

Two basic functions: coalesc for coalescent-based simulation, and forward for forward-in-time simulation

Author(s)

François Munoz [aut, cre], Matthias Grenié [aut], Pierre Denelle [aut], Adrien Taudière [ctb], Fabien Laroche [ctb], Caroline Tucker [ctb], Cyrille Violle [ctb]

Maintainer: François Munoz <francois.munoz@hotmail.fr>

References

Hurtt, G. C. and S. W. Pacala (1995). "The consequences of recruitment limitation: reconciling chance, history and competitive differences between plants." Journal of Theoretical Biology 176(1): 1-12.

Hubbell, S. P. (2001). "The Unified Neutral Theory of Biodiversity". Princeton University Press.

Gravel, D., C. D. Canham, M. Beaudet and C. Messier (2006). "Reconciling niche and neutrality: the continuum hypothesis." Ecology Letters 9(4): 399-409.

Munoz, F., P. Couteron, B. R. Ramesh and R. S. Etienne (2007). "Estimating parameters of neutral communities: from one Single Large to Several Small samples." Ecology 88(10): 2482-2488.

Munoz, F., B. R. Ramesh and P. Couteron (2014). "How do habitat filtering and niche conservatism affect community composition at different taxonomic resolutions?" Ecology 95(8): 2179-2191.

Examples

## Coalescent-based simulation of stabilizing habitat filtering around 
## t = 0.5
J <- 100; theta <- 50; m <- 0.5;
comm <- coalesc(J, m, theta, filt = function(x) 0.5 - abs(0.5 - x))
plot_comm(comm)

## Forward-in-time simulation of stabilizing habitat filtering around 
## t = 0.5, over 100 time steps

# A regional pool including 100 species each including 10 individuals
pool <- sort(rep(as.character(1:100), 10))

# Initial community composed of 10 species each including 10 individuals, 
# with trait information for niche-based dynamics
initial <- data.frame(sp = sort(rep(as.character(1:10), 10)),
                      trait = runif(100))
final <- forward(initial = initial, prob = 0.5, gens = 100, pool = pool,
                 filt = function(x) 0.5 - abs(0.5 - x))
plot_comm(final)

Compute absolute and relative abundances in the local community and the reference pool

Description

Compute the abundances and relative abundances of species in simulated communities and in the corresponding species pools. The input must be an output of either coalesc or the forward functions.

Usage

abund(x)

Arguments

Value

Author(s)

F. Munoz, P. Denelle and M. Grenie

Examples

# Simulation of a neutral community including 500 individuals
J <- 500; theta <- 50; m <- 0.05;
comm1a <- coalesc(J, m, theta) 
abund1a <- abund(comm1a) 
  
# Log-series distribution of regional abundances
fit <- vegan::fisherfit(abund1a$pool$ab)
freq <- as.numeric(names(fit$fisher))
plot(log(freq), fit$fisher,
     xlab = "Frequency (log)",
     ylab = "Species", type = "n")
rect(log(freq - 0.5), 0, log(freq + 0.5), fit$fisher, col="skyblue")

alpha <- fit$estimate
k <- fit$nuisance

curve(alpha * k^exp(x) / exp(x), log(0.5), max(log(freq)),
      col = "red", lwd = 2, add = TRUE)

# Relationship between local and regional abundances
par(mfrow=c(1, 2))
plot(abund1a$pool[rownames(abund1a$com), "relab"],
  abund1a$com$relab,
  main = "m = 0.05",
  xlab = "Regional abundance",
  ylab = "Local abundance",
  log = "xy")
abline(0,1)

# With higher immigration rate
m <- 0.95
comm1b <- coalesc(J, m, theta) 
abund1b <- abund(comm1b) 
plot(abund1b$pool[rownames(abund1b$com),"relab"],
  abund1b$com$relab,
  main = "m = 0.95",
  xlab = "Regional abundance",
  ylab = "Local abundance",
  log = "xy")
abline(0,1)

Coalescent-based simulation of ecological communities undergoing both neutral and niche-base dynamics

Description

Simulates the composition of a community based on immigration from a regional pool, habitat filtering depending on local environment and species traits, and local birth-death stochastic dynamics.

Usage

coalesc(J, m = 1, theta = NULL, filt = NULL, pool = NULL, traits = NULL, 
        Jpool = 50 * J, verbose = FALSE)

Arguments

Details

Coalescent-based simulation of a community of size J. This generic function can simulate a neutral community (if filt = NULL) or a community undergoing both neutral and niche-based dynamics. In the latter case, filt(t) represents the relative ability of immigrants with trait values t in the regional pool to enter the community.

Value

Author(s)

F. Munoz

References

Hurtt, G. C. and S. W. Pacala (1995). "The consequences of recruitment limitation: reconciling chance, history and competitive differences between plants." Journal of Theoretical Biology 176(1): 1-12.

Gravel, D., C. D. Canham, M. Beaudet and C. Messier (2006). "Reconciling niche and neutrality: the continuum hypothesis." Ecology Letters 9(4): 399-409.

Munoz, F., P. Couteron, B. R. Ramesh and R. S. Etienne (2007). "Estimating parameters of neutral communities: from one Single Large to Several Small samples." Ecology 88(10): 2482-2488.

Munoz, F., B. R. Ramesh and P. Couteron (2014). "How do habitat filtering and niche conservatism affect community composition at different taxonomic resolutions?" Ecology 95(8): 2179-2191.

Examples

# Simulation of a neutral community including 100 individuals
J <- 500; theta <- 50; m <- 0.1
comm1 <- coalesc(J, m, theta) 
# Regional and local trait distributions
plot_comm(comm1)
 
# Define a regional pool of species with equal abundances
pool <- cbind(1:10000, rep(1:500, 20), rep(NA, 10000))
# Uniform distribution of trait values
t.sp <- runif(500)
# No intraspecific variation
pool[,3] <- t.sp[pool[,2]]
# Generate a neutral community drawn from the pool
comm2<- coalesc(J, m, pool = pool) 
plot_comm(comm2)

# Directional habitat filtering toward t = 0
comm3 <- coalesc(J, m, filt = function(x) 1 - x, pool = pool)
# Regional and local trait distributions
plot_comm(comm3)

# Function for environmental filtering
sigma <- 0.1
filt_gaussian <- function(t, x) exp(-(x - t)^2/(2*sigma^2))

# Stabilizing habitat filtering around t = 0.1
comm4a <- coalesc(J, m, filt = function(x) filt_gaussian(0.1, x), pool = pool)
plot_comm(comm4a)
# Stabilizing habitat filtering around t = 0.5
comm4b <- coalesc(J, m, theta, filt = function(x) filt_gaussian(0.5, x),
                  pool = pool)
plot_comm(comm4b)
# Stabilizing habitat filtering around t = 0.9
comm4c <- coalesc(J, m, theta, filt = function(x) filt_gaussian(0.9, x),
                  pool = pool)
plot_comm(comm4c)

# Mean trait values in communities reflect the influence of habitat filtering
mean(comm4a$com[, 3])
mean(comm4b$com[, 3])
mean(comm4c$com[, 3])

# Disruptive habitat filtering around t = 0.5
comm5 <- coalesc(J, m, filt = function(x) abs(0.5 - x), pool = pool)
plot_comm(comm5)

# Multi-modal habitat filtering 
t.sp <- rnorm(500)
pool[, 3] <- t.sp[pool[,2]]
comm6 <- coalesc(J, m, filt = function(x) sin(3*x) + 1, pool = pool)
plot_comm(comm6)

Estimation of neutral and non-neutral parameters of community assembly using Approximate Bayesian Computation (ABC)

Description

Estimates parameters of neutral migration-drift dynamics (through migration rate m and parameters of environmental filtering (through a filtering function filt.abc()) from the composition of a local community and the related regional pool.

Usage

  coalesc_abc(comm.obs, pool = NULL, multi = "single", traits = NULL,
              f.sumstats, filt.abc = NULL, params = NULL, 
              theta.max = NULL, nb.samp = 10^6, parallel = TRUE,
              tol = NULL, pkg = NULL, method="rejection")
  do.simul(J, pool = NULL, multi = "single", nb.com = NULL, 
              traits = NULL, f.sumstats = NULL, filt.abc = NULL,
              params, theta.max = NULL, nb.samp = 10^6, 
              parallel = TRUE, tol = NULL, pkg = NULL,
              method = "rejection")

Arguments

Details

coalesc_abc() performs ABC estimation for one (if multi = FALSE, default) or several communities (if multi = TRUE) related to the same regional pool.

do.simul() provides the simulated communities used in ABC estimation, and is not intended to be used directly.

Value

Author(s)

F. Munoz

References

Jabot, F., and J. Chave. 2009. Inferring the parameters of the neutral theory of biodiversity using phylogenetic information and implications for tropical forests. Ecology Letters 12:239-248.

Csillery, K., M. G. B. Blum, O. E. Gaggiotti, and O. Francois. 2010. Approximate Bayesian computation (ABC) in practice. Trends in Ecology & Evolution 25:410-418.

Csillery, K., O. Francois, and M. G. Blum. 2012. abc: an R package for Approximate Bayesian Computation (ABC). Methods in Ecology and Evolution 3:475-479.

See Also

abc() in abc package, parLapply() in parallel package.

Examples

# Trait-dependent filtering function
filt_gaussian <- function(t, params) exp(-(t-params[1])^2/(2*params[2]^2))

# Definition of parameters and their range
params <- data.frame(rbind(c(0, 1), c(0.05, 1)))
row.names(params) <- c("topt", "sigmaopt")
# Number of values to sample in prior distributions
nb.samp <- 10^6 # Should be large

## Not run: 
# Basic summary statistics
f.sumstats <- function(com) array(dimnames=list(c("cwm", "cwv", "cws",
                                                  "cwk", "S", "Es")),
                                  c(mean(com[,3]), var(com[,3]), 
                                    e1071::skewness(com[,3]),  
                                    e1071::kurtosis(com[,3]),
                                    vegan::specnumber(table(com[,2])),
                                    vegan::diversity(table(com[,2]))))

# An observed community is here simulated (known parameters)
comm <- coalesc(J = 400, m = 0.5, theta = 50,
                filt = function(x) filt_gaussian(x, c(0.2, 0.1)))

# ABC estimation of the parameters based on observed community composition 
## Warning: this function may take a while
res <- coalesc_abc(comm$com, comm$pool, f.sumstats = f.sumstats,
                   filt.abc = filt_gaussian, params = params, 
                   nb.samp = nb.samp, parallel = TRUE, 
                   pkg = c("e1071","vegan"), method = "neuralnet")
plot(res$abc, param = res$par)
hist(res$abc)

# Cross validation
## Warning: this function is slow
res$cv <- abc::cv4abc(param = res$par, sumstat = res$ss, nval = 1000,
                      tols = c(0.01, 0.1, 1), method = "neuralnet")
plot(res$cv)

# Multiple community option
# When the input is a site-species matrix, use argument multi="tab"
# See vignette Barro_Colorado for more details

# When the input is a list of communities, use argument multi="seqcom"
comm.obs <- list()

comm.obs[[1]] <- cbind(rep(1,400), coalesc(J = 400, m = 0.5, filt = function(x)
                                           filt_gaussian(x, c(0.2, 0.1)),
                                           pool = comm$pool)$com))
comm.obs[[2]] <- cbind(rep(2,400), coalesc(J = 400, m = 0.5, filt = function(x) 
                                           filt_gaussian(x, c(0.5, 0.1)),
                                           pool = comm$pool)$com))
comm.obs[[3]] <- cbind(rep(3,400), coalesc(J = 400, m = 0.5, filt = function(x) 
                                           filt_gaussian(x, c(0.8, 0.1)),
                                           pool = comm$pool)$com))

comm.obs <- lapply(comm.obs, as.matrix)

res <- coalesc_abc(comm.obs, comm$pool, multi="seqcom", f.sumstats=f.sumstats,
                   filt.abc = filt_gaussian, params = params, nb.samp = nb.samp,
                   parallel = TRUE, pkg = c("e1071","vegan"), tol = 0.1,
                   method = "neuralnet")
                   
lapply(res$abc, summary)


## End(Not run)

Simulation of neutral and niche-based community dynamics forward in time

Description

Simulates niche-based (habitat filtering and/or limiting similarity) and neutral community dynamics from a given initial composition, over a given number of generations.

Usage

forward(initial, prob = 0, d = 1, gens = 150, keep = FALSE,
        pool = NULL, limit.sim = FALSE, coeff.lim.sim = 1,
        sigm = 0.1, filt = NULL, prob.death = NULL,
        method.dist = "euclidean", plot_gens = FALSE)
get_number_of_gens(given_size, pool, nbrep = 5, prob = 1, d = 1, 
                   gens = NULL, limit.sim = FALSE, 
                   coeff.lim.sim = 1, sigm = 0.1, filt = NULL, 
                   prob.death = NULL, method.dist = "euclidean",
                   plot_gens = FALSE)
pick(com, d = 1, prob = 0, pool = NULL, prob.death = prob.death,
     limit.sim = NULL, coeff.lim.sim = 1, sigm = 0.1, filt = NULL, 
     new.index = new.index, method.dist = "euclidean")
pick.mutate(com, d = 1, prob.of.mutate = 0, new.index = 0)
pick.immigrate(com, d = 1, prob.of.immigrate = 0, pool, 
               prob.death = NULL, limit.sim = NULL, coeff.lim.sim = 1,
               sigm = 0.1, filt = NULL, method.dist = "euclidean")

Arguments

Details

It is a zero-sum game, so that the number of individuals of the community is fixed to the number of individuals in initial community.

When niche-based dynamics are simulated, the niche-based constraints influence both immigration and mortality.

Function get_number_of_gen() allows determining the number of generations needed to reach stationary richness for given parameterization of forward(). The target number of generation is based on assessing the change point in species richness change over time for replicate simulated communities with random initial composition. A conservative measure is proposed as the maximum time to reach stationary richness over the replicate simulated communities.

Functions pick.immigrate() and pick.mutate() are used to simulate immigration and speciation events within a time step. They are embedded in forward and are not really intended for the end user.

Value

Author(s)

F. Munoz, derived from the untb function of R. Hankin.

References

For neutral dynamics, S. P. Hubbell 2001. "The Unified Neutral Theory of Biodiversity". Princeton University Press.

Examples

## Not run: 
# Initial community composed of 10 species each including 10 individuals
initial1 <- rep(as.character(1:10), each = 10)

# Simulation of speciation and drift dynamics over 100 time steps
final1 <- forward(initial = initial1, prob = 0.1, gens = 1000)
# The final community includes new species (by default names begins with "new")
final1$com$sp # includes new species generated by speciation events

# A regional pool including 100 species each including 10 individuals
pool <- rep(as.character(1:100), each = 10)

# Simulation of migration and drift dynamics over 1000 time steps
final2 <- forward(initial = initial1, prob = 0.1, gens = 1000, pool = pool)
# The final community includes species that have immigrated from the pool
final2$com$sp # includes new species that immigrated from the pool

# Initial community composed of 10 species each including 10 individuals, 
# with trait information for niche-based dynamics
initial2 <- data.frame(sp = rep(as.character(1:10), each = 10), 
                      trait = runif(100))

# Simulation of stabilizing hab. filtering around t = 0.5, over 1000 time steps
sigm <- 0.1
filt_gaussian <- function(t,x) exp(-(x - t)^2/(2*sigm^2))
final3 <- forward(initial = initial2, prob = 0.1, gens = 1000, pool = pool, 
                 filt = function(x) filt_gaussian(0.5,x))
plot_comm(final3) # trait distribution in final community

# With higher immigration
final4 <- forward(initial = initial2, prob = 0.8, gens = 1000, pool = pool, 
                 filt = function(x) filt_gaussian(0.5,x))
plot_comm(final4) # should be closer to 0.5

# Simulation of limiting similarity, over 1000 time steps
final5 <- forward(initial = initial2, prob = 0.1, gens = 1000, pool = pool, 
                 limit.sim = TRUE)
plot_comm(final5)

# Stronger limiting similarity
final6 <- forward(initial = initial2, prob = 0.1, gens = 1000, pool = pool, 
                 limit.sim = TRUE, coeff.lim.sim = 20)
plot_comm(final6) # the distribution will be more even

# Variation of community richness with time
final7 <- forward(initial = initial2, prob = 0.1, gens = 1000, pool = pool, 
                 limit.sim = TRUE, keep = TRUE, plot_gens = TRUE)

# Check stationarity
plot(unlist(lapply(final7$com_t, function(x) length(unique(x[, 2])))), 
     xlab = "Time step", ylab = "Community richness") 

# Index of limiting similarity over time
plot(final7$dist.t, xlab = "Time step", ylab = "Limiting similarity")

## End(Not run)

Regional vs. Local trait distributions of abundances

Description

Graphical function to used on the output of coalesc() or forward() functions.It aims at plotting links between regional and local trait/abundance distributions.

Usage

plot_comm(x, type = "trait", seltrait = 1, main = NULL)

Arguments

Details

If type = "trait", the function provides density plots of the trait or abundance distributions in the regional pool and in a local community. If type = "abund", the function displays the relationship between regional and local species relative abundances. By default type = "trait". To be used on the output of coalesc() or forward() functions.

Value

Return two stacked ggplot2 density plots if type = "trait" and a biplot if type = "abund".

Author(s)

F. Munoz; P. Denelle

Examples

# Simulation of a neutral community including 100 individuals
J <- 500; theta <- 50; m <- 0.1;
comm1 <- coalesc(J, m, theta) 
plot_comm(comm1)
plot_comm(comm1, type = "abund")

# Stabilizing habitat filtering around t = 0.5
comm2 <- coalesc(J, m, theta, filt = function(x) 0.5 - abs(0.5 - x))
plot_comm(comm2)
plot_comm(comm2, type = "abund")

Generates Correlated Traits

Description

Create two random vectors of traits correlated between each other or a vector of traits correlated to an existing one. The linear correlation is defined by the parameter rho.

Usage

tcor(n, rho = 0.5, mar.fun = rnorm, x = NULL, ...)

Arguments

Details

rho parameter is set to 0.5 by default. x = NULL by default. Code adapted from: http://stats.stackexchange.com/questions/15011/generate-a-random-variable-with-a-defined-correlation-to-an-existing-variable

Value

Return a data.frame with two numeric columns, each column defining a trait.

Author(s)

P. Denelle F. Munoz

Examples

# With no predefined trait
traits <- tcor(n = 10000, rho = 0.8)
plot(traits[, 1], traits[, 2])
cor(traits[, 1], traits[, 2])

# With existing trait
existing_trait <- rnorm(10000, 10, 1)
traits <- tcor(n = 10000, rho = 0.8, x = existing_trait)
plot(traits[, 1], traits[, 2])
cor(traits[, 1], traits[, 2])