Topic modelling without Python

A BERTopic-class pipeline — embed documents, reduce with UMAP, cluster with HDBSCAN, name the clusters — is usually a Python affair. With relm it is a handful of base-R lines over CRAN packages the R ecosystem already has, running fully local with zero Python.

The steps: llm_embed() the abstracts, lay them out with uwot::umap(), find clusters with dbscan::hdbscan(), and name each cluster by asking the model itself with llm_generate().

The corpus

The package ships a small synthetic sample so this runs offline; for the real ~5,000-abstract arXiv corpus, use tests/demos/fetch-abstracts.R.

csv <- system.file("extdata", "abstracts-sample.csv", package = "relm")
abstracts <- if (nzchar(csv)) {
  utils::read.csv(csv, stringsAsFactors = FALSE)
} else {
  # fallback when the vignette is knitted from the source tree
  utils::read.csv(file.path("..", "inst", "extdata", "abstracts-sample.csv"),
                  stringsAsFactors = FALSE)
}
nrow(abstracts)
#> [1] 1200
table(abstracts$label)
#> 
#>        Climate and atmospheric science Condensed matter and quantum materials 
#>                                    150                                    150 
#>             Cosmology and astrophysics              Cryptography and security 
#>                                    150                                    150 
#>                          Deep learning         Epidemiology and public health 
#>                                    150                                    150 
#>         Genomics and molecular biology        Number theory and combinatorics 
#>                                    150                                    150

Embed, reduce, cluster, name

set.seed() plus n_sgd_threads = 1 makes the UMAP layout reproducible; HDBSCAN is deterministic given the layout. Cluster names come from the model itself.

library(relm)
m <- llm(model_path)

emb <- llm_embed(m, abstracts$text, pooling = "mean", normalize = TRUE)

set.seed(20240707)
coords <- uwot::umap(emb, n_neighbors = 15, min_dist = 0.1, n_components = 2,
                     metric = "cosine", n_sgd_threads = 1)
cluster <- dbscan::hdbscan(coords, minPts = 15)$cluster # 0 = noise

# Name each cluster from its most central abstracts.
ids <- sort(setdiff(unique(cluster), 0))
labels <- vapply(ids, function(g) {
  members <- which(cluster == g)
  cen <- colMeans(coords[members, , drop = FALSE])
  medoid <- members[which.min(rowSums(sweep(coords[members, , drop = FALSE], 2, cen)^2))]
  prompt <- paste0(
    "Reply with a SHORT topic label of 2 to 4 words for these abstracts, ",
    "and nothing else.\n\n",
    paste0("- ", abstracts$text[medoid], collapse = "\n"), "\n\nTopic label:"
  )
  trimws(gsub("[\"'`.]", "", llm_generate(m, prompt, max_tokens = 12, temperature = 0, seed = 1)))
}, character(1))
names(labels) <- ids
close(m)
# Illustrative layout (no model/uwot/dbscan at render time): eight Gaussian blobs
# keyed off the sample's known topics, so the map matches what the real pipeline
# produces. Run the chunk above on your own model for the genuine embedding map.
set.seed(20240707)
labs <- unique(abstracts$label)
ang <- setNames(2 * pi * (seq_along(labs) - 1) / length(labs), labs)
cx <- 6 * cos(ang[abstracts$label]); cy <- 6 * sin(ang[abstracts$label])
coords <- cbind(cx + rnorm(nrow(abstracts), 0, 0.8),
                cy + rnorm(nrow(abstracts), 0, 0.8))
cluster <- as.integer(factor(abstracts$label, levels = labs))
labels <- setNames(labs, seq_along(labs))
# a synthetic embedding (one direction per topic + noise) so the quality and
# structure metrics below have something to read without a model
D <- 24
dirs <- matrix(rnorm(length(labs) * D), length(labs), D)
dirs <- dirs / sqrt(rowSums(dirs^2))
emb <- dirs[cluster, ] + matrix(rnorm(nrow(abstracts) * D, sd = 0.5), nrow(abstracts), D)
emb <- emb / sqrt(rowSums(emb^2))

The labelled map

ids <- sort(setdiff(unique(cluster), 0))
pal <- adjustcolor(grDevices::hcl.colors(length(ids), "Dark 3"), 0.6)
col <- rep(adjustcolor("grey70", 0.45), length(cluster))
for (j in seq_along(ids)) col[cluster == ids[j]] <- pal[j]

