Nothing
rspecies <-
function(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)
{
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
type <- match.arg(type)
if (is.function(family))
family <- family()
if (!(family$family %in% c("gaussian", "poisson", "binomial"))) {
print(family)
stop("'family' not allowed")
}
if (is.matrix(b)) {
b <- lapply(1:spp, function(z) b[, z])
}
if (type == "glm.nb") {
if (length(theta == 1))
theta <- rep(theta, spp)
lambda <- lapply(1:spp, function(i) pmax(exp(data.matrix(x) %*% b[[i]]), .Machine$double.eps))
shape <- lapply(1:spp, function(i) lambda[[i]] / theta[i])
y <- sapply(shape, function(z) rpois(n, rgamma(n, z, scale=theta)))
} else {
mu <- lapply(1:spp, function(i) data.matrix(x) %*% b[[i]])
if (length(sigma == 1))
sigma <- rep(sigma, spp)
dia <- sigma^2
if (type == "glmm") {
SIGMA <- lapply(1:spp, function(i) dia[i] * cor)
require(MASS)
mu <- lapply(1:spp, function(i) MASS::mvrnorm(1, mu[[i]], SIGMA[[i]]))
}
if (family$family == "gaussian" && type == "glm") {
y <- sapply(1:spp, function(i) rnorm(n, family$linkinv(mu[[i]]), dia[[i]]))
} else {
RFUN <- switch(family$family,
"gaussian" = function(z) return(z),
"poisson" = function(z) rpois(n, family$linkinv(z)),
"binomial" = function(z) rbinom(n, 1, family$linkinv(z)))
y <- sapply(mu, RFUN)
}
}
out <- matrix(y, n, spp)
if (attrib) {
class(out) <- c("rspecies", "matrix")
attr(out, "call") <- match.call()
attr(out, "seed") <- seed
}
out
}
print.rspecies <- function(x, ...) {
class(x) <- class(x)[-1]
attr(x, "call") <- NULL
attr(x, "seed") <- NULL
print(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.