R/simulate.rda.R

In vegan: Community Ecology Package

`simulate.rda` <-
    function(object, nsim = 1, seed = NULL, indx = NULL, rank = "full",
             correlated = FALSE, ...) 
{
    ## Fail if there is no constrained component (it could be possible
    ## to change the function to handle unconstrained ordination, too,
    ## when rank < "full", but that would require redesign)
    if (is.null(object$CCA))
        stop("function can be used only with constrained ordination")
    
    ## Handle RNG: code directly from stats::simulate.lm
    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))
    }
    ## indx can be an output of permute::shuffleSet in which case it
    ## is a nsim x nrow matrix, or it may be a single vector, in which
    ## case it will changed to shuffleSet
    if (!is.null(indx))
        if (is.vector(indx))
            dim(indx) <- c(1, length(indx))
    ## If nsim is missing, take it from indx (if given)
    if (missing(nsim) && !is.null(indx))
        nsim <- nrow(indx)
    ## Check that dims match
    if (!is.null(indx))
        if(nrow(indx) != nsim)
            stop(gettextf("'nsim' (%d) and no. of 'indx' rows (%d) do not match",
                          nsim, nrow(indx)))
    ## Proper simulation: very similar for simulate.lm, but produces
    ## an array of response matrices

    ftd <- predict(object, type = "response", rank = rank)
    ## pRDA: add partial Fit to the constrained
    if (!is.null(object$pCCA))
        ftd <- ftd + object$pCCA$Fit
    ## if(is.null(indx)), we have parametric Gaussian simulation and
    ## need to generate sd matrices. The residuals sd is always taken
    ## from the unconstrained (residual) component $CA$Xbar. If
    ## species are uncorrelated, we need only species sd's, but if
    ## correlated, we also need species covariances.
    if (!correlated)
        dev <- outer(rep(1, nrow(ftd)), apply(object$CA$Xbar, 2, sd))
    else
        dev <- cov(object$CA$Xbar)
    ## Generate an array
    ans <- array(0, c(dim(ftd), nsim))
    for (i in seq_len(nsim)) {
        if (!is.null(indx))
            ans[,,i] <- as.matrix(ftd + object$CA$Xbar[indx[i,],])
        else if (!correlated)
            ans[,,i] <- as.matrix(ftd + matrix(rnorm(length(ftd), sd = dev),
                                               nrow = nrow(ftd)))
        else {
            ans[,,i] <- t(apply(ftd, 1,
                                function(x) mvrnorm(1, mu = x, Sigma = dev))) 
        }
    }
    ## set RNG attributes
    if (is.null(indx))
        attr(ans, "seed") <- RNGstate
    else 
        attr(ans, "seed") <- "index"
    ## set commsim attributes if nsim > 1, else return a 2-dim matrix
    if (nsim == 1) {
        ans <- ans[,,1]
        attributes(ans) <- attributes(ftd)
    } else {
        dimnames(ans) <- list(rownames(ftd), colnames(ftd),
                              paste("sim", seq_len(nsim), sep = "_"))
        attr(ans, "data") <- round(ftd + object$CA$Xbar, 12)
        attr(ans, "method") <- paste("simulate", ifelse(is.null(indx),
                                                        "parametric", "index"))
        attr(ans, "binary") <- FALSE
        attr(ans, "isSeq") <- FALSE
        attr(ans, "mode") <- "double"
        attr(ans, "start") <- 1L
        attr(ans, "end") <- as.integer(nsim)
        attr(ans, "thin") <- 1L
        class(ans) <- c("simulate.rda", "simmat", "array")
    }
    ans
}

### simulate. cca was cloned from simulate.rda.  Works with internal
### Chi-square standardized form, and at the end back-standardizes
### with row and column totals and matrix grand totals. This does not
### still guarantee that all marginal totals are positive.

