tab=X             # = X raw table of occurrences
row.covariates=Y  # = Y environmental variables
ndim=3            # number of eigenvec/eigenvals to be extracted

# symSqrt() makes the square root of matrix cc or its inverse
`symSqrt` <-
function(cc,inv=FALSE)
{
  e<-eigen(cc)
  ev<-e$values
  ind<-which(ev>0)
  lbd<-rep(0,length(ev))
  if (inv) lbd[ind]<-1/sqrt(ev[ind]) else lbd[ind]<-sqrt(ev[ind])
  kk<-e$vectors
  return(kk%*%(lbd*t(kk)))
}



> anacor
function (tab, ndim = 2, row.covariates, col.covariates, scaling = c("Benzecri", 
    "Benzecri"), eps = 1e-06) 
{
    tab <- as.matrix(tab)
    scaling[1] <- match.arg(scaling[1], c("standard", "centroid", 
        "Benzecri", "Goodman"))
    scaling[2] <- match.arg(scaling[2], c("standard", "centroid", 
        "Benzecri", "Goodman"))
    if (missing(row.covariates)) {
        row.covariates <- NULL
    }
    else {
        row.covariates <- as.matrix(row.covariates)
        if (nrow(row.covariates) != nrow(tab)) 
            stop("Matrix with row covariates does not match table dimensions!")
    }
    if (missing(col.covariates)) {
        col.covariates <- NULL
    }
    else {
        col.covariates <- as.matrix(col.covariates)
        if (nrow(col.covariates) != ncol(tab)) 
            stop("Matrix with column covariates does not match table dimensions!")
    }
    if (ndim > min(dim(tab)) - 1) 
        stop("Too many dimensions!")
    name <- deparse(substitute(tab))
    if (any(is.na(tab))) 
        tab <- reconstitute(tab, eps = eps)

    n <- dim(tab)[1]             # = r rows
    m <- dim(tab)[2]             # = c columns
    N <- sum(tab)                # = T total
    tab <- as.matrix(tab)
    prop <- as.vector(t(tab))/N  # = P correspondence matrix of frequencies
    r <- rowSums(tab)            # = R row sums
    c <- colSums(tab)            # = C column sums
    qdim <- ndim + 1
    r <- ifelse(r == 0, 1, r)
    c <- ifelse(c == 0, 1, c)
    ROW <- FALSE
    COL <- FALSE

    z <- tab/sqrt(outer(r, c))   # = Ps matrix of centered & scaled frequencies for CA
    chisq <- N * (sum(z^2) - 1)  # ??

    if (is.matrix(row.covariates)) {
        xr <- as.matrix(cbind(rep(1, n), row.covariates)) # = Y in regression form
        nr <- dim(xr)[2]               # number of regression components
        cx <- crossprod(xr, r * xr)    # weighted squarred Y weighted by R
        or <- symSqrt(cx, inv = TRUE)  # inverse of matrix square root of cx
        ROW <- TRUE
    }
    if (is.matrix(col.covariates)) {
        yr <- cbind(rep(1, m), col.covariates)
        mr <- dim(yr)[2]
        cy <- crossprod(yr, c * yr)
        oc <- symSqrt(cy, inv = TRUE)
        COL <- TRUE
    }

    if (ROW) {z <- or %*% crossprod(xr, tab)} else { z <- tab/sqrt(r)}                   

# ^ cross product CCA z for svd

    if (COL) {z <- z %*% yr %*% oc} else {z <- z/outer(rep(1, dim(z)[1]), sqrt(c))}

# ^ else option utilized here further modifies z

    sv <- svd(z, nu = qdim, nv = qdim)   # svd decomposition
    sval <- N * ((sv$d[-1])^2)           # eigenvals scaled by N dropping first column
    sigmavec <- (sv$d)[2:qdim]           # singular values dropping first column

    dimlab <- paste("D", 1:ndim, sep = "")
    Tr <- Xstar <- isetCorRow <- NULL
    Tc <- Ystar <- isetCorCol <- NULL


    if (ROW) {
        x <- (xr %*% or %*% (sv$u))[, -1]     # 
        U <- or %*% (sv$u)[, -1]              # = u 
        Tr <- U %*% (diag(sigmavec))          #
        dimnames(Tr) <- list(paste("beta", 1:(dim(Tr)[1]), sep = ""), 
            dimlab)

        if (COL) 
            y <- (yr %*% oc %*% (sv$v))[, -1]            
        else y <- ((sv$v)/sqrt(c))[, -1]        #
        Xstar <- (diag(1/r) %*% tab %*% y)      #        
        dimnames(Xstar) <- list(rownames(tab), dimlab)
        isetCorRow <- cor(Xstar, row.covariates)
    }
    else {
        x <- ((sv$u)/sqrt(r))[, -1]


    }
    if (COL) {
        y <- (yr %*% oc %*% (sv$v))[, -1]
        V <- oc %*% (sv$v)[, -1]
        Tc <- V %*% (diag(sigmavec))
        dimnames(Tc) <- list(paste("beta", 1:(dim(Tc)[1]), sep = ""), 
            dimlab)
        Ystar <- (diag(1/c) %*% t(tab) %*% x)
        dimnames(Ystar) <- list(colnames(tab), dimlab)
        isetCorCol <- cor(Ystar, col.covariates)
    }
    else {
        y <- ((sv$v)/sqrt(c))[, -1]
    }


    x <- x * sqrt(N)         #
    y <- y * sqrt(N)         #


    if (scaling[2] == "centroid") 
        y <- y * outer(rep(1, m), sigmavec)  #
    if (scaling[1] == "centroid") 
        x <- x * outer(rep(1, n), sigmavec)  #


    if (scaling[2] == "Goodman") 
        y <- y * outer(rep(1, m), sqrt(sigmavec))
    if (scaling[1] == "Goodman") 
        x <- x * outer(rep(1, n), sqrt(sigmavec))
    if (scaling[2] == "Benzecri") 
        y <- y * outer(rep(1, m), sigmavec)
    if (scaling[1] == "Benzecri") 
        x <- x * outer(rep(1, n), sigmavec)


    benzres <- benzdist(scaling, x, y, z, tab, n, m, r, c, row.covariates, 
        col.covariates)

    prob <- sval/chisq
    pcum <- cumsum(prob)

    cs.mat <- cbind(sval, prob, pcum)

    res.deriv <- gsvdDer(tab, ndim)
    res.deriv.scaling <- gsvdScal(res.deriv, scaling)
    res.acov <- acovuv(res.deriv.scaling, prop, n, m, N, ndim)
    se.sigma <- sqrt(diag(res.acov$acovd))

    colnames(x) <- colnames(y) <- dimlab
    rownames(x) <- rownames(tab)
    rownames(y) <- colnames(tab)
    rownames(cs.mat) <- paste("Component", 1:dim(cs.mat)[1])
    colnames(cs.mat) <- c("Chisq", "Proportion", "Cumulative Proportion")


    result <- list(datname = name, tab = tab, ndim = ndim, row.covariates = row.covariates, 
        col.covariates = col.covariates, row.scores = x, col.scores = y, 
        chisq.decomp = cs.mat, chisq = chisq, singular.values = sigmavec, 
        se.singular.values = se.sigma, left.singvec = sv$u[, 
            -1], right.singvec = sv$v[, -1], eigen.values = sigmavec^2, 
        scaling = scaling, bdmat = benzres[1:4], rmse = benzres[5:6], 
        row.acov = res.acov$acovu, col.acov = res.acov$acovv, 
        cancoef = list(rows = Tr, columns = Tc), sitescores = list(rows = Xstar, 
            columns = Ystar), isetcor = list(rows = isetCorRow, 
            columns = isetCorCol))

    class(result) <- "anacor"
    result
}
<environment: namespace:anacor>