From a14f1a166db1f15e454712d9bb52ed895533bd0c Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 11 May 2017 10:58:11 +0300 Subject: [PATCH 01/37] anova.cca(..., by="term") uses directly permutest.cca we used to manipulate formula within anova.ccabyterm, refit models and call anova.ccalist --- R/anova.ccabyterm.R | 36 ++++++++++++------------------------ 1 file changed, 12 insertions(+), 24 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 7355327b..425b5eb3 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. +### 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 `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", From 7f3b961823b81de3d62c19b299ce0dae3a80d0e2 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 11 May 2017 11:04:35 +0300 Subject: [PATCH 02/37] fix F.perm attribute --- R/anova.ccabyterm.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 425b5eb3..803792e3 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -33,7 +33,7 @@ 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 } From 03a03d78efb531c69effb4dd0001a1e8e8679137 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 12 May 2017 10:53:05 +0300 Subject: [PATCH 03/37] add sequential test for axes based on 1-df constrasts in permutest This is similar to the forward selection test in Legendre, Oksanen, ter Braak, MEE 2, 269-277 (2011) with following differences: - the F statistic uses always the same denominator from the residual unconstrained model like all vegan anova.cca tests. However, the numerator should be the same. - The previous axes are permuted similarly as tested axes instead of being permuted like Conditions. However, Conditions prior to constrained analysis are still permuted like defined by the model argument. The method needs testing and verification before incorporation into master. --- R/anova.cca.R | 7 +++++-- R/anova.ccabyterm.R | 40 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 2 deletions(-) diff --git a/R/anova.cca.R b/R/anova.cca.R index d78b838b..695e21f8 100644 --- a/R/anova.cca.R +++ b/R/anova.cca.R @@ -37,7 +37,7 @@ return(anova.ccanull(object)) ## by cases if (!is.null(by)) { - by <- match.arg(by, c("terms", "margin", "axis")) + by <- match.arg(by, c("terms", "margin", "axis", "seqaxis")) if (is.null(object$terms)) stop("model must be fitted with formula interface") sol <- switch(by, @@ -51,7 +51,10 @@ "axis" = anova.ccabyaxis(object, permutations = permutations, model = model, parallel = parallel, - cutoff = cutoff)) + cutoff = cutoff), + "seqaxis" = anova.ccabysaxis(object, + permutations = permutations, + model = model, parallel=parallel)) attr(sol, "Random.seed") <- seed attr(sol, "control") <- control return(sol) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 803792e3..38ceb510 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -189,3 +189,43 @@ class(out) <- c("anova.cca", "anova", "data.frame") out } + +### Sequential test for axes is based on permustest by="onedf": we +### replace the original QR decomposition with the QR decomposition of +### linear component axes (u) and asd for 1-df contrasts. This is +### essentially similar as partialling out previous constrained axes +### except that it implies permutation model "direct", i.e., previous +### axes (Conditions) and residuals use the same permutation. This differs from + +anova.ccabysaxis <- + function(object, permutations, model, parallel) +{ + ## replace original QR decomposition of constraints with QR + ## decomposition of LC of constraints. + object$CCA$QR <- qr(object$CCA$u) + sol <- permutest(object, permutations = permutations, + model = model, by = "onedf", parallel = parallel) + ## Reformat + 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(sol$termlabels, "Residual"), + c("Df", varname, "F", "Pr(>F)")) + head <- paste0("Permutation test for ", object$method, " under ", + model, " model\n", + "Sequential test for axes\n", + howHead(attr(permutations, "control"))) + mod <- paste("Model:", c(object$call)) + attr(out, "heading") <- c(head, mod) + attr(out, "F.perm") <- sol$F.perm + class(out) <- c("anova.cca", "anova","data.frame") + out +} From 0a49fc10e11e1ec19de1f07bf036ef397e540029 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 12 May 2017 11:03:36 +0300 Subject: [PATCH 04/37] not yet implemented for partial models --- R/anova.ccabyterm.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 38ceb510..6787e4ce 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -202,6 +202,8 @@ anova.ccabysaxis <- { ## replace original QR decomposition of constraints with QR ## decomposition of LC of constraints. + if (!is.null(object$pCCA) && object$pCCA$rank) + stop("by = 'seqaxis' is not implemented yet for partial models") object$CCA$QR <- qr(object$CCA$u) sol <- permutest(object, permutations = permutations, model = model, by = "onedf", parallel = parallel) From c5c53143372ca958ae483b3c0c6b2df1475f55e8 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 12 May 2017 12:54:19 +0300 Subject: [PATCH 05/37] new anova.cca/permustest.cca by="terms" needs assign attribute of terms --- R/adonis2.R | 2 ++ 1 file changed, 2 insertions(+) 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) From 2e1881cae4190669dbc0db8f0533033254c5058b Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 22 May 2017 10:50:56 +0300 Subject: [PATCH 06/37] implement sequential test for axes in partial ordination --- R/anova.ccabyterm.R | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 6787e4ce..7a57838f 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -197,14 +197,15 @@ ### except that it implies permutation model "direct", i.e., previous ### axes (Conditions) and residuals use the same permutation. This differs from -anova.ccabysaxis <- +`anova.ccabysaxis` <- function(object, permutations, model, parallel) { ## replace original QR decomposition of constraints with QR ## decomposition of LC of constraints. - if (!is.null(object$pCCA) && object$pCCA$rank) - stop("by = 'seqaxis' is not implemented yet for partial models") - object$CCA$QR <- qr(object$CCA$u) + LC <- object$CCA$u + if (!is.null(object$pCCA) && object$pCCA$rank > 0) + LC <- cbind(qr.X(object$pCCA$QR), LC) + object$CCA$QR <- qr(LC) sol <- permutest(object, permutations = permutations, model = model, by = "onedf", parallel = parallel) ## Reformat From b47e2ccaf7d6dbcdc502177cbbb200a88908df43 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 22 May 2017 11:07:49 +0300 Subject: [PATCH 07/37] somewhat clumsy fix to bug of handling partial models in by = "axis" anova(, by = "axis") did not include partial terms in marginal tests. As the result, the P-values were usually too high (close to 1). --- R/anova.ccabyterm.R | 3 +++ 1 file changed, 3 insertions(+) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 7a57838f..426259c0 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -154,6 +154,9 @@ for (i in seq_along(eig)) { part <- paste("~ . +Condition(", paste(names(LC)[-i], collapse = "+"), ")") + ## handle partial models (clumsily) + if (!is.null(object$pCCA) && object$pCCA$rank > 0) + part <- paste(part, "Condition(qr.X(object$pCCA$QR))", sep="+") upfla <- update(fla, part) ## only one axis, and cannot partial out? if (length(eig) == 1) From fe1e6570a0ded61a89069018f2487657567576a4 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Tue, 23 May 2017 13:19:39 +0300 Subject: [PATCH 08/37] document by = "seqaxis" in anova.cca --- man/anova.cca.Rd | 62 ++++++++++++++++++++++++++++---------------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/man/anova.cca.Rd b/man/anova.cca.Rd index 7500bbec..e52b3682 100644 --- a/man/anova.cca.Rd +++ b/man/anova.cca.Rd @@ -44,14 +44,15 @@ number of permutations required, or a permutation matrix where each row gives the permuted indices.} - \item{by}{Setting \code{by = "axis"} will assess significance for - each constrained axis, and setting \code{by = "terms"} will assess - significance for each term (sequentially from first to last), and - setting \code{by = "margin"} will assess the marginal effects of - the terms (each marginal term analysed in a model with all other - variables). Function \code{permutest} accepts choices - \code{"terms"} for sequential test of terms, and \code{"onedf"} - for sequential test of one-degree-of-freedom constrasts.} + \item{by}{Setting \code{by = "terms"} will assess significance for + each term (sequentially from first to last), \code{by = "margin"} + will assess the marginal effects of the terms (each marginal term + analysed in a model with all other variables), \code{by = "axis"} + will assess marginal effects of axes, and \code{by = "seqaxis"} the + sequential effects of axes. Function \code{permutest} accepts + choices \code{"terms"} for sequential test of terms, and + \code{"onedf"} for sequential test of one-degree-of-freedom + constrasts.} \item{model}{Permutation model: \code{model="direct"} permutes community data, \code{model="reduced"} permutes residuals @@ -104,29 +105,28 @@ Setting \code{first = TRUE} will perform a test for the first constrained eigenvalue. Argument \code{first} can be set either in \code{anova.cca} or in \code{permutest.cca}. It is also possible to - perform significance tests for each axis or for each term - (constraining variable) using argument \code{by} in - \code{anova.cca}. Setting \code{by = "axis"} will perform separate - significance tests for each constrained axis. All previous - constrained axes will be used as conditions (\dQuote{partialled - 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 - 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 - 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 - labels that are analysed, or a fitted model of lower scope. The - marginal effects are also known as \dQuote{Type III} effects, but - the current function only evaluates marginal terms. It will, for - instance, ignore main effects that are included in interaction - terms. In calculating pseudo-\eqn{F}, all terms are compared to the - same residual of the full model. + perform significance tests for each term (constraining variable) or + for each axis using argument \code{by}. 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 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 labels that are analysed, or a + fitted model of lower scope. The marginal effects are also known as + \dQuote{Type III} effects, but the current function only evaluates + marginal terms. It will, for instance, ignore main effects that are + included in interaction terms. Setting \code{by = "axis"} will + perform marginal test for each constrained axis. All previous + constrained axes will be used as conditions (\dQuote{partialled 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 calculations in large models. Setting \code{by = "seqaxis"} + will perform sequential test for each constrained axis. In all cases, + the terms are compared to the same residual of the full model when + calculating the pseudo-\eqn{F}. Community data are permuted with choice \code{model="direct"}, and residuals after partial CCA/ RDA/ dbRDA with choice From c4562d72f2ac8ce6b4f6e9857c3c59028d97ba18 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 29 May 2017 16:12:11 +0300 Subject: [PATCH 09/37] Revert "document by = "seqaxis" in anova.cca" This reverts commit fe1e6570a0ded61a89069018f2487657567576a4. --- man/anova.cca.Rd | 62 ++++++++++++++++++++++++++++---------------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/man/anova.cca.Rd b/man/anova.cca.Rd index e52b3682..7500bbec 100644 --- a/man/anova.cca.Rd +++ b/man/anova.cca.Rd @@ -44,15 +44,14 @@ number of permutations required, or a permutation matrix where each row gives the permuted indices.} - \item{by}{Setting \code{by = "terms"} will assess significance for - each term (sequentially from first to last), \code{by = "margin"} - will assess the marginal effects of the terms (each marginal term - analysed in a model with all other variables), \code{by = "axis"} - will assess marginal effects of axes, and \code{by = "seqaxis"} the - sequential effects of axes. Function \code{permutest} accepts - choices \code{"terms"} for sequential test of terms, and - \code{"onedf"} for sequential test of one-degree-of-freedom - constrasts.} + \item{by}{Setting \code{by = "axis"} will assess significance for + each constrained axis, and setting \code{by = "terms"} will assess + significance for each term (sequentially from first to last), and + setting \code{by = "margin"} will assess the marginal effects of + the terms (each marginal term analysed in a model with all other + variables). Function \code{permutest} accepts choices + \code{"terms"} for sequential test of terms, and \code{"onedf"} + for sequential test of one-degree-of-freedom constrasts.} \item{model}{Permutation model: \code{model="direct"} permutes community data, \code{model="reduced"} permutes residuals @@ -105,28 +104,29 @@ Setting \code{first = TRUE} will perform a test for the first constrained eigenvalue. Argument \code{first} can be set either in \code{anova.cca} or in \code{permutest.cca}. It is also possible to - perform significance tests for each term (constraining variable) or - for each axis using argument \code{by}. 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 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 labels that are analysed, or a - fitted model of lower scope. The marginal effects are also known as - \dQuote{Type III} effects, but the current function only evaluates - marginal terms. It will, for instance, ignore main effects that are - included in interaction terms. Setting \code{by = "axis"} will - perform marginal test for each constrained axis. All previous - constrained axes will be used as conditions (\dQuote{partialled 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 calculations in large models. Setting \code{by = "seqaxis"} - will perform sequential test for each constrained axis. In all cases, - the terms are compared to the same residual of the full model when - calculating the pseudo-\eqn{F}. + perform significance tests for each axis or for each term + (constraining variable) using argument \code{by} in + \code{anova.cca}. Setting \code{by = "axis"} will perform separate + significance tests for each constrained axis. All previous + constrained axes will be used as conditions (\dQuote{partialled + 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 + 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 + 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 + labels that are analysed, or a fitted model of lower scope. The + marginal effects are also known as \dQuote{Type III} effects, but + the current function only evaluates marginal terms. It will, for + instance, ignore main effects that are included in interaction + 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 From 13e9e6cdda6e0c27f954aa128862543b479efeca Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 29 May 2017 16:14:20 +0300 Subject: [PATCH 10/37] Revert "implement sequential test for axes in partial ordination" This reverts commit 2e1881cae4190669dbc0db8f0533033254c5058b. --- R/anova.ccabyterm.R | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 426259c0..25669cc9 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -200,15 +200,14 @@ ### except that it implies permutation model "direct", i.e., previous ### axes (Conditions) and residuals use the same permutation. This differs from -`anova.ccabysaxis` <- +anova.ccabysaxis <- function(object, permutations, model, parallel) { ## replace original QR decomposition of constraints with QR ## decomposition of LC of constraints. - LC <- object$CCA$u - if (!is.null(object$pCCA) && object$pCCA$rank > 0) - LC <- cbind(qr.X(object$pCCA$QR), LC) - object$CCA$QR <- qr(LC) + if (!is.null(object$pCCA) && object$pCCA$rank) + stop("by = 'seqaxis' is not implemented yet for partial models") + object$CCA$QR <- qr(object$CCA$u) sol <- permutest(object, permutations = permutations, model = model, by = "onedf", parallel = parallel) ## Reformat From e68c4fdd0cc9065aa8693784c17c2784e3e10c9e Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 29 May 2017 16:15:58 +0300 Subject: [PATCH 11/37] Revert "not yet implemented for partial models" This reverts commit 0a49fc10e11e1ec19de1f07bf036ef397e540029. --- R/anova.ccabyterm.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 25669cc9..8ab8b2e8 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -205,8 +205,6 @@ anova.ccabysaxis <- { ## replace original QR decomposition of constraints with QR ## decomposition of LC of constraints. - if (!is.null(object$pCCA) && object$pCCA$rank) - stop("by = 'seqaxis' is not implemented yet for partial models") object$CCA$QR <- qr(object$CCA$u) sol <- permutest(object, permutations = permutations, model = model, by = "onedf", parallel = parallel) From 5977db9aa772504d6953d20f3b160be4f8eaf8b6 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 29 May 2017 16:16:26 +0300 Subject: [PATCH 12/37] Revert "add sequential test for axes based on 1-df constrasts in permutest" This reverts commit 03a03d78efb531c69effb4dd0001a1e8e8679137. --- R/anova.cca.R | 7 ++----- R/anova.ccabyterm.R | 40 ---------------------------------------- 2 files changed, 2 insertions(+), 45 deletions(-) diff --git a/R/anova.cca.R b/R/anova.cca.R index 695e21f8..d78b838b 100644 --- a/R/anova.cca.R +++ b/R/anova.cca.R @@ -37,7 +37,7 @@ return(anova.ccanull(object)) ## by cases if (!is.null(by)) { - by <- match.arg(by, c("terms", "margin", "axis", "seqaxis")) + by <- match.arg(by, c("terms", "margin", "axis")) if (is.null(object$terms)) stop("model must be fitted with formula interface") sol <- switch(by, @@ -51,10 +51,7 @@ "axis" = anova.ccabyaxis(object, permutations = permutations, model = model, parallel = parallel, - cutoff = cutoff), - "seqaxis" = anova.ccabysaxis(object, - permutations = permutations, - model = model, parallel=parallel)) + cutoff = cutoff)) attr(sol, "Random.seed") <- seed attr(sol, "control") <- control return(sol) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 8ab8b2e8..129e56c6 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -192,43 +192,3 @@ class(out) <- c("anova.cca", "anova", "data.frame") out } - -### Sequential test for axes is based on permustest by="onedf": we -### replace the original QR decomposition with the QR decomposition of -### linear component axes (u) and asd for 1-df contrasts. This is -### essentially similar as partialling out previous constrained axes -### except that it implies permutation model "direct", i.e., previous -### axes (Conditions) and residuals use the same permutation. This differs from - -anova.ccabysaxis <- - function(object, permutations, model, parallel) -{ - ## replace original QR decomposition of constraints with QR - ## decomposition of LC of constraints. - object$CCA$QR <- qr(object$CCA$u) - sol <- permutest(object, permutations = permutations, - model = model, by = "onedf", parallel = parallel) - ## Reformat - 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(sol$termlabels, "Residual"), - c("Df", varname, "F", "Pr(>F)")) - head <- paste0("Permutation test for ", object$method, " under ", - model, " model\n", - "Sequential test for axes\n", - howHead(attr(permutations, "control"))) - mod <- paste("Model:", c(object$call)) - attr(out, "heading") <- c(head, mod) - attr(out, "F.perm") <- sol$F.perm - class(out) <- c("anova.cca", "anova","data.frame") - out -} From b6c4335a25d8c8f9c932fc02fd69c1ce86800c9d Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 29 May 2017 17:13:28 +0300 Subject: [PATCH 13/37] only partial out previous axes in anova.cca(..., by = "axes") --- R/anova.ccabyterm.R | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 129e56c6..c2187e3d 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -152,12 +152,13 @@ Fstat <- Fstat[seq_len(ncol(LC))] } for (i in seq_along(eig)) { - part <- paste("~ . +Condition(", - paste(names(LC)[-i], collapse = "+"), ")") - ## handle partial models (clumsily) - if (!is.null(object$pCCA) && object$pCCA$rank > 0) - part <- paste(part, "Condition(qr.X(object$pCCA$QR))", sep="+") - upfla <- update(fla, part) + if (i > 1) { + part <- paste("~ . +Condition(", + paste(names(LC)[seq_len(i)], collapse = "+"), ")") + upfla <- update(fla, part) + } else { + upfla <- fla + } ## only one axis, and cannot partial out? if (length(eig) == 1) mod <- permutest(object, permutations, model = model, @@ -166,7 +167,7 @@ mod <- permutest(update(object, upfla, data = LC), permutations, model = model, - parallel = parallel) + 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) From 8fd6730118769b46e7316fcb741ad8f66c4de8bf Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Wed, 31 May 2017 11:17:12 +0300 Subject: [PATCH 14/37] fix a bug, handle Conditions & use original data in anova.ccabyaxis - bug: we partialled out with current RDA-axis instead of all previous after the first one so that permutations were off by one axis - Conditions: we use now original formula with its Condition(), whereas previously we dropped them - original data: we use the original formula for constraints and only use LC scores for Conditions Still FAILS with subset in cca-object-tests --- R/anova.ccabyterm.R | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index c2187e3d..b1def162 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -140,7 +140,7 @@ object <- update(object, subset = object$subset) } LC <- as.data.frame(LC) - fla <- reformulate(names(LC)) + fla <- formula(object) Pvals <- rep(NA, ncol(LC)) F.perm <- matrix(ncol = ncol(LC), nrow = nperm) environment(object$terms) <- environment() @@ -151,10 +151,13 @@ Df <- Df[seq_len(ncol(LC))] Fstat <- Fstat[seq_len(ncol(LC))] } + axnams <- colnames(LC) + mf <- model.frame(object) + LC <- cbind(mf, LC) for (i in seq_along(eig)) { if (i > 1) { part <- paste("~ . +Condition(", - paste(names(LC)[seq_len(i)], collapse = "+"), ")") + paste(axnams[seq_len(i-1)], collapse = "+"), ")") upfla <- update(fla, part) } else { upfla <- fla From 7503aa3bc31dbfa8e9bdf6ccc2de00943c46121a Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Wed, 31 May 2017 12:10:56 +0300 Subject: [PATCH 15/37] more robust model.frame.cca & subset handling in anova.ccabyterm FAILS still in tests/vegan-tests. Now it seems to me that the problem is the old one: if 'data' name is equal to a function name (in vegan-tests it data name is 'df'), the function is found first in ordiParseFormula. It seems that anova by = "terms" now works with this, and perhaps also by = "margin". We got to get rid of update() within nested anova functions! --- R/anova.ccabyterm.R | 11 +++++++++-- R/model.frame.cca.R | 9 ++++++--- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index b1def162..55847141 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -125,18 +125,26 @@ 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 + ## constraints and model frame LC <- object$CCA$u + mf <- model.frame(object) ## 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 + tmp <- matrix(NA, nrow = length(object$subset), + ncol = ncol(mf)) + tmp <- as.data.frame(tmp) + colnames(tmp) <- colnames(mf) + tmp[object$subset,] <- mf + mf <- tmp object <- update(object, subset = object$subset) } LC <- as.data.frame(LC) @@ -152,7 +160,6 @@ Fstat <- Fstat[seq_len(ncol(LC))] } axnams <- colnames(LC) - mf <- model.frame(object) LC <- cbind(mf, LC) for (i in seq_along(eig)) { if (i > 1) { 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 } From fa9f9fd945bd4868dc853489f96adb396bdf459f Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Wed, 31 May 2017 12:52:59 +0300 Subject: [PATCH 16/37] update tests for working & failing anova.cca cases some anova.cca tests were commented out because they used to fail, but now they work (by = "terms"), whereas some that used to work now fail (by = "axis") and some are unchanged (by = "margin") --- tests/vegan-tests.R | 25 ++++----- tests/vegan-tests.Rout.save | 133 ++++++++++++++++++++------------------------ 2 files changed, 70 insertions(+), 88 deletions(-) diff --git a/tests/vegan-tests.R b/tests/vegan-tests.R index b83f18d7..179cc095 100644 --- a/tests/vegan-tests.R +++ b/tests/vegan-tests.R @@ -38,34 +38,31 @@ 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="axis", permutations=99) +anova(m, by="term", permutations=99) # failed before 2.5-0 +##anova(m, by="margin", permutations=99) # does not work with missing data +##anova(m, by="axis", permutations=99) # fails in 2.5-0 ## 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="axis", permutations=99) +anova(p, by="term", permutations=99) # failed before 2.5-0 +##anova(p, by="margin", permutations=99) # does not accept NA data +## anova(p, by="axis", permutations=99) # fails in 2.5-0 ## 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="axis", permutations=99) +anova(p, by="term", permutations=99) # failed before 2.5-0 +##anova(p, by="margin", permutations=99) # does not work with missing data +##anova(p, by="axis", permutations=99) # fails in 2.5-0 ### attach()ed data frame instead of data= attach(df) 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="term", permutations=99) # failed before 2.5-0 ##anova(q, by="margin", permutations=99) -anova(q, by="axis", permutations=99) +##anova(q, by="axis", permutations=99) # fails in 2.5-0 ### Check that constrained ordination functions can be embedded. ### The data.frame 'df' is still attach()ed. foo <- function(bar, Y, X, ...) diff --git a/tests/vegan-tests.Rout.save b/tests/vegan-tests.Rout.save index b0373ab0..1b7bff04 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-05-29 r72746) -- "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,26 +65,22 @@ 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="axis", permutations=99) +> anova(m, by="term", permutations=99) # failed before 2.5-0 Permutation test for cca under reduced model -Marginal tests for axes +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 = 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) +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) # does not work with missing data +> ##anova(m, by="axis", permutations=99) # fails in 2.5-0 > ## capscale > p <- capscale(fla, data=df, na.action=na.exclude, subset = Use != "Pasture" & spno > 7) > anova(p, permutations=99) @@ -98,26 +94,22 @@ 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="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 = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = object$subset) - 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 -Residual 5 30.1605 +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.59 +spno 1 3.167 0.5250 0.91 +Residual 5 30.160 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +> ##anova(p, by="margin", permutations=99) # does not accept NA data +> ## anova(p, by="axis", permutations=99) # fails in 2.5-0 > ## see that capscale can be updated and also works with 'dist' input > dis <- vegdist(dune) > p <- update(p, dis ~ .) @@ -130,26 +122,22 @@ Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.a Df SumOfSqs F Pr(>F) Model 6 1.54840 1.6423 0.11 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.03 * +poly(A1, 2) 2 0.30602 0.9737 0.52 +spno 1 0.07121 0.4532 0.74 +Residual 5 0.78568 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +> ##anova(p, by="margin", permutations=99) # does not work with missing data +> ##anova(p, by="axis", permutations=99) # fails in 2.5-0 > ### attach()ed data frame instead of data= > attach(df) > q <- cca(fla, na.action = na.omit, subset = Use != "Pasture" & spno > 7) @@ -164,25 +152,22 @@ Model 6 1.25838 1.3106 0.16 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="axis", permutations=99) +> anova(q, by="term", permutations=99) # failed before 2.5-0 Permutation test for cca under reduced model -Marginal tests for axes +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 = object$subset) - 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 -Residual 5 0.80011 +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.38 +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(q, by="margin", permutations=99) +> ##anova(q, by="axis", permutations=99) # fails in 2.5-0 > ### Check that constrained ordination functions can be embedded. > ### The data.frame 'df' is still attach()ed. > foo <- function(bar, Y, X, ...) @@ -322,12 +307,12 @@ Permutation: free Number of permutations: 99 Model: cca(formula = dune ~ A1 + Moisture + Condition(Management), data = dune.env, subset = object$subset) - 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 -Residual 10 0.91705 + Df ChiSquare F Pr(>F) +CCA1 1 0.27109 2.9561 0.04 * +CCA2 1 0.14057 1.5329 0.02 * +CCA3 1 0.08761 0.9553 0.02 * +CCA4 1 0.05624 0.6132 0.01 ** +Residual 10 0.91705 --- 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) @@ -362,9 +347,9 @@ 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 @@ -390,8 +375,8 @@ 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 @@ -418,9 +403,9 @@ 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 @@ -446,8 +431,8 @@ 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 +892,4 @@ Cumulative Proportion 0.999490 1.0000000 > > proc.time() user system elapsed - 6.044 0.056 6.096 + 4.012 0.028 4.038 From 595253e527453a6dd30d11eedf62d58c577fdfce Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Wed, 31 May 2017 18:58:54 +0300 Subject: [PATCH 17/37] do not update() formula, but refit with ordConstrained anova.cca by="axes" refits models with ordConstrained instead of updating formula and refitting. This should be more robust in nested functions. The ordConstrained() function gained new option method="pass" which will pass the old dependent data unmodified to allow various models. --- R/anova.ccabyterm.R | 65 +++++++++++++++++------------------------------------ R/ordConstrained.R | 8 ++++--- 2 files changed, 26 insertions(+), 47 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 55847141..1e2bfdbb 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -125,59 +125,36 @@ resdf <- nobs(object) - length(eig) - max(object$pCCA$rank, 0) - 1 Fstat <- eig/object$CA$tot.chi*resdf Df <- rep(1, length(eig)) - ## constraints and model frame + call <- 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) + X <- qr.X(object$CCA$QR) LC <- object$CCA$u - mf <- model.frame(object) - ## 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 - tmp <- matrix(NA, nrow = length(object$subset), - ncol = ncol(mf)) - tmp <- as.data.frame(tmp) - colnames(tmp) <- colnames(mf) - tmp[object$subset,] <- mf - mf <- tmp - object <- update(object, subset = object$subset) + ## In CA we need to de-weight X and Z + if (attr(Y, "METHOD") == "CA") { + invw <- 1/sqrt(attr(Y, "RW")) + if (!is.null(Z)) + Z <- invw * Z + X <- invw * X } - LC <- as.data.frame(LC) - fla <- formula(object) 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) - LC <- cbind(mf, LC) for (i in seq_along(eig)) { if (i > 1) { - part <- paste("~ . +Condition(", - paste(axnams[seq_len(i-1)], collapse = "+"), ")") - upfla <- update(fla, part) - } else { - upfla <- fla + object <- ordConstrained(Y, X, cbind(Z, LC[, seq_len(i-1)]), "pass") } - ## only one axis, and cannot partial out? - if (length(eig) == 1) + 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, first = TRUE) + } 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) @@ -197,7 +174,7 @@ model, " model\n", "Marginal tests for axes\n", howHead(attr(permutations, "control"))) - mod <- paste("Model:", c(object$call)) + mod <- paste("Model:", c(call)) attr(out, "heading") <- c(head, mod) attr(out, "F.perm") <- F.perm class(out) <- c("anova.cca", "anova", "data.frame") diff --git a/R/ordConstrained.R b/R/ordConstrained.R index 8602367c..6c185647 100644 --- a/R/ordConstrained.R +++ b/R/ordConstrained.R @@ -357,17 +357,18 @@ `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) ## Partial @@ -388,5 +389,6 @@ out <- c(head, call = match.call(), list("pCCA" = partial, "CCA" = constraint, "CA" = resid)) + class(out) <- "cca" out } From 365703dbab18534a525714a2e6923b5c036ee8f3 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Wed, 31 May 2017 19:20:03 +0300 Subject: [PATCH 18/37] conditions must be NULL if they were not used --- R/anova.ccabyterm.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 1e2bfdbb..b90da7c8 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -132,6 +132,8 @@ 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 <- qr.X(object$CCA$QR) LC <- object$CCA$u ## In CA we need to de-weight X and Z From 861f6b1526405b5a1dd4c62d982854c22b623435 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Wed, 31 May 2017 19:37:35 +0300 Subject: [PATCH 19/37] get back old anova headers for anova.cca(..., by="axis") --- R/anova.ccabyterm.R | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index b90da7c8..cc07148a 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -125,7 +125,8 @@ resdf <- nobs(object) - length(eig) - max(object$pCCA$rank, 0) - 1 Fstat <- eig/object$CA$tot.chi*resdf Df <- rep(1, length(eig)) - call <- object$call + ## save object: will be modified later + origobj <- object ## constraints and model matrices Y <- object$Ybar if (is.null(Y)) @@ -162,21 +163,21 @@ if (Pvals[i] > cutoff) break } - out <- data.frame(c(Df, resdf), c(eig, object$CA$tot.chi), + out <- data.frame(c(Df, resdf), c(eig, origobj$CA$tot.chi), c(Fstat, NA), c(Pvals,NA)) rownames(out) <- c(names(eig), "Residual") - if (inherits(object, c("capscale", "dbrda")) && object$adjust == 1) + if (inherits(origobj, c("capscale", "dbrda")) && origobj$adjust == 1) varname <- "SumOfSqs" - else if (inherits(object, "rda")) + else if (inherits(origobj, "rda")) varname <- "Variance" else varname <- "ChiSquare" colnames(out) <- c("Df", varname, "F", "Pr(>F)") - head <- paste0("Permutation test for ", object$method, " under ", + head <- paste0("Permutation test for ", origobj$method, " under ", model, " model\n", "Marginal tests for axes\n", howHead(attr(permutations, "control"))) - mod <- paste("Model:", c(call)) + mod <- paste("Model:", c(origobj$call)) attr(out, "heading") <- c(head, mod) attr(out, "F.perm") <- F.perm class(out) <- c("anova.cca", "anova", "data.frame") From f3e9aa420ccde2127c24960912f1ef00b9c5cb57 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Wed, 31 May 2017 19:38:47 +0300 Subject: [PATCH 20/37] anova.cca(..., by="axis") passes now more tests --- tests/vegan-tests.R | 10 ++-- tests/vegan-tests.Rout.save | 116 +++++++++++++++++++++++++++++++++----------- 2 files changed, 92 insertions(+), 34 deletions(-) diff --git a/tests/vegan-tests.R b/tests/vegan-tests.R index 179cc095..ac7d2216 100644 --- a/tests/vegan-tests.R +++ b/tests/vegan-tests.R @@ -40,29 +40,29 @@ m <- cca(fla, data=df, na.action=na.exclude, subset = Use != "Pasture" & spno anova(m, permutations=99) anova(m, by="term", permutations=99) # failed before 2.5-0 ##anova(m, by="margin", permutations=99) # does not work with missing data -##anova(m, by="axis", permutations=99) # fails in 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) anova(p, by="term", permutations=99) # failed before 2.5-0 ##anova(p, by="margin", permutations=99) # does not accept NA data -## anova(p, by="axis", permutations=99) # fails in 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) anova(p, by="term", permutations=99) # failed before 2.5-0 ##anova(p, by="margin", permutations=99) # does not work with missing data -##anova(p, by="axis", permutations=99) # fails in 2.5-0 +anova(p, by="axis", permutations=99) ### attach()ed data frame instead of data= attach(df) 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) # failed before 2.5-0 +anova(q, by="term", permutations=99) # failed before 2.5-0 ##anova(q, by="margin", permutations=99) -##anova(q, by="axis", permutations=99) # fails in 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. foo <- function(bar, Y, X, ...) diff --git a/tests/vegan-tests.Rout.save b/tests/vegan-tests.Rout.save index 1b7bff04..a4c43ad9 100644 --- a/tests/vegan-tests.Rout.save +++ b/tests/vegan-tests.Rout.save @@ -1,7 +1,7 @@ -R Under development (unstable) (2017-05-29 r72746) -- "Unsuffered Consequences" -Copyright (C) 2017 The R Foundation for Statistical Computing -Platform: x86_64-pc-linux-gnu (64-bit) +R version 3.3.1 (2016-06-21) -- "Bug in Your Hair" +Copyright (C) 2016 The R Foundation for Statistical Computing +Platform: x86_64-apple-darwin13.4.0 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. @@ -80,7 +80,23 @@ Residual 5 0.80011 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ##anova(m, by="margin", permutations=99) # does not work with missing data -> ##anova(m, by="axis", permutations=99) # fails in 2.5-0 +> anova(m, by="axis", permutations=99) +Permutation test for cca under reduced model +Marginal 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 = Use != "Pasture" & spno > 7) + Df ChiSquare F Pr(>F) +CCA1 1 0.46993 2.9366 0.06 . +CCA2 1 0.26217 1.6384 0.68 +CCA3 1 0.19308 1.2066 0.86 +CCA4 1 0.18345 1.1464 0.77 +CCA5 1 0.08871 0.5544 0.97 +CCA6 1 0.06104 0.3815 0.91 +Residual 5 0.80011 +--- +Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ## capscale > p <- capscale(fla, data=df, na.action=na.exclude, subset = Use != "Pasture" & spno > 7) > anova(p, permutations=99) @@ -90,7 +106,7 @@ 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) -Model 6 59.582 1.6462 0.04 * +Model 6 59.582 1.6462 0.03 * Residual 5 30.160 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 @@ -103,13 +119,29 @@ 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.59 -spno 1 3.167 0.5250 0.91 +poly(A1, 2) 2 10.150 0.8413 0.62 +spno 1 3.167 0.5250 0.92 Residual 5 30.160 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ##anova(p, by="margin", permutations=99) # does not accept NA data -> ## anova(p, by="axis", permutations=99) # fails in 2.5-0 +> anova(p, by="axis", permutations=99) +Permutation test for capscale under reduced model +Marginal 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 = Use != "Pasture" & spno > 7) + Df Variance F Pr(>F) +CAP1 1 25.0252 4.1487 0.02 * +CAP2 1 15.8759 2.6319 0.33 +CAP3 1 8.0942 1.3419 0.86 +CAP4 1 5.0675 0.8401 0.96 +CAP5 1 3.5671 0.5914 0.98 +CAP6 1 1.9520 0.3236 0.96 +Residual 5 30.1605 +--- +Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ## see that capscale can be updated and also works with 'dist' input > dis <- vegdist(dune) > p <- update(p, dis ~ .) @@ -120,7 +152,7 @@ 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.1 Residual 5 0.78568 > anova(p, by="term", permutations=99) # failed before 2.5-0 Permutation test for capscale under reduced model @@ -130,14 +162,28 @@ 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 1.17117 2.4844 0.03 * -poly(A1, 2) 2 0.30602 0.9737 0.52 -spno 1 0.07121 0.4532 0.74 +Management 3 1.17117 2.4844 0.02 * +poly(A1, 2) 2 0.30602 0.9737 0.55 +spno 1 0.07121 0.4532 0.90 Residual 5 0.78568 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ##anova(p, by="margin", permutations=99) # does not work with missing data -> ##anova(p, by="axis", permutations=99) # fails in 2.5-0 +> anova(p, by="axis", permutations=99) +Permutation test for capscale under reduced model +Marginal 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.52 +CAP3 1 0.14701 0.9355 0.98 +CAP4 1 0.11879 0.7560 0.98 +CAP5 1 0.04213 0.2681 0.99 +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) @@ -148,11 +194,11 @@ 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.15 Residual 5 0.80011 > ## commented tests below fail in vegan 2.1-40 because number of > ## observations changes -> anova(q, by="term", permutations=99) # failed before 2.5-0 +> 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 @@ -161,13 +207,27 @@ 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.38 -spno 1 0.08318 0.5198 0.96 +poly(A1, 2) 2 0.35127 1.0976 0.40 +spno 1 0.08318 0.5198 0.98 Residual 5 0.80011 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ##anova(q, by="margin", permutations=99) -> ##anova(q, by="axis", permutations=99) # fails in 2.5-0 +> anova(q, by="axis", permutations=99) +Permutation test for cca under reduced model +Marginal tests for axes +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) +CCA1 1 0.46993 2.9366 0.15 +CCA2 1 0.26217 1.6384 0.75 +CCA3 1 0.19308 1.2066 0.90 +CCA4 1 0.18345 1.1464 0.81 +CCA5 1 0.08871 0.5544 0.98 +CCA6 1 0.06104 0.3815 0.88 +Residual 5 0.80011 > ### Check that constrained ordination functions can be embedded. > ### The data.frame 'df' is still attach()ed. > foo <- function(bar, Y, X, ...) @@ -306,15 +366,13 @@ Marginal tests for axes Permutation: free Number of permutations: 99 -Model: cca(formula = dune ~ A1 + Moisture + Condition(Management), data = dune.env, subset = object$subset) - Df ChiSquare F Pr(>F) -CCA1 1 0.27109 2.9561 0.04 * -CCA2 1 0.14057 1.5329 0.02 * -CCA3 1 0.08761 0.9553 0.02 * -CCA4 1 0.05624 0.6132 0.01 ** -Residual 10 0.91705 ---- -Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +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.10 +CCA2 1 0.14057 1.5329 0.49 +CCA3 1 0.08761 0.9553 0.80 +CCA4 1 0.05624 0.6132 0.80 +Residual 10 0.91705 > all.equal(tab[,2], c(m$CCA$eig, m$CA$tot.chi), check.attributes=FALSE) [1] TRUE > tab[nrow(tab),1] == m$CA$rank @@ -586,7 +644,7 @@ Run 11 stress 0.1192686 Run 12 stress 0.1192683 Run 13 stress 0.1886532 Run 14 stress 0.1183186 -... Procrustes: rmse 2.260742e-06 max resid 5.293932e-06 +... Procrustes: rmse 2.260745e-06 max resid 5.293932e-06 ... Similar to previous best Run 15 stress 0.1192678 Run 16 stress 0.1183186 @@ -892,4 +950,4 @@ Cumulative Proportion 0.999490 1.0000000 > > proc.time() user system elapsed - 4.012 0.028 4.038 + 3.327 0.098 3.413 From 2fd4c5f1c39a308036955f85f4d2405291802473 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 1 Jun 2017 12:01:31 +0300 Subject: [PATCH 21/37] model.matrix.cca is based on qr.X This is more robust than the old way of deeply nested function using model.frame based on ordiParseFormula and often failing when embedding these calls. CHANGE: old version returned weighted model matrix for cca, but this returns unweighted, and provides a separate method for rda. FIXME: model matrix for constraints also contains partial terms. --- NAMESPACE | 1 + R/model.matrix.cca.R | 42 +++++++++++++++++++++++++++++------------- man/model.matrix.cca.Rd | 13 ++++++++++--- 3 files changed, 40 insertions(+), 16 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 4852e05d..f537925f 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/model.matrix.cca.R b/R/model.matrix.cca.R index e2083fbc..b0a048eb 100644 --- a/R/model.matrix.cca.R +++ b/R/model.matrix.cca.R @@ -1,19 +1,35 @@ `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 <- w * qr.X(object$CCA$QR) 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) + 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/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) From 93b6afa916000fe2b5cbd600d17ec78125ac8b42 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 1 Jun 2017 13:02:02 +0300 Subject: [PATCH 22/37] remove conditions from the model matrix of constraints --- R/model.matrix.cca.R | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/R/model.matrix.cca.R b/R/model.matrix.cca.R index b0a048eb..4cd0d91b 100644 --- a/R/model.matrix.cca.R +++ b/R/model.matrix.cca.R @@ -5,8 +5,13 @@ w <- 1/sqrt(object$rowsum) if (!is.null(object$pCCA)) Z <- w * qr.X(object$pCCA$QR) - if (!is.null(object$CCA)) - X <- w * qr.X(object$CCA$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(Z)) m$Conditions <- Z @@ -23,8 +28,11 @@ X <- Z <- NULL if (!is.null(object$pCCA)) Z <- qr.X(object$pCCA$QR) - if (!is.null(object$CCA)) + 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 From c9d1bf0248201f56a5369b44c781b50a36bafd62 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 1 Jun 2017 13:31:02 +0300 Subject: [PATCH 23/37] anova.cca by="margin" no longer uses update() Function is now more robust and embeddable, because it does not need to search data for updated models in nested functions. After this, no anova.cca by= case needs update(). Function also works now with missing data with global deletion of observations with any missing value. --- R/anova.ccabyterm.R | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index cc07148a..cd5f5013 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -67,8 +67,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 From 1080f9ee14729f31732a8086007d0d30beb6811a Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 1 Jun 2017 13:38:39 +0300 Subject: [PATCH 24/37] use model.matrix in anova.ccabyaxis --- R/anova.ccabyterm.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index cd5f5013..8c8b62f2 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -148,15 +148,15 @@ Z <- qr.X(object$pCCA$QR) else Z <- NULL - X <- qr.X(object$CCA$QR) - LC <- object$CCA$u - ## In CA we need to de-weight X and Z - if (attr(Y, "METHOD") == "CA") { - invw <- 1/sqrt(attr(Y, "RW")) - if (!is.null(Z)) - Z <- invw * Z - X <- invw * X + X <- model.matrix(object) + if (!is.null(object$pCCA)) { + Z <- X$Conditions + X <- X$Constraints + } else { + Z <- NULL } + LC <- object$CCA$u + Pvals <- rep(NA, ncol(LC)) F.perm <- matrix(ncol = ncol(LC), nrow = nperm) axnams <- colnames(LC) From b292ec030fbeec9685a368fd048847af1b94b553 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 1 Jun 2017 13:51:00 +0300 Subject: [PATCH 25/37] anova.cca by="margin" handles now missing data with listwise deletion --- R/anova.ccabyterm.R | 3 --- 1 file changed, 3 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 8c8b62f2..41231149 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -46,9 +46,6 @@ { 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 From adb475324598ba05e33e2298e49a7f1435daadb7 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Thu, 1 Jun 2017 13:51:37 +0300 Subject: [PATCH 26/37] anova.cca by= methods pass now all vegan-tests --- tests/vegan-tests.R | 8 +-- tests/vegan-tests.Rout.save | 140 +++++++++++++++++++++++++++++--------------- 2 files changed, 98 insertions(+), 50 deletions(-) diff --git a/tests/vegan-tests.R b/tests/vegan-tests.R index ac7d2216..f5eb5075 100644 --- a/tests/vegan-tests.R +++ b/tests/vegan-tests.R @@ -39,20 +39,20 @@ spno <- specnumber(dune) m <- cca(fla, data=df, na.action=na.exclude, subset = Use != "Pasture" & spno > 7) anova(m, permutations=99) anova(m, by="term", permutations=99) # failed before 2.5-0 -##anova(m, by="margin", permutations=99) # does not work with missing data +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) anova(p, by="term", permutations=99) # failed before 2.5-0 -##anova(p, by="margin", permutations=99) # does not accept NA data +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) anova(p, by="term", permutations=99) # failed before 2.5-0 -##anova(p, by="margin", permutations=99) # does not work with missing data +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) @@ -61,7 +61,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) # failed before 2.5-0 -##anova(q, by="margin", permutations=99) +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 a4c43ad9..f5806f2a 100644 --- a/tests/vegan-tests.Rout.save +++ b/tests/vegan-tests.Rout.save @@ -1,7 +1,7 @@ -R version 3.3.1 (2016-06-21) -- "Bug in Your Hair" -Copyright (C) 2016 The R Foundation for Statistical Computing -Platform: x86_64-apple-darwin13.4.0 (64-bit) +R Under development (unstable) (2017-05-31 r72750) -- "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. You are welcome to redistribute it under certain conditions. @@ -79,7 +79,18 @@ 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) # does not work with missing data +> 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 @@ -88,12 +99,12 @@ 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) -CCA1 1 0.46993 2.9366 0.06 . -CCA2 1 0.26217 1.6384 0.68 -CCA3 1 0.19308 1.2066 0.86 -CCA4 1 0.18345 1.1464 0.77 -CCA5 1 0.08871 0.5544 0.97 -CCA6 1 0.06104 0.3815 0.91 +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 @@ -106,7 +117,7 @@ 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) -Model 6 59.582 1.6462 0.03 * +Model 6 59.582 1.6462 0.04 * Residual 5 30.160 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 @@ -119,12 +130,23 @@ 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.62 -spno 1 3.167 0.5250 0.92 +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) # does not accept NA data +> 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 @@ -133,12 +155,12 @@ 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) -CAP1 1 25.0252 4.1487 0.02 * -CAP2 1 15.8759 2.6319 0.33 -CAP3 1 8.0942 1.3419 0.86 -CAP4 1 5.0675 0.8401 0.96 -CAP5 1 3.5671 0.5914 0.98 -CAP6 1 1.9520 0.3236 0.96 +CAP1 1 25.0252 4.1487 0.03 * +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 @@ -152,7 +174,7 @@ 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.1 +Model 6 1.54840 1.6423 0.14 Residual 5 0.78568 > anova(p, by="term", permutations=99) # failed before 2.5-0 Permutation test for capscale under reduced model @@ -163,12 +185,23 @@ 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 1.17117 2.4844 0.02 * -poly(A1, 2) 2 0.30602 0.9737 0.55 -spno 1 0.07121 0.4532 0.90 +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) # does not work with missing data +> 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 Marginal tests for axes @@ -178,10 +211,10 @@ 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.52 +CAP2 1 0.45691 2.9078 0.48 CAP3 1 0.14701 0.9355 0.98 -CAP4 1 0.11879 0.7560 0.98 -CAP5 1 0.04213 0.2681 0.99 +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= @@ -194,7 +227,7 @@ 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.15 +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 @@ -208,11 +241,22 @@ Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit 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.98 +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) +> 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 @@ -220,14 +264,16 @@ 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) -CCA1 1 0.46993 2.9366 0.15 -CCA2 1 0.26217 1.6384 0.75 -CCA3 1 0.19308 1.2066 0.90 -CCA4 1 0.18345 1.1464 0.81 -CCA5 1 0.08871 0.5544 0.98 -CCA6 1 0.06104 0.3815 0.88 -Residual 5 0.80011 + Df ChiSquare F Pr(>F) +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 > ### Check that constrained ordination functions can be embedded. > ### The data.frame 'df' is still attach()ed. > foo <- function(bar, Y, X, ...) @@ -367,12 +413,14 @@ Permutation: free Number of permutations: 99 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.10 -CCA2 1 0.14057 1.5329 0.49 -CCA3 1 0.08761 0.9553 0.80 -CCA4 1 0.05624 0.6132 0.80 -Residual 10 0.91705 + Df ChiSquare F Pr(>F) +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 > all.equal(tab[,2], c(m$CCA$eig, m$CA$tot.chi), check.attributes=FALSE) [1] TRUE > tab[nrow(tab),1] == m$CA$rank @@ -644,7 +692,7 @@ Run 11 stress 0.1192686 Run 12 stress 0.1192683 Run 13 stress 0.1886532 Run 14 stress 0.1183186 -... Procrustes: rmse 2.260745e-06 max resid 5.293932e-06 +... Procrustes: rmse 2.260742e-06 max resid 5.293932e-06 ... Similar to previous best Run 15 stress 0.1192678 Run 16 stress 0.1183186 @@ -950,4 +998,4 @@ Cumulative Proportion 0.999490 1.0000000 > > proc.time() user system elapsed - 3.327 0.098 3.413 + 4.144 0.044 4.184 From 5b5dc764525ce3ebaa173602409dc57677d64d05 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 09:20:19 +0300 Subject: [PATCH 27/37] expand comments on implementation of anova.cca by= cases --- R/anova.ccabyterm.R | 42 ++++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 41231149..e6b86eff 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -1,11 +1,11 @@ -### 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. 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 calls directly permutest.cca +### by = "terms" calls directly permutest.cca which decomposes the +### inertia between successive terms within compiled C code. `anova.ccabyterm` <- function(object, permutations, model, parallel) @@ -38,8 +38,17 @@ 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, ...) @@ -117,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) From 7394670b50aedc87d6d9bc304f59bcb175e9c783 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 10:02:46 +0300 Subject: [PATCH 28/37] reorder code --- R/anova.ccabyterm.R | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index e6b86eff..0f9a8202 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -157,8 +157,21 @@ resdf <- nobs(object) - length(eig) - max(object$pCCA$rank, 0) - 1 Fstat <- eig/object$CA$tot.chi*resdf Df <- rep(1, length(eig)) - ## save object: will be modified later - origobj <- object + + ## 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", + "Marginal 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)) @@ -195,22 +208,11 @@ if (Pvals[i] > cutoff) break } - out <- data.frame(c(Df, resdf), c(eig, origobj$CA$tot.chi), + 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(origobj, c("capscale", "dbrda")) && origobj$adjust == 1) - varname <- "SumOfSqs" - else if (inherits(origobj, "rda")) - varname <- "Variance" - else - varname <- "ChiSquare" colnames(out) <- c("Df", varname, "F", "Pr(>F)") - head <- paste0("Permutation test for ", origobj$method, " under ", - model, " model\n", - "Marginal tests for axes\n", - howHead(attr(permutations, "control"))) - mod <- paste("Model:", c(origobj$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 From 7d16c7f742b661cfe1104c46364856b2ad6d1d32 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 10:19:09 +0300 Subject: [PATCH 29/37] update anova.cca help --- man/anova.cca.Rd | 59 +++++++++++++++++++++++++------------------------------- 1 file changed, 26 insertions(+), 33 deletions(-) diff --git a/man/anova.cca.Rd b/man/anova.cca.Rd index 7500bbec..97596129 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. From c057adaf47cfdf3b946cc2180b095c02aefb4b04 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 10:39:19 +0300 Subject: [PATCH 30/37] more examples on anova.cca() usage --- man/anova.cca.Rd | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/man/anova.cca.Rd b/man/anova.cca.Rd index 97596129..004c1b45 100644 --- a/man/anova.cca.Rd +++ b/man/anova.cca.Rd @@ -171,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 } From cf6e3baf5bc1f4c2916a45f1128496964e686433 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 10:39:32 +0300 Subject: [PATCH 31/37] anova.ccalist could scramble long model formulae --- R/anova.ccalist.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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"))) From c218adffa7dfc7a2105d700be05aecc2975918db Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 11:59:00 +0300 Subject: [PATCH 32/37] ordConstrained gets method information from initialized Y earlier we used the supplied method argument, but after allowing method = "pass" this was insufficient and made distance-based methods wrongly analysed. Consequently anova by = "margin" and by = "axis" were wrong. It was also necessary to make permutest.cca identify distance-based methods from the x$method item as an alternative to class. This is somewhat hacky, and more fundamental change may be needed if ordConstrained(..., method="pass") is used more widely. --- R/ordConstrained.R | 18 ++++++++++++------ R/permutest.cca.R | 2 +- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/R/ordConstrained.R b/R/ordConstrained.R index 6c185647..360a3643 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 @@ -370,7 +376,7 @@ "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) diff --git a/R/permutest.cca.R b/R/permutest.cca.R index 5603ae75..050f08c7 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")) || x$method == "dbrda" ## C function to get the statististics in one loop getF <- function(indx, ...) { From c46c4306de1b9967c5bfd02910d894b47e6a54bc Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 14:16:31 +0300 Subject: [PATCH 33/37] ordConstrained() must set class to make all work with method="pass" Partial models gave inconsistent results in dbrda() because ordiYbar(..., "CA") did not know input was distance-based. This also appeared in anova(..., by="axis") which uses partial models. --- R/ordConstrained.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/ordConstrained.R b/R/ordConstrained.R index 360a3643..3307e5b5 100644 --- a/R/ordConstrained.R +++ b/R/ordConstrained.R @@ -395,6 +395,10 @@ out <- c(head, call = match.call(), list("pCCA" = partial, "CCA" = constraint, "CA" = resid)) - class(out) <- "cca" + class(out) <- switch(attr(Y, "METHOD"), + "CA" = "cca", + "PCA" = c("rda", "cca"), + "CAPSCALE" = c("capscale", "rda", "cca"), + "DISTBASED" = c("dbrda", "rda", "cca")) out } From 1a1fafccdea9aefc814d3971a2cdae55cce2bae8 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 14:34:47 +0300 Subject: [PATCH 34/37] inherits(x, "dbrda") works also with ordConstrained(..., method="pass") --- R/permutest.cca.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/permutest.cca.R b/R/permutest.cca.R index 050f08c7..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")) || x$method == "dbrda" + isDB <- inherits(x, c("dbrda")) # only dbrda is distance-based ## C function to get the statististics in one loop getF <- function(indx, ...) { From 757b45840c32a5fd1e2742189940f933f9a92d02 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 14:47:29 +0300 Subject: [PATCH 35/37] we now have forward tests for axes insteaad of marginal --- R/anova.ccabyterm.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R index 0f9a8202..e146e031 100644 --- a/R/anova.ccabyterm.R +++ b/R/anova.ccabyterm.R @@ -168,7 +168,7 @@ head <- paste0("Permutation test for ", object$method, " under ", model, " model\n", - "Marginal tests for axes\n", + "Forward tests for axes\n", howHead(attr(permutations, "control"))) head <- c(head, paste("Model:", c(object$call))) From 32e8d884a539d156aabd9bd9e9c0ab9d2e266788 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 2 Jun 2017 14:49:19 +0300 Subject: [PATCH 36/37] update vegan and cca-object tests reference output --- tests/cca-object-tests.Rout.save | 104 +++++++++++++++++++-------------------- tests/vegan-tests.Rout.save | 22 ++++----- 2 files changed, 61 insertions(+), 65 deletions(-) diff --git a/tests/cca-object-tests.Rout.save b/tests/cca-object-tests.Rout.save index b6877865..2079e8aa 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-01 r72753) -- "Unsuffered Consequences" Copyright (C) 2017 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) @@ -194,13 +194,13 @@ Eigenvalues for unconstrained axes: > ## names > sort(names(mcca)) [1] "CA" "CCA" "Ybar" "call" "call" - [6] "colsum" "grand.total" "inertia" "method" "na.action" -[11] "pCCA" "rowsum" "subset" "terminfo" "terms" -[16] "tot.chi" + [6] "colsum" "grand.total" "inertia" "method" "method" +[11] "na.action" "pCCA" "rowsum" "subset" "terminfo" +[16] "terms" "tot.chi" > sort(names(mrda)) [1] "CA" "CCA" "Ybar" "call" "colsum" "inertia" - [7] "method" "na.action" "pCCA" "subset" "terminfo" "terms" -[13] "tot.chi" + [7] "method" "method" "na.action" "pCCA" "subset" "terminfo" +[13] "terms" "tot.chi" > sort(names(mcap)) [1] "CA" "CCA" "Ybar" "adjust" "call" "colsum" [7] "inertia" "method" "pCCA" "sqrt.dist" "terminfo" "terms" @@ -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.156 0.060 3.210 diff --git a/tests/vegan-tests.Rout.save b/tests/vegan-tests.Rout.save index f5806f2a..aa81ed22 100644 --- a/tests/vegan-tests.Rout.save +++ b/tests/vegan-tests.Rout.save @@ -1,5 +1,5 @@ -R Under development (unstable) (2017-05-31 r72750) -- "Unsuffered Consequences" +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) @@ -93,7 +93,7 @@ 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 @@ -149,7 +149,7 @@ 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 @@ -204,7 +204,7 @@ 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 -Marginal tests for axes +Forward tests for axes Permutation: free Number of permutations: 99 @@ -259,7 +259,7 @@ 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 @@ -408,7 +408,7 @@ 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 @@ -446,7 +446,7 @@ 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 @@ -474,7 +474,7 @@ 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 @@ -502,7 +502,7 @@ 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 @@ -530,7 +530,7 @@ 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 @@ -998,4 +998,4 @@ Cumulative Proportion 0.999490 1.0000000 > > proc.time() user system elapsed - 4.144 0.044 4.184 + 4.132 0.044 4.171 From c9e055705eff45139e917bf40b6f2a0407e58f64 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Tue, 6 Jun 2017 11:11:37 +0300 Subject: [PATCH 37/37] ordConstrained sets method: do not do it again --- R/capscale.R | 1 - R/cca.default.R | 1 - R/dbrda.R | 1 - R/rda.default.R | 2 +- tests/cca-object-tests.Rout.save | 14 +++++++------- 5 files changed, 8 insertions(+), 11 deletions(-) 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/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/tests/cca-object-tests.Rout.save b/tests/cca-object-tests.Rout.save index 2079e8aa..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-06-01 r72753) -- "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) @@ -194,13 +194,13 @@ Eigenvalues for unconstrained axes: > ## names > sort(names(mcca)) [1] "CA" "CCA" "Ybar" "call" "call" - [6] "colsum" "grand.total" "inertia" "method" "method" -[11] "na.action" "pCCA" "rowsum" "subset" "terminfo" -[16] "terms" "tot.chi" + [6] "colsum" "grand.total" "inertia" "method" "na.action" +[11] "pCCA" "rowsum" "subset" "terminfo" "terms" +[16] "tot.chi" > sort(names(mrda)) [1] "CA" "CCA" "Ybar" "call" "colsum" "inertia" - [7] "method" "method" "na.action" "pCCA" "subset" "terminfo" -[13] "terms" "tot.chi" + [7] "method" "na.action" "pCCA" "subset" "terminfo" "terms" +[13] "tot.chi" > sort(names(mcap)) [1] "CA" "CCA" "Ybar" "adjust" "call" "colsum" [7] "inertia" "method" "pCCA" "sqrt.dist" "terminfo" "terms" @@ -2539,4 +2539,4 @@ Alopgeni -0.197400341 0.18131362 > > proc.time() user system elapsed - 3.156 0.060 3.210 + 3.092 0.072 3.163