op <- par(mar = c(4.2, 4.2, 3, 1), mgp = c(2.5, 0.7, 0), las = 1)
xr <- range(coords[, 1]); yr <- range(coords[, 2]) # pad so edge labels fit
plot(coords, col = col, pch = 19, cex = 0.5, xlab = "UMAP 1", ylab = "UMAP 2",
     xlim = xr + c(-1, 1) * 0.10 * diff(xr), ylim = yr + c(-1, 1) * 0.06 * diff(yr),
     main = sprintf("Topic modelling without Python (%s)",
                    if (run_real) "real embeddings" else "illustrative"))
halo <- function(x, y, s) {
  o <- 0.006 * diff(par("usr")[1:2])
  for (dx in c(-1, 1)) for (dy in c(-1, 1)) text(x + dx * o, y + dy * o, s, col = "white", font = 2)
  text(x, y, s, font = 2)
}
for (j in seq_along(ids)) {
  members <- which(cluster == ids[j])
  cen <- colMeans(coords[members, , drop = FALSE])
  med <- members[which.min(rowSums(sweep(coords[members, , drop = FALSE], 2, cen)^2))]
  halo(coords[med, 1], coords[med, 2], labels[[as.character(ids[j])]])
}

UMAP scatter of abstracts coloured into topic clusters with topic labels.

par(op)

Each colour is a discovered topic; grey points are HDBSCAN’s noise class. These are the model’s suggested groupings, not ground truth — so the next three steps give the reader ways to judge them.

How good are the topics?

A map invites the question how good are these clusters, really? Two base-R metrics answer it, and they disagree in an informative way. Simplified silhouette (on the UMAP coordinates) scores each point by how much closer it sits to its own topic’s centre than to the nearest other topic — separation in the map you see. Cohesion is the mean cosine of a topic’s documents to their centroid in the original embedding space, where UMAP’s density distortion cannot flatter the result.

simplified_silhouette <- function(coords, cluster) {
  ids <- sort(setdiff(unique(cluster), 0))
  cen <- t(vapply(ids, function(g) colMeans(coords[cluster == g, , drop = FALSE]),
                  numeric(ncol(coords))))
  d <- vapply(seq_along(ids), function(i) sqrt(rowSums(sweep(coords, 2, cen[i, ])^2)),
              numeric(nrow(coords)))
  sil <- rep(NA_real_, nrow(coords))
  for (i in seq_along(ids)) {
    m <- which(cluster == ids[i]); a <- d[m, i]; b <- apply(d[m, -i, drop = FALSE], 1, min)
    sil[m] <- (b - a) / pmax(a, b)
  }
  sil
}
cohesion <- function(emb, members) {           # = mean cosine of members to the centroid
  U <- emb[members, , drop = FALSE]; U <- U / sqrt(rowSums(U^2))
  sqrt(sum(colMeans(U)^2))
}
ids <- sort(setdiff(unique(cluster), 0))
sil <- simplified_silhouette(coords, cluster)
quality <- data.frame(
  topic = labels[as.character(ids)],
  n = as.integer(table(cluster)[as.character(ids)]),
  silhouette = round(vapply(ids, function(g) mean(sil[cluster == g]), numeric(1)), 2),
  cohesion = round(vapply(ids, function(g) cohesion(emb, which(cluster == g)), numeric(1)), 2)
)
quality[order(-quality$silhouette), ]
#>                                    topic   n silhouette cohesion
#> 3         Epidemiology and public health 150       0.76     0.37
#> 6             Cosmology and astrophysics 150       0.75     0.40
#> 2 Condensed matter and quantum materials 150       0.74     0.39
#> 4        Number theory and combinatorics 150       0.74     0.39
#> 7        Climate and atmospheric science 150       0.74     0.39
#> 8                          Deep learning 150       0.74     0.39
#> 5         Genomics and molecular biology 150       0.73     0.40
#> 1              Cryptography and security 150       0.72     0.40

On real decoder-LLM embeddings the two metrics tell different stories: silhouette varies from topic to topic (weak, overlapping clusters score lower — honestly), while cohesion is uniformly high because these embeddings are anisotropic (every document points into a narrow cone of the space). That gap is exactly why both are reported: the map can look cleaner than the raw geometry is.

What is each topic about?

