sdray / adespatial Goto Github PK

I just found a question on Researchgate, why adespatial's beta.div.comp(..., coef="BJ", quant=FALSE) and betapart's beta.multi(..., index.family="jaccard") return different output, despite one would assume the functions run similar calculations, given the index family is set to "jaccard" and the method of Baselga et al is used.

So I tested the following:

require("ade4")
data(doubs)
A = doubs$fish[-8,]
A <- ifelse(A > 0, 1, 0)

require("adespatial")
beta.div.comp(A, coef="BJ", quant=FALSE)$part

Resulting in:

     BDtotal         Repl          Nes Repl/BDtotal  Nes/BDtotal 
   0.3258676    0.1674413    0.1584263    0.5138323    0.4861677

Whereas when I now ran

require("betapart")
beta.multi(A, index.family="jaccard")

it returns

$beta.JTU
[1] 0.7885784

$beta.JNE
[1] 0.1470249

$beta.JAC
[1] 0.9356033

Obviously, the values are not in the same order, but also they are completely different.

So I went to github and copied the code from the relevant functions of both packages. I changed variable names where variables had a pendant in the other package, to indicate which values or parts of the script are basically the same for both packages.
This is the result:

require("adespatial")
require("betapart")
require("ade4")
# load some example data
data(doubs)
A = doubs$fish[-8,]
A <- ifelse(A > 0, 1, 0)

#-----------------------------------------------------------------------------
#### adespatial
### beta.div.comp
# first, run the function
beta.div.comp(A, coef="BJ", quant=FALSE)

## what it does:
n <- nrow(A)
a <- A %*% t(A)
b <- A %*% (1 - t(A))
c <- (1 - A) %*% t(A)
min.bc <- pmin(b, c)
D <- (b + c) / (a + b + c)               # Jaccard dissimilarity
repl <- 2 * min.bc / (a + 2 * min.bc)    # replacement, turnover
rich <- D - repl

D <- as.dist(D)
repl <- as.dist(repl)
rich <- as.dist(rich)

## output values:
# turnover/replacement
total.div <- sum(D) / (n * (n - 1))     # == mean(D) / 2; btw, why is it not "just" the mean?
# nestedness
repl.div <- sum(repl) / (n * (n - 1))   # == mean(repl) / 2
# total
rich.div <- sum(rich) / (n * (n - 1))   # == mean(rich) / 2

## the following produces the same values using betapart:
mean(beta.pair(A, index.family = "jaccard")$beta.jac) / 2
mean(beta.pair(A, index.family = "jaccard")$beta.jtu) / 2
mean(beta.pair(A, index.family = "jaccard")$beta.jne) / 2

#-----------------------------------------------------------------------------
#### betapart
### beta.multi
# first, run the function
beta.multi(A, index.family="jaccard")

## what it does:
a <- A %*% t(A)
c <-  abs(sweep(a, 2, diag(a)))
sumSi <- sum(diag(a))                   # species by site richness
St <- sum(colSums(A) > 0)               # regional species richness; or ncol(A), if all columns contain values > 0
ms.a <- sumSi - St                      # multi site shared species term

max.bc <- pmax(c, t(c))
min.bc <- pmin(c, t(c))

sum.max.bc <- sum(max.bc[lower.tri(max.bc)]) # == sum(as.dist(max.bc))
sum.min.bc <- sum(min.bc[lower.tri(min.bc)]) # == sum(as.dist(min.bc))

## output values:
# turnover/replacement
beta.jtu <- (2 * sum.min.bc) / (ms.a + (2 * sum.min.bc))
# nestedness
beta.jne <- (ms.a / (ms.a + (2 * sum.min.bc))) * ((sum.max.bc - sum.min.bc) / ((ms.a) + sum.max.bc + sum.min.bc))
# total
beta.jac <- (sum.min.bc + sum.max.bc) / (ms.a + sum.min.bc + sum.max.bc)

As you can see, there are some basic equations (the ones described in the relevant papers on the partitioning of beta diversity), which are similar for both approaches. However, the adespatial function first calculates some diversity matrices and then sums them up while the betapart approach first summarises the input matrix to obtain single values and then applies the equations for betadiversity decomposition.

Now, my question would be: Why are there different outputs? Are there errors in the code, or are the functions supposed to behave differently? In case the differences are intended, could you please elaborate why or add to the documentation?

I've been reading through the adespatial docs and the cited literature (Dray 2006, and the 3 Bauman 2018 papers) and I continue to struggle identifying what the "global" MEM selection method is doing.

In Bauman 2018 "Optimizing the choice of a spatial weighting matrix in eigenvector-based methods" writes

... our method consists in (1) performing a global test (based on the R2 of the model considering all MEM variables as explanatory variables)

