sdray / adespatial Goto Github PK

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.

PLiiiize

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)
  }

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

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)

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?

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!

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

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.