{"id":1880,"date":"2025-04-30T12:51:39","date_gmt":"2025-04-30T10:51:39","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=1880"},"modified":"2025-05-02T09:12:48","modified_gmt":"2025-05-02T07:12:48","slug":"fast-grouped-counts-and-means-in-r","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2025\/04\/30\/fast-grouped-counts-and-means-in-r\/","title":{"rendered":"Fast Grouped Counts and Means in R"},"content":{"rendered":"\n

Edited on 2025-05-01: Multiple improvements by Christian, especially on making Polars neater, DuckDB faster, and the plot easier to read.<\/em><\/p>\n\n\n\n

From time to time, the following questions pop up: <\/p>\n\n\n\n

    \n
  1. How to calculate grouped counts and (weighted) means?<\/li>\n\n\n\n
  2. What are fast ways to do it in R?<\/li>\n<\/ol>\n\n\n\n

    This blog post presents a couple of approaches and then compares their speed with a naive (non-scientific!) benchmark.<\/p>\n\n\n\n

    \"\"<\/figure>\n\n\n\n

    Base R<\/h2>\n\n\n\n

    There are many ways to calculate grouped counts and means in base R, e.g., aggregate()<\/code>, tapply()<\/code>, by()<\/code>, split()<\/code> + lapply()<\/code>. In my experience, the fastest way is a combination of tabulate()<\/code> and rowsum()<\/code>.<\/p>\n\n\n

    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    # Make data\nset.seed(1)\n\nn <- 1e6\n\ny <- rexp(n)\nw <- runif(n)\ng <- factor(sample(LETTERS[1:3], n, TRUE))\ndf <- data.frame(y = y, g = g, w = w)\n\n# Grouped counts\ntabulate(g)\n# 333469 333569 332962\n\n# Grouped means\nrowsum(y, g) \/ tabulate(g)\n      [,1]\n# A 1.000869\n# B 1.001043\n# C 1.000445\n\n# Grouped weighted mean\nws <- rowsum(data.frame(y = y * w, w), g)\nws[, 1L] \/ ws[, 2L]\n# 1.0022749 1.0017816 0.9997058<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    But: tabulate()<\/code> ignores missing values. To avoid problems, create an explicit missing level via factor(x, exclude = NULL<\/code>).<\/p>\n\n\n\n
    Let’s turn to some other approaches.<\/p>\n\n\n\n
    dplyr<\/h2>\n\n\n\n
    Not optimized for speed or memory, but the de-facto standard in data processing with R. I love its syntax.<\/p>\n\n\n
    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    library(tidyverse)\n\ndf <- tibble(df)\n\n# Grouped counts\ndplyr::count(df, g)\n\n# Grouped means\ndf |>\n  group_by(g) |>\n  summarize(mean(y))\n\n# Grouped weighted means\ndf |>\n  group_by(g) |>\n  summarize(sum(w * y) \/ sum(w))<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    data.table<\/a><\/h2>\n\n\n\n
    Does not need an introduction. Since 2006 the<\/em> package for fast data manipulation written in C.<\/p>\n\n\n
    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    library(data.table)\n\ndt <- data.table(df)\n\n# Grouped counts (use keyby for sorted output)\ndt[, .N, by = g]\n#         g      N\n#    <fctr>  <int>\n# 1:      C 332962\n# 2:      B 333569\n# 3:      A 333469\n\n# Grouped means\ndt[, mean(y), by = g]\n\n# Grouped weighted means\ndt[, sum(w * y) \/ sum(w), by = g]\ndt[, weighted.mean(y, w), by = g]<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    DuckDB<\/h2>\n\n\n\n
    Extremely powerful query engine \/ database system written in C++, with initial release in 2019, and R bindings since 2020. Allows larger-than-RAM calculations. <\/p>\n\n\n
    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    library(duckdb)\n\ncon <- dbConnect(duckdb())\n\n# only registers: duckdb_register(con, name = \"df\", df = df)\ndbWriteTable(con, name = \"df\", value = df)\n\ndbGetQuery(con, \"SELECT g, COUNT(*) N FROM df GROUP BY g\")\ndbGetQuery(con, \"SELECT g, AVG(y) AS mean FROM df GROUP BY g\")\ncon |> \n  dbGetQuery(\n  \"\n  SELECT g, SUM(y * w) \/ sum(w) as wmean\n  FROM df\n  GROUP BY g\n  \"\n  )\n#   g     wmean\n# 1 A 1.0022749\n# 2 B 1.0017816\n# 3 C 0.9997058<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    collapse<\/h2>\n\n\n\n
    C\/C++-based package for data transformation and statistical computing. {collapse} was initially released on CRAN in 2020. It can do much more than grouped calculations, check it out!<\/p>\n\n\n
    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    library(collapse)\n\nfcount(g)\nfnobs(g, g) # Faster and does not need memory, but ignores missing values\nfmean(y, g = g)\nfmean(y, g = g, w = w)\n#         A         B         C\n# 1.0022749 1.0017816 0.9997058<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    Polars<\/h2>\n\n\n\n
    R bindings of the fantastic Polars project<\/a> that started in 2020. First R release in 2022. Currently under heavy revision.<\/p>\n\n\n\n
    The current package is not up-to-date with the main project, thus we expect the revised version (available in this branch<\/a>)  to be faster.<\/p>\n\n\n
    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    # Sys.setenv(NOT_CRAN = \"true\")\n# install.packages(\"polars\", repos = \"https:\/\/community.r-multiverse.org\")\nlibrary(polars)\n\ndfp <- as_polars_df(df)\n\n# Grouped counts\ndfp$get_column(\"g\")$value_counts()\n# Faster, but eats more memory\ndfp$select(\"g\")$with_columns(pl$lit(1L))$group_by(\"g\")$sum()\n\n# Grouped means\n(\n  dfp\n  $group_by(\"g\")\n  $agg(pl$col(\"y\")$mean())\n)\n\n# Grouped weighted means\n(\n  dfp\n  $group_by(\"g\")\n  $agg((pl$col(\"y\") * pl$col(\"w\"))$sum() \/ pl$col(\"w\")$sum())\n)\n# shape: (3, 2)\n# \u250c\u2500\u2500\u2500\u2500\u2500\u252c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510\n# \u2502 g   \u2506 y        \u2502\n# \u2502 --- \u2506 ---      \u2502\n# \u2502 cat \u2506 f64      \u2502\n# \u255e\u2550\u2550\u2550\u2550\u2550\u256a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2561\n# \u2502 C   \u2506 0.999706 \u2502\n# \u2502 B   \u2506 1.001782 \u2502\n# \u2502 A   \u2506 1.002275 \u2502\n# \u2514\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    Naive Benchmark<\/h2>\n\n\n\n
    Let’s compare the speed of these approaches for sample sizes up to 10^8 using a Windows system with an Intel i7-13700H CPU.<\/p>\n\n\n\n
    This is not at all meant as a scientific benchmark!<\/p>\n\n\n
    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    # We run the code in a fresh session\nlibrary(tidyverse) # 2.0.0\nlibrary(duckdb) # 1.2.1\nlibrary(data.table) # 1.16.0\nlibrary(collapse) # 2.0.19\nlibrary(polars) # 0.22.3\n\npolars_info() # 8 threads\nsetDTthreads(8)\ncon <- dbConnect(duckdb(config = list(threads = \"8\")))\n\nset.seed(1)\n\nN <- 10^(5:8)\nm_queries <- 3\nresults <- vector(\"list\", length(N) * m_queries)\n\nfor (i in seq_along(N)) {\n  n <- N[i]\n\n  # Create data\n  y <- rexp(n)\n  w <- runif(n)\n  g <- factor(sample(LETTERS, n, TRUE))\n\n  df <- tibble(y = y, g = g, w = w)\n  dt <- data.table(df)\n  dfp <- as_polars_df(df)\n  dbWriteTable(con, name = \"df\", value = df, overwrite = TRUE)\n\n  # Grouped counts\n  results[[1 + (i - 1) * m_queries]] <- bench::mark(\n    base = tabulate(g),\n    dplyr = dplyr::count(df, g),\n    data.table = dt[, .N, by = g],\n    polars = dfp$get_column(\"g\")$value_counts(),\n    collapse = fcount(g),\n    duckdb = dbGetQuery(con, \"SELECT g, COUNT(*) N FROM df GROUP BY g\"),\n    check = FALSE,\n    min_iterations = 3,\n  ) |>\n    bind_cols(n = n, query = \"counts\")\n\n  results[[2 + (i - 1) * m_queries]] <- bench::mark(\n    base = rowsum(y, g) \/ tabulate(g),\n    dplyr = df |> group_by(g) |> summarize(mean(y)),\n    data.table = dt[, mean(y), by = g],\n    polars = dfp$group_by(\"g\")$agg(pl$col(\"y\")$mean()),\n    collapse = fmean(y, g = g),\n    duckdb = dbGetQuery(con, \"SELECT g, AVG(y) AS mean FROM df GROUP BY g\"),\n    check = FALSE,\n    min_iterations = 3\n  ) |>\n    bind_cols(n = n, query = \"means\")\n\n  results[[3 + (i - 1) * m_queries]] <- bench::mark(\n    base = {\n      ws <- rowsum(data.frame(y = y * w, w), g)\n      ws[, 1L] \/ ws[, 2L]\n    },\n    dplyr = df |> group_by(g) |> summarize(sum(w * y) \/ sum(w)),\n    data.table = dt[, sum(w * y) \/ sum(w), by = g],\n    polars = (\n      dfp\n      $group_by(\"g\")\n      $agg((pl$col(\"y\") * pl$col(\"w\"))$sum() \/ pl$col(\"w\")$sum())\n    ),\n    collapse = fmean(y, g = g, w = w),\n    duckdb = dbGetQuery(\n      con,\n      \"SELECT g, SUM(y * w) \/ sum(w) as wmean FROM df GROUP BY g\"\n    ),\n    check = FALSE,\n    min_iterations = 3\n  ) |>\n    bind_cols(n = n, query = \"weighted means\")\n}\n\nresults_df <- bind_rows(results) |>\n  group_by(n, query) |>\n  mutate(\n    time = as.numeric(median) * 1000, # ms\n    n = as.factor(n),\n    approach = as.character(expression),\n    relative = as.numeric(time \/ min(time))\n  ) |>\n  ungroup()\n\nggplot(\n  results_df, aes(y = time, x = n, group = approach, color = approach)\n) +\n  geom_point() +\n  geom_line() +\n  scale_y_log10(labels = scales::label_number()) +\n  facet_wrap(\"query\") +\n  labs(x = \"n\", y = \"Time [ms]\", color = element_blank()) +\n  theme_gray(base_size = 14)\n\n<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    \"\"<\/figure>\n\n\n\n
    Memory<\/h2>\n\n\n\n
    What about memory? {dplyr}, {data.table}, and rowsum()<\/code> require a lot of it, as does collapse::fcount()<\/code>. For the other approaches, almost no memory is required, or profmem can’ t measure it.<\/p>\n\n\n\n
    Final words<\/h2>\n\n\n\n