计算 R 中数据框中每对分类单元之间的差异答案

【问题标题】：Compute dissimilarity between each pairs of categorical unit in a data frame in R计算 R 中数据框中每对分类单元之间的差异
【发布时间】：2020-05-23 09:19:06
【问题描述】：

从偶然事件matrix，我们可以计算每对行之间的差异，然后将输出转换为data.frame。

例如，使用 Bray-Curtis 距离，我们可以：

# 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

现在，我很想知道是否有可能从长格式 data frame 中获得相同的输出 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 解决方案存在。

【问题讨论】：

标签： r dataframe dplyr data.table

【解决方案1】：

这会实现您的目标吗？基本上，它只是将长格式数据重新排列成类似矩阵的数据帧并从中计算 BC。我想你的实际数据集是长格式的。

library(tidyverse)

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

【讨论】：

【解决方案2】：

ecodist::bcdist 调用了 Bray Curtis 距离的 C 实现，这在时间上很难被击败。但是，它是单线程的，因此，一种可能的方法是通过 Rcpp 使用 OpenMP 并行化计算：

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

【讨论】：