Вычислить различие между каждой парой категориальной единицы в кадре данных в R

Question

Из непредвиденных обстоятельств matrix мы можем вычислить различие между каждой парой строк и затем преобразовать результат в data.frame.

Например, с расстоянием Брея-Кертиса мы можем получить :

# 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 (т.е. data frame, имеющий 3 столбца) из длинного формата data frame в качестве входных данных. Например, если мы преобразуем приведенный выше пример matrix как:

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

, возможно ли получить тот же вывод bray, содержащий различие Брея-Кертиса между каждой парой строк? Бьюсь об заклад, существуют эффективные dplyr или data.table решения.

Answer 1

ecodist::bcdist вызывает C реализацию расстояния Брея Кертиса, которое довольно сложно преодолеть с точки зрения времени. Однако он является однопоточным и, следовательно, возможный подход заключается в распараллеливании вычислений с использованием OpenMP через R cpp:

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

Answer 2

Будет ли это достичь того, что вы после. По сути, он просто переставляет данные длинного формата в матричный массив данных и вычисляет из него B C. Я представляю, что ваш фактический набор данных представлен в длинном формате.

library(tidyverse)

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