使用 beta分散器计算组外样本与质心之间的多元距离

问题描述 投票：0回答：1

我在生态数据上使用 

vegan
包和 
betadisper
来计算样本与其质心之间的距离。但是，我不确定如何计算不同组中样本与质心之间的距离。我有跨多个站点的几年数据，分为两组：之前（2013 年之前采样的年份）和之后（2013 年及之后采样的年份）。对于每个站点，我想确定之前质心的位置，然后计算每年采样的质心与该质心之间的距离。例如，2007 年与“之前”、2008 年与“之前”、...、2020 年与“之前”。

使用 beta分散器中返回的主坐标，我可以计算质心之间的距离（在确保正定特征值之后），但我不确定如何找到属于不同组的样本之间的距离（例如 2013 年及以后的样本）属于“之后”，但我想知道它们到“之前”质心的距离）。

在下面的示例代码中，我过滤掉了一个站点并分配了一个名为“year_period”的分组变量来指导

betadisper
分组。最终我想对所有站点的整个数据集执行此操作。为了澄清，我的最终目标是找到每个采样年份与给定站点的“之前”质心之间的距离。

这是数据示例，其中行是给定地点的年份，列是物种计数。

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

计算质心之间距离的脚本（目前仅适用于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")

非常感谢！

r vegan mds multivariate-time-series
1个回答
1
投票

这被交叉发布到 github 中的 vegan，我在那里给出了这个函数，它将执行所要求的操作：

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

这是一个概念验证函数，在边缘情况下可能会失败（例如，这会降低维度）。这也适用于半度量索引（具有非半定相异矩阵）。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.