R/factorfit.R

In vegan: Community Ecology Package

`factorfit` <-
    function (X, P, permutations = 0, strata = NULL, w,  ...)
{
    EPS <- sqrt(.Machine$double.eps)
    P <- as.data.frame(P)
    ## Check that all variables are factors, and coerce if necessary
    if(any(!sapply(P, is.factor)))
        P <- data.frame(lapply(P, function(x)
                        if (is.factor(x)) x else factor(x)))
    P <- droplevels(P) ## make sure only the used levels are present
    if (any(!sapply(P, is.factor)))
        stop("All non-numeric variables must be factors")
    NR <- nrow(X)
    NC <- ncol(X)
    NF <- ncol(P)
    if (missing(w) || is.null(w))
        w <- 1
    if (length(w) == 1)
        w <- rep(w, NR)
    r <- NULL
    pval <- NULL
    totvar <- .C("goffactor", as.double(X), as.integer(rep(0, NR)),
                 as.double(w), as.integer(NR), as.integer(NC), as.integer(1),
                 double(1), double(1), double(1), var = double(1), PACKAGE = "vegan")$var
    sol <- centroids.cca(X, P, w)
    var.id <- rep(names(P), sapply(P, nlevels))
    ## make permutation matrix for all variables handled in the next loop
    permat <- getPermuteMatrix(permutations, NR, strata = strata)
    permutations <- nrow(permat)

    for (i in seq_along(P)) {
        A <- as.integer(P[[i]])
        NL <- nlevels(P[[i]])
        invar <- .C("goffactor", as.double(X), as.integer(A - 1), as.double(w),
                    as.integer(NR), as.integer(NC),
                    as.integer(NL), double(NL), double(NL), double(NL),
                    var = double(1), PACKAGE = "vegan")$var
        r.this <- 1 - invar/totvar
        r <- c(r, r.this)
        if (permutations) {
            A <- as.integer(P[[i]])
            NL <- nlevels(P[[i]])
            ptest <- function(indx, ...) {
                take <- A[indx]
                invar <- .C("goffactor", as.double(X),
                            as.integer(take -  1), as.double(w),
                            as.integer(NR), as.integer(NC),
                            as.integer(NL), double(NL), double(NL), double(NL),
                            var = double(1), PACKAGE = "vegan")$var
                1 - invar/totvar
            }
            tmp <- sapply(seq_len(permutations),
                          function(indx,...) ptest(permat[indx,], ...))
            pval.this <- (sum(tmp >= r.this - EPS) + 1)/(permutations + 1)
            pval <- c(pval, pval.this)
        }
    }
    if (is.null(colnames(X)))
        colnames(sol) <- paste("Dim", 1:ncol(sol), sep = "")
    else colnames(sol) <- colnames(X)
    names(r) <- names(P)
    if (!is.null(pval))
        names(pval) <- names(P)
    out <- list(centroids = sol, r = r, permutations = permutations,
                pvals = pval, var.id = var.id)
    out$control <- attr(permat, "control")
    class(out) <- "factorfit"
    out
}