diff --git a/NAMESPACE b/NAMESPACE
index 00998520..fc01fcdf 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -230,6 +230,7 @@ S3method(logLik, radfit.frame)
 # model.frame, model.matrix: stats
 S3method(model.frame, cca)
 S3method(model.matrix, cca)
+S3method(model.matrix, rda)
 # multipart: vegan
 S3method(multipart, default)
 S3method(multipart, formula)
diff --git a/R/adonis2.R b/R/adonis2.R
index 159545d2..d61e005d 100644
--- a/R/adonis2.R
+++ b/R/adonis2.R
@@ -78,6 +78,8 @@
     formula[[2]] <- NULL                # to remove the lhs
     rhs.frame <- model.frame(formula, data, drop.unused.levels = TRUE) # to get the data frame of rhs
     rhs <- model.matrix(formula, rhs.frame) # and finally the model.matrix
+    assign <- attr(rhs, "assign") ## assign attribute
+    sol$terminfo$assign <- assign[assign > 0]
     rhs <- rhs[,-1, drop=FALSE] # remove the (Intercept) to get rank right
     rhs <- scale(rhs, scale = FALSE, center = TRUE) # center
     qrhs <- qr(rhs)
diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R
index 7355327b..e146e031 100644
--- a/R/anova.ccabyterm.R
+++ b/R/anova.ccabyterm.R
@@ -1,43 +1,31 @@
-### Implementation of by-cases for vegan 2.2 versions of
-### anova.cca. These are all internal functions that are not intended
-### to be called by users in normal sessions, but they should be
-### called from anova.cca (2.2). Therefore the user interface is rigid
-### and input is not checked. The 'permutations' should be a
-### permutation matrix.
+### Implementation of by-cases for anova.cca. These are all internal
+### functions that are not intended to be called by users in normal
+### sessions, but they should be called from anova.cca. Therefore the
+### user interface is rigid and input is not checked. The
+### 'permutations' should be a permutation matrix.
 
