Calculate multivariate distance between samples outside of group and centroid using betadisper

Question

I am using the vegan package with betadisper on ecological data to calculate the distance between samples and their centroid. However, I am not sure how to calculate the distance between samples and a centroid in a different group. I have several years of data spanning multiple sites, with two groupings: before (years sampled pre-2013) and after (years sampled 2013 and beyond). For each site, I want to determine the location of the before centroid, then calculate the distance between each year sampled and that centroid. E.g., 2007 vs. 'before', 2008 vs. 'before', ..., 2020 vs. 'before'.

Using the principal coordinates returned in betadisper, I can calculate the distances (after ensuring positive-definite eigenvalues) between centroids, but I am not sure how to find the distances between samples that belong to a different group (e.g., samples 2013 and beyond belong to 'after', but I want to know their distance to the 'before' centroid).

In the sample code below I filtered out a single site and assigned a grouping variable called 'year_period' to guide the betadisper grouping. Eventually I want to perform this for the entire dataset across all sites. To clarify, my ultimate goal is to find the distance between each sampled year and the 'before' centroid for a given site.

Here is a sample of the data, where rows are years for a given site, and columns are counts of species.

dist_dat <- structure(list(year_period = c("2007 before", "2008 before", 
"2009 before", "2010 before", "2011 before", "2012 before", "2013 after", 
"2014 after", "2015 after", "2016 after", "2017 after", "2018 after", 
"2019 after", "2020 after"), year = 2007:2020, site = c("HOPKINS_UC", 
"HOPKINS_UC", "HOPKINS_UC", "HOPKINS_UC", "HOPKINS_UC", "HOPKINS_UC", 
"HOPKINS_UC", "HOPKINS_UC", "HOPKINS_UC", "HOPKINS_UC", "HOPKINS_UC", 
"HOPKINS_UC", "HOPKINS_UC", "HOPKINS_UC"), species1 = c(0.166666666666667, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species2 = c(3.83333333333333, 
8.5, 9.66666666666667, 6.16666666666667, 8.66666666666667, 11.6666666666667, 
16.3333333333333, 1, 0.833333333333333, 0.833333333333333, 0.333333333333333, 
3.33333333333333, 1, 1.66666666666667), species3 = c(1.16666666666667, 
7.66666666666667, 6.5, 1.83333333333333, 2, 2.16666666666667, 
4.16666666666667, 24, 43.5, 319.333333333333, 188.5, 324, 239.833333333333, 
3402.5), species4 = c(0.166666666666667, 1.83333333333333, 0.333333333333333, 
0, 0.166666666666667, 1, 0, 1.66666666666667, 5.5, 14.1666666666667, 
20, 66.3333333333333, 15.5, 99), species5 = c(93.3333333333333, 
65.5, 72.1666666666667, 39.6666666666667, 68.6666666666667, 58.5, 
48.6666666666667, 29, 17.3333333333333, 6.33333333333333, 23.6666666666667, 
35.8333333333333, 17.5, 23.6666666666667), species6 = c(3.5, 
3.83333333333333, 3.33333333333333, 4.16666666666667, 3.83333333333333, 
3, 4.16666666666667, 3.33333333333333, 4.66666666666667, 3, 2.83333333333333, 
4, 5, 4.66666666666667), species7 = c(0, 0.5, 1.33333333333333, 
0.333333333333333, 0.5, 0, 0.5, 0.5, 1.33333333333333, 0, 5, 
11.8333333333333, 17.3333333333333, 120.333333333333), species8 = c(0.666666666666667, 
1.16666666666667, 0.833333333333333, 0.666666666666667, 1.16666666666667, 
0.5, 1.16666666666667, 0, 0, 0, 0, 0.166666666666667, 0, 0), 
    species9 = c(0, 0.333333333333333, 0.166666666666667, 0, 
    0, 0.333333333333333, 0.166666666666667, 0, 0, 0, 0, 0, 0, 
    0), species10 = c(1.33333333333333, 0.5, 0.666666666666667, 
    0.666666666666667, 0, 0.666666666666667, 0.666666666666667, 
    0, 0, 0, 0, 0, 0, 0), species11 = c(0, 0, 0, 0, 0, 0.166666666666667, 
    0.166666666666667, 0, 0.333333333333333, 0, 0, 0.166666666666667, 
    0.166666666666667, 0), species12 = c(0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 1.5), species13 = c(0, 0.5, 0.166666666666667, 
    0, 0, 0.166666666666667, 0.5, 0.5, 0, 0, 0.166666666666667, 
    0.166666666666667, 0, 0.333333333333333), species14 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species15 = c(0.166666666666667, 
    0.166666666666667, 0, 0.5, 0.166666666666667, 0, 0.166666666666667, 
    0, 0, 0.166666666666667, 0, 0.333333333333333, 0.333333333333333, 
    1.16666666666667), species16 = c(0.333333333333333, 0.166666666666667, 
    0.5, 0.166666666666667, 0.5, 0.166666666666667, 0.166666666666667, 
    0.5, 0.666666666666667, 0, 0.5, 2.5, 1.83333333333333, 3), 
    species17 = c(0, 0, 0, 0, 0, 0.333333333333333, 0, 0, 0, 
    0, 0, 0, 0, 0), species18 = c(0, 0.166666666666667, 0.333333333333333, 
    0.5, 0, 0.333333333333333, 0.833333333333333, 0, 0.833333333333333, 
    0, 1.5, 2.66666666666667, 1.5, 2.5), species19 = c(0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.166666666666667, 0), species20 = c(2.33333333333333, 
    2.16666666666667, 4.16666666666667, 3.83333333333333, 5.33333333333333, 
    6.16666666666667, 3.33333333333333, 4.33333333333333, 2.33333333333333, 
    4.5, 5, 6.66666666666667, 3.33333333333333, 8.66666666666667
    ), species21 = c(0, 0, 0, 0.166666666666667, 0, 0, 0.166666666666667, 
    0, 0, 0, 0, 0, 0, 0), species22 = c(0.5, 2.66666666666667, 
    0.833333333333333, 0.5, 0.833333333333333, 3.33333333333333, 
    0, 2.16666666666667, 0, 0.333333333333333, 1.16666666666667, 
    6.66666666666667, 0.333333333333333, 2), species23 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species24 = c(1.16666666666667, 
    0.833333333333333, 1, 1.5, 0.5, 3, 0.666666666666667, 3, 
    4.16666666666667, 2.16666666666667, 2, 4.33333333333333, 
    2, 3.66666666666667), species25 = c(0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0), species26 = c(0, 0, 0.333333333333333, 
    0.666666666666667, 0.166666666666667, 0.166666666666667, 
    0, 0.333333333333333, 0.333333333333333, 0.5, 0.666666666666667, 
    0.5, 1.5, 0.666666666666667), species27 = c(0.166666666666667, 
    0.5, 0.166666666666667, 0.333333333333333, 0.5, 1.33333333333333, 
    0, 1.66666666666667, 0.166666666666667, 0.333333333333333, 
    0.5, 1.83333333333333, 0.166666666666667, 1.16666666666667
    ), species28 = c(0, 0, 0, 0, 0, 0.166666666666667, 0, 0, 
    0.166666666666667, 0, 0, 0, 0, 0), species29 = c(0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species30 = c(0, 0.333333333333333, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species31 = c(0, 0, 
    0.333333333333333, 0.166666666666667, 0.666666666666667, 
    1, 1, 0.5, 0.166666666666667, 0.166666666666667, 0, 0.666666666666667, 
    0.333333333333333, 0.833333333333333), species32 = c(2.66666666666667, 
    0, 3.33333333333333, 0.5, 1.5, 0.666666666666667, 0.166666666666667, 
    0, 0.166666666666667, 2.83333333333333, 9.33333333333333, 
    0.5, 0.666666666666667, 13), species33 = c(0.333333333333333, 
    0, 0.5, 0.333333333333333, 0.333333333333333, 0.5, 0, 0.166666666666667, 
    0.166666666666667, 0.166666666666667, 0.166666666666667, 
    0.5, 0.166666666666667, 0), species34 = c(0.5, 0.166666666666667, 
    0, 0.5, 0.666666666666667, 2.16666666666667, 3.83333333333333, 
    7.83333333333333, 1.16666666666667, 1, 5.5, 11.6666666666667, 
    6, 5.5), species35 = c(0.166666666666667, 1, 0, 0, 0.666666666666667, 
    1.83333333333333, 0.333333333333333, 0.166666666666667, 0.166666666666667, 
    0.166666666666667, 0, 0.666666666666667, 0, 0), species36 = c(0.5, 
    2.33333333333333, 3.66666666666667, 3.66666666666667, 5.16666666666667, 
    5.83333333333333, 5, 3.66666666666667, 0, 0.333333333333333, 
    1.16666666666667, 8.66666666666667, 1, 2.5), species37 = c(0.166666666666667, 
    0, 0.166666666666667, 0, 0, 0.166666666666667, 0, 0, 0.166666666666667, 
    0, 0, 2.66666666666667, 0.666666666666667, 0.5), species38 = c(0.5, 
    0, 0.166666666666667, 0, 0, 0.166666666666667, 0, 0.166666666666667, 
    0.166666666666667, 0, 0.833333333333333, 1.16666666666667, 
    0, 0.333333333333333), species39 = c(0.333333333333333, 0.333333333333333, 
    0.5, 0, 1.83333333333333, 8.83333333333333, 77.1666666666667, 
    14, 51.6666666666667, 178.166666666667, 241, 164.666666666667, 
    31.5, 269.833333333333), species40 = c(0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0), species41 = c(0, 0, 1, 0, 0, 0, 
    0, 0, 0.166666666666667, 0, 0, 0, 0, 0), species42 = c(0, 
    0, 0, 0, 0, 0, 0, 0.166666666666667, 0, 0, 0.166666666666667, 
    0.166666666666667, 0, 0.333333333333333), species43 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species44 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species45 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species46 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5), species47 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species48 = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), species49 = c(1.33333333333333, 
    0.833333333333333, 0.333333333333333, 0.666666666666667, 
    0.666666666666667, 0.833333333333333, 0.166666666666667, 
    0.5, 0.666666666666667, 0, 1.5, 3.5, 0.833333333333333, 1
    )), row.names = c(NA, -14L), class = c("tbl_df", "tbl", "data.frame"
))

Script for calculating distance between centroids (currently only works for year_period)

#create a distance matrix
stan_max_distmat <- vegdist(dist_dat[4:52], method = "bray", na.rm = T)

#use betadisper to reduce vegdist to principal coords
disper_mat <- betadisper(stan_max_distmat, type="centroid", 
                           group=dist_dat$year_period)


shift_dist <- reshape2::melt(as.matrix(sqrt(dist(disper_mat$centroids[,disper_mat$eig>0]^2)-
                                        dist(disper_mat$centroids[,disper_mat$eig<0]^2))))%>%
            tibble::rownames_to_column("distance")

Many thanks!

Answer 1

This was crossposted to vegan in github where I gave this function that will do what was asked:

### Function to find distances from each sampling unit to each centroid
### for vegan::betadisper result.
### x (input): result object from vegan::betadisper
`betadistances` <-
    function(x)
 {
     cnt <- x$centroids
     coord <- x$vectors
     pos <- which(x$eig >= 0)
     neg <- which(x$eig < 0)
     d <- apply(cnt[,pos], 1,
                function(z) rowSums(sweep(coord[,pos], 2, z)^2))
     if (length(neg))
         d <- d - apply(cnt[, neg], 1,
                        function(z) rowSums(sweep(coord[,neg], 2, z)^2))
     d <- as.data.frame(sqrt(d))
     cbind("group" = x$group, d)
 }

This is a proof-of-the-concept function that can fail in marginal cases (for instance, this will drop dimensions). This will also work with semimetric indices (with non-semidefinite dissimilarity matrices).