Select the optimal number of topics from a pool of topic models

Identify the number of topics that best describes the corpus with the Test 1 statistic of Lewis and Grossetti (2022), a Pearson chi-square goodness-of-fit test on each document's word distribution. The routine evaluates each model of the grid, optionally calibrates the p-values under the fitted-model null, and selects the optimal topic count with one of two rules; see Details.

Usage

optop_select(
  topic_models,
  weighted_dfm,
  q = 0.95,
  alpha = 0.05,
  selection = c("sequential", "min"),
  calibrate = c("none", "bootstrap", "moment"),
  n_boot = 200L,
  doc_lengths = NULL,
  seed = NULL,
  n_threads = 1L,
  do_plot = TRUE,
  verbose = TRUE
)

optimal_topic(
  topic_models,
  weighted_dfm,
  q = 0.95,
  alpha = 0.05,
  selection = c("sequential", "min"),
  calibrate = c("none", "bootstrap", "moment"),
  n_boot = 200L,
  doc_lengths = NULL,
  seed = NULL,
  n_threads = 1L,
  do_plot = TRUE,
  verbose = TRUE
)

Arguments

topic_models: A list of fitted topic models spanning the candidate values of \(K\), one model per \(K\) (an unordered grid is sorted by increasing \(K\)). Supported classes: topicmodels::LDA() fits (VEM or Gibbs), topicmodels::CTM() fits, the seededlda models (seededlda::textmodel_lda(), seededlda::textmodel_seededlda(), seededlda::textmodel_seqlda()) and NLPstudio fits (nlp_topic_fit). Classes can be mixed within a grid as long as every model was fitted on the same corpus and vocabulary; the test consumes the models only through their fitted word probabilities (see Details). Elements may also be loader functions returning a fit: each is materialized on demand and released, so very large grids never sit in memory at once (extracted weights are still cached while the grid fits a budget of 256 MB by default, configurable via options(optop.cache_mb = )).
weighted_dfm: A weighted quanteda::dfm containing word proportions for each document, built from the same counts dfm the models were fitted on; it is recommended that document ids are available via quanteda::docid(). May also be an optop_corpus() of proportion shards for corpora past the size of a single dfm: shards are evaluated one at a time and combine exactly (a one-shard corpus reproduces the plain call bit for bit, calibrated p-values included; multi-shard runs agree up to floating-point summation order). With a corpus, doc_lengths must be named.
q: Numeric in \((0, 1]\). Cumulative probability mass retained as "relatively important" words; the remaining words are collapsed into a single bin whose mass stays strictly below \(1 - q\). Equals \(1 - I^K\) in the paper's notation; the default 0.95 matches the paper's numerical setting \(I^K = 0.05\).
alpha: Numeric in \([0, 1]\). Significance level used by the "sequential" selection rule (default 0.05); ignored by "min".
selection: Character; how the optimal K is chosen (see Details): "sequential" (default) or "min".
calibrate: Character; how the p-values are computed (see the Calibration part of Details): "none" (default, the asymptotic chi-square of Equation 8), "bootstrap" (parametric bootstrap under the fitted-model null), or "moment" (exact multinomial moments matched to a scaled chi-square).
n_boot: Integer; number of bootstrap replicates per model when calibrate = "bootstrap" (default 200L, minimum 20).
doc_lengths: Numeric vector with the length (total token count) of each document; required when calibrate != "none", typically quanteda::ntoken() of the counts dfm the models were fitted on. If named, names are matched to the document ids of weighted_dfm; if unnamed, the vector must follow the row order of weighted_dfm.
seed: Optional integer passed to set.seed() before the bootstrap, for reproducible calibrated p-values (default NULL).
n_threads: Integer; number of OpenMP threads used by the compiled cores (the per-document statistic loop and the bootstrap), default 1L. Results are bit-identical for any value: every document owns a deterministic RNG stream and all reductions run in a fixed order, so n_threads only affects wall time. On builds without OpenMP (e.g. the default macOS toolchain) the value is ignored and the cores run single-threaded. For thread and BLAS configuration on server machines (including the mclapply oversubscription trap), see the package website article https://contefranz.github.io/OpTop/articles/very-large-corpora.html.
do_plot: Logical; if TRUE, plot the standardized statistic versus topics with vertical and horizontal guides at the selected optimum and a subtitle reporting the selection method, the selected \(K\), and the calibration if any (default TRUE).
verbose: Logical; if TRUE (default), report progress (a cli progress bar across the model grid) and a selection summary. Regardless of verbose, corrective actions (dropping documents the models never saw, reordering models or features, the sequential rule's fallback) are always signalled.

Value

A data.table with columns:

topic: integer number of topics (\(K\)).
OpTop: standardized Test 1 statistic (raw statistic divided by df).
df: degrees of freedom \(\sum_j P_j\) of the raw statistic.
pval: the p-value the selection rules use: asymptotic upper tail for calibrate = "none", calibrated otherwise.
pval_chisq: only when calibrate != "none": the uncalibrated asymptotic chi-square p-value, for comparison.

Details

optop_select() was called optimal_topic() before OpTop 0.19.0; the old name keeps working as an alias that notes the rename once per session and will not be removed before a 2.0. New code should use optop_select(), which follows the package-wide optop_ naming.

For each document \(j\), the fitted word probabilities are sorted in decreasing order and the smallest head whose cumulative mass exceeds \(q\) is kept (\(P_j\) words); the remaining words are collapsed into a single "min" bin. The document contributes \(P_j + 1\) times its Pearson term over the \(P_j + 1\) bins, and the corpus statistic is the sum over documents (Equation 8 of the paper), asymptotically chi-square with \(\sum_j P_j\) degrees of freedom. The returned OpTop column reports the standardized statistic (raw statistic divided by its degrees of freedom, the version plotted in the paper's Figure 2), while pval is the p-value the selection rules consume: the upper tail of the raw statistic on its full degrees of freedom by default, or the calibrated value when calibrate != "none".

The statistic touches each model only through its fitted word probabilities \(I_j = \theta_j^\top \Phi\), where \(\theta_j\) is the document's topic distribution and \(\Phi\) the topic-word matrix. Any implementation that provides those two (row-stochastic) matrices is therefore admissible, whatever its estimation method; that is what the supported classes listed under topic_models have in common.

Selection rules.

"sequential" (default): scan \(K\) upward and select the smallest \(K\) the test fails to reject (pval > alpha), the classical sequential scheme for model order. If every model is rejected, the rule falls back to the global minimum with a warning.
"min": select the \(K\) with the minimum standardized statistic, the rule used in the published case study.

The pre-0.9.9 "legacy" rule, deprecated since 0.9.9, was removed in 0.16.0.

Calibration. The chi-square reference of Equation 8 is a yardstick rather than an exact null law, for two reasons. First, the classical Pearson asymptotics hold for the count statistic, whose scale factor is the document length \(N_j\); Equation 8 works on proportions scaled by the bin count \(P_j + 1\) instead, so the statistic's null magnitude is off by roughly \(N_j / (P_j + 1)\) per document. Second, the expected probabilities are estimated from the same data they are tested against. The practical consequence is that with \(\sum_j P_j\) degrees of freedom in the thousands the chi-square quantiles are razor-thin and upper-tail p-values saturate at 0 or 1 unless the fit is genuinely borderline; alpha is not a true Type-I error rate.

Calibration replaces the chi-square reference with the distribution of the statistic under the conditional fitted-model null: document \(j\) is \(\mathrm{Multinomial}(N_j, I_j)\), where \(I_j\) are the \(K\)-model's fitted word probabilities and \(N_j\) the observed document lengths (hence doc_lengths). Two properties make this exact and fast:

the per-document envelope (sorted fitted probabilities, cutoff at \(q\), collapsed min bin) depends only on the model, never on the data, and the statistic touches the data only through sums over those fixed bins;
a multinomial collapsed over bins is multinomial on the collapsed probabilities, so the null can be simulated directly on the \(P_j + 1\) bins, exactly equivalent to simulating whole documents over the vocabulary, at a tiny fraction of the cost.

calibrate = "bootstrap" draws \(B\) = n_boot null replicates \(T^\ast\) of the statistic \(T\) this way and reports the empirical upper-tail p-value with the standard finite-sample correction, \((1 + \#\{T^\ast \ge T\}) / (B + 1)\). This is the reference method: alpha becomes a genuine Type-I error rate with respect to the conditional null, at bootstrap resolution \(1 / (B + 1)\).

calibrate = "moment" uses the exact multinomial moments of the per-document Pearson term (Haldane, 1937): over \(k\) bins with probabilities \(p_b\) and length \(n\), the count statistic \(X^2\) has \(E[X^2] = k - 1\) and \(\mathrm{Var}[X^2] = 2(k - 1) + (\sum_b 1/p_b - k^2 - 2k + 2)/n\); scaling by the statistic's \(k/n\) factor and summing over independent documents gives the null mean \(\mu\) and variance \(\sigma^2\) of the corpus statistic, which are matched to a scaled chi-square \(a\,\chi^2_\nu\) (Satterthwaite: \(a = \sigma^2 / (2\mu)\), \(\nu = 2\mu^2 / \sigma^2\)). Closed form, no simulation: a fast approximation that corrects the location and scale of the reference but not its higher moments.

One caveat applies to both: the null holds the estimated \(\theta\) and \(\phi\) fixed (no per-replicate re-fitting of the model; the "double bootstrap" would be exact but computationally prohibitive). Calibrated p-values are therefore conditional on the fitted parameters and do not account for estimation noise in \(\theta\) and \(\phi\).

Input alignment. weighted_dfm must be a quanteda::dfm of word proportions (row-wise), built from the same counts dfm the models were fitted on. Its features must match the models' common vocabulary: when they are the same set in a different order the columns are reordered automatically (signalled), any other mismatch is an error. Document identifiers are taken from quanteda::docid(weighted_dfm) and matched against each model's own document identifiers: documents not present in every model are dropped (with a warning), and each retained dfm row is paired with the corresponding row of \(\theta\) per model. Neither the row order of weighted_dfm nor the order in which each model saw the documents matters: alignment is always by identifier.

Performance note. Both hot paths run in C++ compiled code and parallelize over documents with OpenMP (see n_threads): the statistic evaluates each block of documents with one BLAS product and sorts the fitted probabilities per document, and the bootstrap draws the null replicates directly on the collapsed bins, fused with the Pearson reduction, at about \(B \times \mathrm{df}\) floating-point operations per model. Extracted model weights are cached across the evaluation while the grid fits options(optop.cache_mb = ) (default 256): raising it on large servers keeps big grids resident, and results never depend on the budget. See the "Computational efficiency" section of the vignette for measurements.

References

Lewis, C. M. and Grossetti, F. (2022). A statistical approach for optimal topic model identification. Journal of Machine Learning Research, 23(58), 1–20. https://jmlr.org/papers/v23/19-297.html

Haldane, J. B. S. (1937). The exact value of the moments of the distribution of chi-square. Biometrika, 29, 133–143.

Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2(6), 110–114.

Davison, A. C. and Hinkley, D. V. (1997). Bootstrap Methods and their Application. Cambridge University Press, Cambridge.

Examples

# simulate a corpus whose truth is known, then let Test 1 recover it:
# the candidate grid holds two random models and the true K = 3 model
rdirich <- function(n, k) {
  g <- matrix(stats::rgamma(n * k, shape = 1), n, k)
  g / rowSums(g)
}
set.seed(1)
theta <- rdirich(40, 3)
rownames(theta) <- sprintf("d%02d", 1:40)
phi <- rdirich(3, 60)
colnames(phi) <- sprintf("w%02d", 1:60)
counts <- sim_dfm(theta, phi, doc_length = 400, seed = 2)
wprop <- quanteda::dfm_weight(counts, scheme = "prop")

models <- lapply(c(2, 3, 4), function(k) {
  if (k == 3) return(optop_model(theta, phi))
  th <- rdirich(40, k); rownames(th) <- rownames(theta)
  ph <- rdirich(k, 60); colnames(ph) <- colnames(phi)
  optop_model(th, ph)
})

test1 <- optop_select(models, wprop, do_plot = FALSE)
#> 
#> ── Optimal topic selection ──
#> 
#> ℹ Evaluating 3 models on 40 documents and 60 features (q = 0.95, alpha = 0.05, selection = sequential)
#> ✔ Optimal model has 3 topics (selected by sequential adequacy scan at alpha = 0.05)
#> ℹ Completed in 0.02s
test1
#>    topic     OpTop    df         pval
#>    <num>     <num> <num>        <num>
#> 1:     2 1.6555559  1845 4.324752e-63
#> 2:     3 0.1317413  2001 1.000000e+00
#> 3:     4 1.0828900  2110 4.249976e-03

# bootstrap-calibrated p-values: document lengths come from the counts
test1_cal <- optop_select(models, wprop,
                          calibrate = "bootstrap", n_boot = 50,
                          doc_lengths = quanteda::ntoken(counts),
                          seed = 42, do_plot = FALSE)
#> 
#> ── Optimal topic selection ──
#> 
#> ℹ Evaluating 3 models on 40 documents and 60 features (q = 0.95, alpha = 0.05, selection = sequential)
#> ℹ Calibrating p-values by parametric bootstrap (B = 50, seed = 42) under the fitted-model null
#> ✔ Optimal model has 3 topics (selected by sequential adequacy scan at alpha = 0.05, bootstrap-calibrated)
#> ℹ Completed in 0.07s

# \donttest{
# the same flow with real engine fits (topicmodels VEM)
fits <- lapply(2:4, function(k) {
  topicmodels::LDA(quanteda::convert(counts, to = "topicmodels"),
                   k = k, method = "VEM", control = list(seed = 500 + k))
})
optop_select(fits, wprop, do_plot = FALSE)
#> 
#> ── Optimal topic selection ──
#> 
#> ℹ Evaluating 3 models on 40 documents and 60 features (q = 0.95, alpha = 0.05, selection = sequential)
#> ✔ Optimal model has 2 topics (selected by sequential adequacy scan at alpha = 0.05)
#> ℹ Completed in 0.02s
#>    topic     OpTop    df  pval
#>    <num>     <num> <num> <num>
#> 1:     2 0.1686357  2029     1
#> 2:     3 0.1224388  1998     1
#> 3:     4 0.2431973  2040     1
# }