-### by = terms builds models as a sequence of adding terms and submits
-### this to anova.ccalist
+### by = "terms" calls directly permutest.cca which decomposes the
+### inertia between successive terms within compiled C code.
 
 `anova.ccabyterm` <-
     function(object, permutations, model, parallel)
 {
-    ## We need term labels but without Condition() terms
-    trms <- terms(object)
-    trmlab <- attr(trms, "term.labels")
-    trmlab <- trmlab[trmlab %in% attr(terms(object$terminfo),
-                                      "term.labels")]
-    ntrm <- length(trmlab)
-    m0 <- update(object, paste(".~.-", paste(trmlab, collapse = "-")))
-    mods <- list(m0)
-    for (i in seq_along(trmlab)) {
-        fla <- paste(". ~ . + ", trmlab[i])
-        mods[[i+1]] <- update(mods[[i]], fla)
-    }
     ## The result
-    sol <- anova.ccalist(mods, permutations = permutations,
-                         model = model, parallel = parallel)
+    sol <- permutest(object, permutations = permutations,
+                     model = model, by = "terms", parallel = parallel)
     ## Reformat
-    out <- data.frame(c(sol[-1, 3], sol[ntrm+1, 1]),
-                      c(sol[-1, 4], sol[ntrm+1, 2]),
-                      c(sol[-1, 5], NA),
-                      c(sol[-1, 6], NA))
+    EPS <- sqrt(.Machine$double.eps)
+    Pval <- (colSums(sweep(sol$F.perm, 2, sol$F.0 - EPS, ">=")) + 1) /
+        (sol$nperm + 1)
+    out <- data.frame(sol$df, sol$chi, c(sol$F.0, NA), c(Pval, NA))
+
     if (inherits(object, c("capscale", "dbrda")) && object$adjust == 1)
         varname <- "SumOfSqs"
     else if (inherits(object, "rda"))
         varname <- "Variance"
     else
         varname <- "ChiSquare"
-    dimnames(out) <- list(c(trmlab, "Residual"),
+    dimnames(out) <- list(c(sol$termlabels, "Residual"),
                           c("Df", varname, "F", "Pr(>F)"))
     head <- paste0("Permutation test for ", object$method, " under ",
                    model, " model\n",
@@ -45,22 +33,28 @@
                    howHead(attr(permutations, "control")))
     mod <- paste("Model:", c(object$call))
     attr(out, "heading") <- c(head, mod)
-    attr(out, "F.perm") <- attr(sol, "F.perm")
+    attr(out, "F.perm") <- sol$F.perm
     class(out) <- c("anova.cca", "anova","data.frame")
     out
 }
 
-## by = margin: this is not a anova.ccalist case, but we omit each
-## term in turn and compare against the complete model.
+## by = "margin": we omit each term in turn and compare against the
+## complete model. This does not involve partial terms (Conditions) on
+## other variables, but the permutations remain similar in "direct"
+## and "reduced" (default) models (perhaps this model should not be
+## used with "full" models?). This is basically similar decomposition
+## as by="term", but compares models without each term in turn against
+## the complete model in separate calls to permutest.cca. From vegan
+## 2.5-0 this does not update model formula -- this avoids scoping
+## issues and makes the function more robust when embedded in other
+## functions. Instead, we call ordConstrained with method="pass" with
+## modified constraint matrix.
 
 `anova.ccabymargin` <-
     function(object, permutations, scope, ...)
 {
     EPS <- sqrt(.Machine$double.eps)
     nperm <- nrow(permutations)
-    ## Refuse to handle models with missing data
-    if (!is.null(object$na.action))
-        stop("by = 'margin' models cannot handle missing data")
     ## We need term labels but without Condition() terms
     if (!is.null(scope) && is.character(scope))
         trms <- scope
@@ -79,8 +73,21 @@
     ## (vegan 2.0) where other but 'nm' were partialled out within
     ## Condition(). Now we only fit the model without 'nm' and compare
     ## the difference against the complete model.
-    mods <- lapply(trmlab, function(nm, ...)
-           permutest(update(object, paste(".~.-", nm)),
+    Y <- ordiYbar(object, "init")
+    X <- model.matrix(object)
+    ## we must have Constraints to get here, but we may also have
+    ## Conditions
+    if (!is.null(object$pCCA)) {
+        Z <- X$Conditions
+        X <- X$Constraints
+    } else {
+        Z <- NULL
+    }
+    ass <- object$terminfo$assign
+    if (is.null(ass))
+        stop("old time result object: update() your model")
+    mods <- lapply(unique(ass), function(i, ...)
+           permutest(ordConstrained(Y, X[, ass != i, drop=FALSE], Z, "pass"),
                      permutations, ...), ...)
     ## Chande in df
     Df <- sapply(mods, function(x) x$df[2]) - dfbig
@@ -119,7 +126,20 @@
     out
 }
 
-### Marginal test for axes
+### by = "axis" uses partial model: we use the original constraints,
+### but add previous axes 1..(k-1) to Conditions when evaluating the
+### significance of axis k which is compared against the first
+### eigenvalue of the permutations. To avoid scoping issues, this
+### calls directly ordConstrained() with modified Conditions (Z) and
+### original Constraints (X) instead of updating formula. This
+### corresponds to "forward" model in Legendre, Oksanen, ter Braak
+### (2011).
+
+### In 2.2-x to 2.4-3 we used "marginal model" where original
+### Constraints were replaced with LC scores axes (object$CCA$u), and
+### all but axis k were used as Conditions when evaluating the
+### significance of axis k. My (J.Oksanen) simulations showed that
+### this gave somewhat biased results.
 
 `anova.ccabyaxis` <-
     function(object, permutations, model, parallel, cutoff = 1)
@@ -137,45 +157,52 @@
     resdf <- nobs(object) - length(eig) - max(object$pCCA$rank, 0) - 1
     Fstat <- eig/object$CA$tot.chi*resdf
     Df <- rep(1, length(eig))
-    ## Marginal P-values
-    LC <- object$CCA$u
-    ## missing values?
-    if (!is.null(object$na.action))
-        LC <- napredict(structure(object$na.action,
-                                  class = "exclude"), LC)
-    ## subset?
-    if (!is.null(object$subset)) {
-        tmp <- matrix(NA, nrow = length(object$subset),
-                      ncol = ncol(LC))
-        tmp[object$subset,] <- LC
-        LC <- tmp
-        object <- update(object, subset = object$subset)
+
+    ## collect header and varname here: 'object' is modified later
+    if (inherits(object, c("capscale", "dbrda")) && object$adjust == 1)
+        varname <- "SumOfSqs"
+    else if (inherits(object, "rda"))
+        varname <- "Variance"
+    else
+        varname <- "ChiSquare"
+
+    head <- paste0("Permutation test for ", object$method, " under ",
+                   model, " model\n",
+                   "Forward tests for axes\n",
+                   howHead(attr(permutations, "control")))
+    head <- c(head, paste("Model:", c(object$call)))
+
+    ## constraints and model matrices
+    Y <- object$Ybar
+    if (is.null(Y))
+        stop("old style result object does not work: update() your model")
+    if (!is.null(object$pCCA))
+        Z <- qr.X(object$pCCA$QR)
+    else
+        Z <- NULL
+    X <- model.matrix(object)
+    if (!is.null(object$pCCA)) {
+        Z <- X$Conditions
+        X <- X$Constraints
+    } else {
+        Z <- NULL
     }
-    LC <- as.data.frame(LC)
-    fla <- reformulate(names(LC))
+    LC <- object$CCA$u
+
     Pvals <- rep(NA, ncol(LC))
     F.perm <- matrix(ncol = ncol(LC), nrow = nperm)
-    environment(object$terms) <- environment()
-    ## in dbrda, some axes can be imaginary, but we only want to have
-    ## an analysis of real-valued dimensions, and we must adjust data
-    if (ncol(LC) < length(eig)) {
-        eig <- eig[seq_len(ncol(LC))]
-        Df <- Df[seq_len(ncol(LC))]
-        Fstat <- Fstat[seq_len(ncol(LC))]
-    }
+    axnams <- colnames(LC)
     for (i in seq_along(eig)) {
-        part <- paste("~ . +Condition(",
-                      paste(names(LC)[-i], collapse = "+"), ")")
-        upfla <- update(fla, part)
-        ## only one axis, and cannot partial out?
-        if (length(eig) == 1)
+        if (i > 1) {
+            object <- ordConstrained(Y, X, cbind(Z, LC[, seq_len(i-1)]), "pass")
+        }
+        if (length(eig) == i) {
             mod <- permutest(object, permutations, model = model,
                              parallel = parallel)
-        else
-            mod <-
-                permutest(update(object, upfla, data = LC),
-                          permutations, model = model,
-                          parallel = parallel)
+        } else {
+            mod <- permutest(object, permutations, model = model,
+                             parallel = parallel, first = TRUE)
+        }
         Pvals[i] <- (sum(mod$F.perm >= mod$F.0 - EPS) + 1) / (nperm + 1)
         F.perm[ , i] <- mod$F.perm
         if (Pvals[i] > cutoff)
@@ -184,19 +211,8 @@
     out <- data.frame(c(Df, resdf), c(eig, object$CA$tot.chi),
                       c(Fstat, NA), c(Pvals,NA))
     rownames(out) <- c(names(eig), "Residual")
-    if (inherits(object, c("capscale", "dbrda")) && object$adjust == 1)
-        varname <- "SumOfSqs"
-    else if (inherits(object, "rda"))
-        varname <- "Variance"
-    else
-        varname <- "ChiSquare"
     colnames(out) <- c("Df", varname, "F", "Pr(>F)")
-    head <- paste0("Permutation test for ", object$method, " under ",
-                   model, " model\n",
-                   "Marginal tests for axes\n",
-                   howHead(attr(permutations, "control")))
-    mod <- paste("Model:", c(object$call))
-    attr(out, "heading") <- c(head, mod)
+    attr(out, "heading") <- head
     attr(out, "F.perm") <- F.perm
     class(out) <- c("anova.cca", "anova", "data.frame")
     out
diff --git a/R/anova.ccalist.R b/R/anova.ccalist.R
index bf98358e..0da55626 100644
--- a/R/anova.ccalist.R
+++ b/R/anova.ccalist.R
@@ -86,7 +86,8 @@
                             c("ResDf", paste0("Res", varname), "Df",
                               varname, "F", "Pr(>F)"))
     ## Collect header information
-    formulae <- sapply(object, function(z) deparse(formula(z)))
+    formulae <- sapply(object,
+                       function(z) deparse(formula(z), width.cutoff = 500))
     head <- paste0("Permutation tests for ", method, " under ",
                   mods[[big]]$model, " model\n",
                    howHead(attr(permutations, "control")))
diff --git a/R/capscale.R b/R/capscale.R
index 2dad29b0..895a00dd 100644
--- a/R/capscale.R
+++ b/R/capscale.R
@@ -143,7 +143,6 @@
     sol$terminfo <- ordiTerminfo(d, data)
     sol$call$formula <- formula(d$terms, width.cutoff = 500)
     sol$call$formula[[2]] <- formula[[2]]
-    sol$method <- "capscale"
     sol$sqrt.dist <- sqrt.dist
     if (!is.na(X$ac) && X$ac > 0) {
         sol$ac <- X$ac
diff --git a/R/cca.default.R b/R/cca.default.R
index 3cad2a21..611e0a21 100644
--- a/R/cca.default.R
+++ b/R/cca.default.R
@@ -32,7 +32,6 @@
     if (!is.null(sol$pCCA) && sol$pCCA$tot.chi == 0)
         pCCA$rank <- 0
     sol <- c(list(call = call,
-                  method = "cca",
                   inertia =  "mean squared contingency coefficient"),
              sol)
     class(sol) <- "cca"
diff --git a/R/dbrda.R b/R/dbrda.R
index 32976049..2c6b067f 100644
--- a/R/dbrda.R
+++ b/R/dbrda.R
@@ -134,7 +134,6 @@
     sol$terminfo <- ordiTerminfo(d, data)
     sol$call$formula <- formula(d$terms, width.cutoff = 500)
     sol$call$formula[[2]] <- formula[[2]]
-    sol$method <- "dbrda"
     sol$sqrt.dist <- sqrt.dist
     if (!is.na(ac) && ac > 0) {
         sol$ac <- ac
diff --git a/R/model.frame.cca.R b/R/model.frame.cca.R
index ecabd493..21e78ec6 100644
--- a/R/model.frame.cca.R
+++ b/R/model.frame.cca.R
@@ -1,17 +1,20 @@
 `model.frame.cca` <-
-    function (formula, ...) 
+    function (formula, ...)
 {
     if (inherits(formula, "prc"))
         stop("model.frame does not work with 'prc' results")
     call <- formula$call
-    m <- match(c("formula", "data", "na.action", "subset"), names(call), 
+    m <- match(c("formula", "data", "na.action", "subset"), names(call),
         0)
     call <- call[c(1, m)]
+    ## subset must be evaluated before ordiParseFormula
+    if (!is.null(call$subset))
+        call$subset <- formula$subset
     call[[1]] <- as.name("ordiParseFormula")
     out <- eval(call, environment(), parent.frame())
     mf <- out$modelframe
     attr(mf, "terms") <- out$terms.expand
-    if (!is.null(out$na.action)) 
+    if (!is.null(out$na.action))
         attr(mf, "na.action") <- out$na.action
     mf
 }
diff --git a/R/model.matrix.cca.R b/R/model.matrix.cca.R
index e2083fbc..4cd0d91b 100644
--- a/R/model.matrix.cca.R
+++ b/R/model.matrix.cca.R
@@ -1,19 +1,43 @@
 `model.matrix.cca` <-
-    function (object, ...) 
+    function(object, ...)
 {
-    if (inherits(object, "prc"))
-        stop("model.matrix does not work with 'prc' results")
-    call <- object$call
-    m <- match(c("formula", "data", "na.action", "subset"), names(call), 
-        0)
-    call <- call[c(1, m)]
-    call[[1]] <- as.name("ordiParseFormula")
-    out <- eval(call, environment(), parent.frame())[c("Z", "Y")]
+    X <- Z <- NULL
+    w <- 1/sqrt(object$rowsum)
+    if (!is.null(object$pCCA))
+        Z <- w * qr.X(object$pCCA$QR)
+    if (!is.null(object$CCA)) {
+        X <- qr.X(object$CCA$QR)
+        ## First columns come from Z
+        if (!is.null(Z))
+            X <- X[, -seq_len(ncol(Z)), drop = FALSE]
+        X <- w * X
+    }
     m <- list()
-    if (!is.null(out$Z))
-        m$Conditions <- out$Z
-    if (!is.null(out$Y))
-        m$Constraints <- out$Y
+    if (!is.null(Z))
+        m$Conditions <- Z
+    if (!is.null(X))
+        m$Constraints <- X
+    if (length(m) == 1)
+        m <- m[[1]]
+    m
+}
+
+`model.matrix.rda` <-
+    function(object, ...)
+{
+    X <- Z <- NULL
+    if (!is.null(object$pCCA))
+        Z <- qr.X(object$pCCA$QR)
+    if (!is.null(object$CCA)) {
+        X <- qr.X(object$CCA$QR)
+        if (!is.null(Z))
+            X <- X[, -seq_len(ncol(Z)), drop=FALSE]
+    }
+    m <- list()
+    if (!is.null(Z))
+        m$Conditions <- Z
+    if (!is.null(X))
+        m$Constraints <- X
     if (length(m) == 1)
         m <- m[[1]]
     m
diff --git a/R/ordConstrained.R b/R/ordConstrained.R
index 8602367c..3307e5b5 100644
--- a/R/ordConstrained.R
+++ b/R/ordConstrained.R
@@ -97,19 +97,25 @@
 
 ### COMMON HEADER INFORMATION FOR ORDINATION MODELS
 
-`ordHead`<- function(Y, method)
+`ordHead`<- function(Y)
 {
-    if (method == "dbrda")
+    method <- attr(Y, "METHOD")
+    headmethod <- switch(method,
+                         "CA" = "cca",
+                         "PCA" = "rda",
+                         "CAPSCALE" = "capscale",
+                         "DISTBASED" = "dbrda")
+    if (method == "DISTBASED")
         totvar <- sum(diag(Y))
     else
         totvar <- sum(Y^2)
-    head <- list("tot.chi" = totvar, "Ybar" = Y)
-    if (method == "cca")
+    head <- list("tot.chi" = totvar, "Ybar" = Y, "method" = headmethod)
+    if (method == "CA")
         head <- c(list("grand.total" = attr(Y, "tot"),
                        "rowsum" = attr(Y, "RW"),
                        "colsum" = attr(Y, "CW")),
                   head)
-    else if (method == "rda")
+    else if (method == "PCA")
         head <- c(list("colsum" = sqrt(colSums(Y^2))),
                   head)
     head
@@ -357,19 +363,20 @@
 
 `ordConstrained` <-
     function(Y, X = NULL, Z = NULL,
-             method = c("cca", "rda", "capscale", "dbrda"),
+             method = c("cca", "rda", "capscale", "dbrda", "pass"),
              arg = FALSE)
 {
     method = match.arg(method)
     partial <- constraint <- resid <- NULL
-    ## init
+    ## init; "pass" returns unchanged Y, presumably from previous init
     Y <- switch(method,
                 "cca" = initCA(Y),
                 "rda" = initPCA(Y, scale = arg),
                 "capscale" = initCAP(Y),
-                "dbrda" = initDBRDA(Y))
+                "dbrda" = initDBRDA(Y),
+                "pass" = Y)
     ## header info for the model
-    head <- ordHead(Y, method)
+    head <- ordHead(Y)
     ## Partial
     if (!is.null(Z)) {
         out <- ordPartial(Y, Z)
@@ -388,5 +395,10 @@
     out <- c(head,
              call = match.call(),
              list("pCCA" = partial, "CCA" = constraint, "CA" = resid))
+    class(out) <- switch(attr(Y, "METHOD"),
+                         "CA" = "cca",
+                         "PCA" = c("rda", "cca"),
+                         "CAPSCALE" = c("capscale", "rda", "cca"),
+                         "DISTBASED" = c("dbrda", "rda", "cca"))
     out
 }
diff --git a/R/permutest.cca.R b/R/permutest.cca.R
index 5603ae75..1d676b03 100644
--- a/R/permutest.cca.R
+++ b/R/permutest.cca.R
@@ -32,7 +32,7 @@ permutest.default <- function(x, ...)
     ## special cases
     isCCA <- !inherits(x, "rda")    # weighting
     isPartial <- !is.null(x$pCCA)   # handle conditions
-    isDB <- inherits(x, c("dbrda"))
+    isDB <- inherits(x, c("dbrda")) # only dbrda is distance-based
     ## C function to get the statististics in one loop
     getF <- function(indx, ...)
     {
diff --git a/R/rda.default.R b/R/rda.default.R
index f1adf60e..b7d061a3 100644
--- a/R/rda.default.R
+++ b/R/rda.default.R
@@ -15,7 +15,7 @@
     sol$call <- call
     inertia <- if (scale) "correlations" else "variance"
     sol <- c(sol,
-             list(method = "rda", "inertia" = inertia))
+             list("inertia" = inertia))
     class(sol) <- c("rda", "cca")
     sol
 }
diff --git a/man/anova.cca.Rd b/man/anova.cca.Rd
index 7500bbec..004c1b45 100644
--- a/man/anova.cca.Rd
+++ b/man/anova.cca.Rd
@@ -32,10 +32,10 @@
 \arguments{
 
   \item{object}{One or several result objects from \code{\link{cca}},
-    \code{\link{rda}} or \code{\link{capscale}}. If there are several
-    result objects, they are compared against each other in the order
-    they were supplied. For a single object, a test specified in
-    \code{by} or an overall test is given.}
+    \code{\link{rda}}, \code{\link{dbrda}} or \code{\link{capscale}}. If
+    there are several result objects, they are compared against each
+    other in the order they were supplied. For a single object, a test
+    specified in \code{by} or an overall test is given.}
 
   \item{x}{A single ordination result object.}
  
@@ -71,7 +71,7 @@
     \code{permutations = how(\dots, blocks)} instead. }
 
   \item{cutoff}{Only effective with \code{by="axis"} where stops
-    permutations after an axis exceeds the \code{cutoff}.}
+    permutations after an axis exceeds the \code{cutoff} \eqn{p}-value.}
 
   \item{scope}{Only effective with \code{by="margin"} where it can be
     used to select the marginal terms for testing. The default is to
@@ -88,13 +88,12 @@
 
 \details{
 
-  Functions \code{anova.cca} and \code{permutest.cca} implement 
-  ANOVA like permutation tests for the joint effect of constraints in
-  \code{\link{cca}}, \code{\link{rda}} or \code{\link{capscale}}.
-  Functions \code{anova.cca} and \code{permutest.cca} differ in
-  printout style and in interface.  Function \code{permutest.cca} is
-  the proper workhorse, but \code{anova.cca} passes all parameters to
-  \code{permutest.cca}.
+  Functions \code{anova.cca} and \code{permutest.cca} implement ANOVA
+  like permutation tests for the joint effect of constraints in
+  \code{\link{cca}}, \code{\link{rda}}, \code{\link{dbrda}} or
+  \code{\link{capscale}}. Function \code{anova} is intended as a more
+  user-friendly alternative to \code{permutest} (that is the real
+  workhorse). 
 
   Function \code{anova} can analyse a sequence of constrained
   ordination models. The analysis is based on the differences in
@@ -112,12 +111,12 @@
   out}) and a test for the first constrained eigenvalues is
   performed (Legendre et al. 2011). 
   You can stop permutation tests after exceeding a given
-   significance level with argument \code{cutoff} to speed up
+  significance level with argument \code{cutoff} to speed up
   calculations in large models. Setting \code{by = "terms"} will
   perform separate significance test for each term (constraining
   variable). The terms are assessed sequentially from first to last,
   and the order of the terms will influence their
-  significance. Setting \code{by = "margin"} will perform separate
+  significances. Setting \code{by = "margin"} will perform separate
   significance test for each marginal term in a model with all other
   terms. The marginal test also accepts a \code{scope} argument for
   the \code{\link{drop.scope}} which can be a character vector of term
@@ -128,24 +127,23 @@
   terms. In calculating pseudo-\eqn{F}, all terms are compared to the
   same residual of the full model. 
   
-  Community data are permuted with choice \code{model="direct"},
-  and residuals after partial CCA/ RDA/ dbRDA with choice
-  \code{model="reduced"} (default).   If there is no partial CCA/
-  RDA/ dbRDA stage, \code{model="reduced"} simply permutes the data
-  and is equivalent to \code{model="direct"}.  The test statistic is
+  Community data are permuted with choice \code{model="direct"}, and
+  residuals after partial CCA/ RDA/ dbRDA with choice
+  \code{model="reduced"} (default).  If there is no partial CCA/ RDA/
+  dbRDA stage, \code{model="reduced"} simply permutes the data and is
+  equivalent to \code{model="direct"}.  The test statistic is
   \dQuote{pseudo-\eqn{F}}, which is the ratio of constrained and
   unconstrained total Inertia (Chi-squares, variances or something
-  similar), each divided by their respective ranks.  If there are no
-  conditions (\dQuote{partial} terms), the sum of all eigenvalues
-  remains constant, so that pseudo-\eqn{F} and eigenvalues would give
-  equal results.  In partial CCA/ RDA/ dbRDA, the effect of
+  similar), each divided by their respective degrees of freedom.  If
+  there are no conditions (\dQuote{partial} terms), the sum of all
+  eigenvalues remains constant, so that pseudo-\eqn{F} and eigenvalues
+  would give equal results.  In partial CCA/ RDA/ dbRDA, the effect of
   conditioning variables (\dQuote{covariables}) is removed before
-  permutation, and these residuals are added to the non-permuted
-  fitted values of partial CCA (fitted values of \code{X ~ Z}).
-  Consequently, the total Chi-square is not fixed, and test based on
+  permutation, and the total Chi-square is not fixed, and test based on
   pseudo-\eqn{F} would differ from the test based on plain
-  eigenvalues. CCA is a weighted method, and environmental data are
-  re-weighted at each permutation step using permuted weights.  }
+  eigenvalues.
+
+}
 
 \value{ 
   The function \code{anova.cca} calls \code{permutest.cca} and fills an
@@ -155,11 +153,6 @@
   \code{F.perm} (the permuted test statistics).
 }
 
-\note{
-  Some cases of \code{anova} need access to the original data on
-  constraints (at least \code{by = "term"} and \code{by = "margin"}),
-  and they may fail if data are unavailable.
-}
 \references{
   Legendre, P. and Legendre, L. (2012). \emph{Numerical Ecology}. 3rd
   English ed. Elsevier.
@@ -178,11 +171,15 @@
   semiautomatic model building (see \code{\link{deviance.cca}}). }
 
 \examples{
-data(varespec)
-data(varechem)
-vare.cca <- cca(varespec ~ Al + P + K, varechem)
+data(varespec, varechem)
+mod <- cca(varespec ~ Al + P + K, varechem)
 ## overall test
-anova(vare.cca)
+anova(mod)
+## tests for individual terms
+anova(mod, by="term")
+anova(mod, by="margin")
+## test for adding all environmental variables
+anova(mod, cca(varespec ~ ., varechem))
 }
 \keyword{ multivariate }
 \keyword{ htest }
diff --git a/man/model.matrix.cca.Rd b/man/model.matrix.cca.Rd
index edb5040c..b5e1e79a 100644
--- a/man/model.matrix.cca.Rd
+++ b/man/model.matrix.cca.Rd
@@ -1,6 +1,7 @@
 \name{model.matrix.cca}
 \Rdversion{1.1}
 \alias{model.matrix.cca}
+\alias{model.matrix.rda}
 \alias{model.frame.cca}
 %- Also NEED an '\alias' for EACH other topic documented here.
 \title{
@@ -33,19 +34,25 @@
 \details{
   The constrained ordination method objects do not save data on model 
   frame or design matrix, and the functions must reconstruct the
-  information in the session. This will fail if the data sets and variables
-  of the original model are unavailable.
+  information from the result object.
 }
 \value{
   Returns a data frame (\code{model.frame}) or an unnamed matrix or a 
   list of two matrices called \code{Constraints} and \code{Conditions} 
   (\code{model.matrix}). 
 }
+
+\note{The \code{model.matrix} returns the unweighted model matrix also
+  for \code{\link{cca}}. Prior to \pkg{vegan} version 2.5-0 it returned
+  the weighted model matrix.}
+
 \author{
  Jari Oksanen
 }
 \seealso{
- \code{\link{model.frame}}, \code{\link{model.matrix}}.
+ \code{\link{model.frame}}, \code{\link{model.matrix}}. The model matrix
+ is built from the QR decomposition component of \code{\link{cca.object}}
+ using \code{\link{qr.X}}.
 }
 \examples{
 data(dune)
diff --git a/tests/cca-object-tests.Rout.save b/tests/cca-object-tests.Rout.save
index b6877865..f515f512 100644
--- a/tests/cca-object-tests.Rout.save
+++ b/tests/cca-object-tests.Rout.save
@@ -1,5 +1,5 @@
 
-R Under development (unstable) (2017-05-02 r72640) -- "Unsuffered Consequences"
+R Under development (unstable) (2017-06-05 r72769) -- "Unsuffered Consequences"
 Copyright (C) 2017 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu (64-bit)
 
@@ -2059,100 +2059,96 @@ Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 > 
 > anova(mcca, permutations = per, by="axis")
 Permutation test for cca under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 49
 
 Model: cca(formula = dune ~ Condition(Management) + Manure + A1, data = dune.env)
-         Df ChiSquare      F Pr(>F)  
-CCA1      1   0.22257 2.5154   0.02 *
-CCA2      1   0.12260 1.3855   0.32  
-CCA3      1   0.06390 0.7222   0.92  
-CCA4      1   0.04055 0.4583   0.98  
-Residual 12   1.06181                
----
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+         Df ChiSquare      F Pr(>F)
+CCA1      1   0.22257 2.5154   0.14
+CCA2      1   0.12260 1.3855   0.70
+CCA3      1   0.06390 0.7222   1.00
+CCA4      1   0.04055 0.4583   0.90
+Residual 12   1.06181              
 > anova(mrda, permutations = per, by="axis")
 Permutation test for rda under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 49
 
 Model: rda(formula = dune ~ Condition(Management) + Manure + A1, data = dune.env)
-         Df Variance      F Pr(>F)
-RDA1      1    8.539 2.6895   0.12
-RDA2      1    4.330 1.3636   0.34
-RDA3      1    2.206 0.6948   0.92
-RDA4      1    1.716 0.5405   1.00
-Residual 12   38.102              
+         Df Variance      F Pr(>F)  
+RDA1      1    8.539 2.6895   0.08 .
+RDA2      1    4.330 1.3636   0.72  
+RDA3      1    2.206 0.6948   0.96  
+RDA4      1    1.716 0.5405   0.86  
+Residual 12   38.102                
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 > anova(mrda1, permutations = per, by="axis")
 Permutation test for rda under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 49
 
 Model: rda(formula = dune ~ Condition(Management) + Manure + A1, data = dune.env, scale = TRUE)
          Df Variance      F Pr(>F)
-RDA1      1   2.6117 2.0784   0.14
-RDA2      1   1.6991 1.3521   0.34
-RDA3      1   1.1039 0.8785   0.80
-RDA4      1   0.7261 0.5778   0.96
+RDA1      1   2.6117 2.0784   0.18
+RDA2      1   1.6991 1.3521   0.76
+RDA3      1   1.1039 0.8785   0.96
+RDA4      1   0.7261 0.5778   0.88
 Residual 12  15.0789              
 > anova(mcap, permutations = per, by="axis")
 Permutation test for capscale under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 49
 
 Model: capscale(formula = dune ~ Condition(Management) + Manure + A1, data = dune.env, distance = "bray")
-         Df SumOfSqs      F Pr(>F)  
-CAP1      1  0.67370 3.9998   0.02 *
-CAP2      1  0.24300 1.4427   0.34  
-CAP3      1  0.10457 0.6209   0.96  
-CAP4      1  0.05136 0.3049   1.00  
-Residual 12  2.02118                
----
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+         Df SumOfSqs      F Pr(>F)
+CAP1      1  0.67370 3.9998   0.12
+CAP2      1  0.24300 1.4427   0.78
+CAP3      1  0.10457 0.6209   0.94
+CAP4      1  0.05136 0.3049   0.98
+Residual 12  2.02118              
 > anova(mdb, permutations = per, by="axis")
 Permutation test for dbrda under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 49
 
 Model: dbrda(formula = dune ~ Condition(Management) + Manure + A1, data = dune.env, distance = "bray")
-         Df SumOfSqs      F Pr(>F)  
-dbRDA1    1  0.66660 4.4542   0.02 *
-dbRDA2    1  0.23850 1.5936   0.34  
-dbRDA3    1  0.09102 0.6082   0.98  
-dbRDA4    1  0.03844 0.2569   0.98  
-Residual 12  1.79587                
----
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+         Df SumOfSqs      F Pr(>F)
+dbRDA1    1  0.66660 4.4542   0.10
+dbRDA2    1  0.23850 1.5936   0.76
+dbRDA3    1  0.09102 0.6082   0.94
+dbRDA4    1  0.03844 0.2569   0.96
+Residual 12  1.79587              
 > anova(mancap, permutations = per, by="axis")
 Permutation test for capscale under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 49
 
 Model: capscale(formula = dune ~ Condition(Management) + Manure + A1, data = dune.env, distance = "manhattan")
          Df Variance      F Pr(>F)
 CAP1      1   125.66 3.3341   0.14
-CAP2      1    49.31 1.3083   0.46
-CAP3      1    24.24 0.6431   0.90
-CAP4      1    13.87 0.3681   1.00
+CAP2      1    49.31 1.3083   0.80
+CAP3      1    24.24 0.6431   0.98
+CAP4      1    13.87 0.3681   0.92
 Residual 12   452.26              
 > anova(mandb, permutations = per, by="axis")
 Permutation test for dbrda under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 49
 
 Model: dbrda(formula = dune ~ Condition(Management) + Manure + A1, data = dune.env, distance = "manhattan")
          Df Variance      F Pr(>F)
 dbRDA1    1   123.16 3.8146   0.12
-dbRDA2    1    47.37 1.4671   0.42
-dbRDA3    1    21.91 0.6787   0.86
-dbRDA4    1    10.13 0.3137   0.96
+dbRDA2    1    47.37 1.4671   0.74
+dbRDA3    1    21.91 0.6787   0.96
+dbRDA4    1    10.13 0.3137   0.88
 Residual 12   387.45              
 > 
 > ## the following do not all work with partial models
@@ -2543,4 +2539,4 @@ Alopgeni -0.197400341  0.18131362
 > 
 > proc.time()
    user  system elapsed 
-  3.656   0.052   3.706 
+  3.092   0.072   3.163 
diff --git a/tests/vegan-tests.R b/tests/vegan-tests.R
index b83f18d7..f5eb5075 100644
--- a/tests/vegan-tests.R
+++ b/tests/vegan-tests.R
@@ -38,24 +38,21 @@ spno <- specnumber(dune)
 ## cca/rda
 m <-  cca(fla, data=df,  na.action=na.exclude,  subset = Use != "Pasture" & spno > 7)
 anova(m, permutations=99)
-## vegan 2.1-40 cannot handle missing data in next two
-##anova(m, by="term", permutations=99)
-##anova(m, by="margin", permutations=99)
+anova(m, by="term", permutations=99) # failed before 2.5-0
+anova(m, by="margin", permutations=99) # works since 2.5-0
 anova(m, by="axis", permutations=99)
 ## capscale
 p <- capscale(fla, data=df, na.action=na.exclude, subset = Use != "Pasture" & spno > 7)
 anova(p, permutations=99)
-## vegan 2.1-40 cannot handle missing data in next two
-##anova(p, by="term", permutations=99)
-##anova(p, by="margin", permutations=99)
+anova(p, by="term", permutations=99) # failed before 2.5-0
+anova(p, by="margin", permutations=99) # works since 2.5-0
 anova(p, by="axis", permutations=99)
 ## see that capscale can be updated and also works with 'dist' input
 dis <- vegdist(dune)
 p <- update(p, dis ~ .)
 anova(p, permutations=99)
-## vegan 2.1-40 cannot handle missing data in next two
-##anova(p, by="term", permutations=99)
-##anova(p, by="margin", permutations=99)
+anova(p, by="term", permutations=99) # failed before 2.5-0
+anova(p, by="margin", permutations=99) # works since 2.5-0
 anova(p, by="axis", permutations=99)
 ### attach()ed data frame instead of data=
 attach(df)
@@ -63,8 +60,8 @@ q <- cca(fla, na.action = na.omit, subset = Use != "Pasture" & spno > 7)
 anova(q, permutations=99)
 ## commented tests below fail in vegan 2.1-40 because number of
 ## observations changes
-##anova(q, by="term", permutations=99) 
-##anova(q, by="margin", permutations=99)
+anova(q, by="term", permutations=99) # failed before 2.5-0
+anova(q, by="margin", permutations=99) # works since 2.5-0
 anova(q, by="axis", permutations=99)
 ### Check that constrained ordination functions can be embedded.
 ### The data.frame 'df' is still attach()ed.
diff --git a/tests/vegan-tests.Rout.save b/tests/vegan-tests.Rout.save
index b0373ab0..aa81ed22 100644
--- a/tests/vegan-tests.Rout.save
+++ b/tests/vegan-tests.Rout.save
@@ -1,6 +1,6 @@
 
-R Under development (unstable) (2016-11-22 r71678) -- "Unsuffered Consequences"
-Copyright (C) 2016 The R Foundation for Statistical Computing
+R Under development (unstable) (2017-06-01 r72753) -- "Unsuffered Consequences"
+Copyright (C) 2017 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu (64-bit)
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -65,24 +65,47 @@ Model     6   1.25838 1.3106   0.09 .
 Residual  5   0.80011                
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-> ## vegan 2.1-40 cannot handle missing data in next two
-> ##anova(m, by="term", permutations=99)
-> ##anova(m, by="margin", permutations=99)
+> anova(m, by="term", permutations=99) # failed before 2.5-0
+Permutation test for cca under reduced model
+Terms added sequentially (first to last)
+Permutation: free
+Number of permutations: 99
+
+Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+            Df ChiSquare      F Pr(>F)  
+Management   3   0.82392 1.7163   0.03 *
+poly(A1, 2)  2   0.35127 1.0976   0.39  
+spno         1   0.08318 0.5198   0.96  
+Residual     5   0.80011                
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+> anova(m, by="margin", permutations=99) # works since 2.5-0
+Permutation test for cca under reduced model
+Marginal effects of terms
+Permutation: free
+Number of permutations: 99
+
+Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+            Df ChiSquare      F Pr(>F)
+Management   3   0.55418 1.1544   0.21
+poly(A1, 2)  2   0.32940 1.0292   0.44
+spno         1   0.08318 0.5198   0.93
+Residual     5   0.80011              
 > anova(m, by="axis", permutations=99)
 Permutation test for cca under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 99
 
-Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = object$subset)
-         Df ChiSquare      F Pr(>F)   
-CCA1      1   0.46993 2.9366   0.01 **
-CCA2      1   0.26217 1.6384   0.18   
-CCA3      1   0.19308 1.2066   0.29   
-CCA4      1   0.18345 1.1464   0.38   
-CCA5      1   0.08871 0.5544   0.76   
-CCA6      1   0.06104 0.3815   0.88   
-Residual  5   0.80011                 
+Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+         Df ChiSquare      F Pr(>F)  
+CCA1      1   0.46993 2.9366   0.07 .
+CCA2      1   0.26217 1.6384   0.75  
+CCA3      1   0.19308 1.2066   0.97  
+CCA4      1   0.18345 1.1464   0.90  
+CCA5      1   0.08871 0.5544   0.98  
+CCA6      1   0.06104 0.3815   0.93  
+Residual  5   0.80011                
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 > ## capscale
@@ -98,23 +121,46 @@ Model     6   59.582 1.6462   0.04 *
 Residual  5   30.160                
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-> ## vegan 2.1-40 cannot handle missing data in next two
-> ##anova(p, by="term", permutations=99)
-> ##anova(p, by="margin", permutations=99)
+> anova(p, by="term", permutations=99) # failed before 2.5-0
+Permutation test for capscale under reduced model
+Terms added sequentially (first to last)
+Permutation: free
+Number of permutations: 99
+
+Model: capscale(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+            Df Variance      F Pr(>F)   
+Management   3   46.265 2.5566   0.01 **
+poly(A1, 2)  2   10.150 0.8413   0.63   
+spno         1    3.167 0.5250   0.79   
+Residual     5   30.160                 
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+> anova(p, by="margin", permutations=99) # works since 2.5-0
+Permutation test for capscale under reduced model
+Marginal effects of terms
+Permutation: free
+Number of permutations: 99
+
+Model: capscale(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+            Df Variance      F Pr(>F)
+Management   3  28.7515 1.5888   0.13
+poly(A1, 2)  2   8.0847 0.6701   0.88
+spno         1   3.1669 0.5250   0.90
+Residual     5  30.1605              
 > anova(p, by="axis", permutations=99)
 Permutation test for capscale under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 99
 
-Model: capscale(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = object$subset)
+Model: capscale(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
          Df Variance      F Pr(>F)  
 CAP1      1  25.0252 4.1487   0.03 *
-CAP2      1  15.8759 2.6319   0.07 .
-CAP3      1   8.0942 1.3419   0.25  
-CAP4      1   5.0675 0.8401   0.64  
-CAP5      1   3.5671 0.5914   0.85  
-CAP6      1   1.9520 0.3236   0.96  
+CAP2      1  15.8759 2.6319   0.39  
+CAP3      1   8.0942 1.3419   0.78  
+CAP4      1   5.0675 0.8401   0.95  
+CAP5      1   3.5671 0.5914   0.96  
+CAP6      1   1.9520 0.3236   0.92  
 Residual  5  30.1605                
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
@@ -128,28 +174,49 @@ Number of permutations: 99
 
 Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
          Df SumOfSqs      F Pr(>F)
-Model     6  1.54840 1.6423   0.11
+Model     6  1.54840 1.6423   0.14
 Residual  5  0.78568              
-> ## vegan 2.1-40 cannot handle missing data in next two
-> ##anova(p, by="term", permutations=99)
-> ##anova(p, by="margin", permutations=99)
-> anova(p, by="axis", permutations=99)
+> anova(p, by="term", permutations=99) # failed before 2.5-0
 Permutation test for capscale under reduced model
-Marginal tests for axes
+Terms added sequentially (first to last)
 Permutation: free
 Number of permutations: 99
 
-Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = object$subset)
-         Df SumOfSqs      F Pr(>F)  
-CAP1      1  0.77834 4.9533   0.02 *
-CAP2      1  0.45691 2.9078   0.03 *
-CAP3      1  0.14701 0.9355   0.51  
-CAP4      1  0.11879 0.7560   0.65  
-CAP5      1  0.04213 0.2681   0.94  
-CAP6      1  0.00522 0.0332   1.00  
-Residual  5  0.78568                
+Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+            Df SumOfSqs      F Pr(>F)  
+Management   3  1.17117 2.4844   0.02 *
+poly(A1, 2)  2  0.30602 0.9737   0.53  
+spno         1  0.07121 0.4532   0.84  
+Residual     5  0.78568                
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+> anova(p, by="margin", permutations=99) # works since 2.5-0
+Permutation test for capscale under reduced model
+Marginal effects of terms
+Permutation: free
+Number of permutations: 99
+
+Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+            Df SumOfSqs      F Pr(>F)
+Management   3  0.64888 1.3765   0.22
+poly(A1, 2)  2  0.26160 0.8324   0.69
+spno         1  0.07121 0.4532   0.91
+Residual     5  0.78568              
+> anova(p, by="axis", permutations=99)
+Permutation test for capscale under reduced model
+Forward tests for axes
+Permutation: free
+Number of permutations: 99
+
+Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7)
+         Df SumOfSqs      F Pr(>F)
+CAP1      1  0.77834 4.9533   0.13
+CAP2      1  0.45691 2.9078   0.48
+CAP3      1  0.14701 0.9355   0.98
+CAP4      1  0.11879 0.7560   0.96
+CAP5      1  0.04213 0.2681   1.00
+CAP6      1  0.00522 0.0332   1.00
+Residual  5  0.78568              
 > ### attach()ed data frame instead of data=
 > attach(df)
 > q <- cca(fla, na.action = na.omit, subset = Use != "Pasture" & spno > 7)
@@ -160,26 +227,50 @@ Number of permutations: 99
 
 Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = Use != "Pasture" & spno > 7)
          Df ChiSquare      F Pr(>F)
-Model     6   1.25838 1.3106   0.16
+Model     6   1.25838 1.3106   0.11
 Residual  5   0.80011              
 > ## commented tests below fail in vegan 2.1-40 because number of
 > ## observations changes
-> ##anova(q, by="term", permutations=99) 
-> ##anova(q, by="margin", permutations=99)
+> anova(q, by="term", permutations=99) # failed before 2.5-0
+Permutation test for cca under reduced model
+Terms added sequentially (first to last)
+Permutation: free
+Number of permutations: 99
+
+Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = Use != "Pasture" & spno > 7)
+            Df ChiSquare      F Pr(>F)  
+Management   3   0.82392 1.7163   0.02 *
+poly(A1, 2)  2   0.35127 1.0976   0.40  
+spno         1   0.08318 0.5198   0.95  
+Residual     5   0.80011                
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+> anova(q, by="margin", permutations=99) # works since 2.5-0
+Permutation test for cca under reduced model
+Marginal effects of terms
+Permutation: free
+Number of permutations: 99
+
+Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = Use != "Pasture" & spno > 7)
+            Df ChiSquare      F Pr(>F)
+Management   3   0.55418 1.1544   0.33
+poly(A1, 2)  2   0.32940 1.0292   0.44
+spno         1   0.08318 0.5198   0.94
+Residual     5   0.80011              
 > anova(q, by="axis", permutations=99)
 Permutation test for cca under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 99
 
-Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = object$subset)
+Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = Use != "Pasture" & spno > 7)
          Df ChiSquare      F Pr(>F)  