This makes it sounds like the Global test method is only possible in a univariate case. It does so by calculating lm(y ~ all_mems).

In Bauman et al 2018 "Disentangling good from bad practices in the selection of spatial or phylogenetic
eigenvectors" the first approach is referred to as the "AIC approach." Is this the same as the "global"? It does this by sorting MEMs based on their R2 (which, again, sounds only possible in the univariate case). Then, MEMs are added in a stepwise manner to a regression. Then the subset of mems that minimizes AICc is returned.

However, this figure

when coupled with the source code's .testGlobal

This function performs "Redundancy Analysis" from {vegan}. Looking into this is a bit hard as the only source cited is from Legendre & Legendre p. 650 which is a textbook that can not readily be viewed. But based on Legendre & Anderson 1999—which contains math notation beyond my comprehension—RDA is a type of ANOVA that derives dummy variables from PCoA.

All that to say, it is unclear as to what the global method is actually doing under the hood. Further, all my reading suggests that it can only be used in a univariate case. However adespatial supports multivariate data with the global method.

Can you clarify how this method works? Perhaps there is a paper I've missed that spells it out clearly.

@sdray @aursiber rgeos was archived two days ago, so please do move to submit an updated version of adespatial as soon as possible. See https://github.com/r-spatial/evolution https://r-spatial.github.io/evolution/ for details.

Dear adespatial admins,

Thank you for the great package! I would like to perform a distance based Moran's Eigenvector Map (db-MEM) for further use in RDA to detect local adapations. Following the guidelines for spatial weighting I truncated the distance matrix based on the
minimum spanning tree (graph-based MEM) and want to weight my neighbours based on (highly accurate) cumulative landscape resistance values. However, I am unable to get my landscape resistances in the right formatting for weighting.

Neighbour list object:
Number of regions: 29
Number of nonzero links: 164
Percentage nonzero weights: 19.50059
Average number of links: 5.655172

mem.gab<- scores.listw(listw_NN, wt = resistance_weights)), MEM.autocor = c("positive"))

Error in bicenter.wt(w, row.wt = wt, col.wt = wt) :
length of row.wt must equal the number of rows in x

I can directly integrate the resistance matrix with the mat2listw but then I can't truncate the matrix as I understand is advised.

Characteristics of weights list object:
Neighbour list object:
Number of regions: 29
Number of nonzero links: 812
Percentage nonzero weights: 96.55172
Average number of links: 28

Weights style: M
Weights constants summary:
n nn S0 S1 S2
M 29 841 13149.41 902809.7 33894518

Would you have any advise on how to format the resistance matrix for weighting?

Thanks a lot!
With kind regards,
Frederik Van Daele

Some time ago, soi.graph() stopped using compiled code because of questions with regard to whether that code was free or not. It then used rgeos and RANN. RANN is being dropped as Suggests: in this release of spdep, and dbscan is replacing it. RANN did fast k-nearest neighbours, but its porting of ANN did not permit its use for distance threshold neighbours, while dbscan ports ANN in a slightly better way and can search for distance neighbours. This means that your l. 187 should have checked that RANN was available, and must now check that dbscan is available - a simple conditioning through require("dbscan") is enough.

