Computing and filtering pairwise differences in R and duckdb
We have listeners transcribe children’s speech, and to filter out potentially unreliable listeners, we will select the first 2 listeners that are within 10 percentage points (.10 proportion units) of each other.
library(dplyr) withr::local_seed(20260828) data <- tibble::tibble( id = 1:20, mu = rnorm(20) |> plogis(), phi = round(rgamma(20, shape = 40)), a = floor(mu * phi), b = phi - a, y = Map(rbeta, 4, a, b) ) |> tidyr::unnest(y) |> group_by(id) |> mutate(listener_num = seq_along(y)) |> ungroup() |> select(id, listener_num, y) data #> # A tibble: 80 × 3 #> id listener_num y #> <int> <int> <dbl> #> 1 1 1 0.226 #> 2 1 2 0.214 #> 3 1 3 0.132 #> 4 1 4 0.168 #> 5 2 1 0.911 #> 6 2 2 0.888 #> 7 2 3 0.891 #> 8 2 4 0.819 #> 9 3 1 0.147 #> 10 3 2 0.217 #> # ℹ 70 more rows
Let’s work through the problem for a single vector of observations.
outer() will populate an outer product matrix from each combination of
values from two vectors. Here is a visualization of how elements from
vectors are paired off:
x1 <- letters[1:5] x2 <- 1:5 outer(x1, x2, FUN = paste0) #> [,1] [,2] [,3] [,4] [,5] #> [1,] "a1" "a2" "a3" "a4" "a5" #> [2,] "b1" "b2" "b3" "b4" "b5" #> [3,] "c1" "c2" "c3" "c4" "c5" #> [4,] "d1" "d2" "d3" "d4" "d5" #> [5,] "e1" "e2" "e3" "e4" "e5"
For pairwise differences, here are the steps in a very base R style:
- Compute pairwise differences into matrix
- Keep just the lower triangle of the matrix
- Get the row and column indices of pairs that meet the criteria
- Find the first such pair
Row and column indices correspond to data collection order (row 1 and column 1 are listener 1). So, the first such pair is the first one with the smallest maximum row/column index. For example, if pairs 1-4 and 2-3 both satisfy the criteria, 2-3 has to be first because 4 would not have been collected yet by the time 2-3 satisfied the criteria.
limit <- .1 xs <- data |> filter(id == 1) |> pull(y) xs #> [1] 0.2264930 0.2141150 0.1319205 0.1681183 diffs <- outer(xs, xs, FUN = "-") diffs #> [,1] [,2] [,3] [,4] #> [1,] 0.00000000 0.01237806 0.09457258 0.05837471 #> [2,] -0.01237806 0.00000000 0.08219452 0.04599665 #> [3,] -0.09457258 -0.08219452 0.00000000 -0.03619787 #> [4,] -0.05837471 -0.04599665 0.03619787 0.00000000 pairs <- which(abs(diffs) <= limit, arr.ind = TRUE) # Keep the lower triangle of the matrix pairs <- pairs[pairs[, "col"] < pairs[, "row"], 1:2, drop = FALSE] pairs #> row col #> [1,] 2 1 #> [2,] 3 1 #> [3,] 4 1 #> [4,] 3 2 #> [5,] 4 2 #> [6,] 4 3 # Bc we have a lower triangle matrix, row index > col index, # so keep smallest row index row_pair <- which.min(pairs[, "row"]) pairs[row_pair, ] #> row col #> 2 1 pair <- pairs[row_pair, , drop = TRUE] |> unname() pair #> [1] 2 1
Wrapping these steps into a function, we get:
find_first_consistent_pair <- function(xs, limit = .1) { diffs <- outer(xs, xs, FUN = "-") pairs <- which(abs(diffs) <= limit, arr.ind = TRUE) # Keep the lower triangle of the matrix pairs <- pairs[pairs[, "col"] < pairs[, "row"], 1:2, drop = FALSE] if (nrow(pairs) == 0) return(c(NA_integer_, NA_integer_)) row_pair <- which.min(pairs[, "row"]) pair <- pairs[row_pair, , drop = TRUE] |> unname() pair } data_with_pairs <- data |> group_by(id) |> mutate( in_pair = listener_num %in% find_first_consistent_pair(y) ) |> ungroup() data_with_pairs #> # A tibble: 80 × 4 #> id listener_num y in_pair #> <int> <int> <dbl> <lgl> #> 1 1 1 0.226 TRUE #> 2 1 2 0.214 TRUE #> 3 1 3 0.132 FALSE #> 4 1 4 0.168 FALSE #> 5 2 1 0.911 TRUE #> 6 2 2 0.888 TRUE #> 7 2 3 0.891 FALSE #> 8 2 4 0.819 FALSE #> 9 3 1 0.147 TRUE #> 10 3 2 0.217 TRUE #> # ℹ 70 more rows
(Aside: I know I have couple of other base R approaches sitting around my computer. I can’t remember where I stashed them though. The point of this note, in fact, is to put this bit of code somewhere more permanent than a random R file on my machine.)
Database version
In production, I have these intelligibility values computed from a duckdb database so I would like to make this computation using R code that can be converted to the duckdb dialect of SQL.
First, let’s spin up a duckdb database.
db <- withr::local_db_connection(DBI::dbConnect(duckdb::duckdb())) #> duckdb is storing downloaded extensions and secrets under ~/.duckdb: #> ℹ C:\Users\Tristan/.duckdb #> This persists across sessions and is shared with the DuckDB CLI and other clients. #> ℹ Run duckdb(shared_home = FALSE) to use a temporary directory instead. #> ℹ See ?duckdb_storage for details and alternatives. dplyr::copy_to(db, data, "data") tbl(db, "data") #> # A query: ?? x 3 #> # Database: DuckDB 1.5.5 [Tristan@Windows 10 x64:R 4.6.0/:memory:] #> id listener_num y #> <int> <int> <dbl> #> 1 1 1 0.226 #> 2 1 2 0.214 #> 3 1 3 0.132 #> 4 1 4 0.168 #> 5 2 1 0.911 #> 6 2 2 0.888 #> 7 2 3 0.891 #> 8 2 4 0.819 #> 9 3 1 0.147 #> 10 3 2 0.217 #> # ℹ more rows
To do the outer product type of pairing, we do a self-join on the
tables, but we can tweak the joining criteria to keep just the
lower-triangle of pairs. I create a self_left_join() function for
friendly |> piping.
self_left_join <- function(x, ...) left_join(x, x, ...) group_vars <- rlang::syms(c("id")) tbl_pairs <- tbl(db, "data") |> self_left_join( join_by( # want A:B comparison. ordering removes redundant A:A and B:A comparisons !!! group_vars, x$listener_num < y$listener_num, ), suffix = c("_left", "_right") ) tbl_pairs #> # A query: ?? x 5 #> # Database: DuckDB 1.5.5 [Tristan@Windows 10 x64:R 4.6.0/:memory:] #> id listener_num_left y_left listener_num_right y_right #> <int> <int> <dbl> <int> <dbl> #> 1 1 1 0.226 4 0.168 #> 2 1 2 0.214 4 0.168 #> 3 1 3 0.132 4 0.168 #> 4 2 1 0.911 4 0.819 #> 5 2 2 0.888 4 0.819 #> 6 2 3 0.891 4 0.819 #> 7 3 1 0.147 4 0.209 #> 8 3 2 0.217 4 0.209 #> 9 3 3 0.208 4 0.209 #> 10 4 1 0.232 4 0.236 #> # ℹ more rows
Because listener_num_left < listener_num_right as a result of the table join,
the first such pair is the first one with the smallest listener_num_right:
limit <- .1 tbl_selected_pairs <- tbl_pairs |> group_by(!!! group_vars) |> mutate(diffs = abs(y_left - y_right)) |> filter(diffs <= limit) |> filter(listener_num_right == min(listener_num_right, na.rm = TRUE)) |> filter(listener_num_left == min(listener_num_left, na.rm = TRUE)) |> select( !!! group_vars, listener_num_left, listener_num_right ) tbl_selected_pairs #> # A query: ?? x 3 #> # Database: DuckDB 1.5.5 [Tristan@Windows 10 x64:R 4.6.0/:memory:] #> # Groups: id #> id listener_num_left listener_num_right #> <int> <int> <int> #> 1 15 2 3 #> 2 11 2 3 #> 3 16 1 2 #> 4 10 1 3 #> 5 18 1 2 #> 6 1 1 2 #> 7 3 1 2 #> 8 6 1 3 #> 9 2 1 2 #> 10 4 1 2 #> 11 12 1 2 #> 12 20 1 3 #> 13 8 1 2 #> 14 13 1 2 #> 15 19 1 2 #> 16 7 1 3 #> 17 5 1 2 #> 18 9 2 3 #> 19 14 1 2 #> 20 17 1 3
The data are in a wide format right now, so we need to pivot them into a longer shape. I am going to use duckdb’s own functions (written in all caps) to accomplish this task:
tbl_pairs_to_keep <- tbl_selected_pairs |> mutate( # unpivoting by nesting and unnesting values listener_num = LIST_VALUE(listener_num_left, listener_num_right) |> UNNEST() ) |> ungroup() |> select(!!! group_vars, listener_num) tbl_pairs_to_keep #> # A query: ?? x 2 #> # Database: DuckDB 1.5.5 [Tristan@Windows 10 x64:R 4.6.0/:memory:] #> id listener_num #> <int> <int> #> 1 4 1 #> 2 4 2 #> 3 12 1 #> 4 12 2 #> 5 20 1 #> 6 20 3 #> 7 10 1 #> 8 10 3 #> 9 18 1 #> 10 18 2 #> # ℹ more rows
Finally, we can filter down to the desired listener ids.
tbl_keep <- tbl(db, "data") |> inner_join(tbl_pairs_to_keep, by = join_by(id, listener_num)) tbl_keep #> # A query: ?? x 3 #> # Database: DuckDB 1.5.5 [Tristan@Windows 10 x64:R 4.6.0/:memory:] #> id listener_num y #> <int> <int> <dbl> #> 1 1 1 0.226 #> 2 1 2 0.214 #> 3 2 1 0.911 #> 4 2 2 0.888 #> 5 3 1 0.147 #> 6 3 2 0.217 #> 7 4 1 0.232 #> 8 4 2 0.317 #> 9 5 1 0.295 #> 10 5 2 0.267 #> # ℹ more rows
Each tbl_ here is a SQL query until we finally collect() the data
into R, so we can marvel at the SQL we generated:
dplyr::show_query(tbl_keep) #> <SQL> #> SELECT "data".* #> FROM "data" #> INNER JOIN ( #> SELECT #> id, #> UNNEST(LIST_VALUE(listener_num_left, listener_num_right)) AS listener_num #> FROM ( #> SELECT #> id, #> listener_num_left, #> y_left, #> listener_num_right, #> y_right, #> diffs, #> MIN(listener_num_left) OVER (PARTITION BY id) AS col02 #> FROM ( #> SELECT *, MIN(listener_num_right) OVER (PARTITION BY id) AS col01 #> FROM ( #> SELECT *, ABS(y_left - y_right) AS diffs #> FROM ( #> SELECT #> data_LHS.id AS id, #> data_LHS.listener_num AS listener_num_left, #> data_LHS.y AS y_left, #> data_RHS.listener_num AS listener_num_right, #> data_RHS.y AS y_right #> FROM "data" AS data_LHS #> LEFT JOIN "data" AS data_RHS #> ON ( #> data_LHS.id = data_RHS.id AND #> data_LHS.listener_num < data_RHS.listener_num #> ) #> ) AS q01 #> ) AS q01 #> WHERE (diffs <= 0.1) #> ) AS q01 #> WHERE (listener_num_right = col01) #> ) AS q01 #> WHERE (listener_num_left = col02) #> ) AS RHS #> ON ("data".id = RHS.id AND "data".listener_num = RHS.listener_num)
Finally, we can check that the two versions agree:
a <- data_with_pairs |> filter(in_pair) |> select(1, 2, 3) |> arrange(id, listener_num) b <- tbl_keep |> collect() |> arrange(id, listener_num) all(a == b) #> [1] TRUE
Leave a comment