-CCA1      1   0.46993 2.9366   0.03 *
-CCA2      1   0.26217 1.6384   0.22  
-CCA3      1   0.19308 1.2066   0.31  
-CCA4      1   0.18345 1.1464   0.32  
-CCA5      1   0.08871 0.5544   0.70  
-CCA6      1   0.06104 0.3815   0.88  
+CCA1      1   0.46993 2.9366   0.05 *
+CCA2      1   0.26217 1.6384   0.82  
+CCA3      1   0.19308 1.2066   0.93  
+CCA4      1   0.18345 1.1464   0.88  
+CCA5      1   0.08871 0.5544   0.97  
+CCA6      1   0.06104 0.3815   0.89  
 Residual  5   0.80011                
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
@@ -317,16 +408,16 @@ Eigenvalues for unconstrained axes:
 
 > tab
 Permutation test for cca under reduced model
-Marginal tests for axes
+Forward tests for axes
 Permutation: free
 Number of permutations: 99
 
-Model: cca(formula = dune ~ A1 + Moisture + Condition(Management), data = dune.env, subset = object$subset)
+Model: cca(formula = dune ~ A1 + Moisture + Condition(Management), data = dune.env, subset = A1 > 3)
          Df ChiSquare      F Pr(>F)  
-CCA1      1   0.27109 2.9561   0.04 *
-CCA2      1   0.14057 1.5329   0.22  
-CCA3      1   0.08761 0.9553   0.71  
-CCA4      1   0.05624 0.6132   0.98  
+CCA1      1   0.27109 2.9561   0.09 .
+CCA2      1   0.14057 1.5329   0.56  
+CCA3      1   0.08761 0.9553   0.78  
+CCA4      1   0.05624 0.6132   0.76  
 Residual 10   0.91705                
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
@@ -355,16 +446,16 @@ Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 > cap.model.cond <- capscale(X ~ A + B + Condition(C))
 > anova(cap.model.cond, by="axis", strata=C)  # -> error pre r2287
 Permutation test for capscale under reduced model