Some functions (e.g. mem.select) could take an orthobasis directly as input (e.g., generated by adephylo to allows to deal with spatial and phylogenetic autocor). For mem.slect, this would require to update the MIR procedure (as it requires the listw object to copute Moran's I)

This is a slightly tricky problem because, if I run the first line of the example code of the aem.time function :

out <- aem.time(20, moran = TRUE)

There is the following warning :

Warning messages:
1: In mat2listw(mat01) :
style is M (missing); style should be set to a valid value
2: In mat2listw(mat01) :
no-neighbour observations found, set zero.policy to TRUE;
this warning will soon become an error

But there are no error. However, when calling the object out, there is the following error message :

Error in print.listw(x) :
regions with no neighbours found, use zero.policy=TRUE

After a bit of digging and some clarification from r-spatial/spdep#136, found that the warnings and error can be solved by changing L134 of the aem.time function from

lw <- mat2listw(mat01)

to

lw <- mat2listw(mat01, style = "W", zero.policy = TRUE)

Please see the below reprex:

library(spdep)
library(sfdep)
library(adespatial)

geo <- guerry$geometry

nb <- st_knn(geo, 10)
#> ! Polygon provided. Using point on surface.

wt <- st_kernel_weights(nb, geo, "gaussian", adaptive = TRUE)
#> ! Polygon provided. Using point on surface.


knn_lw <- nb2listw(nb, wt, "B")
contiguity_lw <- nb2listw(st_contiguity(geo))
del_lw <- nb2listw(st_nb_delaunay(geo), style = "B")
#> ! Polygon provided. Using point on surface.

x <- guerry$literacy

select_res <- listw.select(x, list(knn_lw, contiguity_lw, del_lw))
#> Procedure stopped (alpha criteria): pvalue for variable 10 is 0.054000 (> 0.050000)
#> Procedure stopped (adjR2thresh criteria) adjR2cum = 0.792070 with 10 variables (> 0.780129)
#> Procedure stopped (adjR2thresh criteria) adjR2cum = 0.706226 with 5 variables (> 0.703867)
#> Error in `[<-.data.frame`(`*tmp*`, , 1, value = c(0.763253300428723, 0.780129191853915, : replacement has 3 rows, data has 0

varipart function in ade4 should be improved to facilitate computation. One basic function for computation (weighted, unweighted cases / dudi, df matrix or vector as entries) and one for tests. Unweighted/parametric case by QR ?This should speed up the computations

Anyone have a solution?

ERROR: package installation failed
Error: Failed to install 'adespatial' from GitHub:
System command 'Rcmd.exe' failed, exit status: 1, stdout & stderr were printed
In addition: There were 28 warnings (use warnings() to see them)

Hey Stéphane,

I've been unable to return Moran's I values from the aem.time function. It looks like the code is missing from the function installed from CRAN while your code has it commented out.

EDIT: I just read in your pull requests that you've merged most of the functions from AEM to adespatial but haven't quite finished the aem.time function. I've managed to code my own workaround by downloading the AEM package separately and copying your code in (with minor modifications). I wanted to leave this issue up in case anyone else runs into the same problem.

  if (moran) {
    require(AEM, quietly=TRUE)
    nb <- spdep::cell2nb(n,1)
    links <- spdep::listw2sn(spdep::nb2listw(nb))[,1:2]
    fr.to.aem <- links[(links[,2] - links[,1]) == 1, 1:2] #remove link duplicates
    res <- AEM::moran.I.multi(E.svd$u[,1:k], link=fr.to.aem, weight=w, plot.res=plot.moran)
    Moran <- res$res.mat[,1:2]
    positive <- rep(FALSE,k)
    positive[which(Moran[,1] > res$expected)] <- TRUE
    Moran <- cbind(as.data.frame(Moran), positive)
    colnames(Moran) <- c("Moran","p.value","Positive")
    out <- list(E=E, values=E.svd$d[1:k]^2/(n-1), aem=E.svd$u[,1:k],
                Moran=Moran, expected_Moran=res$expected)
  }

Dear @sdray ,
I'm trying to run the example of the new function msr.4thcorner, but I keep getting this error:

Error in hist.default(vec.sim, plot = FALSE, nclass = 10) : character(0)
In addition: Warning messages:
1: In max(vec.sim) : no non-missing arguments to max; returning -Inf
2: In min(vec.sim) : no non-missing arguments to min; returning Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf

I'm using

> R.Version()$version.string; R.Version()$platform
[1] "R version 3.4.0 (2017-04-21)"
[1] "x86_64-apple-darwin15.6.0"

and installed from github

> packageVersion("adespatial")
[1] ‘0.0.8’

I also tried to run it with my own data, but I get the exact same error. Que pensez vous?

Merci boucoup!

PLiiiize

Hello,
While I was following this sPCA tutorial, I found that the multivariate tests to detect local structuring (local.rtest) does not seem to work properly in the adespatial package. However, it works fine in the adegenet package.

library(adegenet) ; library(adespatial)
data(spcaIllus)
obj <- spcaIllus$dat4
mySpca <- spca(obj,type=1,scannf=FALSE,plot.nb=FALSE,nfposi=1,nfnega=0)

myGtest <- adegenet::global.rtest(obj$tab,mySpca$lw,nperm=99) ; myGtest # Works just fine
myGtest <- adespatial::global.rtest(obj$tab,mySpca$lw,nperm=99) ; myGtest # Works just fine

myLtest <- adegenet::local.rtest(X = obj$tab, listw = mySpca$lw, nperm = 99) ; myLtest # Works just fine
myLtest <- adespatial::local.rtest(X = obj$tab, listw = mySpca$lw, nperm = 99) ; myLtest # Observation: NA -> not working

I tried to understand why, and the only difference between adegenet:::local.rtest and adespatial:::local.rtest is the call to the orthobasis.listw function.

• adegenet:::local.rtest calls adegenet:::.orthobasis.listw
• adespatial:::local.rtest calls adespatial:::orthobasis.listw

The issue could be problematic since using local.rtest without precision on the package calls the adespatial::local.rtest which doesn't work properly.

Outputs should be randtest and krandtest objects. To do that, we should be able to deal with extra information in these objects (modify print and plot functions)