R Under development (unstable) (2024-05-12 r86534 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ### vegan-tests: unit tests for vegan functions > > ### This file contains unit tests for vegan functions. This file is > ### run in R CMD check and its results are compared against previously > ### saved results in vegan-tests.Rout.save. If you change tests, you > ### must generate new vegan-tests.Rout.save in this directory. > > ### The current plan is that tests/ are not included in the CRAN > ### release, but only in the development versin of vegan in R-Forge. > > ### The tests here are not intended for human reading. The tests need > ### not be ecological or biologically meaningful, but they are only > ### intended for testing strange arguments, protect against > ### regressions and test correctness of results. > > ### The tests are in a single file and you should clean work space > ### after the unit test. You should set random number seed (if needed) > ### for comparison against vegan-tests.Rout.save, and you should > ### remove the seed after the unit test. If possible, avoid very long > ### lasting tests. > > ###<--- BEGIN TESTS ---> > suppressPackageStartupMessages(require(vegan)) > ###<--- BEGIN anova.cca test ---> > ### anova.cca tests: should work with (1) subset, (2) missing values, > ### (3) with functions of variables poly(A1,2), (4) variables in data > ### frame attached or in data=, or (5) in working environment > set.seed(4711) > data(dune) > data(dune.env) > df <- dune.env > df$Management[c(1,5)] <- NA > ## formula > fla <- as.formula("dune ~ Management + poly(A1, 2) + spno") > ### variable in the .GlobalEnv > spno <- specnumber(dune) > ### data= argument > ## cca/rda > m <- cca(fla, data=df, na.action=na.exclude, subset = Use != "Pasture" & spno > 7) > anova(m, permutations=99) Permutation test for cca under reduced model 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) Model 6 1.25838 1.3106 0.06 . Residual 5 0.80011 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(m, by="term", permutations=99) # failed before 2.5-0 Permutation test for cca under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 99 Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df ChiSquare F Pr(>F) Management 3 0.82392 1.7163 0.02 * poly(A1, 2) 2 0.35127 1.0976 0.37 spno 1 0.08318 0.5198 0.96 Residual 5 0.80011 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(m, by="margin", permutations=99) # works since 2.5-0 Permutation test for cca under reduced model Marginal effects of terms Permutation: free Number of permutations: 99 Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df ChiSquare F Pr(>F) Management 3 0.55418 1.1544 0.20 poly(A1, 2) 2 0.32940 1.0292 0.30 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 Forward tests for axes Permutation: free Number of permutations: 99 Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df ChiSquare F Pr(>F) CCA1 1 0.46993 2.9366 0.08 . CCA2 1 0.26217 1.6384 0.96 CCA3 1 0.19308 1.2066 1.00 CCA4 1 0.18345 1.1464 CCA5 1 0.08871 0.5544 CCA6 1 0.06104 0.3815 Residual 5 0.80011 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ## adonis > adonis2(fla, data = dune.env) Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 adonis2(formula = fla, data = dune.env) Df SumOfSqs R2 F Pr(>F) Management 3 1.4686 0.34161 3.0480 0.003 ** poly(A1, 2) 2 0.5829 0.13559 1.8146 0.057 . spno 1 0.1596 0.03713 0.9940 0.420 Residual 13 2.0879 0.48567 Total 19 4.2990 1.00000 --- 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) Permutation test for capscale under reduced model 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) Model 6 59.582 1.6462 0.04 * Residual 5 30.160 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(p, by="term", permutations=99) # failed before 2.5-0 Permutation test for capscale under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 99 Model: capscale(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df Variance F Pr(>F) Management 3 46.265 2.5566 0.01 ** poly(A1, 2) 2 10.150 0.8413 0.66 spno 1 3.167 0.5250 0.90 Residual 5 30.160 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(p, by="margin", permutations=99) # works since 2.5-0 Permutation test for capscale under reduced model Marginal effects of terms Permutation: free Number of permutations: 99 Model: capscale(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df Variance F Pr(>F) Management 3 28.7515 1.5888 0.06 . poly(A1, 2) 2 8.0847 0.6701 0.85 spno 1 3.1669 0.5250 0.88 Residual 5 30.1605 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(p, by="axis", permutations=99) Permutation test for capscale under reduced model Forward tests for axes Permutation: free Number of permutations: 99 Model: capscale(formula = dune ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df Variance F Pr(>F) CAP1 1 25.0252 4.1487 0.03 * CAP2 1 15.8759 2.6319 0.37 CAP3 1 8.0942 1.3419 0.86 CAP4 1 5.0675 0.8401 0.95 CAP5 1 3.5671 0.5914 0.96 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 ~ .) > anova(p, permutations=99) Permutation test for capscale under reduced model 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) Model 6 1.54840 1.6423 0.11 Residual 5 0.78568 > anova(p, by="term", permutations=99) # failed before 2.5-0 Permutation test for capscale under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 99 Model: capscale(formula = 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.04 * poly(A1, 2) 2 0.30602 0.9737 0.52 spno 1 0.07121 0.4532 0.86 Residual 5 0.78568 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(p, by="margin", permutations=99) # works since 2.5-0 Permutation test for capscale under reduced model Marginal effects of terms Permutation: free Number of permutations: 99 Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df SumOfSqs F Pr(>F) Management 3 0.64888 1.3765 0.25 poly(A1, 2) 2 0.26160 0.8324 0.67 spno 1 0.07121 0.4532 0.81 Residual 5 0.78568 > anova(p, by="axis", permutations=99) Permutation test for capscale under reduced model Forward tests for axes Permutation: free Number of permutations: 99 Model: capscale(formula = dis ~ Management + poly(A1, 2) + spno, data = df, na.action = na.exclude, subset = Use != "Pasture" & spno > 7) Df SumOfSqs F Pr(>F) CAP1 1 0.77834 4.9533 0.08 . CAP2 1 0.45691 2.9078 0.48 CAP3 1 0.14701 0.9355 0.99 CAP4 1 0.11879 0.7560 0.99 CAP5 1 0.04213 0.2681 1.00 CAP6 1 0.00522 0.0332 Residual 5 0.78568 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ### 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) Permutation test for cca under reduced model 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) Model 6 1.25838 1.3106 0.06 . Residual 5 0.80011 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ## 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 Permutation test for cca under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 99 Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = Use != "Pasture" & spno > 7) Df ChiSquare F Pr(>F) Management 3 0.82392 1.7163 0.04 * poly(A1, 2) 2 0.35127 1.0976 0.33 spno 1 0.08318 0.5198 0.94 Residual 5 0.80011 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(q, by="margin", permutations=99) # works since 2.5-0 Permutation test for cca under reduced model Marginal effects of terms Permutation: free Number of permutations: 99 Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = Use != "Pasture" & spno > 7) Df ChiSquare F Pr(>F) Management 3 0.55418 1.1544 0.29 poly(A1, 2) 2 0.32940 1.0292 0.42 spno 1 0.08318 0.5198 0.93 Residual 5 0.80011 > anova(q, by="axis", permutations=99) Permutation test for cca under reduced model Forward 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.09 . CCA2 1 0.26217 1.6384 0.98 CCA3 1 0.19308 1.2066 1.00 CCA4 1 0.18345 1.1464 CCA5 1 0.08871 0.5544 CCA6 1 0.06104 0.3815 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, ...) + { + bar <- match.fun(bar) + bar(Y ~ X, ...) + } > foo("cca", dune, Management, na.action = na.omit) Call: cca(formula = Y ~ X, na.action = ..1) Inertia Proportion Rank Total 2.0949 1.0000 Constrained 0.6236 0.2977 3 Unconstrained 1.4713 0.7023 14 Inertia is scaled Chi-square 2 observations deleted due to missingness Eigenvalues for constrained axes: CCA1 CCA2 CCA3 0.31573 0.20203 0.10584 Eigenvalues for unconstrained axes: CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10 CA11 0.4478 0.1910 0.1788 0.1409 0.1202 0.0949 0.0732 0.0570 0.0531 0.0448 0.0312 CA12 CA13 CA14 0.0181 0.0104 0.0098 > foo("rda", dune, Management, na.action = na.omit) Call: rda(formula = Y ~ X, na.action = ..1) Inertia Proportion Rank Total 85.1176 1.0000 Constrained 32.7765 0.3851 3 Unconstrained 52.3412 0.6149 14 Inertia is variance 2 observations deleted due to missingness Eigenvalues for constrained axes: RDA1 RDA2 RDA3 15.956 13.621 3.199 Eigenvalues for unconstrained axes: PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 14.893 9.136 6.042 5.674 3.638 2.865 2.504 1.968 1.888 1.239 0.959 PC12 PC13 PC14 0.779 0.501 0.255 > foo("capscale", dune, Management, dist="jaccard", na.action = na.omit) Call: bar(formula = Y ~ X, distance = "jaccard", na.action = ..2) Inertia Proportion Rank Total 5.1432 1.0000 Constrained 1.6453 0.3199 3 Unconstrained 3.4979 0.6801 14 Inertia is squared Jaccard distance Species scores projected from 'Y' 2 observations deleted due to missingness Eigenvalues for constrained axes: CAP1 CAP2 CAP3 0.8504 0.6047 0.1902 Eigenvalues for unconstrained axes: MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 1.0474 0.4406 0.4386 0.4054 0.2847 0.1947 0.1546 0.1506 0.0957 0.0935 0.0761 MDS12 MDS13 MDS14 0.0603 0.0436 0.0120 > foo("capscale", vegdist(dune), Management, na.action = na.omit) Call: bar(formula = Y ~ X, na.action = ..1) Inertia Proportion Rank Total 3.7652 RealTotal 3.9493 1.0000 Constrained 1.4456 0.3660 3 Unconstrained 2.5037 0.6340 13 Imaginary -0.1841 Inertia is squared Bray distance 2 observations deleted due to missingness Eigenvalues for constrained axes: CAP1 CAP2 CAP3 0.7910 0.5497 0.1050 Eigenvalues for unconstrained axes: MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 1.0756 0.3691 0.3349 0.2695 0.1651 0.0931 0.0726 0.0673 0.0286 0.0174 0.0093 MDS12 MDS13 0.0011 0.0001 > foo("capscale", dune, Management, na.action = na.omit) # fails in 2.2-1 Call: bar(formula = Y ~ X, na.action = ..1) Inertia Proportion Rank Total 85.1176 1.0000 Constrained 32.7765 0.3851 3 Unconstrained 52.3412 0.6149 14 Inertia is mean squared Euclidean distance Species scores projected from 'Y' 2 observations deleted due to missingness Eigenvalues for constrained axes: CAP1 CAP2 CAP3 15.956 13.621 3.199 Eigenvalues for unconstrained axes: MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8 MDS9 MDS10 MDS11 14.893 9.136 6.042 5.674 3.638 2.865 2.504 1.968 1.888 1.239 0.959 MDS12 MDS13 MDS14 0.779 0.501 0.255 > ## adonis must be done with detached 'df' or it will be used instead > ## of with(dune.env, ...) > detach(df) > with(dune.env, foo("adonis2", dune, Management)) Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 bar(formula = Y ~ X) Df SumOfSqs R2 F Pr(>F) X 3 1.4686 0.34161 2.7672 0.006 ** Residual 16 2.8304 0.65839 Total 19 4.2990 1.00000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ## the test case reported in github issue #285 by @ktmbiome > var <- "Moisture" > adonis2(dune ~ dune.env[, var]) Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 adonis2(formula = dune ~ dune.env[, var]) Df SumOfSqs R2 F Pr(>F) dune.env[, var] 3 1.7282 0.40199 3.5851 0.001 *** Residual 16 2.5709 0.59801 Total 19 4.2990 1.00000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > rm(var) > ### > > ### Check that statistics match in partial constrained ordination > m <- cca(dune ~ A1 + Moisture + Condition(Management), dune.env, subset = A1 > 3) > tab <- anova(m, by = "axis", permutations = 99) > m Call: cca(formula = dune ~ A1 + Moisture + Condition(Management), data = dune.env, subset = A1 > 3) Inertia Proportion Rank Total 2.0976 1.0000 Conditional 0.6251 0.2980 3 Constrained 0.5555 0.2648 4 Unconstrained 0.9170 0.4372 10 Inertia is scaled Chi-square Eigenvalues for constrained axes: CCA1 CCA2 CCA3 CCA4 0.27109 0.14057 0.08761 0.05624 Eigenvalues for unconstrained axes: CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10 0.31042 0.13634 0.11974 0.09408 0.07763 0.06425 0.04449 0.02925 0.02785 0.01299 > tab Permutation test for cca under reduced model Forward tests for axes 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.05 * CCA2 1 0.14057 1.5329 0.65 CCA3 1 0.08761 0.9553 1.00 CCA4 1 0.05624 0.6132 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 [1] TRUE > ## clean-up > rm(df, spno, fla, m, p, q, tab, dis, foo, .Random.seed) > ### <--- END anova.cca test ---> > > ### Sven Neulinger reported failures in > ### partial analysis which (mostly) were fixed in r2087. Below his test. > > set.seed(4711) > X <- matrix(rnorm(30*6), 30, 6) > > A <- factor(rep(rep(c("a","b"), each=3),5)) > B <- factor(rep(c("a","b","c"), 10)) > ## Sven Neulinger's tests failed still in 2.2-1, now due to look-up > ## order: function stats::C was found before matrix 'C'. The test was > ## OK when non-function name was used ('CC'). > C <- factor(rep(c(1:5), each=6)) > > ## partial db-RDA > 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 Forward tests for axes Blocks: strata Permutation: free Number of permutations: 999 Model: capscale(formula = X ~ A + B + Condition(C)) Df Variance F Pr(>F) CAP1 1 0.2682 1.3075 0.779 CAP2 1 0.0685 0.3339 0.997 CAP3 1 0.0455 0.2217 0.997 Residual 22 4.5130 > anova(cap.model.cond, by="terms", strata=C) # -> error pre r2287 Permutation test for capscale under reduced model Terms added sequentially (first to last) Blocks: strata Permutation: free Number of permutations: 999 Model: capscale(formula = X ~ A + B + Condition(C)) Df Variance F Pr(>F) A 1 0.1316 0.6415 0.696 B 2 0.2506 0.6108 0.821 Residual 22 4.5130 > > ## db-RDA without conditional factor > cap.model <- capscale(X ~ A + B) > anova(cap.model, by="axis", strata=C) # -> no error Permutation test for capscale under reduced model Forward tests for axes Blocks: strata Permutation: free Number of permutations: 999 Model: capscale(formula = X ~ A + B) Df Variance F Pr(>F) CAP1 1 0.2682 1.3267 0.785 CAP2 1 0.0685 0.3388 0.995 CAP3 1 0.0455 0.2249 0.995 Residual 26 5.2565 > anova(cap.model, by="terms", strata=C) # -> no error Permutation test for capscale under reduced model Terms added sequentially (first to last) Blocks: strata Permutation: free Number of permutations: 999 Model: capscale(formula = X ~ A + B) Df Variance F Pr(>F) A 1 0.1316 0.6509 0.665 B 2 0.2506 0.6198 0.834 Residual 26 5.2565 > > # partial RDA > 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 Forward tests for axes Blocks: strata Permutation: free Number of permutations: 999 Model: rda(formula = X ~ A + B + Condition(C)) Df Variance F Pr(>F) RDA1 1 0.2682 1.3075 0.770 RDA2 1 0.0685 0.3339 0.993 RDA3 1 0.0455 0.2217 0.993 Residual 22 4.5130 > anova(rda.model.cond, by="terms", strata=C) # -> error pre r2287 Permutation test for rda under reduced model Terms added sequentially (first to last) Blocks: strata Permutation: free Number of permutations: 999 Model: rda(formula = X ~ A + B + Condition(C)) Df Variance F Pr(>F) A 1 0.1316 0.6415 0.656 B 2 0.2506 0.6108 0.825 Residual 22 4.5130 > > # RDA without conditional factor > rda.model <- rda(X ~ A + B) > anova(rda.model, by="axis", strata=C) # -> no error Permutation test for rda under reduced model Forward tests for axes Blocks: strata Permutation: free Number of permutations: 999 Model: rda(formula = X ~ A + B) Df Variance F Pr(>F) RDA1 1 0.2682 1.3267 0.765 RDA2 1 0.0685 0.3388 0.996 RDA3 1 0.0455 0.2249 0.996 Residual 26 5.2565 > anova(rda.model, by="terms", strata=C) # -> no error Permutation test for rda under reduced model Terms added sequentially (first to last) Blocks: strata Permutation: free Number of permutations: 999 Model: rda(formula = X ~ A + B) Df Variance F Pr(>F) A 1 0.1316 0.6509 0.690 B 2 0.2506 0.6198 0.823 Residual 26 5.2565 > ## clean.up > rm(X, A, B, C, cap.model.cond, cap.model, rda.model.cond, rda.model) > ### end Sven Neulinger's tests > > ### Benedicte Bachelot informed us that several anova.cca* functions > ### failed if community data name was the same as a function name: the > ### function name was found first, and used instead ofa data. This > ### seems to be related to the same problem that Sven Neulinger > ### communicated, and his examples still faile if Condition or strata > ### are function names. However, the following examples that failed > ### should work now: > > set.seed(4711) > cca <- dune > m <- cca(cca ~ ., dune.env) > anova(m, by="term") Permutation test for cca under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 Model: cca(formula = cca ~ A1 + Moisture + Management + Use + Manure, data = dune.env) Df ChiSquare F Pr(>F) A1 1 0.22476 2.5704 0.012 * Moisture 3 0.51898 1.9783 0.006 ** Management 3 0.39543 1.5074 0.031 * Use 2 0.10910 0.6238 0.898 Manure 3 0.25490 0.9717 0.486 Residual 7 0.61210 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > m <- capscale(cca ~ ., dune.env) > anova(m, by="term") Permutation test for capscale under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 Model: capscale(formula = cca ~ A1 + Moisture + Management + Use + Manure, data = dune.env) Df Variance F Pr(>F) A1 1 8.1148 2.7156 0.009 ** Moisture 3 21.6497 2.4150 0.004 ** Management 3 19.1153 2.1323 0.007 ** Use 2 4.7007 0.7865 0.724 Manure 3 9.6257 1.0737 0.360 Residual 7 20.9175 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > rm(m, cca) > > ### end Benedicte Bachelot tests > > ### Richard Telford tweeted this example on 23/2/2015. Fails in 2.2-1, > ### but should work now. Also issue #100 in github.com/vegandevs/vegan. > set.seed(4711) > data(dune, dune.env) > foo <- function(x, env) { + m <- rda(x ~ Manure + A1, data = env) + anova(m, by = "margin") + } > out <- lapply(dune, foo, env = dune.env) > out$Poatriv Permutation test for rda under reduced model Marginal effects of terms Permutation: free Number of permutations: 999 Model: rda(formula = x ~ Manure + A1, data = env) Df Variance F Pr(>F) Manure 4 4.7257 5.7006 0.011 * A1 1 0.0153 0.0736 0.783 Residual 14 2.9014 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > rm(foo, out) > ### end Richard Telford test > > ### github issue #291 reported that anova(mod, by="margin") gave wrong > ### results in vegan 2.5-2 when 'mod' had only one constraining > ### variable. In such corner case, all the following models should be > ### equal > > set.seed(1046) > z <- runif(20) > p <- shuffleSet(20, 99) > mod <- rda(dune ~ z) > (a0 <- anova(mod, permutations=p)) Permutation test for rda under reduced model Permutation: free Number of permutations: 99 Model: rda(formula = dune ~ z) Df Variance F Pr(>F) Model 1 7.376 1.7298 0.08 . Residual 18 76.748 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (at <- anova(mod, permutations=p, by="term")) Permutation test for rda under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 99 Model: rda(formula = dune ~ z) Df Variance F Pr(>F) z 1 7.376 1.7298 0.08 . Residual 18 76.748 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (am <- anova(mod, permutations=p, by="margin")) Permutation test for rda under reduced model Marginal effects of terms Permutation: free Number of permutations: 99 Model: rda(formula = dune ~ z) Df Variance F Pr(>F) z 1 7.376 1.7298 0.08 . Residual 18 76.748 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (aa <- anova(mod, permutations=p, by="axis")) Permutation test for rda under reduced model Forward tests for axes Permutation: free Number of permutations: 99 Model: rda(formula = dune ~ z) Df Variance F Pr(>F) RDA1 1 7.376 1.7298 0.08 . Residual 18 76.748 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (p1 <- permutest(mod, permutations=p, by="onedf")) Permutation test for rda under reduced model Permutation: free Number of permutations: 99 Model: rda(formula = dune ~ z) Permutation test for all constrained eigenvalues Df Inertia F Pr(>F) z 1 7.376 1.7298 0.08 . Residual 18 76.748 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > all.equal(permustats(a0)$permutations, permustats(at)$permutations) [1] TRUE > all.equal(permustats(a0)$permutations, permustats(am)$permutations) [1] TRUE > all.equal(permustats(a0)$permutations, permustats(aa)$permutations) [1] TRUE > all.equal(permustats(a0)$permutations, permustats(p1)$permutations) [1] TRUE > rm(z,p,mod,a0,at,am,aa,p1) > > ### nestednodf: test case by Daniel Spitale in a comment to News on > ### the release of vegan 1.17-6 in vegan.r-forge.r-project.org. > x <- c(1,0,1,1,1,1,1,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,0,0,0) > m1 <- matrix(x, nrow=5, ncol=5, byrow=FALSE)# as in Fig 2 Almeida-Neto et al 2008. > (nodf1 <- nestednodf(m1, order = FALSE, weighted = FALSE)) N columns : 53.33333 N rows : 63.33333 NODF : 58.33333 Matrix fill: 0.56 > ## Now the same matrix but with abundance data > x <- c(5,0,2,1,1,4,1,1,0,0,0,7,1,1,0,3,1,0,0,0,9,1,0,0,0) > m <- matrix(x, nrow=5, ncol=5, byrow=FALSE) > (nodfq <- nestednodf(m, order = FALSE, weighted = FALSE)) N columns : 53.33333 N rows : 63.33333 NODF : 58.33333 Matrix fill: 0.56 > identical(nodf1, nodfq) [1] TRUE > rm(x, m, m1, nodfq, nodf1) > ### end nestednodf > > ### envfit & plot.envfit: latter failed if na.action resulted in only > ### observation with a given factor level was removed. plot.envfit would > ### fail with error about too long subscript > ### fixed case where data presented to envfit also has extraneous levels > data(dune) > data(dune.env) > ## add a new level to one of the factors > levels(dune.env$Management) <- c(levels(dune.env$Management), "foo") > ## fit nMDS and envfit > set.seed(1) > mod <- metaMDS(dune, trace = 0) > ef <- envfit(mod, dune.env, permutations = 99) > plot(mod) > plot(ef, p.max = 0.1) > rm(mod, ef) > ### end envfit & plot.envfit > > ### protest (& Procrustes analysis): Stability of the permutations and > ### other results. > data(mite) > mod <- rda(mite) > x <- scores(mod, display = "si", choices=1:6) > set.seed(4711) > xp <- x[sample(nrow(x)),] > pro <- protest(x, xp, permutations = how(nperm = 99)) > pro Call: protest(X = x, Y = xp, permutations = how(nperm = 99)) Procrustes Sum of Squares (m12 squared): 0.9389 Correlation in a symmetric Procrustes rotation: 0.2471 Significance: 0.46 Permutation: free Number of permutations: 99 > pro$t [1] 0.2142734 0.2317691 0.2325867 0.2318738 0.2297942 0.2547687 0.2150992 [8] 0.2856655 0.2323035 0.2430091 0.2607687 0.2653226 0.2127120 0.2482851 [15] 0.3044707 0.2797892 0.2644764 0.2560291 0.2916665 0.2330199 0.2434141 [22] 0.2891063 0.2871878 0.2609612 0.2158827 0.2297220 0.2417230 0.2244848 [29] 0.2177450 0.2528551 0.2616422 0.2584395 0.2412174 0.2244347 0.2264776 [36] 0.2244727 0.1973865 0.2837535 0.2634514 0.2823068 0.2699507 0.2750179 [43] 0.2628074 0.2458801 0.2635264 0.2801017 0.2296173 0.1947667 0.2929171 [50] 0.2422339 0.1831137 0.2688563 0.2555202 0.2043385 0.2287129 0.2366098 [57] 0.2141830 0.1858646 0.2616734 0.1974915 0.3197986 0.2245312 0.2712518 [64] 0.2171044 0.2241602 0.2100868 0.2183290 0.2143512 0.2642714 0.2572791 [71] 0.2806032 0.2401609 0.2063171 0.2016924 0.2684739 0.2484569 0.2200211 [78] 0.2487178 0.1925942 0.2957915 0.2991842 0.2359773 0.2440686 0.2504427 [85] 0.2919524 0.2324999 0.2086198 0.2298386 0.2727169 0.3253668 0.2292110 [92] 0.2193247 0.3517265 0.2777563 0.3713275 0.2330743 0.1986322 0.2027622 [99] 0.2281227 > rm(x, xp, pro) > ### end protest > > ### Check that functions related to predict.rda work correctly for all > ### constrained ordination methods. > > ### simulate.rda/cca/capscale: based on predict.* and the following > ### should get back the data > data(dune, dune.env) > ind <- seq_len(nrow(dune)) > target <- as.matrix(dune) > ## rda > mod <- rda(dune ~ Condition(Moisture) + Management + A1, dune.env) > dat <- simulate(mod, indx = ind) > all.equal(dat, target, check.attributes = FALSE) [1] TRUE > ## cca > mod <- cca(dune ~ Condition(Moisture) + Management + A1, dune.env) > dat <- simulate(mod, indx = ind) > all.equal(dat, target, check.attributes = FALSE) [1] TRUE > ## capscale: Euclidean distances -- non-Euclidean distances have an > ## imaginary component and will not give back the data. > d <- dist(dune) > mod <- capscale(d ~ Condition(Moisture) + Management + A1, dune.env) > dat <- simulate(mod, indx = ind) > all.equal(dat, d, check.attributes = FALSE) [1] TRUE > ## clean up > rm(ind, target, mod, dat, d) > ### end simulate.* > > ### test metaMDS works with long expression for comm > ### originally reported to GLS by Richard Telford > data(varespec) > set.seed(1) > mod <- metaMDS(subset(varespec, select = colSums(varespec) > 0, subset = rowSums(varespec) > 0), trace=0) > mod Call: metaMDS(comm = subset(varespec, select = colSums(varespec) > 0, subset = rowSums(varespec) > 0), trace = 0) global Multidimensional Scaling using monoMDS Data: wisconsin(sqrt(subset(varespec, select = colSums(varespec) > 0, subset = rowSums(varespec) > 0))) Distance: bray Dimensions: 2 Stress: 0.1825658 Stress type 1, weak ties Best solution was not repeated after 20 tries The best solution was from try 17 (random start) Scaling: centring, PC rotation, halfchange scaling Species: expanded scores based on 'wisconsin(sqrt(subset(varespec, select = colSums(varespec) > 0, subset = rowSums(varespec) > 0)))' > rm(mod) > ### The above should run without error & last lines tests changes to the > ### printed output > > ## dbrda tests > > ## the following three should be all equal > data(varespec, varechem) > (mr <- rda(varespec ~ Al + P + Condition(pH), varechem)) Call: rda(formula = varespec ~ Al + P + Condition(pH), data = varechem) Inertia Proportion Rank Total 1825.6594 1.0000 Conditional 234.0961 0.1282 1 Constrained 424.4375 0.2325 2 Unconstrained 1167.1258 0.6393 20 Inertia is variance Eigenvalues for constrained axes: RDA1 RDA2 314.60 109.84 Eigenvalues for unconstrained axes: PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 692.3 171.3 106.1 59.9 35.1 30.1 22.4 16.8 (Showing 8 of 20 unconstrained eigenvalues) > (md <- dbrda(varespec ~ Al + P + Condition(pH), varechem)) Call: dbrda(formula = varespec ~ Al + P + Condition(pH), data = varechem) Inertia Proportion Rank Total 1825.6594 1.0000 Conditional 234.0961 0.1282 1 Constrained 424.4375 0.2325 2 Unconstrained 1167.1258 0.6393 20 Inertia is mean squared Euclidean distance Eigenvalues for constrained axes: dbRDA1 dbRDA2 314.60 109.84 Eigenvalues for unconstrained axes: MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8 692.3 171.3 106.1 59.9 35.1 30.1 22.4 16.8 (Showing 8 of 20 unconstrained eigenvalues) > (mc <- capscale(varespec ~ Al + P + Condition(pH), varechem)) Call: capscale(formula = varespec ~ Al + P + Condition(pH), data = varechem) Inertia Proportion Rank Total 1825.6594 1.0000 Conditional 234.0961 0.1282 1 Constrained 424.4375 0.2325 2 Unconstrained 1167.1258 0.6393 20 Inertia is mean squared Euclidean distance Species scores projected from 'varespec' Eigenvalues for constrained axes: CAP1 CAP2 314.60 109.84 Eigenvalues for unconstrained axes: MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8 692.3 171.3 106.1 59.9 35.1 30.1 22.4 16.8 (Showing 8 of 20 unconstrained eigenvalues) > ## the following two should be zero (within 1e-15) > p <- shuffleSet(nrow(varespec), 999) > all(abs(permustats(anova(mr, permutations=p))$permutations - + permustats(anova(md, permutations=p))$permutations) + < sqrt(.Machine$double.eps)) [1] TRUE > > all(abs(permustats(anova(mr, permutations=p))$permutations - + permustats(anova(mc, permutations=p))$permutations) + < sqrt(.Machine$double.eps)) [1] TRUE > ## eigenvals returns a list now (>= 2.5-0) > data(varespec, varechem) > mod <- cca(varespec ~ Al + P + Condition(pH), varechem) > ev <- summary(eigenvals(mod)) > stopifnot(inherits(ev, "matrix")) > stopifnot(!is.list(ev)) > ev Importance of components: CCA1 CCA2 CA1 CA2 CA3 CA4 CA5 Eigenvalue 0.2756 0.13713 0.3504 0.2220 0.2010 0.17366 0.13219 Proportion Explained 0.1423 0.07078 0.1809 0.1146 0.1038 0.08964 0.06823 Cumulative Proportion 0.1423 0.21304 0.3939 0.5085 0.6123 0.70193 0.77016 CA6 CA7 CA8 CA9 CA10 CA11 CA12 Eigenvalue 0.09979 0.09005 0.07377 0.05072 0.03499 0.02761 0.02337 Proportion Explained 0.05151 0.04648 0.03808 0.02618 0.01806 0.01425 0.01206 Cumulative Proportion 0.82167 0.86815 0.90623 0.93241 0.95047 0.96472 0.97678 CA13 CA14 CA15 CA16 CA17 CA18 Eigenvalue 0.014048 0.009628 0.007234 0.005155 0.003573 0.002545 Proportion Explained 0.007251 0.004970 0.003734 0.002661 0.001845 0.001314 Cumulative Proportion 0.984031 0.989001 0.992735 0.995396 0.997240 0.998554 CA19 CA20 Eigenvalue 0.0018068 0.0009945 Proportion Explained 0.0009326 0.0005133 Cumulative Proportion 0.9994867 1.0000000 > > proc.time() user system elapsed 2.93 0.23 3.11