[Vegan-commits] r322 - in branches/1.11-0: . R inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Apr 20 17:44:59 CEST 2008
Author: jarioksa
Date: 2008-04-20 17:44:59 +0200 (Sun, 20 Apr 2008)
New Revision: 322
Modified:
branches/1.11-0/DESCRIPTION
branches/1.11-0/R/adonis.R
branches/1.11-0/inst/ChangeLog
branches/1.11-0/inst/NEWS
Log:
merged adonis fixes to the 1.11 branch
Modified: branches/1.11-0/DESCRIPTION
===================================================================
--- branches/1.11-0/DESCRIPTION 2008-04-20 06:26:51 UTC (rev 321)
+++ branches/1.11-0/DESCRIPTION 2008-04-20 15:44:59 UTC (rev 322)
@@ -1,7 +1,7 @@
Package: vegan
Title: Community Ecology Package
Version: 1.11-4
-Date: April 10, 2008
+Date: April 20, 2008
Author: Jari Oksanen, Roeland Kindt, Pierre Legendre, Bob O'Hara, Gavin L. Simpson,
M. Henry H. Stevens
Maintainer: Jari Oksanen <jari.oksanen at oulu.fi>
Modified: branches/1.11-0/R/adonis.R
===================================================================
--- branches/1.11-0/R/adonis.R 2008-04-20 06:26:51 UTC (rev 321)
+++ branches/1.11-0/R/adonis.R 2008-04-20 15:44:59 UTC (rev 322)
@@ -16,9 +16,13 @@
rhs <- model.matrix(formula, rhs.frame) # and finally the model.matrix
options(contrasts=op.c)
grps <- attr(rhs, "assign")
+ qrhs <- qr(rhs)
+ ## Take care of aliased variables and pivoting in rhs
+ rhs <- rhs[, qrhs$pivot, drop=FALSE]
+ rhs <- rhs[, 1:qrhs$rank, drop=FALSE]
+ grps <- grps[qrhs$pivot][1:qrhs$rank]
u.grps <- unique(grps)
- nterms <- max(u.grps)
-
+ nterms <- length(u.grps) - 1
H.s <- lapply(2:length(u.grps),
function(j) {Xj <- rhs[, grps %in% u.grps[1:j] ]
qrX <- qr(Xj, tol=1e-7)
@@ -33,22 +37,23 @@
I <- diag(n)
ones <- matrix(1,nrow=n)
A <- -(dmat)/2
-
G <- -.5 * dmat %*% (I - ones%*%t(ones)/n)
-
SS.Exp.comb <- sapply(H.s, function(hat) sum( diag(G %*% hat) ) )
-
SS.Exp.each <- c(SS.Exp.comb - c(0,SS.Exp.comb[-nterms]) )
- SS.Res <- sum(diag( ( G %*% (I-H.s[[nterms]]) )))
+ H.snterm <- H.s[[nterms]]
+ if (length(H.s) > 1)
+ for (i in length(H.s):2)
+ H.s[[i]] <- H.s[[i]] - H.s[[i-1]]
+ SS.Res <- sum(diag( ( G %*% (I-H.snterm))))
df.Exp <- sapply(u.grps[-1], function(i) sum(grps==i) )
- df.Res <- n - dim(rhs)[2]
- beta <- qr.coef( qr(rhs), as.matrix(lhs) )
+ df.Res <- n - qrhs$rank
+ beta <- qr.coef(qrhs, as.matrix(lhs) )
colnames(beta) <- colnames(lhs)
F.Mod <- (SS.Exp.each/df.Exp) / (SS.Res/df.Res)
- f.test <- function(H, G, I, df.Exp, df.Res){
+ f.test <- function(H, G, I, df.Exp, df.Res, H.snterm){
(sum( diag(G %*% H) )/df.Exp) /
- ( sum(diag( G %*% (I-H) ))/df.Res) }
+ ( sum(diag( G %*% (I-H.snterm) ))/df.Res) }
SS.perms <- function(H, G, I){
c(SS.Exp.p = sum( diag(G%*%H) ),
@@ -65,7 +70,7 @@
## SS.s <- sapply(G.p, function(Gs) { SS.perms(H, Gs, I) } )
f.perms <- sapply(1:nterms, function(i) {
sapply(1:permutations, function(j) {
- f.test(H.s[[i]], G.p[[j]], I, df.Exp[i], df.Res)
+ f.test(H.s[[i]], G.p[[j]], I, df.Exp[i], df.Res, H.snterm)
} )
})
SumsOfSqs = c(SS.Exp.each, SS.Res, sum(SS.Exp.each) + SS.Res)
@@ -75,7 +80,7 @@
F.Model = c(F.Mod, NA,NA),
R2 = SumsOfSqs/SumsOfSqs[length(SumsOfSqs)],
P = c(rowSums(t(f.perms) > F.Mod)/permutations, NA, NA))
- rownames(tab) <- c(attr(attr(rhs.frame, "terms"), "term.labels"),
+ rownames(tab) <- c(attr(attr(rhs.frame, "terms"), "term.labels")[u.grps],
"Residuals", "Total")
colnames(tab)[ncol(tab)] <- "Pr(>F)"
out <- list(aov.tab = tab, call = match.call(),
Modified: branches/1.11-0/inst/ChangeLog
===================================================================
--- branches/1.11-0/inst/ChangeLog 2008-04-20 06:26:51 UTC (rev 321)
+++ branches/1.11-0/inst/ChangeLog 2008-04-20 15:44:59 UTC (rev 322)
@@ -2,9 +2,11 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
-Version 1.11-4 (Opened April 10th, 2008)
+Version 1.11-4 (Opened April 10, 2008)
- * adonis: merged Rev 312 (not dropping unused factor levels)
+ * adonis: merged Rev 312 (not dropping unused factor levels), Rev
+ 315 (wrong df in deficit rank models), Rev 320 (wrong statistics
+ in permutations: a critical bug).
Version 1.11-3 (April 8, 2008)
Modified: branches/1.11-0/inst/NEWS
===================================================================
--- branches/1.11-0/inst/NEWS 2008-04-20 06:26:51 UTC (rev 321)
+++ branches/1.11-0/inst/NEWS 2008-04-20 15:44:59 UTC (rev 322)
@@ -2,6 +2,14 @@
VEGAN RELEASE VERSIONS
+CHANGES IN VEGAN VERSION 1.11-4
+
+ - adonis: there was a critical bug in adonis code, and wrong (and
+ different) statistics were used in permutations so that P-values
+ were grossly wrong in multi-variables models (single variable
+ models were OK). In addition, df were wrong in deficit rank
+ models, and unused factor levels were not dropped.
+
CHANGES IN VEGAN VERSION 1.11-3
GENERAL
More information about the Vegan-commits
mailing list