-Marginal tests for axes
+Forward tests for axes
 Blocks:  strata 
 Permutation: free
 Number of permutations: 999
 
 Model: capscale(formula = X ~ A + B + Condition(C))
          Df Variance      F Pr(>F)
-CAP1      1   0.2682 1.3075  0.242
-CAP2      1   0.0685 0.3339  0.921
-CAP3      1   0.0455 0.2217  0.966
+CAP1      1   0.2682 1.3075  0.783
+CAP2      1   0.0685 0.3339  0.998
+CAP3      1   0.0455 0.2217  0.967
 Residual 22   4.5130              
 > anova(cap.model.cond, by="terms", strata=C)  # -> error pre r2287
 Permutation test for capscale under reduced model
@@ -383,15 +474,15 @@ Residual 22   4.5130
 > cap.model <- capscale(X ~ A + B)
 > anova(cap.model, by="axis", strata=C)  # -> no error
 Permutation test for capscale under reduced model
-Marginal tests for axes
+Forward tests for axes
 Blocks:  strata 
 Permutation: free
 Number of permutations: 999
 
 Model: capscale(formula = X ~ A + B)
          Df Variance      F Pr(>F)
-CAP1      1   0.2682 1.3267  0.240
-CAP2      1   0.0685 0.3388  0.913
+CAP1      1   0.2682 1.3267  0.764
+CAP2      1   0.0685 0.3388  0.993
 CAP3      1   0.0455 0.2249  0.964
 Residual 26   5.2565              
 > anova(cap.model, by="terms", strata=C)  # -> no error
