[Vegan-commits] r1552 - in pkg/vegan: . R inst tests tests/Examples
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 23 12:09:21 CET 2011
Author: jarioksa
Date: 2011-03-23 12:09:21 +0100 (Wed, 23 Mar 2011)
New Revision: 1552
Modified:
pkg/vegan/DESCRIPTION
pkg/vegan/R/anova.ccabyaxis.R
pkg/vegan/inst/ChangeLog
pkg/vegan/tests/Examples/vegan-Ex.Rout.save
pkg/vegan/tests/vegan-tests.R
pkg/vegan/tests/vegan-tests.Rout.save
Log:
bug fix: anova.ccabyaxis gave wrong results in partial models (and version to 1.18-26)
Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION 2011-03-23 05:39:33 UTC (rev 1551)
+++ pkg/vegan/DESCRIPTION 2011-03-23 11:09:21 UTC (rev 1552)
@@ -1,7 +1,7 @@
Package: vegan
Title: Community Ecology Package
-Version: 1.18-25
-Date: March 10, 2011
+Version: 1.18-26
+Date: March 23, 2011
Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre,
R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens,
Helene Wagner
Modified: pkg/vegan/R/anova.ccabyaxis.R
===================================================================
--- pkg/vegan/R/anova.ccabyaxis.R 2011-03-23 05:39:33 UTC (rev 1551)
+++ pkg/vegan/R/anova.ccabyaxis.R 2011-03-23 11:09:21 UTC (rev 1552)
@@ -1,10 +1,6 @@
"anova.ccabyaxis" <-
function (object, cutoff = 1, ...)
{
- ## FIXME: wrong results in partial model
- ## BUG -- NEEDS AN URGENT FIX
- if (!is.null(object$pCCA))
- stop("anova by 'axis' does not work in partial models: call the developers")
cutoff <- cutoff + sqrt(.Machine$double.eps)
rnk <- object$CCA$rank
if (!max(rnk, 0))
@@ -18,11 +14,26 @@
} else {
u <- object$CCA$u
}
+ ## Get conditions
+ if (!is.null(object$pCCA)) {
+ CondMat <- qr.X(object$pCCA$QR)
+ ## deweight if CCA
+ if (!inherits(object, "rda"))
+ CondMat <- sweep(CondMat, 1, sqrt(object$rowsum), "/")
+ }
+ else
+ CondMat <- NULL
## pad with NA rows if there is a subset
if (!is.null(object$subset)) {
lc <- matrix(NA, nrow=length(object$subset),
- ncol = ncol(u))
+ ncol = NCOL(u))
lc[object$subset,] <- u
+ if (!is.null(CondMat)) {
+ tmp <- matrix(NA, nrow=length(object$subset),
+ ncol = NCOL(CondMat))
+ tmp[object$subset,] <- CondMat
+ CondMat <- tmp
+ }
object$call$subset <- object$subset
} else {
lc <- u
@@ -38,6 +49,8 @@
environment(object$terms) <- environment()
fla <- paste(". ~ ", axnam[1], "+ Condition(",
paste(axnam[-1], collapse="+"),")")
+ if (!is.null(CondMat))
+ fla <- paste(fla, " + Condition(CondMat)")
fla <- update(formula(object), fla)
sol <- anova(update(object, fla, data=lc), ...)
out[c(1, rnk + 1), ] <- sol
@@ -52,11 +65,12 @@
for (.ITRM in 2:rnk) {
fla <- paste(".~", axnam[.ITRM], "+Condition(",
paste(axnam[-(.ITRM)], collapse="+"),")")
+ if (!is.null(CondMat))
+ fla <- paste(fla, "+ Condition(CondMat)")
fla <- update(formula(object), fla)
sol <- update(object, fla, data = lc)
assign(".Random.seed", seed, envir = .GlobalEnv)
- out[.ITRM, ] <- as.matrix(anova(sol, ...))[1,
- ]
+ out[.ITRM, ] <- as.matrix(anova(sol, ...))[1,]
if (out[.ITRM, "N.Perm"] > bigperm) {
bigperm <- out[.ITRM, "N.Perm"]
bigseed <- get(".Random.seed", envir = .GlobalEnv,
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2011-03-23 05:39:33 UTC (rev 1551)
+++ pkg/vegan/inst/ChangeLog 2011-03-23 11:09:21 UTC (rev 1552)
@@ -2,14 +2,17 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
-Version 1.18-25 (opened March 10, 2011)
+Version 1.18-26 (opened March 23, 2011)
* anova.ccabyaxis: I (JO) noticed in an ORDNEWS message that
- anova(..., by = "axis") drops the original conditions in partial
- models and so gives wrong results. The first thing here prevents
- the analysis, and a final fix brings along a new release (this
- also breaks the tests).
+ anova(..., by = "axis") ignores the original Conditions in partial
+ models and so gives wrong results in partial models. The fix finds
+ the conditioning matrix as qr.X(object$pCCA$QR), and uses this as
+ a Condition in updated models. Added a regression test that checks
+ that statistics and residual df match.
+Version 1.18-25 (closed March 23, 2011)
+
* ordilabel: gained argument 'xpd' to draw labels outside the plot
region.
Modified: pkg/vegan/tests/Examples/vegan-Ex.Rout.save
===================================================================
--- pkg/vegan/tests/Examples/vegan-Ex.Rout.save 2011-03-23 05:39:33 UTC (rev 1551)
+++ pkg/vegan/tests/Examples/vegan-Ex.Rout.save 2011-03-23 11:09:21 UTC (rev 1552)
@@ -22,7 +22,7 @@
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('vegan')
-This is vegan 1.18-25
+This is vegan 1.18-26
>
> assign(".oldSearch", search(), pos = 'CheckExEnv')
> cleanEx()
@@ -4668,7 +4668,7 @@
Formula:
y ~ s(x1, x2, k = knots)
-<environment: 0x103156c00>
+<environment: 0x10234bed0>
Estimated degrees of freedom:
6.2955 total = 7.295494
@@ -4684,7 +4684,7 @@
Formula:
y ~ s(x1, x2, k = knots)
-<environment: 0x103286148>
+<environment: 0x10538e238>
Estimated degrees of freedom:
4.9207 total = 5.920718
@@ -4840,7 +4840,7 @@
Formula:
y ~ s(x1, x2, k = knots)
-<environment: 0x104eeb520>
+<environment: 0x1064ab778>
Estimated degrees of freedom:
8.9275 total = 9.927492
@@ -4853,7 +4853,7 @@
Formula:
y ~ s(x1, x2, k = knots)
-<environment: 0x107318ed0>
+<environment: 0x1067618b0>
Estimated degrees of freedom:
7.753 total = 8.75294
@@ -4866,7 +4866,7 @@
Formula:
y ~ s(x1, x2, k = knots)
-<environment: 0x1071ce4a0>
+<environment: 0x106033698>
Estimated degrees of freedom:
8.8962 total = 9.89616
@@ -7495,7 +7495,7 @@
Formula:
y ~ s(x1, x2, k = knots)
-<environment: 0x10456e8b0>
+<environment: 0x10718d600>
Estimated degrees of freedom:
2 total = 3
@@ -7972,7 +7972,7 @@
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
-Time elapsed: 118.697 1.564 126.298 0 0
+Time elapsed: 115.054 1.293 119.552 0 0
> grDevices::dev.off()
null device
1
Modified: pkg/vegan/tests/vegan-tests.R
===================================================================
--- pkg/vegan/tests/vegan-tests.R 2011-03-23 05:39:33 UTC (rev 1551)
+++ pkg/vegan/tests/vegan-tests.R 2011-03-23 11:09:21 UTC (rev 1552)
@@ -55,8 +55,16 @@
anova(q, by="margin", perm=100)
anova(q, by="axis", perm=100)
detach(df)
+### Check that statistics match in partial constrained ordination
+m <- cca(dune ~ A1 + Moisture + Condition(Management), dune.env, subset = A1 > 3)
+tab <- anova(m, by = "axis", perm.max = 100)
+m
+tab
+all.equal(tab[,2], c(m$CCA$eig, m$CA$tot.chi), check.attributes=FALSE)
+tab[nrow(tab),1] == m$CA$rank
+
## clean-up
-rm(df, spno, fla, m, p, q, .Random.seed)
+rm(df, spno, fla, m, p, q, tab, .Random.seed)
### <--- END anova.cca test --->
### nestednodf: test case by Daniel Spitale in a comment to News on
Modified: pkg/vegan/tests/vegan-tests.Rout.save
===================================================================
--- pkg/vegan/tests/vegan-tests.Rout.save 2011-03-23 05:39:33 UTC (rev 1551)
+++ pkg/vegan/tests/vegan-tests.Rout.save 2011-03-23 11:09:21 UTC (rev 1552)
@@ -1,6 +1,6 @@
-R version 2.12.1 (2010-12-16)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
@@ -190,8 +190,45 @@
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> detach(df)
+> ### Check that statistics match in partial constrained ordination
+> m <- cca(dune ~ A1 + Moisture + Condition(Management), dune.env, subset = A1 > 3)
+> tab <- anova(m, by = "axis", perm.max = 100)
+> m
+Call: cca(formula = dune ~ A1 + Moisture + Condition(Management), data
+= dune.env, subset = A1 > 3)
+
+ Inertia Proportion Rank
+Total 2.0976 1.0000
+Conditional 0.6251 0.2980 3
+Constrained 0.5555 0.2648 4
+Unconstrained 0.9170 0.4372 10
+Inertia is mean squared contingency coefficient
+
+Eigenvalues for constrained axes:
+ CCA1 CCA2 CCA3 CCA4
+0.27109 0.14057 0.08761 0.05624
+
+Eigenvalues for unconstrained axes:
+ CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10
+0.31042 0.13634 0.11974 0.09408 0.07763 0.06425 0.04449 0.02925 0.02785 0.01299
+
+> tab
+Model: cca(formula = dune ~ A1 + Moisture + Condition(Management), data = dune.env, subset = c(TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE))
+ Df Chisq F N.Perm Pr(>F)
+CCA1 1 0.2711 2.9561 99 0.02 *
+CCA2 1 0.1406 1.5329 99 0.10 .
+CCA3 1 0.0876 0.9553 99 0.46
+CCA4 1 0.0562 0.6132 99 0.78
+Residual 10 0.9170
+---
+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+> all.equal(tab[,2], c(m$CCA$eig, m$CA$tot.chi), check.attributes=FALSE)
+[1] TRUE
+> tab[nrow(tab),1] == m$CA$rank
+[1] TRUE
+>
> ## clean-up
-> rm(df, spno, fla, m, p, q, .Random.seed)
+> rm(df, spno, fla, m, p, q, tab, .Random.seed)
> ### <--- END anova.cca test --->
>
> ### nestednodf: test case by Daniel Spitale in a comment to News on
More information about the Vegan-commits
mailing list