[Vegan-commits] r1388 - in branches/1.17: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 2 10:38:50 CET 2010
Author: jarioksa
Date: 2010-12-02 10:38:50 +0100 (Thu, 02 Dec 2010)
New Revision: 1388
Modified:
branches/1.17/R/as.mlm.cca.R
branches/1.17/R/as.mlm.rda.R
branches/1.17/R/capscale.R
branches/1.17/R/eigenvals.R
branches/1.17/R/fitted.procrustes.R
branches/1.17/R/intersetcor.R
branches/1.17/R/isomap.R
branches/1.17/R/mantel.correlog.R
branches/1.17/R/ordiParseFormula.R
branches/1.17/R/pcnm.R
branches/1.17/R/plot.prc.R
branches/1.17/R/plot.procrustes.R
branches/1.17/R/prc.R
branches/1.17/R/print.cca.R
branches/1.17/R/print.summary.prc.R
branches/1.17/R/procrustes.R
branches/1.17/R/scores.rda.R
branches/1.17/R/summary.prc.R
branches/1.17/R/treedist.R
branches/1.17/R/vif.cca.R
branches/1.17/R/wcmdscale.R
branches/1.17/inst/ChangeLog
branches/1.17/man/capscale.Rd
branches/1.17/man/eigenvals.Rd
branches/1.17/man/pcnm.Rd
branches/1.17/man/plot.cca.Rd
branches/1.17/man/prc.Rd
branches/1.17/man/procrustes.Rd
branches/1.17/man/treedive.Rd
branches/1.17/man/wcmdscale.Rd
Log:
merge most things up to r1385 except fitspecaccum & friends and ordiR2step to branches/1.17
Modified: branches/1.17/R/as.mlm.cca.R
===================================================================
--- branches/1.17/R/as.mlm.cca.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/as.mlm.cca.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -5,7 +5,10 @@
wa <- x$CCA$wa
wa <- sweep(wa, 1, sqrt(w), "*")
X <- qr.X(x$CCA$QR)
- colnames(X) <- colnames(X)[x$CCA$QR$pivot]
+ ## qr.X gives wrong column names now, and they are fixed here
+ ## (hopefully fixed in 2.12.1, but that's only a wish).
+ if (getRversion() <= "2.12.0")
+ colnames(X)[x$CCA$QR$pivot] <- colnames(X)
lm(wa ~ . - 1, data = as.data.frame(X))
}
Modified: branches/1.17/R/as.mlm.rda.R
===================================================================
--- branches/1.17/R/as.mlm.rda.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/as.mlm.rda.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -2,7 +2,10 @@
function (x)
{
X <- qr.X(x$CCA$QR)
- colnames(X) <- colnames(X)[x$CCA$QR$pivot]
+ ## We hope that column names will be fixed in R 2.12.1 (but
+ ## perhaps in vain).
+ if (getRversion() <= "2.12.0")
+ colnames(X)[x$CCA$QR$pivot] <- colnames(X)
lm(x$CCA$wa ~ . - 1, data = as.data.frame(X))
}
Modified: branches/1.17/R/capscale.R
===================================================================
--- branches/1.17/R/capscale.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/capscale.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -94,7 +94,7 @@
names(negax) <- paste("NEG", seq_along(negax), sep="")
sol$CA$eig <- c(sol$CA$eig, negax)
sol$CA$imaginary.chi <- sum(negax)
- sol$tot.chi <- sol$tot.chi + sol$CA$imaginary.chi
+ sol$tot.chi <- sol$tot.chi + abs(sol$CA$imaginary.chi)
sol$CA$imaginary.rank <- length(negax)
sol$CA$imaginary.u.eig <- X$negaxes
}
Modified: branches/1.17/R/eigenvals.R
===================================================================
--- branches/1.17/R/eigenvals.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/eigenvals.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -10,15 +10,19 @@
function(x, ...)
{
## svd and eigen return unspecified 'list', see if this could be
- ## either of them
+ ## either of them (like does cmdscale)
out <- NA
if (is.list(x)) {
## eigen
if (length(x) == 2 && all(names(x) %in% c("values", "vectors")))
out <- x$values
## svd: return squares of singular values
- if (length(x) == 3 && all(names(x) %in% c("d", "u", "v")))
+ else if (length(x) == 3 && all(names(x) %in% c("d", "u", "v")))
out <- x$d^2
+ ## cmdscale() will return all eigenvalues from R 2.12.1
+ else if (getRversion() > "2.12.0" &&
+ all(c("points","eig","GOF") %in% names(x)))
+ out <- x$eig
}
class(out) <- "eigenvals"
out
@@ -52,7 +56,7 @@
else
out <- c(x$CCA$eig, x$CA$eig)
if (!is.null(out))
- class(out) <- c("eigenvals")
+ class(out) <- "eigenvals"
out
}
@@ -61,24 +65,35 @@
function(x, ...)
{
out <- x$eig
- class(out) <-"eigenvals"
+ class(out) <- "eigenvals"
out
}
+## pcnm (in vegan)
+`eigenvals.pcnm` <-
+ function(x, ...)
+{
+ out <- x$values
+ class(out) <- "eigenvals"
+ out
+}
+
`print.eigenvals` <-
function(x, ...)
{
- print(unclass(x), ...)
+ print(zapsmall(unclass(x), ...))
invisible(x)
}
`summary.eigenvals` <-
function(object, ...)
{
- vars <- object/sum(object)
+ ## abs(object) is to handle neg eigenvalues of wcmdscale and
+ ## capscale
+ vars <- object/sum(abs(object))
importance <- rbind(`Eigenvalue` = object,
`Proportion Explained` = round(vars, 5),
- `Cumulative Proportion`= round(cumsum(vars), 5))
+ `Cumulative Proportion`= round(cumsum(abs(vars)), 5))
out <- list(importance = importance)
class(out) <- c("summary.eigenvals", "summary.prcomp")
out
Modified: branches/1.17/R/fitted.procrustes.R
===================================================================
--- branches/1.17/R/fitted.procrustes.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/fitted.procrustes.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,8 +1,27 @@
-"fitted.procrustes" <-
+`fitted.procrustes` <-
function(object, truemean = TRUE, ...)
{
fit <- object$Yrot
if (truemean)
- fit <- sweep(fit, 2, object$translation, "+")
+ fit <- sweep(fit, 2, object$xmean, "+")
fit
}
+
+## Like above, except when takes 'newata'
+
+`predict.procrustes` <-
+ function(object, newdata, truemean = TRUE, ...)
+{
+ if (missing(newdata))
+ return(fitted(object, truemean = truemean))
+ if (object$symmetric)
+ stop(gettextf("'predict' not available for symmetric procrustes analysis with 'newdata'"))
+ Y <- as.matrix(newdata)
+ ## scaling and rotation
+ Y <- object$scale * Y %*% object$rotation
+ ## translation: always
+ Y <- sweep(Y, 2, object$translation, "+")
+ if (!truemean)
+ Y <- sweep(Y, 2, object$xmean, "-")
+ Y
+}
Modified: branches/1.17/R/intersetcor.R
===================================================================
--- branches/1.17/R/intersetcor.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/intersetcor.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -7,5 +7,9 @@
w <- object$rowsum
wa <- sweep(object$CCA$wa, 1, sqrt(w), "*")
}
- cor(qr.X(object$CCA$QR), wa)
+ X <- qr.X(object$CCA$QR)
+ ## current R (2.12.0) uses wrong column names in pivoted qr.X()
+ if (getRversion() <= "2.12.0")
+ colnames(X)[object$CCA$QR$pivot] <- colnames(X)
+ cor(X, wa)
}
Modified: branches/1.17/R/isomap.R
===================================================================
--- branches/1.17/R/isomap.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/isomap.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -3,6 +3,13 @@
{
dist <- isomapdist(dist, ...)
out <- cmdscale(dist, k=ndim, eig=TRUE)
+ ## some versions of cmdscale may return NaN points corresponding
+ ## to negative eigenvalues.
+ if ((naxes <- sum(out$eig > 0)) < ndim) {
+ out$points <- out$points[, seq(naxes), drop = FALSE]
+ warning(gettextf("isomap returns only %d axes with positive eigenvalues",
+ naxes))
+ }
npoints <- nrow(out$points)
net <- matrix(FALSE, nrow=npoints, ncol=npoints)
net[lower.tri(net)][attr(dist, "net")] <- TRUE
Modified: branches/1.17/R/mantel.correlog.R
===================================================================
--- branches/1.17/R/mantel.correlog.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/mantel.correlog.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,155 +1,157 @@
-`mantel.correlog` <-
- function(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL,
- cutoff=TRUE, r.type="pearson", nperm=999, mult="holm",
- progressive=TRUE)
-{
- r.type <- match.arg(r.type, c("pearson", "spearman", "kendall"))
- mult <- match.arg(mult, c("sidak", p.adjust.methods))
-
- epsilon <- .Machine$double.eps
- D.eco <- as.matrix(D.eco)
-
- ## Geographic distance matrix
- if(!is.null(D.geo)) {
- if(!is.null(XY))
- stop("You provided both a geographic distance matrix and a list of site coordinates. Which one should the function use?")
- D.geo <- as.matrix(D.geo)
- } else {
- if(is.null(XY)) {
+`mantel.correlog` <-
+ function(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL,
+ cutoff=TRUE, r.type="pearson", nperm=999, mult="holm",
+ progressive=TRUE)
+{
+ r.type <- match.arg(r.type, c("pearson", "spearman", "kendall"))
+ mult <- match.arg(mult, c("sidak", p.adjust.methods))
+
+ epsilon <- .Machine$double.eps
+ D.eco <- as.matrix(D.eco)
+
+ ## Geographic distance matrix
+ if(!is.null(D.geo)) {
+ if(!is.null(XY))
+ stop("You provided both a geographic distance matrix and a list of site coordinates. Which one should the function use?")
+ D.geo <- as.matrix(D.geo)
+ } else {
+ if(is.null(XY)) {
stop("You did not provide a geographic distance matrix nor a list of site coordinates")
- } else {
- D.geo <- as.matrix(dist(XY))
- }
- }
-
- n <- nrow(D.geo)
- if(n != nrow(D.eco))
- stop("Numbers of objects in D.eco and D.geo are not equal")
- n.dist <- n*(n-1)/2
- vec.D <- as.vector(as.dist(D.geo))
- vec.DD <- as.vector(D.geo)
-
- ## Number of classes and breakpoints
-
- if(!is.null(break.pts)) {
- ## Use the list of break points
- if(n.class > 0)
- stop("You provided both a number of classes and a list of break points. Which one should the function use?")
- n.class = length(break.pts) - 1
-
- } else {
- ## No breakpoints have been provided: equal-width classes
- if(n.class == 0) {
+ } else {
+ D.geo <- as.matrix(dist(XY))
+ }
+ }
+
+ n <- nrow(D.geo)
+ if(n != nrow(D.eco))
+ stop("Numbers of objects in D.eco and D.geo are not equal")
+ n.dist <- n*(n-1)/2
+ vec.D <- as.vector(as.dist(D.geo))
+ vec.DD <- as.vector(D.geo)
+
+ ## Number of classes and breakpoints
+
+ if(!is.null(break.pts)) {
+ ## Use the list of break points
+ if(n.class > 0)
+ stop("You provided both a number of classes and a list of break points. Which one should the function use?")
+ n.class = length(break.pts) - 1
+
+ } else {
+ ## No breakpoints have been provided: equal-width classes
+ if(n.class == 0) {
## Use Sturges rule to determine the number of classes
n.class <- ceiling(1 + log(n.dist, base=2))
- }
- ## Compute the breakpoints from n.class
- start.pt <- min(vec.D)
- end.pt <- max(vec.D)
- width <- (end.pt - start.pt)/n.class
- break.pts <- vector(length=n.class+1)
- break.pts[n.class+1] <- end.pt
- for(i in 1:n.class)
- break.pts[i] <- start.pt + width*(i-1)
- }
-
- half.cl <- n.class %/% 2
-
- ## Move the first breakpoint a little bit to the left
- break.pts[1] <- break.pts[1] - epsilon
-
- ## Find the break points and the class indices
- class.ind <- break.pts[1:n.class] +
- (0.5*(break.pts[2:(n.class+1)]-break.pts[1:n.class]))
-
- ## Create the matrix of distance classes
- vec2 <- vector(length=n^2)
- for(i in 1:n^2)
- vec2[i] <- min( which(break.pts >= vec.DD[i]) ) - 1
-
- ## Start assembling the vectors of results
- class.index <- NA
- n.dist <- NA
- mantel.r <- NA
- mantel.p <- NA
- ## check.sums = matrix(NA,n.class,1)
-
- ## Create a model-matrix for each distance class, then compute a Mantel test
- for(k in 1:n.class) {
- class.index <- c(class.index, class.ind[k])
- vec3 <- rep(0, n*n)
- sel <- which(vec2 == k)
- vec3[sel] <- 1
- mat.D2 <- matrix(vec3,n,n)
- diag(mat.D2) <- 0
- n.dis <- sum(mat.D2)
- n.dist <- c(n.dist, n.dis)
- if(n.dis == 0) {
- mantel.r <- c(mantel.r, NA)
- mantel.p <- c(mantel.p, NA)
- } else {
- row.sums <- apply(mat.D2, 1, sum)
- ## check.sums[k,1] = length(which(row.sums == 0))
- if((cutoff==FALSE) ||
- !(cutoff==TRUE && k > half.cl && any(row.sums == 0))) {
- temp <- mantel(mat.D2, D.eco, method=r.type, permutations=nperm)
- mantel.r <- c(mantel.r, -temp$statistic)
- temp.p <- temp$signif
- if(temp$statistic >= 0) {
- temp.p <- ((temp.p*nperm)+1)/(nperm+1)
- } else {
- temp.p <- (((1-temp.p)*nperm)+1)/(nperm+1)
- }
- mantel.p <- c(mantel.p, temp.p)
- } else {
- mantel.r <- c(mantel.r, NA)
- mantel.p <- c(mantel.p, NA)
- }
- }
- }
-
- mantel.res <- cbind(class.index, n.dist, mantel.r, mantel.p)
- mantel.res <- mantel.res[-1,]
-
- ## Note: vector 'mantel.p' starts with a NA value
- mantel.p <- mantel.p[-1]
- n.tests <- length(which(!is.na(mantel.p)))
-
- if(mult == "none") {
- colnames(mantel.res) <-
- c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)")
- } else {
- ## Correct P-values for multiple testing
- if(progressive) {
- p.corr <- mantel.p[1]
- if(mult == "sidak") {
- for(j in 2:n.tests)
- p.corr <- c(p.corr, 1-(1-mantel.p[j])^j)
- } else {
- for(j in 2:n.tests) {
- temp <- p.adjust(mantel.p[1:j], method=mult)
- p.corr <- c(p.corr, temp[j])
- }
- }
- } else {
- ## Correct all p-values for 'n.tests' simultaneous tests
- if(mult == "sidak") {
- p.corr <- 1 - (1 - mantel.p[1:n.tests])^n.tests
- } else {
- p.corr <- p.adjust(mantel.p[1:n.tests], method=mult)
- }
- }
- temp <- c(p.corr, rep(NA,(n.class-n.tests)))
- mantel.res <- cbind(mantel.res, temp)
- colnames(mantel.res) <-
- c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)", "Pr(corrected)")
- }
- rownames(mantel.res) <-
- rownames(mantel.res,do.NULL = FALSE, prefix = "D.cl.")
-
- ## Output the results
- res <- list(mantel.res=mantel.res, n.class=n.class, break.pts=break.pts,
- mult=mult, n.tests=n.tests, call=match.call() )
- class(res) <- "mantel.correlog"
- res
-}
+ }
+ ## Compute the breakpoints from n.class
+ start.pt <- min(vec.D)
+ end.pt <- max(vec.D)
+ width <- (end.pt - start.pt)/n.class
+ break.pts <- vector(length=n.class+1)
+ break.pts[n.class+1] <- end.pt
+ for(i in 1:n.class)
+ break.pts[i] <- start.pt + width*(i-1)
+ }
+
+ half.cl <- n.class %/% 2
+
+ ## Move the first breakpoint a little bit to the left
+ break.pts[1] <- break.pts[1] - epsilon
+
+ ## Find the break points and the class indices
+ class.ind <- break.pts[1:n.class] +
+ (0.5*(break.pts[2:(n.class+1)]-break.pts[1:n.class]))
+
+ ## Create the matrix of distance classes
+ vec2 <- vector(length=n^2)
+ for(i in 1:n^2)
+ vec2[i] <- min( which(break.pts >= vec.DD[i]) ) - 1
+
+ ## Start assembling the vectors of results
+ class.index <- NA
+ n.dist <- NA
+ mantel.r <- NA
+ mantel.p <- NA
+ ## check.sums = matrix(NA,n.class,1)
+
+ ## Create a model-matrix for each distance class, then compute a Mantel test
+ for(k in 1:n.class) {
+ class.index <- c(class.index, class.ind[k])
+ vec3 <- rep(0, n*n)
+ sel <- which(vec2 == k)
+ vec3[sel] <- 1
+ mat.D2 <- matrix(vec3,n,n)
+ diag(mat.D2) <- 0
+ n.dis <- sum(mat.D2)
+ n.dist <- c(n.dist, n.dis)
+ if(n.dis == 0) {
+ mantel.r <- c(mantel.r, NA)
+ mantel.p <- c(mantel.p, NA)
+ } else {
+ row.sums <- apply(mat.D2, 1, sum)
+ ## check.sums[k,1] = length(which(row.sums == 0))
+ if((cutoff==FALSE) ||
+ !(cutoff==TRUE && k > half.cl && any(row.sums == 0))) {
+ temp <- mantel(mat.D2, D.eco, method=r.type, permutations=nperm)
+ mantel.r <- c(mantel.r, -temp$statistic)
+ temp.p <- temp$signif
+
+ ## The mantel() function produces a one-tailed p-value
+ ## (H1: r>0) Here, compute a one-tailed p-value in
+ ## direction of the sign
+ if(temp$statistic < 0) {
+ temp.p <- ((1-temp.p)*(nperm+1))/(nperm+1)
+ }
+ mantel.p <- c(mantel.p, temp.p)
+ } else {
+ mantel.r <- c(mantel.r, NA)
+ mantel.p <- c(mantel.p, NA)
+ }
+ }
+ }
+
+ mantel.res <- cbind(class.index, n.dist, mantel.r, mantel.p)
+ mantel.res <- mantel.res[-1,]
+
+ ## Note: vector 'mantel.p' starts with a NA value
+ mantel.p <- mantel.p[-1]
+ n.tests <- length(which(!is.na(mantel.p)))
+
+ if(mult == "none") {
+ colnames(mantel.res) <-
+ c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)")
+ } else {
+ ## Correct P-values for multiple testing
+ if(progressive) {
+ p.corr <- mantel.p[1]
+ if(mult == "sidak") {
+ for(j in 2:n.tests)
+ p.corr <- c(p.corr, 1-(1-mantel.p[j])^j)
+ } else {
+ for(j in 2:n.tests) {
+ temp <- p.adjust(mantel.p[1:j], method=mult)
+ p.corr <- c(p.corr, temp[j])
+ }
+ }
+ } else {
+ ## Correct all p-values for 'n.tests' simultaneous tests
+ if(mult == "sidak") {
+ p.corr <- 1 - (1 - mantel.p[1:n.tests])^n.tests
+ } else {
+ p.corr <- p.adjust(mantel.p[1:n.tests], method=mult)
+ }
+ }
+ temp <- c(p.corr, rep(NA,(n.class-n.tests)))
+ mantel.res <- cbind(mantel.res, temp)
+ colnames(mantel.res) <-
+ c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)", "Pr(corrected)")
+ }
+ rownames(mantel.res) <-
+ rownames(mantel.res,do.NULL = FALSE, prefix = "D.cl.")
+
+ ## Output the results
+ res <- list(mantel.res=mantel.res, n.class=n.class, break.pts=break.pts,
+ mult=mult, n.tests=n.tests, call=match.call() )
+ class(res) <- "mantel.correlog"
+ res
+}
Modified: branches/1.17/R/ordiParseFormula.R
===================================================================
--- branches/1.17/R/ordiParseFormula.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/ordiParseFormula.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -8,6 +8,8 @@
flapart <- fla <- formula <- formula(Terms, width.cutoff = 500)
specdata <- formula[[2]]
X <- eval.parent(specdata, n = envdepth)
+ ## X is usually a matrix, but it is "dist" with capscale():
+ X <- as.matrix(X)
indPartial <- attr(Terms, "specials")$Condition
zmf <- ymf <- Y <- Z <- NULL
formula[[2]] <- NULL
Modified: branches/1.17/R/pcnm.R
===================================================================
--- branches/1.17/R/pcnm.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/pcnm.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,5 +1,5 @@
-"pcnm" <-
- function(dis, threshold, w)
+`pcnm` <-
+ function(dis, threshold, w, dist.ret = FALSE)
{
EPS <- sqrt(.Machine$double.eps)
wa.old <- options(warn = -1)
@@ -19,6 +19,11 @@
res$vectors <- sweep(res$vectors, 2, sqrt(res$values[1:k]), "/")
colnames(res$vectors) <- paste("PCNM", 1:k, sep="")
res$threshold <- threshold
+ if (dist.ret) {
+ attr(dis, "method") <- paste(attr(dis, "method"), "pcnm")
+ attr(dis, "threshold") <- threshold
+ res$dist <- dis
+ }
class(res) <- "pcnm"
res
}
Modified: branches/1.17/R/plot.prc.R
===================================================================
--- branches/1.17/R/plot.prc.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/plot.prc.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,8 +1,10 @@
-"plot.prc" <-
- function (x, species = TRUE, select, scaling = 2, axis = 1, type = "l",
- xlab, ylab, ylim, lty = 1:5, col = 1:6, pch, legpos, cex = 0.8,
- ...)
+`plot.prc` <-
+ function (x, species = TRUE, select, scaling = 3, axis = 1, type = "l",
+ xlab, ylab, ylim, lty = 1:5, col = 1:6, pch, legpos, cex = 0.8,
+ ...)
{
+ ## save level names before getting the summary
+ levs <- x$terminfo$xlev[[2]]
x <- summary(x, scaling = scaling, axis = axis)
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
@@ -10,46 +12,45 @@
xax <- rownames(b)
if (missing(xlab))
xlab <- x$names[1]
- if (missing(ylab))
+ if (missing(ylab))
ylab <- "Effect"
- if (!missing(select))
+ if (!missing(select))
x$sp <- x$sp[select]
- if (missing(ylim))
- ylim <- if (species)
+ if (missing(ylim))
+ ylim <- if (species)
range(b, x$sp, na.rm = TRUE)
else range(b, na.rm = TRUE)
if (species) {
op <- par("mai")
- mrg <- max(strwidth(names(x$sp), cex = cex, units = "in")) +
+ mrg <- max(strwidth(names(x$sp), cex = cex, units = "in")) +
strwidth("mmm", cex = cex, units = "in")
par(mai = c(op[1:3], max(op[4], mrg)))
}
- if (missing(pch))
+ if (missing(pch))
pch <- as.character(1:nrow(b))
- matplot(xax, b, type = type, xlab = xlab, ylab = ylab, ylim = ylim,
+ matplot(xax, b, type = type, xlab = xlab, ylab = ylab, ylim = ylim,
cex = cex, lty = lty, col = col, pch = pch, ...)
abline(h = 0, col = "gray")
if (species) {
- linestack(x$sp, at = par("usr")[2], add = TRUE, hoff = 1,
+ linestack(x$sp, at = par("usr")[2], add = TRUE, hoff = 1,
cex = cex, ...)
rug(x$sp, side = 4)
}
if (missing(legpos)) {
holes <- abs(par("usr")[3:4] - range(b, na.rm = TRUE))
- if (holes[1] > holes[2])
+ if (holes[1] > holes[2])
legpos <- "bottomleft"
else legpos <- "topleft"
}
if (!is.na(legpos)) {
- levs <- rownames(coef(x))
nl <- length(levs)
pp <- type %in% c("b", "p")
pl <- type %in% c("b", "l")
- if (length(lty) == 1)
- lty <- rep(lty, nl)
- legend(legpos, legend = levs, col = col, lty = if (pl)
- lty[1:nl], pch = if (pp)
- pch, cex = cex, title = x$names[2])
+ if (length(lty) == 1)
+ lty <- rep(lty, nl-1)
+ legend(legpos, legend = levs, col = c("gray", col),
+ lty = if (pl) lty[c(1,1:(nl-1))],
+ pch = if (pp) pch, cex = cex, title = x$names[2])
}
invisible()
}
Modified: branches/1.17/R/plot.procrustes.R
===================================================================
--- branches/1.17/R/plot.procrustes.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/plot.procrustes.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,9 +1,16 @@
"plot.procrustes" <-
- function (x, kind = 1, choices = c(1,2), xlab, ylab, main, ar.col = "blue",
- len = 0.05, ...)
+ function (x, kind = 1, choices = c(1,2), to.target = TRUE,
+ type = "p", xlab, ylab, main, ar.col = "blue",
+ len = 0.05, cex = 0.7, ...)
{
- Yrot <- x$Yrot[, choices]
- X <- x$X[, choices]
+ type <- match.arg(type, c("points", "text", "none"))
+ if (to.target) {
+ tails <- x$Yrot[, choices]
+ heads <- x$X[, choices]
+ } else {
+ tails <- x$X[, choices]
+ heads <- x$Yrot[, choices]
+ }
if (missing(main))
main <- "Procrustes errors"
if (kind <= 1) {
@@ -12,12 +19,12 @@
xlab <- paste("Dimension", choices[1])
if (missing(ylab))
ylab <- paste("Dimension", choices[2])
- xrange <- range(Yrot[, 1], X[, 1])
- yrange <- range(Yrot[, 2], X[, 2])
+ xrange <- range(tails[, 1], heads[, 1])
+ yrange <- range(tails[, 2], heads[, 2])
plot(xrange, yrange, xlab = xlab, ylab = ylab, main = main,
- type = "n", asp = 1, ...)
+ type = "n", asp = 1, ...)
if (kind > 0) {
- abline(v = 0, lty = 2)
+ abline(v = 0, lty = 2)
abline(h = 0, lty = 2)
if (ncol(x$rotation) == 2) {
## Draw rotated axes only if they visibly differ from
@@ -41,13 +48,18 @@
text(tmp[2,choices[1]], tmp[2,choices[2]], as.character(k))
}
}
- points(Yrot, ...)
- ow <- options(warn = -1)
- arrows(Yrot[, 1], Yrot[, 2], X[, 1], X[, 2], col = ar.col,
- len = len, ...)
- options(ow)
+ if (type != "none") {
+ ow <- options(warn = -1)
+ arrows(tails[, 1], tails[, 2], heads[, 1], heads[, 2],
+ col = ar.col, len = len, ...)
+ options(ow)
+ if (type == "text" && !is.null(rownames(tails)))
+ ordilabel(tails, cex = cex, ...)
+ else
+ points(tails, cex = cex, ...)
+ }
}
- out <- list(heads = X, points = Yrot)
+ out <- list(heads = heads, points = tails)
class(out) <- "ordiplot"
}
else if (kind == 2) {
@@ -58,7 +70,7 @@
res <- residuals(x)
q <- quantile(res)
plot(res, type = "h", xlab = xlab, ylab = ylab, main = main,
- ...)
+ ...)
abline(h = q[2:4], lty = c(2, 1, 2))
out <- list(sites = cbind(seq(along = res), res))
class(out) <- "ordiplot"
Modified: branches/1.17/R/prc.R
===================================================================
--- branches/1.17/R/prc.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/prc.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,7 +1,7 @@
-"prc" <-
- function (response, treatment, time, ...)
+`prc` <-
+ function (response, treatment, time, ...)
{
- extras <- match.call(expand.dots=FALSE)$...
+ extras <- match.call(expand.dots = FALSE)$...
if (is.null(extras$data))
data <- parent.frame()
else
@@ -9,15 +9,22 @@
y <- deparse(substitute(response))
x <- deparse(substitute(treatment))
z <- deparse(substitute(time))
+ oldcon <- options(contrasts = c("contr.treatment", "contr.poly"))
+ on.exit(options(oldcon))
fla <- as.formula(paste("~", x, "+", z))
- mf <- model.frame(fla, data)
+ mf <- model.frame(fla, data, na.action = na.pass)
if (!all(sapply(mf, is.factor)))
stop(x, " and ", z, " must be factors")
- if (any(sapply(mf, is.ordered)))
+ if (any(sapply(mf, is.ordered)))
stop(x, " or ", z, " cannot be ordered factors")
- fla <- as.formula(paste(y, "~", z, "*", x, "+ Condition(",
- z, ")"))
- mod <- rda(fla, ...)
+ fla.zx <- as.formula(paste("~", z, ":", x))
+ fla.z <- as.formula(paste("~", z))
+ # delete first (control) level from the design matrix
+ X = model.matrix(fla.zx, mf)[,-c(seq_len(nlevels(time)+1))]
+ Z = model.matrix(fla.z, mf)[,-1]
+ mod <- rda(response ~ X + Condition(Z), ...)
+ mod$terminfo$xlev = list(levels(time), levels(treatment))
+ names(mod$terminfo$xlev) = c(paste(z), paste(x))
mod$call <- match.call()
class(mod) <- c("prc", class(mod))
mod
Modified: branches/1.17/R/print.cca.R
===================================================================
--- branches/1.17/R/print.cca.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/print.cca.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -9,9 +9,10 @@
cat("\n")
chi <- c(x$tot.chi, x$pCCA$tot.chi, x$CCA$tot.chi, x$CA$tot.chi,
x$CA$imaginary.chi)
+ props <- abs(chi)/sum(abs(chi[-1]))
rnk <- c(NA, x$pCCA$rank, x$CCA$rank, x$CA$rank, x$CA$imaginary.rank)
- tbl <- cbind(chi, rnk)
- colnames(tbl) <- c("Inertia", "Rank")
+ tbl <- cbind(chi, props, rnk)
+ colnames(tbl) <- c("Inertia", "Proportion", "Rank")
rn <- c("Total", "Conditional", "Constrained", "Unconstrained",
"Imaginary")
rownames(tbl) <- rn[c(TRUE, !is.null(x$pCCA), !is.null(x$CCA),
Modified: branches/1.17/R/print.summary.prc.R
===================================================================
--- branches/1.17/R/print.summary.prc.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/print.summary.prc.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -5,7 +5,9 @@
cat(deparse(x$call), "\n")
cat("Species scores:\n")
print(x$sp, digits=x$digits, ...)
- cat("\nCoefficients for", paste(x$names, collapse=":"), "interaction\n")
+ cat("\nCoefficients for",
+ paste(x$names[2], "+", paste(x$names, collapse=":")),
+ "interaction\n")
cat(paste("which are contrasts to", x$names[2], x$corner, "\n"))
cat(paste(c("rows are",", columns are"), x$names[2:1], collapse=""))
cat("\n")
Modified: branches/1.17/R/procrustes.R
===================================================================
--- branches/1.17/R/procrustes.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/procrustes.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -29,10 +29,13 @@
c <- sum(sol$d)/ctrace(Y)
}
Yrot <- c * Y %*% A
- b <- xmean - t(A %*% ymean)
+ ## Translation (b) needs scale (c) although Mardia et al. do not
+ ## have this. Reported by Christian Dudel.
+ b <- xmean - c * ymean %*% A
R2 <- ctrace(X) + c * c * ctrace(Y) - 2 * c * sum(sol$d)
reslt <- list(Yrot = Yrot, X = X, ss = R2, rotation = A,
- translation = b, scale = c, symmetric = symmetric, call = match.call())
+ translation = b, scale = c, xmean = xmean,
+ symmetric = symmetric, call = match.call())
reslt$svd <- sol
class(reslt) <- "procrustes"
return(reslt)
Modified: branches/1.17/R/scores.rda.R
===================================================================
--- branches/1.17/R/scores.rda.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/scores.rda.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,4 +1,4 @@
-"scores.rda" <-
+`scores.rda` <-
function (x, choices = c(1, 2), display = c("sp", "wa", "cn"),
scaling = 2, const, ...)
{
@@ -25,8 +25,18 @@
nrow(x$CA$u)
else
nrow(x$CCA$u)
+ ## const multiplier of scores
if (missing(const))
const <- sqrt(sqrt((nr-1) * sumev))
+ ## canoco 3 compatibility -- canoco 4 is incompatible
+ ##else if (pmatch(const, "canoco")) {
+ ## const <- (sqrt(nr-1), sqrt(nr))
+ ##}
+ ##
+ ## const[1] for species, const[2] for sites and friends
+ if (length(const) == 1) {
+ const <- c(const, const)
+ }
rnk <- x$CCA$rank
sol <- list()
if ("species" %in% take) {
@@ -38,7 +48,7 @@
v <- sweep(v, 1, x$colsum, "/")
v <- v * sqrt(sumev / (nr - 1))
}
- v <- const * v
+ v <- const[1] * v
}
sol$species <- v
}
@@ -47,7 +57,7 @@
if (scaling) {
scal <- list(slam, 1, sqrt(slam))[[abs(scaling)]]
wa <- sweep(wa, 2, scal, "*")
- wa <- const * wa
+ wa <- const[2] * wa
}
sol$sites <- wa
}
@@ -56,7 +66,7 @@
if (scaling) {
scal <- list(slam, 1, sqrt(slam))[[abs(scaling)]]
u <- sweep(u, 2, scal, "*")
- u <- const * u
+ u <- const[2] * u
}
sol$constraints <- u
}
@@ -78,7 +88,7 @@
if (scaling) {
scal <- list(slam, 1, sqrt(slam))[[abs(scaling)]]
cn <- sweep(cn, 2, scal, "*")
- cn <- const * cn
+ cn <- const[2] * cn
}
sol$centroids <- cn
}
@@ -93,6 +103,9 @@
## Only one type of scores: return a matrix instead of a list
if (length(sol) == 1)
sol <- sol[[1]]
+ ## collapse const if both items identical
+ if (identical(const[1], const[2]))
+ const <- const[1]
attr(sol, "const") <- const
return(sol)
}
Modified: branches/1.17/R/summary.prc.R
===================================================================
--- branches/1.17/R/summary.prc.R 2010-12-01 19:46:21 UTC (rev 1387)
+++ branches/1.17/R/summary.prc.R 2010-12-02 09:38:50 UTC (rev 1388)
@@ -1,20 +1,19 @@
-"summary.prc" <-
- function (object, axis = 1, scaling = 2, digits = 4, ...)
+`summary.prc` <-
+ function (object, axis = 1, scaling = 3, digits = 4, ...)
{
- species <- drop(scores(object, scaling = scaling, display="sp", choices=axis))
- b <- coef(object)[, axis]
+ sc = scores(object, scaling = scaling, display = c("sp", "lc"),
+ choices=axis, ...)
+ ## coef for scaled sites (coef(object) gives for orthonormal)
+ b <- qr.coef(object$CCA$QR, sc$constraints)
prnk <- object$pCCA$rank
lentreat <- length(object$terminfo$xlev[[2]])
- lenb <- length(b)
- b <- b[-(1:(2 * prnk))]
- bx <- b[1:(lentreat - 1)]
- by <- b[lentreat:length(b)]
- b <- cbind(bx, matrix(by, nrow = lentreat - 1, byrow = TRUE))
+ b = matrix(b[-(1:prnk)], nrow = lentreat-1, byrow = TRUE)
rownames(b) <- (object$terminfo$xlev[[2]])[-1]
colnames(b) <- object$terminfo$xlev[[1]]
- out <- list(sp = species, coefficients = b, names = names(object$terminfo$xlev),
- corner = (object$terminfo$xlev[[2]])[1], call = object$call,
- digits = digits)
+ out <- list(sp = drop(sc$species), coefficients = b,
+ names = names(object$terminfo$xlev),
+ corner = (object$terminfo$xlev[[2]])[1],
+ call = object$call, digits = digits)
class(out) <- "summary.prc"
out
}
Modified: branches/1.17/R/treedist.R
===================================================================
--- branches/1.17/R/treedist.R 2010-12-01 19:46:21 UTC (rev 1387)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vegan -r 1388
More information about the Vegan-commits
mailing list