
R version 4.2.1 (2022-06-23) -- "Funny-Looking Kid"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20 (64-bit)

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.09 .
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.39  
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(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.25
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.09 .
CCA2      1   0.26217 1.6384   0.80  
CCA3      1   0.19308 1.2066   0.89  
CCA4      1   0.18345 1.1464   0.84  
CCA5      1   0.08871 0.5544   0.99  
CCA6      1   0.06104 0.3815   0.90  
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.92  
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.98  
CAP5      1  0.04213 0.2681   1.00  
CAP6      1  0.00522 0.0332   1.00  
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.09 .
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.36  
spno         1   0.08318 0.5198   0.95  
Residual     5   0.80011                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> anova(q, by="margin", permutations=99) # works since 2.5-0
Permutation test for cca under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 99

Model: cca(formula = dune ~ Management + poly(A1, 2) + spno, na.action = na.omit, subset = Use != "Pasture" & spno > 7)
            Df ChiSquare      F Pr(>F)
Management   3   0.55418 1.1544   0.30
poly(A1, 2)  2   0.32940 1.0292   0.45
spno         1   0.08318 0.5198   0.96
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.11
CCA2      1   0.26217 1.6384   0.79
CCA3      1   0.19308 1.2066   0.92
CCA4      1   0.18345 1.1464   0.83
CCA5      1   0.08871 0.5544   0.96
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, ...)
+ {
+     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     1.0000     
Constrained    1.4456     0.3840    3
Unconstrained  2.5037     0.6650   13
Imaginary     -0.1841    -0.0489    4
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.10
CCA2      1   0.14057 1.5329   0.56
CCA3      1   0.08761 0.9553   0.82
CCA4      1   0.05624 0.6132   0.81
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
[1] TRUE
> ## clean-up
> rm(df, spno, fla, m, p, q, tab, dis, foo, .Random.seed)
> ### <--- END anova.cca test --->
> 
> ### Sven Neulinger <sneulinger@ifam.uni-kiel.de> 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.964
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.961
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.962
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.969
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.005 **
Management  3   0.39543 1.5074  0.055 . 
Use         2   0.10910 0.6238  0.927   
Manure      3   0.25490 0.9717  0.518   
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 
  1.448   0.051   1.496 
