permutest of dbrda fails with first=TRUE

[jarioksa]

Cleo Tebby found the following code to fail:

library(vegan)
data(pyrifos)
ditch <- gl(12, 1, length=132)
week <- gl(11, 12, labels=c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24))
dose <- factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11))
res_rda_pyrifos_inv_BC<-dbrda(pyrifos~interaction(week,dose) + Condition(week), distance = "bray")
ctrl_pyrifos <- how(plots = Plots(strata = ditch,type = "free"), within = Within(type = "series"), nperm = 99)
permutest(res_rda_pyrifos_inv_BC,  permutations=ctrl_pyrifos, first=TRUE)
#Permutation test for dbrda 
#
#Plots: ditch, plot permutation: free
#Permutation: series
#Number of permutations: 99
# 
#Call: dbrda(formula = pyrifos ~ interaction(week, dose) +
#Condition(week), distance = "bray")
#Error in x$F.perm >= x$F.0 - EPS : invalid comparison with complex values

It seems that the permutation results (x$F.perm) are given in complex format with zero imaginary component (e.g., 10.565131+0i).

[jarioksa]

I have fixes in branch issue-#198. I'll do some testing, and if successful, merge the issue branch to master.

[jarioksa]

I assume I have fixed the bugs. There was some sloppy implementation that we didn't catch in testing when releasing the dbrda and new capscale in vegan-2.4-0.

in RDA-style methods we define the fitted values of Y as Ŷ = HY, where H is the "hat matrix" or projection matrix derived from constraints (independent variables) X. In db-RDA we use Gower-transformed symmetric dissimilarity matrix G, and following Anderson, we find its fitted values as Ĝ = HGH, where H is exactly the same hat matrix as above. In permutest.cca we analyse similarly both RDA-style and distance-based methods, only having some ifs for cases where we need different analysis for db-RDA and RDA. For instance, we do not need such exceptions for sum of all eigenvalues (total inertia) of a component, but we can use HY and HG similarly. However, we need such an exception for the first eigenvalue of db-RDA -- but we didn't have one. So we submitted a non-symmetric HG matrix to eigen analysis instead of symmetric HGH, and that gave complex eigenvalues. I have now fixed those places. (Of course we don't use hat matrix H in vegan, but we use QR decomposition. For db-RDA we may need to use it for both sides of the symmetric matrix G, or two times.)

Here a simple example to show the problem in old code and a verification it works after the fix. Both of the following would fail in permutest before the fix:

mod <- dbrda(dune ~ A1 + Moisture, dune.env, dist="bray")
modc <- dbrda(dune ~ Condition(Management) + A1 + Moisture, dune.env, dist="bray")

They will work after the fix. To see this, it is best to use a bit more controlled procedure:

p <- shuffleSet(nrow(dune), 99) # create permutation matrix
p[1,] <- sort(p[1,]) # make first permutation identical to data
tst <- permutest(mod, first = TRUE, permutations = p)
tst$F.0 - tst$F.perm[1] # compare observed statistic to first permutation

The observed statistic (tst$F.0) was found directly from the result object, and the first permutation is calculated by permutest using non-permuted data. These are identical to floating point precision which verifies that the new code works like expected.

I'll merge this (with usual provision: nothing strange appears), and it will be included in the next release of vegan.

[jarioksa]

Solved and closed with PR #199