`simulate.cca` <-
    function(object, nsim = 1, seed = NULL, indx = NULL, rank = "full",
             correlated = FALSE, ...)
{
    ## Fail if no CCA
    if (is.null(object$CCA))
        stop("function can be used only with constrained ordination")
    ## Handle RNG: code directly from stats::simulate.lm
    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))
    }
    ## Preparations like in simulate.rda()
    if (!is.null(indx))
        if (is.vector(indx))
            dim(indx) <- c(1, length(indx))
    if (missing(nsim) && !is.null(indx))
        nsim <- nrow(indx)
    if (!is.null(indx))
        if(nrow(indx) != nsim)
            stop(gettextf("'nsim' (%d) and no. of 'indx' rows (%d) do not match",
                          nsim, nrow(indx)))   
    ## Need sqrt of rowsums for weighting
    sq.r <- sqrt(object$rowsum)
    ## Fitted value
    ftd <- predict(object, type = "working", rank = rank)
    ## pCCA: add partial Fit to the constrained
    if (!is.null(object$pCCA))
        ftd <- ftd + object$pCCA$Fit
    ## Residual Xbar need weighting and back-weighting
    Xbar <- sweep(object$CA$Xbar, 1, sq.r, "*")
    ## Simulation
    if (correlated)
        dev <- cov(Xbar)
    else
        dev <- outer(rep(1, nrow(ftd)), apply(Xbar, 2, sd))
    ans <- array(0, c(dim(ftd), nsim))
    for (i in seq_len(nsim)) {
        if (is.null(indx)) {
            if (correlated)
                tmp <- mvrnorm(nrow(ftd), numeric(ncol(ftd)), Sigma = dev)
            else
                tmp <- matrix(rnorm(length(ftd), sd = dev), 
                          nrow = nrow(ftd))
            ans[,,i] <- as.matrix(ftd + sweep(tmp, 1, sq.r, "/"))
        }
        else 
            ans[,,i] <- as.matrix(ftd + sweep(Xbar[indx[i,],], 1, sq.r, "/"))
    }
    ## From internal form to the original form with fixed marginal totals
    rc <- object$rowsum %o% object$colsum
    for (i in seq_len(nsim))
        ans[,,i] <- (ans[,,i] * sqrt(rc) + rc) * object$grand.total
    ## RNG attributes
    if (is.null(indx))
        attr(ans, "seed") <- RNGstate
    else
        attr(ans, "seed") <- "index"
    ## set commsim attributes if nsim > 1, else return a 2-dim matrix
    if (nsim == 1) {
        ans <- ans[,,1]
        attributes(ans) <- attributes(ftd)
    } else {
        dimnames(ans) <- list(rownames(ftd), colnames(ftd),
                              paste("sim", seq_len(nsim), sep = "_"))
        obsdata <- ftd + object$CA$Xbar
        obsdata <- (obsdata * sqrt(rc) + rc) * object$grand.total
        attr(ans, "data") <- round(obsdata, 12)
        attr(ans, "method") <- paste("simulate", ifelse(is.null(indx),
                                                        "parametric", "index"))
        attr(ans, "binary") <- FALSE
        attr(ans, "isSeq") <- FALSE
        attr(ans, "mode") <- "double"
        attr(ans, "start") <- 1L
        attr(ans, "end") <- as.integer(nsim)
        attr(ans, "thin") <- 1L
        class(ans) <- c("simulate.cca", "simmat", "array")
    }    
    ans
}


### capscale method: copies simulate.rda as much as possible. Function
### works with the internal metric scaling mapping of fit and error,
### but returns Euclidean distances adjusted to the original scaling
### of input dissimilarities. Only the real components are used, and
### capscale() of simulated dissimilarities have no Imaginary
### component.

`simulate.capscale` <-
    function(object, nsim = 1, seed = NULL, indx = NULL, rank = "full",
             correlated = FALSE, ...) 
{
    ## Fail if no CCA component
    if (is.null(object$CCA))
        stop("function can be used only with constrained ordination")
    if (is.null(indx) && correlated)
        warning("argument 'correlated' does not work and will be ignored")
    ## Handle RNG: code directly from stats::simulate.lm
    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))
    }
    if (nsim > 1)
        .NotYetUsed("nsim")
    ## predict.capscale cannot be used because it returns either
    ## dissimilarities ("response") or scores with the rank of the
    ## constrained solution, and we need rank of the data (not of
    ## constraints).
    if (rank > 0) {
        ftd <- qr.fitted(object$CCA$QR, object$CCA$Xbar)
        ## redo analysis when rank < full
        if (rank < object$CCA$rank) {
            x <- svd(ftd, nu = rank, nv = rank)
            ftd <- x$u %*% diag(x$d[1:rank], nrow=rank) %*% t(x$v)
        }
    } else {
        ftd <- matrix(0, nrow=nrow(object$CA$Xbar),
                      ncol = ncol(object$CA$Xbar))
    }
    ## add partial Fit to the constrained
    if (!is.null(object$pCCA))
        ftd <- ftd + object$pCCA$Fit
    if (is.null(indx))
        ans <- as.data.frame(ftd + matrix(rnorm(length(ftd), 
               sd = outer(rep(1,nrow(ftd)), apply(object$CA$Xbar, 2, sd))), 
               nrow = nrow(ftd)))
    else
        ans <- ftd + object$CA$Xbar[indx,]
    ## return Euclidean distances
    ans <- dist(ans)
    ## remove adjustment done in capscale and put dissimilarities to
    ## (approximately) original scale
    ans <- ans/object$adjust
    if (is.null(indx))
        attr(ans, "seed") <- RNGstate
    else
        attr(ans, "seed") <- indx
    ans
}