在R中计算数据框中每对分类单位之间的差异。

Question

通过一个应急矩阵，我们可以计算每一行之间的相异度并将输出结果转换为一个数据框(data.frame)。

例如，使用布雷-柯蒂斯距离(Bray-Curtis distance)，我们可以得到：

# Generate matrix -------------------------------------------------------------
set.seed(1)
ex <- matrix(data = round(runif(100000), 1), nrow = 1000, ncol = 100)
rownames(ex) <- paste0("row", 1:nrow(ex))
colnames(ex) <- paste0("col", 1:ncol(ex))
ex[1:5, 1:5]
     col1 col2 col3 col4 col5
row1  0.3  0.5  0.9  0.8  0.2
row2  0.4  0.7  1.0  0.5  0.5
row3  0.6  0.4  0.9  0.2  0.0
row4  0.9  1.0  0.4  0.4  0.5
row5  0.2  0.1  0.2  0.8  0.9

# Dissimilarity ---------------------------------------------------------------
# Example of Bray-Curtis
library(ecodist)
bray <- bcdist(ex, rmzero = FALSE)
bray <- as.matrix(bray)
bray[upper.tri(bray)] <- NA
diag(bray) <- NA

# Convert distance matrix into data.frame
bray <- reshape2::melt(bray, varnames = c("id1", "id2"))
# Remove NAs
bray <- bray[complete.cases(bray), ]

head(bray)
   id1  id2     value
2 row2 row1 0.2767599
3 row3 row1 0.3541247
4 row4 row1 0.3588235
5 row5 row1 0.3935618
6 row6 row1 0.2948328
7 row7 row1 0.4045643

现在，我想知道是否可以从长格式的数据框作为输入得到与之前相同的输出 bray（即具有3列的 data frame）。 例如，如果我们将上面提供的示例 matrix 转换为：

# From a data.frame -----------------------------------------------------------
ex_df <- reshape2::melt(ex)
colnames(ex_df) <- c("row", "col", "value")

是否有可能获得相同的bray输出结果，其中包含每一对行之间的Bray-Curtis异质性？ 我敢打赌，存在高效的dplyr或data.table解决方案。

- P. Denelle

Answer 1

这个方案可以实现你的目标吗？基本上，它只是把长格式数据重新排列成类似矩阵的数据框，并从中计算 BC。我想象你的实际数据集是以长格式呈现的。

library(tidyverse)

BC_dist <- ex_df %>% 
  spread(2,3) %>% 
  column_to_rownames("row") %>% 
  bcdist(rmzero = FALSE)

- M_Shimal

Answer 2

ecodist::bcdist 调用了一个C实现的Bray Curtis距离，从时间上来说这是非常难以击败的。然而，它是单线程的，因此一种可能的方法是使用OpenMP通过Rcpp对计算进行并行化：

bcd.cpp:

#include <omp.h>
#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::plugins(openmp)]]

// [[Rcpp::export]]
NumericMatrix bcd(NumericMatrix m) {
    int i, j, k, nr = m.nrow(), nc = m.ncol();
    NumericMatrix res(nr, nr);
    double ms, sum;

    #pragma omp parallel for private(ms, sum, j, k)
    for (i = 0; i < nr - 1; i++) {
        for (j = i + 1; j < nr; j++) {
            ms = 0;
            sum = 0;
            for (k = 0; k < nc; k++) {
                if (m(i, k) < m(j, k)) {
                    ms += m(i, k);
                } else {
                    ms += m(j, k);
                }
                sum += m(i, k) + m(j, k);
            }
            res(j, i) = 1 - 2 * ms / sum;
        }
    }

    return(res);
}

计时代码：

set.seed(0L)
library(ecodist)

nr <- 10000
nc <- 100
m <- matrix(round(runif(nr*nc), 1L), nrow=nr, ncol=nc)

library(Rcpp)
sourceCpp("bcd.cpp")

microbenchmark::microbenchmark(times=3L,
    a1 <- bcdist(m, rmzero = FALSE),
    a2 <- bcd(m))

all.equal(as.vector(a1), a2[lower.tri(a2)])
#[1] TRUE

时间：

Unit: seconds
                            expr       min       lq      mean    median        uq       max neval
 a1 <- bcdist(m, rmzero = FALSE) 24.348883 24.42572 24.496605 24.502548 24.570466 24.638384     3
                    a2 <- bcd(m)  8.365889  8.50686  8.563122  8.647831  8.661739  8.675646     3

- chinsoon12