rspecies: Community Matrix Generation
In rspecies: Multivariate Methods in Ecology with R: Simulating Communities

Description Usage Arguments Details Value Author(s) Examples

View source: R/rspecies.R

These function generate random community matrices, where species' abundances/occurrences and environmental relationships can be explicitly specified. Abundances are modelled by species specific models.

rspecies(n, spp, b=rep(0, spp), x=rep(1, n), type=c("glm", "glm.nb", "glmm"),
    family=gaussian, sigma=1, cor=diag(1, n, n), theta=1, 
    attrib=TRUE, seed=NULL)
## S3 method for class 'rspecies'
print(x, ...)

`n`	The number of sites (rows) in the data matrix.
`spp`	The number of species (columns) in the data matrix.
`b`	Coefficients to be used for determining species' environmental relationships. They can be a vector of length `spp` if only intercepts are provided, or a list with `spp` elements (one for each species, each element must be of the same length p for the same covariates, including the intercept).
`x`	A vector (if only intercepts are provided) or a model matrix with `n` rows and p columns.
`type`	Type of model to be used. `glm` is generalized linear model with error distribution and link function defined by the argument `family`. `glmm.nb` negative binomial (Gamma mixture of Poisson distributions) GLM. `glmm` generalized linear mixed model with random intercept (extra environmental variability). Standard deviation of the random intercept is determided by `sigma`. `cor` can be used to introduce nonindependence among sites (but species are independent).
`family`	Currently `gaussian` `poisson` and `binomial` are supported. If `type="glm.nb"`, the family argument is ignored.
`sigma`	Standard deviation of the Gaussian linear model (`type="glm"`, and `family=gaussian`), or that of the random intercept with zero mean (`type="glmm"`). Can be a vector of length `spp` if species are assumed to have different variances.
`cor`	Correlation matrix for sites. Diagonal should be 1, must be a symmetric square matrix.
`theta`	Scale parameter of the Gamma distribution used when `type="glm.nb"`.
`attrib`	Logical, if attributes (call and seed) should be attached. Can be a vector of length `spp` if species are assumed to have different variances (that is proportional to the Poisson mean controlled by the underlying Gamma distribution).
`seed`	Random seed.
`...`	Further arguments.

The linear model part for all functions is defined as mu = x %*% b. The link function from the family is then used to link this systematic part with the random component of the model.

The GLMM versions use the mean (mu = x %*% b) and the variance covariance matrix (defined by sigma and cor) to generate Multivariate Normal random numbers.

The model matrix should be in the usual format, providing design variables for factors. The formula interface of the model.matrix function can be applied (see Examples).

A community data matrix with n rows and spp columns.

Peter Solymos

## --- random landscape ---

## rnd seed setup
set.seed(1234)
## no. of sites (20x20 grid)
n <- 400
## Random landscape from B Bolker (http://www.zoology.ufl.edu/bolker/landsc.html)
## landsc R package available at http://www.zoology.ufl.edu/bolker
## randomize phase of a 2-D landscape
rland1 <- function(x) {
    n <- nrow(x)   # linear dimension
    l <- length(x) # total number of nodes
    y <- fft(x)
    z <- matrix(complex(modulus=Mod(y),arg=runif(l,-pi,pi)))
    y <- fft(z,inverse=TRUE)/l
    return(matrix(Mod(y),nrow=n))
}
## given a correlation-by-(radial)-distance function r
##  (which can be either a vector or a function)
##  simulate a 2-D landscape with linear size n 
## substitute into a grid of radial distances
simland2 <- function(r,n) {
 g <- outer(1:n,1:n,function(x,y)sqrt(x^2+y^2))
 if (is.function(r)) {
   gr <- apply(g,c(1,2),r)
  }  else if (is.numeric(r)) {
   gr <- matrix(r[round(g)],nrow=n) }
 rland1(gr)  # pass it to a function which FFTs, randomizes phases,
             # and inverse-FFTs
}
## creates binary landscape
simbinland <- function(n, r, prob) {
    r[r>1] <- 1
    g2 <- simland2(r,n)
    gg <- as.numeric(g2 >= quantile(g2, prob))
    dim(gg) <- dim(g2)
    gg
}
## actual landscape generation takes place here
lscape <- simbinland(20, exp(-(1:(2*n))/20)+0.1*rnorm(2*n), 0.5)
hab <- as.factor(array(lscape))
levels(hab) <- c("forest","meadow")
## hab: factor with 2 levels
## temp: temperature, N-S gradient + some rnd noise
## x, y: grid coordinates
x <- data.frame(expand.grid(x=1:20, y=1:20), 
    hab=hab, temp=rep(seq(-1,1,len=20), each=20) + rnorm(n, 0, 0.1))

## --- random species responses ---

## no. of species
m <- 20
## preferences for meadow level of hab (2 groups)
beta.meadow <- c(rnorm(m/2, -1, 0.5), rnorm(m/2, 1, 0.5))
## hate meadows: 1. and 3. quartiles
## love meadows: 2. and 4. quartiles
beta.meadow <- beta.meadow[c(1:5,11:15,6:10,16:20)]
## warm (1. and 2. quartiles) vs. cold requiring species(3. and 4. quartiles)
## unimodal relationship
bu <- c(rnorm(m/2, -0.5, 0.1), rnorm(m/2, 0.5, 0.2))
bh <- rlnorm(m, 0.8, 0.5)
bt <- rlnorm(m, -0.5, 0.2)
uht <- cbind(bu, bh, bt)
b <- t(cbind(t(apply(uht, 1, uht2beta)), beta.meadow))

## --- community matrix generation ---

phi <- runif(m, 0.5, 1.5)
## exponential spatial correlation: rho^d
Cor <- 0.25^as.matrix(dist(x[,1:2]))
diag(Cor) <- 1
## stratified sampling from 400 sites
id <- c(sample(which(x$hab=="forest"), 50), sample(which(x$hab=="meadow"), 50))
X <- model.matrix(~ temp + I(temp^2) + hab, x[id,])
ynb <- rspecies(length(id), m, b, X, type="glm.nb", theta=rnorm(m, 5))
y11 <- rspecies(length(id), m, b, X, type="glm", family=poisson)
y21 <- rspecies(length(id), m, b, X, type="glmm", family=poisson, sigma=1.5*phi)
y22 <- rspecies(length(id), m, b, X, type="glmm", family=poisson, sigma=1.5*phi, cor=Cor[id, id])

## --- exploratory examples ---

plot(0, type="n", ylim=c(0,10), xlim=c(-1,1))
apply(uht, 1, function(i) lines(seq(-1,1,len=100), respgau(seq(-1,1,len=100), i)))

if (require(vegan)) {
plot(cca(y11~temp+hab,x[id,]))
plot(hclust(vegdist(t(y11))))
adonis(y11~temp+hab,x[id,])
}