Naming a cluster with the model is one summary; its distinctive words are another, and they are checkable against the abstracts. Ranking by raw frequency would just return corpus-wide common words. Instead we use the log-odds-ratio with an informative Dirichlet prior (Monroe, Colaresi & Quinn 2008): each word’s over-representation in a topic versus the rest of the corpus, standardised to a z-score, with the prior shrinking common words toward zero. Twenty lines of base R — no Python, no topic-model package.

stopwords <- c("the","a","an","and","or","to","of","in","on","for","with","we","our",
               "that","this","is","are","be","by","as","from","using","study","results",
               "method","approach","paper","show","these","which","can","based","present")
tokenize <- function(t) {
  w <- unlist(regmatches(tolower(t), gregexpr("[[:alpha:]]+", tolower(t))))
  w <- w[nchar(w) >= 3]; w[!w %in% stopwords]
}
toks <- lapply(abstracts$text, tokenize)
vocab <- sort(unique(unlist(toks[cluster %in% ids])))
counts <- vapply(ids, function(g)
  as.integer(table(factor(unlist(toks[cluster == g]), levels = vocab))),
  integer(length(vocab)))
rownames(counts) <- vocab                        # vocabulary x topic counts
yw <- rowSums(counts); alpha <- length(vocab) * yw / sum(yw)  # informative prior
a0 <- sum(alpha); N <- sum(counts)
top_terms <- function(gi, n = 8) {
  y <- counts[, gi]; ni <- sum(y); r <- yw - y
  delta <- log((y + alpha) / (ni + a0 - y - alpha)) -
           log((r + alpha) / (N - ni + a0 - r - alpha))
  z <- delta / sqrt(1 / (y + alpha) + 1 / (r + alpha))
  names(sort(z, decreasing = TRUE))[seq_len(n)]
}
data.frame(
  topic = labels[as.character(ids)],
  distinctive_terms = vapply(seq_along(ids),
    function(gi) paste(top_terms(gi), collapse = ", "), character(1))
)
#>                                    topic
#> 1              Cryptography and security
#> 2 Condensed matter and quantum materials
#> 3         Epidemiology and public health
#> 4        Number theory and combinatorics
#> 5         Genomics and molecular biology
#> 6             Cosmology and astrophysics
#> 7        Climate and atmospheric science
#> 8                          Deep learning
#>                                                                   distinctive_terms
#> 1     surface, knowledge, error, quantum, previously, parameter, security, standard
#> 2              quantum, high, efficient, power, states, temperature, supports, spin
#> 3 case, health, effectiveness, disease, screening, supports, programmes, identifies
#> 4                           theory, number, coding, group, wide, known, under, mild
#> 5            disease, type, precision, identity, resolves, regulatory, cell, burden
#> 6          structure, dark, functions, mass, constant, measurements, mild, spectrum
#> 7   surface, temperature, forecasting, climate, strong, forcing, heterogeneity, sea
#> 8    models, gap, networks, throughput, distribution, efficiency, reduces, training

The terms discriminate the topics rather than echo the corpus — and they let a reader judge each cluster for themselves.

How do the topics relate?

Finally, which topics are near-duplicates and which stand apart? The cosine similarity between topic centroids answers it; hclust() on 1 - cosine turns it into a dendrogram.

centroids <- t(vapply(ids, function(g) colMeans(emb[cluster == g, , drop = FALSE]),
                      numeric(ncol(emb))))
centroids <- centroids / sqrt(rowSums(centroids^2))
cosm <- pmin(pmax(tcrossprod(centroids), -1), 1)
short <- function(s) ifelse(nchar(s) > 22, paste0(substr(s, 1, 19), "..."), s)
hc <- hclust(as.dist(1 - cosm), method = "average")
hc$labels <- short(labels[as.character(ids)])
op <- par(mar = c(2, 4, 3, 9), mgp = c(2.5, 0.7, 0))
plot(as.dendrogram(hc), horiz = TRUE, xlab = "1 - cosine",
     main = "Which topics are near each other")

Dendrogram of topics clustered by centroid cosine similarity.

par(op)

Because decoder-LLM embeddings are anisotropic, every centroid pair sits at a high absolute cosine; the relative structure is what matters, and the demo script draws it as a heatmap on a data-driven scale. The full runnable version — the polished map, the quality panel, per-topic term bars, and the centroid-cosine heatmap — is run_demo_B(extended = TRUE) in tests/demos/demo-B-topics.R (with fetch-abstracts.R for the real ~5,000-abstract arXiv corpus). No Python was involved at any step.