@@ -411,16 +502,16 @@ Residual 26   5.2565
 > rda.model.cond <- rda(X ~ A + B + Condition(C))
 > anova(rda.model.cond, by="axis", strata=C)  # -> no error
 Permutation test for rda under reduced model
-Marginal tests for axes
+Forward tests for axes
 Blocks:  strata 
 Permutation: free
 Number of permutations: 999
 
 Model: rda(formula = X ~ A + B + Condition(C))
          Df Variance      F Pr(>F)
-RDA1      1   0.2682 1.3075  0.286
-RDA2      1   0.0685 0.3339  0.921
-RDA3      1   0.0455 0.2217  0.963
+RDA1      1   0.2682 1.3075  0.774
+RDA2      1   0.0685 0.3339  0.993
+RDA3      1   0.0455 0.2217  0.964
 Residual 22   4.5130              
 > anova(rda.model.cond, by="terms", strata=C)  # -> error pre r2287
 Permutation test for rda under reduced model
@@ -439,15 +530,15 @@ Residual 22   4.5130
 > rda.model <- rda(X ~ A + B)
 > anova(rda.model, by="axis", strata=C)  # -> no error
 Permutation test for rda under reduced model
-Marginal tests for axes
+Forward tests for axes
 Blocks:  strata 
 Permutation: free
 Number of permutations: 999
 
 Model: rda(formula = X ~ A + B)
          Df Variance      F Pr(>F)
-RDA1      1   0.2682 1.3267  0.258
-RDA2      1   0.0685 0.3388  0.898
+RDA1      1   0.2682 1.3267  0.768
+RDA2      1   0.0685 0.3388  0.995
 RDA3      1   0.0455 0.2249  0.971
 Residual 26   5.2565              
 > anova(rda.model, by="terms", strata=C)  # -> no error
@@ -907,4 +998,4 @@ Cumulative Proportion 0.999490 1.0000000
 > 
 > proc.time()
    user  system elapsed 
-  6.044   0.056   6.096 
+  4.132   0.044   4.171