Held-out moment-based specification tests

Answer the question: which parts of the vocabulary does the model still get systematically wrong? A fit index says how much error remains, not where. This function groups the vocabulary into researcher-specified, labeled word sets (the instruments: for example frequency strata, or strata of a training fit score) and tests whether the model consistently over- or under-predicts any group on documents it has never seen, following Section 4 of Lewis and Grossetti (2026).

Usage

optop_moment_test(
  models,
  dtm_eval,
  dtm_train,
  type = c("contrast", "strata", "fit"),
  bins = 5,
  strata = 5,
  min_doc_freq = 5,
  adjust = c("none", "bonferroni", "BH"),
  baseline = NULL,
  n_threads = 1L,
  ...
)

Arguments

models

A list of topic models fitted on the training corpus (any fold-in-capable engine; see optop_index_holdout()).

dtm_eval

A counts document-term matrix of the evaluation documents, or an optop_corpus() of evaluation shards; columns are matched to the training vocabulary by name.

dtm_train

The counts document-term matrix of the training corpus, or an optop_corpus() of training shards (frequencies and document frequencies are pooled), used only to build the instruments (word frequencies and, for type = "fit", the training word-level fit score); never touched by the evaluation residuals.

type

Instrument family:

"contrast" (Test 1): one row contrasting the top against the bottom training-frequency quintile; the scalar t statistic applies.
"strata" (Test 2): bins - 1 rows comparing each training-frequency stratum against the highest-frequency reference stratum.
"fit" (Test 3): strata - 1 rows comparing strata of the training word-level deviance index against the best-fit reference stratum; the instrument matrix depends on \(K\).

bins

Number of frequency strata for type = "strata" (default 5).

strata

Number of fit strata for type = "fit" (default 5).

min_doc_freq

Minimum training document frequency for a word to enter the fit strata (default 5); excluded words carry zero instrument entries, which leaves every row sum at zero.

adjust

Multiple-testing adjustment of the marginal per-stratum t tests: "none" (default), "bonferroni", or "BH".

baseline

Optional training baseline for the fit score (default NULL computes it from dtm_train).

n_threads

Integer; OpenMP threads of the compiled kernels used for the training fit score.

...

Passed to the engine's fold-in routine.

Value

An object of class optop_moment_test: a list with

summary: data.table with one row per \(K\). Reading the columns: q is the number of instrument rows (word-group contrasts); wald the joint test statistic with df = q and pval its p-value, where a small value indicates systematic residual bias along at least one contrast; t and pval_t are the equivalent scalar test when q = 1 (NA otherwise).
moments: per-\(K\) list with g_bar (the mean moments, effect sizes in probability-mass units: how much observed mass exceeds fitted mass along each contrast), sigma (their covariance), n, and marginal (data.table localizing a rejection: one row per stratum with estimate, se, t, pval, and the adjusted pval_adj).
instruments: the instrument matrix and the stratum assignment of every vocabulary word (per \(K\) for type = "fit").
type, adjust, K: the design.

Details

Formally, for each evaluation document the residual probability vector is \(\varepsilon_j = \mathbf d_j - \hat{\mathbf p}^{K,tr}_j\), observed minus fitted word probability under the training-fitted model; the moment vector projects it onto a training-built instrument matrix, \(\mathbf g_j = \mathbf Z\,\varepsilon_j\), and the test asks whether the moments average to zero across evaluation documents. Positive components mean the model under-predicts that direction of the vocabulary.

How it works.

Split upstream. Models are fitted on the training corpus; the evaluation documents stay out of the fitting entirely.
Build the instruments from training data only. Each instrument row contrasts one word group against a reference group, so its moment reads as "average residual mass in this group minus the reference". Because the groups are fixed before any evaluation residual is seen, the test is not fishing in its own evaluation data.
Score every evaluation document. Topic weights are folded in (holding the trained topics fixed), the residual vector is formed, and its projection on the instruments gives one small moment vector per document.
Test. If the model has no systematic bias along the chosen directions, the moments should average to zero; the Wald statistic measures how far from zero the averages are relative to their sampling noise.

Construction. All instruments are built from the training sample only, before any evaluation residual is examined, and every row is exactly mean-zero with \(\|z\|_1 = 2\). The residuals use the evaluation counts on the training vocabulary with document lengths restricted to the aligned tokens, so that observed and fitted vectors are both probability vectors on the training support and \(\mathbf 1^\top \varepsilon_j = 0\) holds identically; out-of-support tokens play no role in the residual moments.

Statistics. With n evaluation documents, the Wald statistic \(\mathcal W = n\,\bar{\mathbf g}^\top \hat\Sigma^{-1} \bar{\mathbf g}\) is asymptotically \(\chi^2_q\) under the null of no systematic residual bias (Proposition 3 of the paper), where \(\hat\Sigma\) is the sample covariance of the document moments. In the scalar case the equivalent t statistic is standard normal and \(\mathcal W = t^2\). Marginal per-stratum t tests localize a rejection, with the requested adjustment.

Interpretation. The null is a statement about the average held-out residual balance of the fitted model, conditional on the training sample; it is not equivalent to correct specification of the model class. Estimation error and the reuse of the evaluation document in folding in its topic weights generically make the moments small but nonzero, so with many evaluation documents the test detects any residual imbalance, whatever its source. Read rejections jointly with the magnitude of the mean moments, which are effect sizes in probability-mass units. When many instruments, topic counts, or families are examined, treat the p-values as exploratory or adjust for multiple testing.

References

Lewis, C. M. and Grossetti, F. (2026). Goodness-of-fit indices and diagnostics for topic models. Working paper.

Hansen, L. P. (1982). Large sample properties of generalized method of moments estimators. Econometrica, 50(4), 1029–1054.

Examples

# \donttest{
# a small synthetic corpus, split into training and evaluation
rdir <- function(n, k, a) {
  g <- matrix(stats::rgamma(n * k, shape = a), nrow = n)
  g / rowSums(g)
}
set.seed(42)
corpus <- sim_dfm(DTW = rdir(120, 6, 0.4), TWW = rdir(6, 300, 0.1),
                  doc_length = rep(400, 120), seed = 1)
dtm_train <- corpus[1:90, ]
dtm_eval <- corpus[91:120, ]
models <- lapply(c(4, 6, 8), function(k) {
  topicmodels::LDA(dtm_train, k = k, control = list(seed = 100 + k))
})

# does residual bias vary across frequency strata?
mt <- optop_moment_test(models, dtm_eval, dtm_train, type = "strata",
                        bins = 4)
mt$summary                  # one Wald test per K
#>        K     q     wald    df         pval     t pval_t
#>    <num> <int>    <num> <int>        <num> <num>  <num>
#> 1:     4     3 23.33836     3 3.432836e-05    NA     NA
#> 2:     6     3 26.47691     3 7.577859e-06    NA     NA
#> 3:     8     3 27.05612     3 5.729999e-06    NA     NA
mt$moments[["6"]]$marginal  # which stratum drives it, at K = 6
#>        stratum     estimate           se        t         pval     pval_adj
#>         <char>        <num>        <num>    <num>        <num>        <num>
#> 1: freq_1_vs_4 0.0003064067 7.082182e-05 4.326445 1.515349e-05 1.515349e-05
#> 2: freq_2_vs_4 0.0002919390 8.655880e-05 3.372725 7.442835e-04 7.442835e-04
#> 3: freq_3_vs_4 0.0004236351 1.310508e-04 3.232602 1.226684e-03 1.226684e-03
# }