0

I wish to carry out a mantel test between a matrix of geographical distances in kilometers (matrix 1) and a corresponding matrix of Bray-curtis distances in which I leave out those values with zero geographical distance. (A distance of 0 km represent neighbouring samples that are not identical and should not be pooled, but are not wholly independent of one another either). Therefore, I wish to keep them in my analysis.

The simple mantel test is easy: I read in the geographical matrix using the basic 

    my_geodist_matrix <-read.table("my_geodist_matrix.txt",header=TRUE,row.names=1)

The matrix looks like this:

        1top    1bottom 2top    2bottom 3top    3bottom
    1top    0   0   20.72   20.72   127.5   127.5
    1bottom 0   0   20.72   20.72   127.5   127.5
    2top    20.72   20.72   0   0   148.137 148.137
    2bottom 20.72   20.72   0   0   148.137 148.137
    3top    127.5   127.5   148.137 148.137 0   0
    3bottom 127.5   127.5   148.137 148.137 0   0

Then, I calculate the BC distances (using the vegan package) thus:

    my_otu_table.dist<-vegdist(my_otu_table,method="bray")

(looks like this):

        1top    1bottom 2top    2bottom 3top
1bottom 0.8341957
  2top 0.9948253 0.9908943 
    2bottom 0.9919492    0.9853757   0.4667466
    3top    0.9482715    0.6923454 0.9926684 0.9882600                                                                                                                                                                                                                                                                                                                                                 3bottom  0.9463127   0.6581957   0.9901074   0.9843602    0.1683397

and do the simple mantel test 

    mantel(my_geodist_matrix ,my_otu_table.dist)

So far, no trouble at all. However, when I try to leave out the pairs with 0 in the geographical matrix and the corresponding Bray-Curtis values, it goes wrong.

    my_geodist_matrix <- lapply(my_geodist_matrix, function(x){replace(x, x == 0, NA)})

and

    my_ otu_table.dist[is.na(my_geodist_matrix)] <- NA

When I do the mantel test now, R gives the following error messages:

    Error in as.data.frame.default(x[[i]], optional = TRUE) : 
      cannot coerce class ""dist"" to a data.frame

or 

    Error in cor(as.vector(xdis), ydis, method = method, use = use) : 
      missing observations in cov/cor

How do I do the desired comparison the best way (and correctly!)?

Christoffer Bugge Harder
  • 151
  • 1
  • 8
  • Can you work out an example that we can reproduce? For instance, using available data sets in the package. However, an obvious problems is that `my_geodist_matrix <- lapply(my_geodist_matrix, function(x){replace(x, x == 0, NA)})` is not a matrix but a simple list: try just `my_geodist_matrix[my_geodist_matrix == 0] <- NA`. Even if correctly defined, you cannot set elements of a `my_otu_table.dist` with indexes of a matrix. – Jari Oksanen Aug 01 '17 at 11:10
  • Just to clarify: the matrix is indexed by rows and columns, but the distances have only the lower diagonal of the matrix. A symmetric matrix has *n x n* items, but the distances only *n x (n-1) / 2* items and so the indices do no match. You must use `as.dist()` to make the matrix match with the distances. So the following may work (provided you fix the matrix part): `my_otu_table.dist[is.na(as.dist(my_geodist_matrix))] <- NA`. (Not tested.) – Jari Oksanen Aug 01 '17 at 18:59
  • Dear Jari, thank you so much for your answer. What a dumb and embarrassing mistake with my list/matrix confusion - I guess I had had my head up somewhere dark for so long that I somehow forgot what "lapply" is short for. I am unfortunately not really able to work out an example using built-in datasets such as varespec/chem, as there appear to be no 0 values in their env/veg.dists. However, your suggestion does work fine - provided, of course, that you include na.rm=TRUE in the mantel test: mantel(my_geodist_matrix ,my_otu_table.dist, na.rm=TRUE) – Christoffer Bugge Harder Aug 06 '17 at 05:58
  • I do notice that you admonish the reader about applying the na.rm=TRUE with care for mantel tests especially for such corresponding positions with missing values in two matrices/distances. However, since I have induced these missing values myself in especially these very matching positions in order to test a specific hypothesis (and correct for the known bias that these matching positions give a falsely inflated association in the mantel statistic because they are not really independent), I allow myself to believe that this is an instance of such appropriate care. – Christoffer Bugge Harder Aug 06 '17 at 06:04
  • Finally, on a side note: I am aware of the difference between distances and matrices. However, it was and is my impression that the mantel test as implemented in the brilliant vegan package does not care whether a, say, BC-distance is tested against a geographical matrix or just the lower triangle of it. In my datasets, mantel tests of all instances along the structure of "my_geodist_matrix" and "my_geodist_matrix[as.dist(my_geodist_matrix)]" yield exactly identical results in all cases. Are there any actual or potential problems here I ought to beware of? – Christoffer Bugge Harder Aug 06 '17 at 06:15

0 Answers0