Moment-Based Inference for Regression with Latent Dirichlet Covariates

Moment-Based Inference for Regression with Latent Dirichlet Covariates

Introduction and Literature Review

Topic models are widely used as preprocessing tools for empirical regression analysis. A corpus of high-dimensional text is first mapped into a low-dimensional representation, typically estimated document-level topic shares, and these topic-derived variables are then used as regressors, controls, or predictors in a downstream model. Examples include newspaper topics in regressions for political violence MuellerRauh2018, news-topic predictors in macroeconomic and financial forecasting regressions LarsenThorsrud2019,Thorsrud2020, topic-based measures of business news in vector autoregressions and stock-market timing models BybeeKellyManelaXiu2024, and topic measures constructed from central-bank transcripts in econometric studies of transparency and deliberation HansenMcMahonPrat2018. In this role, latent Dirichlet allocation and related topic models are best viewed as dimension-reduction steps: they convert sparse word-count data into a small number of latent coordinates intended for subsequent statistical use.

This perspective connects topic-model preprocessing to the use of estimated factors in factor-augmented regressions, where a high-dimensional panel is compressed into a lower-dimensional latent representation before a second-stage analysis StockWatson2002,Bai2003,BaiNg2006. The analogy is useful because it highlights the inferential issue. Once an estimated latent representation is used as a regressor, uncertainty from the first step must be reconciled with inference in the second step. In the empirical topic-modeling workflow, this reconciliation is often absent: researchers estimate topics, construct document-level topic shares, and then condition on those estimated shares as if they were ordinary observed covariates.

We consider a canonical finite-LDA version of this problem. Let document \(i\) contain \(N\) word tokens \(x_{i1},\ldots,x_{iN}\), let \(k\) denote the number of topics, let \(h_i\) denote its latent topic mixture, and let \(O\) denote the topic matrix. Conditional on \(h_i\), the word tokens are generated from the mixture \(Oh_i\). Suppose a scalar response can be written as \[ Y_i=\beta^\top h_i+\varepsilon_i, \] where the residual is orthogonal to the low-order token moments used for identification; a sufficient model-level condition is \[ \mathbb{E}(\varepsilon_i\mid h_i,x_{i1},\ldots,x_{iN})=0 . \] The object of inference is the downstream coefficient \(\beta\). We treat the topic matrix, Dirichlet parameters and regression coefficient as fixed unknown population parameters; the Dirichlet law for \(h_i\) is a random-effects distribution for document heterogeneity, not a prior for a posterior analysis over these unknown parameters. The usual plug-in workflow estimates topic shares \(\hat h_i\) and then regresses \(Y_i\) on \(\hat h_i\). This creates two difficulties.

The first is a first-stage regularity problem. Let \(n\) denote the number of independent documents. For classical downstream inference, the relevant object is not only the first-stage topic estimate but the entire map from the observed corpus to the final regression coefficient. A standard delta-method argument requires this map to be regular in the sense that \[ \sqrt n(\hat\beta-\beta) = n^{-1/2}\sum_{i=1}^n \phi_i+o_p(1), \] for an influence function \(\phi_i\) with consistently estimable variance. Many topic-modeling procedures used in applied work, including Bayesian and approximate-Bayesian implementations of LDA, produce posterior or approximate-posterior summaries under chosen prior, hyperparameter and computational specifications. Such summaries may be useful for exploration, prediction or posterior inference under the specified model, but they do not by themselves provide frequentist confidence intervals for the downstream coefficient \(\beta\). Frequentist validity of the resulting two-stage map must be established separately. In standard plug-in workflows, the second-stage regression instead conditions on estimated topic shares as if they were observed covariates. The resulting standard errors then use only the second-stage regression variation and omit the first-stage component of the influence function; unless that omitted component is asymptotically negligible, the confidence intervals are not valid for the sampling distribution of the full two-step procedure.

The second difficulty is more fundamental. Even if the population topic matrix \(O\) were known, a document-specific topic mixture \(h_i\) cannot be consistently recovered from the words in document \(i\) when document length is fixed. The number of documents may grow, but the information about any particular document's latent mixture remains bounded. Hence the plug-in regressor \(\hat h_i\) contains non-vanishing document-level error. This is the generated-regressor problem emphasized by BattagliaChristensenHansenSacher2024 for regressions using variables extracted from unstructured data. In the present setting it arises even under perfect recovery of the topic matrix: increasing the number of documents improves estimation of population moments, but it does not make each individual \(h_i\) observed. A regression on \(\hat h_i\) can therefore have a wrong probability limit or a persistent finite-document bias, not merely standard errors that condition on a noisy generated regressor.

Integrated latent-variable likelihood procedures or fully Bayesian joint models could also target \(\beta\) under a more complete specification by integrating over the latent document mixtures. The contribution here is complementary: we give a low-order moment route to frequentist confidence intervals for \(\beta\), based on the supervised moment restrictions stated below, without specifying a full likelihood for the joint distribution of responses and words or a prior over the unknown population parameters. The comparison with plug-in regressions is therefore aimed at the common two-step topic-share regression workflow rather than at exhausting all possible latent-variable estimators.

The spectral approach to LDA provides a useful starting point for resolving both difficulties. AnandkumarGeHsuKakadeLiu2012 show that, when the total Dirichlet concentration parameter is known, the topic matrix can be recovered from low-order cross-token word moments after applying Dirichlet-specific corrections. These corrections subtract lower-order terms induced by the common document-level topic mixture, leaving a population object that is diagonal in latent topic coordinates. WangZhu2014 extend this corrected-moment logic to supervised LDA under a related latent-mixture regression structure, introducing response-weighted corrected moments that recover the regression weights. RenWangZhu2018 develop this supervised spectral approach further, giving both a two-stage method that recovers LDA parameters before recovering the regression model and a single-phase method that jointly recovers topic and regression parameters.

We take this supervised spectral moment structure as the starting point. Our supervised identity uses the same response-weighted moment logic, but we formulate it as a conceptually close and operationally different observed-space operator. The response-weighted correction is formed in observed word space, and after projection onto the recovered topic basis the representation \(O^+H^y_{\alpha_0}O\) is diagonal with entries proportional to the components of \(\beta\). The distinction is therefore not the underlying supervised-LDA moment model, but the inferential target and estimation problem: we study frequentist inference for a downstream regression coefficient at fixed document length, without estimating document-level topic shares.

The remaining roadblock, common to the spectral LDA and supervised spectral LDA procedures just discussed, is that the corrected moments are indexed by the total Dirichlet concentration parameter. Write this concentration as \(\alpha_0\). In these procedures, \(\alpha_0\) is supplied as an input when forming the corrected moments. In empirical work this quantity is rarely known, and treating it as a tuning parameter is unsatisfactory for two reasons. First, the downstream procedure then rests on a choice of \(\alpha_0\) rather than on an identification argument for \(\alpha_0\) from the observable word distribution. Practical selection rules are therefore difficult to interpret as estimators of a population parameter. Second, even if a particular rule performs well computationally, its uncertainty has no automatic route into the final confidence interval for \(\beta\). We show that, when the number of topics satisfies \(k\ge3\) and a generic finite-probe condition holds, \(\alpha_0\) is itself identified from observable moments: at the true concentration value, a family of corrected word-moment operators commute, whereas away from the truth they generically fail to commute. This gives a feasible estimator of \(\alpha_0\) and allows its uncertainty to be propagated through the estimator of the downstream coefficient.

We work in a fixed-document-length asymptotic regime: the number of independent documents grows, while the number of words in each document may remain bounded. In this regime, the target is inference on the population coefficient \(\beta\), not recovery of every document-specific topic mixture. Within this setting, the paper makes three contributions.

First, for known total concentration, we give an observed-space operator formulation of the supervised spectral moment identity. This identity starts from the same response-weighted corrected moments as the supervised spectral LDA literature, but casts them in a form tailored to downstream inference: the supervised operator is an observed-space matrix whose topic-basis representation \(O^+H^y_{\alpha_0}O\) is diagonal, and those diagonal entries are proportional to the components of \(\beta\). This identifies the downstream coefficient directly from observable word--response moments, without constructing document-level topic-share regressors.

Second, we identify and estimate the total Dirichlet concentration parameter. The spectral LDA and supervised spectral LDA procedures discussed above form their corrected moments conditional on \(\alpha_0\). We show that, when the number of topics is at least three and a generic probe condition holds, \(\alpha_0\) is the unique value for which a finite family of corrected word-moment operators commute. This gives a feasible estimator of \(\alpha_0\) from the same low-order word moments used for topic recovery.

Third, relative to the supervised spectral recovery results in WangZhu2014 and RenWangZhu2018, we establish frequentist inference for the downstream coefficient in a fixed-dimensional regime. The primitive empirical quantities are distinct-token word moments and response-weighted word moments, averaged across independent documents. With fixed vocabulary dimension \(d\), fixed topic dimension \(k\), and fixed document length \(N\), or after conditioning on a fixed admissible compression dimension, these moments satisfy an ordinary root-\(n\) central limit theorem over documents. The concentration-parameter, spectral, and supervised coefficient maps are smooth under rank and separation conditions. The resulting estimator of \(\beta\) is asymptotically linear with a feasible sandwich variance estimator, so the contribution is a first-order sampling expansion and feasible standard errors rather than only consistency or sample-complexity recovery.

The Monte Carlo experiments are designed to match this inferential logic. They evaluate the finite-sample performance of the proposed confidence intervals and compare them with plug-in regressions using estimated document-topic shares. The proposed intervals have coverage close to nominal in the designs considered. Plug-in regressions undercover substantially, including a known-\(O\) plug-in reconstruction procedure in which the topic matrix is supplied but document-level mixtures are still reconstructed from finite document counts. This isolates the fixed-document generated-regressor problem: increasing the number of documents improves population moment estimation, but it does not make the document-level topic mixtures observed. The simulations are intended to isolate this mechanism under favorable rank, separation and probe conditions, not to exhaust all finite-sample stress tests of the feasible workflow.

The rest of the paper is organized as follows. Section (ref) introduces the finite-LDA model, the downstream regression target, and the observable word and word--response moments. Section (ref) derives the corrected moment operators, proves population identification of the topic matrix and the downstream coefficient, and establishes identification of \(\alpha_0\) through commutativity. Section (ref) defines the proposed moment-based estimators and gives the asymptotic distribution and sandwich variance estimator for the downstream coefficient. Section (ref) reports Monte Carlo evidence on finite-sample performance and on the failure of plug-in topic regressions at fixed document length. Section (ref) gives a real-data illustration using articles from the five general-interest economics journals. Section (ref) discusses limitations and extensions. Proofs and extensions, including fixed-dimensional observed controls and split-sample linear compression for high-dimensional preprocessing, are collected in the Supplement.

Model and target parameter

This section specifies the finite-LDA sampling experiment and the downstream regression target. The point of the formulation is to separate the latent document-level topic mixture, which is useful for defining the model, from the population regression coefficient, which is the object of inference. The corrected spectral moments used for identification are introduced in Section (ref).

Observed data and latent topic structure

We observe independent documents \[ (Y_i,x_{i1},\ldots,x_{iN}),\qquad i=1,\ldots,n, \] where \(Y_i\in\mathbb{R}\) is a scalar document-level response. The vocabulary has \(d\) terms, and each word token is represented by a one-hot vector \(x_{ij}\in\{e_1,\ldots,e_d\}\subset\mathbb{R}^d\). The main theory is stated for a common document length \(N\ge3\), fixed as \(n\to\infty\). The condition \(N\ge3\) is used because the spectral construction relies on third-order cross-token moments. Variable document lengths \(N_i\ge3\) can be handled by normalizing the within-document moment averages document by document, as recorded in Appendix (ref); we maintain the same sampled-document LDA moment structure throughout.

Let \(k\) denote the number of topics, and let \[ O=[O_{:1},\ldots,O_{:k}]\in\mathbb{R}^{d\times k} \] denote the topic matrix. Each column \(O_{:\ell}\) lies in the vocabulary simplex \(\Delta^{d-1}\) and gives the word distribution for topic \(\ell\). The matrix \(O\) and the Dirichlet parameter below are treated as fixed unknown population parameters.

For document \(i\), the latent topic mixture is \[ h_i\sim\operatorname{Dirichlet}(\alpha), \qquad \alpha=(\alpha_1,\ldots,\alpha_k)^\top\in(0,\infty)^k, \] with total concentration

The usual LDA sampling scheme draws latent topic labels \[ z_{ij}\mid h_i \sim \operatorname{categorical}(h_i), \qquad j=1,\ldots,N, \] independently across token positions, and then draws words according to \[ x_{ij}\mid z_{ij}=\ell \sim \operatorname{categorical}(O_{:\ell}) . \] The topic labels \(z_{ij}\) are not observed and will not be used directly. Integrating them out gives the equivalent conditional word distribution

Thus, conditional on \(h_i\), the word tokens are independent and identically distributed with conditional mean

Marginally, tokens in the same document are dependent through the common latent mixture \(h_i\), which is the dependence exploited by the cross-token moment identities below.

The document-specific vector \(h_i\) is not observed. Because \(N\) is fixed, the analysis does not require, and does not assume, that \(h_i\) can be consistently recovered for each document. The inferential target below is instead a finite-dimensional population parameter identified from the joint law of \((Y_i,x_{i1},\ldots,x_{iN})\).

Downstream response model

The response is linked to the latent topic mixture by

for some \(\beta\in\mathbb{R}^k\), with \(\mathbb{E}[Y_i^2]<\infty\). The primitive restriction used by the supervised moment identities is low-order response-token orthogonality over the token positions entering the averaged moments:

A sufficient model-level condition is

or equivalently \(\mathbb{E}(Y_i\mid h_i,x_{i1},\ldots,x_{iN})=\beta^\top h_i\). The weaker condition \(\mathbb{E}(Y_i\mid h_i)=\beta^\top h_i\) alone is not enough for the response-weighted word moments used below if residual variation is systematically related to realized token choices beyond \(h_i\).

It should be noted that, since \(\mathbf{1}^\top h_i=1\), a separate intercept is not identified from a common shift of all topic coefficients. If a model is written as \(c+\tilde\beta^\top h_i\), then it is observationally equivalent to \((\tilde\beta+c\mathbf{1})^\top h_i\). We therefore adopt the normalization that absorbs the intercept into \(\beta\). This is the same normalization issue that arises in a regression with a full set of category indicators: the common level is a convention, while contrasts such as \(\beta_a-\beta_b\) are invariant and often the substantively meaningful objects. Topic labels are arbitrary, so \(\beta\) is always interpreted in the same ordering as the columns of \(O\). Appendix (ref) gives the corresponding fixed-dimensional extension with observed controls.

Observable low-order moments

All population moments below refer to a generic document, and we suppress the document index. By exchangeability, any distinct token positions have the same joint distribution, so we use positions \(1,2,3\) for notation. Write \[ \pi:=\frac{\alpha}{\alpha_0}, \qquad D:=\operatorname{diag}(\alpha). \] The first two Dirichlet moments are

and

The corresponding first and second cross-token word moments are

The use of distinct positions in \(M_2\) is important: \(M_2\) is a cross-token moment, not the second moment of a single multinomial draw.

For a contraction direction \(\eta\in\mathbb{R}^d\), define the flattened third cross-token moment

Let \(\gamma:=O^\top\eta\). Conditional independence and (ref) give

The contracted third Dirichlet moment is

Thus

The corrected moments in the next section subtract lower-order terms induced by the common document mixture; at \(\tau=\alpha_0\), this removes the rank-one term in (ref) and all lower-order third-moment terms in (ref), leaving only the diagonal factor proportional to \(\operatorname{diag}(\alpha\circ\gamma)\).

Finally, define the supervised observable moments

Under (ref) and (ref),

Consequently, whenever \(O\) has full column rank, there exists \(\eta_\beta\in\mathbb{R}^d\) such that

and the supervised moments satisfy the bridge identities

These identities are the population reason why the downstream coefficient can be recovered from observable word--response moments without first constructing document-level estimates of \(h_i\).

Population identification through corrected operators

This section gives the population identities on which the estimators are based. The central point is that the document-level mixture \(h_i\) need not be recovered document by document. Instead, low-order observable moments can be corrected so that the resulting operators are diagonal in the latent topic coordinates. One operator recovers the topic directions; its supervised analogue recovers the downstream coefficient vector. A second consequence of the same diagonalization is an identifying restriction for the unknown concentration mass \(\alpha_0\).

Throughout the section, expectations are under the model of Section (ref). We write \(\tau>0\) for a candidate value of the total concentration and reserve \(\alpha_0\) for the true value. Let \[ \mu=\mathbb{E}[x_1],\qquad M_2=\mathbb{E}[x_1x_2^\top],\qquad T(\eta)=\mathbb{E}\{x_1x_2^\top\langle x_3,\eta\rangle\}, \] where \(\eta\in\mathbb{R}^d\). Define the \(\tau\)-corrected second moment

and the \(\tau\)-corrected contracted third moment

The associated observed-space operator is

where \(B_\tau^+\) is the Moore--Penrose inverse. The operator is generally not symmetric, so eigenvectors below are right eigenvectors.

Corrected LDA moments

The following lemma is the algebraic core of the paper. It states that the correction in (ref)--(ref) removes the off-diagonal Dirichlet terms exactly at \(\tau=\alpha_0\). Let \[ D:=\operatorname{diag}(\alpha),\qquad C_2:=\alpha_0(\alpha_0+1),\qquad C_3:=\alpha_0(\alpha_0+1)(\alpha_0+2). \]

The proof is a direct calculation from the Dirichlet second and third moments and is given in the Supplement. The number of topics is also identified at this point: since \(\mathbb{E}[hh^\top]\) is positive definite and \(O\) has full column rank, \(\operatorname{rank}(M_2)=k\).

Topic directions and supervised coefficients when \texorpdfstring{\(\alpha_0\)}{alpha0} is known

Equation (ref) gives the usual spectral recovery of the topic directions, but in observed coordinates rather than whitened coordinates. If \(\eta\) is chosen so that the coordinates of \(O^\top\eta\) are distinct and nonzero, the right eigenvectors associated with the \(k\) nonzero simple eigenvalues of \(H_{\alpha_0}(\eta)\) are exactly the topic vectors.

The downstream coefficient is identified by a parallel supervised correction. Using the response moments from Section (ref), define for a candidate \(\tau\)

and

The key identity is obtained at \(\tau=\alpha_0\).

The theorem follows from the bridge identities \(m_y=\langle\eta_\beta,\mu\rangle\), \(v_y=M_2\eta_\beta\), and \(T^y=T(\eta_\beta)\) for any \(\eta_\beta\) satisfying \(O^\top\eta_\beta=\beta\). These identities are implied by the low-order orthogonality conditions in (ref); the stronger conditional mean restriction (ref) is only one sufficient primitive condition. Thus \(A_{\alpha_0}^y=A_{\alpha_0}(\eta_\beta)\), and Lemma (ref) applies. Equation (ref) shows that the supervised operator is diagonal in the topic basis. If the entries of \(\beta\) are distinct, this relation may be read as an eigenvalue statement. For inference, however, we use the projection formula (ref), which avoids the complications caused by repeated or zero coefficient values once the topic basis has been fixed.

Identification of the concentration mass

The preceding identities still require the total concentration \(\alpha_0\) as an input. This is a central practical limitation of corrected spectral LDA moments: the correction that diagonalizes the word moments depends on a quantity that is rarely known in applications. We now show that the same operator representation identifies \(\alpha_0\) from observable word moments. The identifying restriction is commutativity.

At the true value, (ref) implies that every corrected operator is diagonal in the same topic basis. Hence any two such operators commute: for every pair \(\eta_1,\eta_2\),

where \([R,S]=RS-SR\). We show that this property is generically unique to the true concentration: away from \(\alpha_0\), mean-orthogonal contractions produce nonzero commutators. Thus commutativity supplies an identifying restriction for the concentration mass. Let

If \(\eta\in\mu^\perp\), then \(w=O^\top\eta\in\alpha^\perp\), because \[ \alpha^\top O^\top\eta=\alpha_0\mu^\top\eta=0. \]

The restriction \(k\ge3\) is used only here. It ensures that \(\alpha^\perp\) has dimension at least two, so non-collinear mean-orthogonal contractions exist. The present commutator argument is therefore formulated for topic models with at least three topics.

For estimation, one does not need to check commutativity over all \(\eta\in\mu^\perp\). Let \(s_1,\ldots,s_L\in\mathbb{R}^d\), define

and let \(\mathcal I\subset\{(\ell,q):1\le \ell<q\le L\}\) be non-empty. Consider the population finite-probe criterion

If there is at least one pair \((\ell,q)\in\mathcal I\) for which \(O^\top v_\ell\) and \(O^\top v_q\) are non-collinear, then

If two raw probes \(s_1,s_2\) are drawn independently from distributions that are absolutely continuous on \(\mathbb{R}^d\), this non-collinearity condition holds with probability one. The sample estimator in Section (ref) is the plug-in version of (ref), using empirical moments and empirical mean-orthogonal projections.

The commutativity result also has a purely unsupervised implication. Combined with the corrected-moment eigenvector identity in Theorem (ref), it implies that the finite-dimensional LDA model is identified from low-order word moments without supplying the total concentration as an input. We record this consequence explicitly.

Estimation and inference

This section turns the population identities of Section (ref) into a feasible estimator. The construction has three steps. First, estimate the observable cross-token moments by within-document averages over distinct token positions. Second, estimate the concentration mass \(\alpha_0\) by minimizing a finite set of empirical commutators. Third, recover the topic matrix from one unsupervised spectral operator and recover the downstream coefficient from the supervised operator. The limit theory treats these steps as one smooth map of empirical moments, after rank selection and topic order have been determined.

Empirical cross-token moments

For document \(i\), define the within-document empirical mean

The second cross-token moment is estimated by averaging over ordered distinct pairs,

Because ordered pairs include both \((a,b)\) and \((b,a)\), \(\hat M_{2,i}\) is symmetric. Equivalently, if \(c_i=\sum_{a=1}^N x_{ia}\) is the word-count vector, then

For a contraction direction \(v\in\mathbb{R}^d\), the third cross-token moment is estimated by

The corresponding sample averages are

The hats in this section refer to the sample size \(n\), which is suppressed when no ambiguity can arise. The estimators in (ref) are unbiased for \(\mu\), \(M_2\), and \(T(v)\), respectively. The use of distinct token positions removes the single-token multinomial contribution that would otherwise appear on the diagonal of the second moment.

The supervised moments are estimated in the same way. Define

Thus \(\hat T^y\) is the response-weighted analogue of the second cross-token moment. The full third-order tensor need not be stored in computation; only the finite number of contractions \(\hat T(v)\) used by the estimator are required.

Large-vocabulary preprocessing. The formulae above are written in the original vocabulary dimension. In large vocabularies, one can instead compress tokens linearly before forming the corrected moments. Appendix (ref) shows that if a fixed-dimensional compression matrix does not lose a topic direction, in the sense that \(\operatorname{rank}(R^\top O)=k\), then the corrected moment identities and the coefficient \(\beta\) are unchanged. When observed controls are included after compression, the control adjustment also requires the latent topic-share scale. Appendix (ref) gives a compressed scale-recovery step based on the compressed first and corrected second moments; this step is used before constructing the control cross-moment \(M_{qh}\). The split-sample PCA version estimates an \(m\)-dimensional compression, with fixed \(m\ge k\), from the uncentered cross-token second moment on an independent split, holds it fixed, and runs the moment estimator on the compressed tokens in the second split. This is intended as a compatibility result for fixed \(k\) and fixed \(m\), not a rigorous growing-vocabulary or growing-rank asymptotic theory.

Plug-in corrected operators

For each candidate \(\tau>0\), define the empirical corrected moments

and

Similarly, the supervised corrected moment is

These are the sample analogues of (ref), (ref), and (ref).

The corrected second moment has rank \(k\) in population, so the inverse in the empirical operator is taken on an estimated rank-\(k\) subspace. Let \[ \hat B_\tau=\sum_{j=1}^d\hat\lambda_j(\tau)\hat u_j(\tau)\hat u_j(\tau)^\top, \qquad \hat\lambda_1(\tau)\ge\cdots\ge\hat\lambda_d(\tau), \] be a spectral decomposition. For a candidate rank \(m\), set

whenever \(\hat\lambda_m(\tau)>0\). If \(m=0\), set \(\hat B_{\tau,0}^+=0\). Values of \(\tau\) for which \(\hat\lambda_m(\tau)\le0\) are excluded from the minimization below. The empirical unsupervised and supervised operators are

If \(k\) is not specified in advance, estimate it from the rank of \(M_2\). Let \(\tilde\lambda_1\ge\cdots\ge\tilde\lambda_d\) be the eigenvalues of \(\hat M_2\), and let \(a_n\) satisfy

Define

When the number of topics is chosen by design or by an external model-selection step, the same formulas apply with \(\hat k\) replaced by that chosen value.

In finite samples, implementing the rank rule in (ref) still requires a numerical choice of threshold, and nearby candidate ranks may have similar empirical second-moment spectra. When a split-sample implementation is used, as in the application below, rank choice can be made on the first split before the response-weighted moment estimation step. We use the commutator identity as a practical response-free rule of thumb: for each candidate \(k\), compute the commutator-based estimate \(\hat\alpha_0(k)\) across several admissible probe draws on the first split, and prefer ranks for which this estimate is interior and stable. This is a finite-sample rule of thumb for rank selection, not an additional rank-consistency theorem. The motivation is that, under the maintained finite-LDA model and the correct rank, admissible finite-probe commutator criteria have the same population minimizer \(\alpha_0\). Thus a working rank for which \(\hat\alpha_0(k)\) moves substantially across reasonable probe draws is less attractive than one for which the concentration estimate is interior and stable. Subsequent estimation and standard errors are computed on the independent second split, conditional on the selected rank, compression, and realized probes.

This split-sample design also allows the commutator calculation to be used a second time as a diagnostic step. Since the working rank is chosen on the first split, the commutator profile can be inspected again on the second split as a holdout diagnostic for the selected finite-dimensional specification. If the same sample were used both to choose the rank and to assess commutator stability, this second calculation would mainly repeat the selection criterion. With an independent second split, boundary behavior, flatness, multimodality, or probe sensitivity provides a separate warning about the stability of the selected spectral diagonalization. We discuss these post-selection diagnostics in Section (ref).

Estimating the concentration mass

Let \(P(u)=I_d-uu^\top/\|u\|^2\) for \(u\ne0\). Choose raw probe vectors \(s_1,\ldots,s_L\in\mathbb{R}^d\), with \(L\ge2\), and define their empirical mean-orthogonal projections

Let \(\mathcal I\subset\{(\ell,q):1\le\ell<q\le L\}\) be a non-empty set of probe pairs. For fixed rank \(m\), define the empirical commutator criterion

where \([R,S]=RS-SR\). If \(\hat\lambda_m(\tau)\le0\), set \(\hat Q_{\mathcal I}^{(m)}(\tau)=+\infty\). The estimator of the total concentration is

where \(\Theta=[\underline\tau,\overline\tau]\subset(0,\infty)\) is compact and contains the true \(\alpha_0\) in its interior. For the true rank \(m=k\), the finite-probe identification result in Section (ref) implies that the population version of (ref) is uniquely minimized at \(\alpha_0\), provided at least one tested pair of projected probes is non-collinear in the latent coordinates. Thus the sample criterion has the correct local target on the event \(\hat k=k\). No uniqueness claim is made here for a misspecified truncation rank.

The probability-one genericity statements are identification results, not finite-sample guarantees. In applied work one should inspect the profile of \(\hat Q_{\mathcal I}^{(m)}(\tau)\), whether \(\hat\alpha_0\) lies near the boundary of \(\Theta\), local curvature around the minimizer, eigenvalue gaps for the ordering operator, and sensitivity to the realized probe and ordering directions. A flat, multimodal or boundary-attaining commutator criterion, or estimates that move materially across admissible probes, should be treated as a substantive warning. It may indicate weak finite-sample identification, an unstable rank choice, ill-conditioned topic separation, an inadmissible compression, or misspecification of the maintained LDA moment structure.

Topic and downstream coefficient estimators

The concentration estimator is then plugged into a single spectral operator to recover the topic columns. Choose an ordering direction \(r\in\mathbb{R}^d\), independent of the data if random, and set

Define

On the high-probability event that \(\hat H^{o}\) has \(\hat k\) eigenvalues separated from the zero cluster, and that these selected eigenvalues are real, simple and mutually separated, let \(\hat o_1(r),\ldots,\hat o_{\hat k}(r)\) be the associated right eigenvectors. The selected nonzero eigenpairs are labeled by decreasing real eigenvalue within this selected cluster, and the columns are normalized by

Off this event, the estimator may be defined by any fixed deterministic convention; this off-event convention is asymptotically irrelevant. Set

The role of \(r\) is only to select a stable ordering of the topic columns. At the population value, the eigenvectors are the columns of \(O\), ordered by the entries of \(O^\top P(\mu)r\).

Finally, form the supervised operator at the estimated concentration,

The downstream coefficient estimator is

where \(\operatorname{diag}(\cdot)\) extracts the diagonal as a vector. This is the direct sample analogue of (ref). It does not estimate document-level topic shares and does not run a second-stage regression on generated regressors.

Asymptotic linearity and standard errors

The asymptotic argument is fixed-dimensional: \(d\), \(k\), and \(N\) are fixed while \(n\to\infty\). Let \(\hat{\mathcal T}_i\) denote the ordered-distinct triple tensor \[ \hat{\mathcal T}_i :=\frac{1}{N(N-1)(N-2)} \sum_{\substack{a,b,c=1\\ a,b,c\;\mathrm{distinct}}}^N x_{ia}\otimes x_{ib}\otimes x_{ic}, \] so that \(\hat T_i(v)=\hat{\mathcal T}_i(I,I,v)\). Collect the document-level moments needed by the estimator in

Let \(Z_0=\mathbb{E} Z_i\) and \(\bar Z=n^{-1}\sum_i Z_i\). Since the word components are bounded, the condition \(\mathbb{E} Y^2<\infty\) is enough for

On the event that the rank estimate is correct and the selected nonzero-cluster ordering eigenvalues are real, simple, mutually separated, and separated from the zero cluster, all estimators above are functions of \(\bar Z\). The commutator minimization is differentiable at the population value because the stacked commutator has a nonzero derivative with respect to \(\tau\) at \(\alpha_0\). The eigenvector and pseudoinverse maps are differentiable on the corresponding fixed-rank, separated-cluster stratum.

The following theorem records the resulting first-order expansion. If the probes or ordering direction are random, the statement is conditional on their realized values.

This inferential expansion is one of the key distinctions from the supervised spectral recovery literature: the goal is not only consistent recovery of supervised-LDA parameters, but a first-order sampling expansion and feasible variance estimator for the downstream coefficient.

The proof of the theorem and corollary is given in the Supplement. Its main steps are standard but useful to separate. The document-level moments satisfy the central limit theorem (ref). The minimizer \(\hat\alpha_0\) has the usual one-dimensional minimum-distance expansion obtained by differentiating the stacked commutator equations at \(\alpha_0\). The topic estimator then follows from first-order perturbation theory for simple right eigenvectors. Finally, (ref) is a smooth map of \((\hat\alpha_0,\hat B_{\hat\alpha_0,\hat k}^+,\hat A_{\hat\alpha_0}^y,\hat O)\); the first-order perturbation of \(\hat O\) contributes only off-diagonal terms inside \(\hat O^+\hat H^y\hat O\), so the diagonal map recovers the coefficient influence function without requiring distinct entries of \(\beta\).

For implementation, let \(G_{\beta,r}\) denote the population map from \(Z_0\) to \(\beta_r\) induced by (ref)--(ref), with the realized probes and ordering direction held fixed. Let

Then the lower block of \(V(r)\) is

All reported standard errors use the influence-function implementation of this derivative. Write the finite-probe commutator map as \(g(\tau,z)\), and let \(\alpha(z)\) be the local minimum-distance solution. At an interior solution,

For fixed \(\alpha_0\), the derivative of the coefficient map is obtained by differentiating the corrected moment operators, the truncated inverse of the corrected second moment, the normalized right eigenvectors used to form \(\hat O(r)\), and the final diagonal map in (ref). The full derivative has the chain-rule form

where \(D_{\beta,r}^{\rm fix}\) is the derivative with \(\alpha_0\) held fixed and \(D_{\beta,\alpha}\) is the derivative of the coefficient map with respect to \(\alpha_0\). The empirical influence values are therefore

Equivalently, with

where \(\hat D_{\beta,r}\) is the sample analogue of (ref). Equations (ref)--(ref) describe the no-control moment vector. If observed controls are included, the same formulas are applied after replacing \(Z_i\), \(\bar Z\), \(D_{\beta,r}\), and \(\Omega\) by the enlarged control moment vector \(Z_i^c\), its mean, the derivative of the control-adjusted map, and \(\operatorname{var}(Z_i^c)\). If the data are first compressed, \(Z_i^c\) is formed from the compressed tokens and includes the compressed scale-recovery moments used by the control adjustment. A pointwise \((1-a)\) confidence interval for the \(j\)-th ordered coefficient is

The same covariance estimate is used for Wald tests and linear contrasts of the ordered coefficient vector. For pointwise confidence intervals we assume \(V_{\beta,jj}(r)>0\). For Wald tests of \(C\beta=c\), we assume \(CV_\beta(r)C^\top\) is nonsingular.

Monte Carlo experiments

This section evaluates the finite-sample behavior of the proposed estimator in the setting for which the theory is designed: downstream regression with a fixed number of words per document. The experiments are organized around two questions. First, does the direct moment-based estimator deliver calibrated inference for the coefficient vector \(\beta\)? Second, what goes wrong if the analyst instead constructs document-level topic estimates and treats them as observed regressors? The number of topics is supplied to all methods, so the exercise isolates the inferential problem rather than the separate rank-selection step.

To make the fixed-document plug-in issue explicit, consider any deterministic reconstruction rule \(\tilde h_i=g_N(c_i;O)\) based on the document counts and the known population topic matrix. If a no-intercept population least-squares regression of \(Y_i\) on \(\tilde h_i\) is nonsingular and \(\mathbb{E}(\tilde h_i\varepsilon_i)=0\), its probability limit is

Thus the plug-in coefficient equals \(\beta\) only under additional restrictions, such as \(\mathbb{E}(\tilde h_i h_i^\top)=\mathbb{E}(\tilde h_i\tilde h_i^\top)\). The correctly specified posterior mean \(\tilde h_i=\mathbb{E}(h_i\mid c_i)\) is a special case satisfying this identity by iterated expectations, but a generic finite-\(N\) reconstruction rule, including the constrained least-squares rule below, need not satisfy it. The simulations therefore use the known-\(O\) plug-in regression as an analytic and numerical benchmark for the persistent generated-regressor problem at fixed document length.

Design and methods

The vocabulary size is \(d=100\), the number of topics is \(k=10\), and each document contains \(N=100\) words. The topic matrix \(O\) is generated once and then held fixed across Monte Carlo replications. Its columns are independently generated from a Dirichlet distribution with concentration parameter \(0.3\), subject to screening restrictions that remove nearly duplicate topics and nearly anchor-word designs. For the realized matrix, the minimum pairwise Hellinger distance between topic columns is \(0.588\), the maximum row dominance is \(0.879\), the minimum Hellinger distance from any row-normalized word profile to a simplex vertex is \(0.250\), and there are no exact zeros.

For document \(i\), \[ h_i\sim\operatorname{Dirichlet}(\alpha),\qquad c_i\mid h_i\sim\operatorname{Multinomial}(N,Oh_i), \] and the response is \[ Y_i=\beta^\top h_i+\varepsilon_i, \qquad \varepsilon_i\sim N(0,\sigma^2), \qquad \beta=(1.0,0.9,\ldots,0.1)^\top . \] The noise variance is calibrated to give population \(R^2=0.35\). We consider a symmetric Dirichlet design, \(\alpha_j=0.5\) for \(j=1,\ldots,10\), so that \(\alpha_0=5\), and an asymmetric design, \(\alpha_j=2j/55\), so that \(\alpha_0=2\). The main downstream comparison uses \(n\in\{3000,5000,15000\}\); each design cell uses \(1000\) Monte Carlo replications. The concentration-parameter table below also reports the original \(n=30000\) design point.

We compare four procedures. The first is an infeasible no-intercept latent-oracle regression that uses the true \(h_i\)'s. The second is the proposed direct spectral estimator, which estimates \(\alpha_0\), recovers the topic matrix, and estimates \(\beta\) from the supervised corrected moments. Its confidence intervals use the influence-function implementation of the first-order variance formula developed above: the derivative of \(\hat\alpha_0\) is obtained from the implicit minimum-distance expansion of the commutator criterion, and the derivative of \(\hat\beta\) is the fixed-\(\alpha_0\) spectral derivative plus the chain-rule contribution from estimating \(\alpha_0\). The third procedure is a known-\(O\) plug-in regression: it uses the true topic matrix \(O\), but estimates document-level topic mixtures by \[ \tilde h_i =\arg\min_{h\in\Delta^{k-1}} \left\|\frac{c_i}{N}-Oh\right\|_2^2, \] and then applies no-intercept ordinary least squares treating \(\tilde h_i\) as observed. The no-intercept specification is natural because \(\mathbf 1^\top h_i=1\). The fourth procedure is the analogous spectral plug-in regression, replacing \(O\) by the unsupervised spectral estimate \(\hat O\). Thus the final two procedures mimic the modular workflow used in many empirical applications. Topic labels are aligned to the true labels by Hungarian matching only when computing simulation errors and coverage.

For each method, Table (ref) reports the root mean squared error \[ \operatorname{RMSE}_\beta = \left\{ \frac{1}{1000}\sum_{r=1}^{1000} \frac{1}{k}\|\hat\beta^{(r)}-\beta\|_2^2 \right\}^{1/2}, \] the average coordinate-wise coverage of nominal \(95\%\) confidence intervals, the minimum coordinate-wise coverage, average interval length, and average Hellinger distance between matched topic columns.

The contrast between the direct estimator and the plug-in regressions is the main finding. In the symmetric design, the direct spectral estimator has average coverage \(0.978\), \(0.969\), and \(0.962\) as \(n\) increases from \(3000\) to \(15000\). Over the same sample sizes, the spectral plug-in regression has coverage \(0.529\), \(0.461\), and \(0.288\). The known-\(O\) plug-in regression, which is given the true topic matrix, also undercovers severely: its coverage falls from \(0.607\) to \(0.272\). Thus the failure is not primarily a topic recovery error. It is the fixed-document generated-regressor error induced by estimating \(h_i\) from only \(N=100\) words and then treating \(\tilde h_i\) as observed.

The asymmetric design gives the same qualitative conclusion. The direct spectral estimator remains close to nominal coverage, whereas the true-topic plug-in coverage falls from \(0.818\) to \(0.506\), and the spectral plug-in coverage falls from \(0.647\) to \(0.460\). The plug-in intervals have lengths comparable to the latent-oracle intervals, so they shrink at the usual \(n^{-1/2}\) rate even though their centers retain finite-document bias. With \(1000\) replications, the Monte Carlo standard error of a single coordinate's nominal \(95\%\) coverage estimate is about \(\{0.95(0.05)/1000\}^{1/2}=0.0069\), far smaller than the reported plug-in undercoverage.

As a further diagnostic, we reran the spectral plug-in regression while giving it the true value of \(\alpha_0\). Across the six design cells in Table (ref), this changed average coverage by at most \(0.039\). The plug-in failure is therefore not explained by uncertainty in the concentration parameter; it comes from using finite-document topic estimates as regressors.

Estimation of the concentration parameter

Table (ref) reports the behavior of the commutator estimator of \(\alpha_0\). Bias and RMSE decrease steadily with \(n\) in both Dirichlet designs. At the two largest sample sizes, the empirical standard deviation and the average estimated standard error are close, supporting the first-order linearization used in the downstream variance formula. Together with Table (ref), this indicates that the proposed method learns \(\alpha_0\) from word moments and delivers calibrated downstream inference without plug-in document-topic regressors.

The simulation designs deliberately hold fixed the topic dimension, topic separation, document length, concentration search interval and probe construction. They are therefore best read as tests of the fixed-document generated-regressor mechanism and of the first-order variance formula in well-separated designs. They do not replace application-specific diagnostics for the commutator criterion, ordering eigengaps, rank selection or compression stability.

Application: abstract topics and citations in economics

We illustrate the method on articles published in the five general-interest economics journals American Economic Review, Econometrica, Journal of Political Economy, Quarterly Journal of Economics, and Review of Economic Studies. Article metadata and citation counts are taken from OpenAlex PriemPiwowarOrr2022,OpenAlexData2026. The response is \[ Y_i=\log\{1+\text{citation count}_i\}, \] and the document consists of the article title and abstract. The analysis is descriptive: a topic coefficient is interpreted as an association between latent abstract content and citation impact, not as a causal effect of writing on that topic. The exercise should also be read under the maintained low-order response-token orthogonality restrictions, which cannot be verified from the observed text alone.

The maintained response equation in this application is \[ Y_i = \beta^\top h_i + \delta^\top q_i + \varepsilon_i, \] where \(h_i\) is the latent topic mixture and \(q_i\) contains publication-year indicators. Since \(\mathbf{1}^\top h_i=1\), a common shift in all topic coefficients is observationally equivalent to an intercept shift. Thus the absolute level of an individual coefficient \(\beta_j\), and a test of \(\beta_j=0\), depends on the normalization implicit in the intercept. The invariant objects are topic contrasts, \[ \beta_a-\beta_b, \] which compare the conditional mean response when latent topic mass is shifted from topic \(b\) to topic \(a\). A reallocation of 10 percentage points of topic mass from \(b\) to \(a\) corresponds to \(0.1(\beta_a-\beta_b)\) log-points in \(Y_i\), holding the year controls fixed. We therefore report the topic coefficient estimates to orient the reader, but interpret the empirical results through joint and pairwise contrast tests.

The raw corpus contains records published between 2000 and 2018. We exclude records whose titles or metadata indicate comments, replies, corrections, editorial material, or very short documents. We lowercase tokens, remove standard English and scholarly stopwords, and keep terms whose document frequency is at least 40. The main specification does not stem words. The resulting analysis sample has 5409 documents and 1771 terms. Because the vocabulary dimension is too large for direct dense third-order moment calculations, we use the split-sample compression described in Appendix (ref). A 10% split of documents is used to estimate a 100-dimensional orthonormal compression from the uncentered cross-token second moment, and the remaining 90% is used for estimation and inference. All reported confidence intervals and standard errors are computed on the second split, conditional on the first-split compression. This follows the split-compression construction in Appendix (ref), under which the first-stage compression is held fixed for second-split inference. The empirical moments use document-specific length normalizations as in Appendix (ref).

Following the split-sample finite-rank diagnostic described in Section (ref), the working topic dimension is selected using only first-split word moments, without using citation outcomes. For each candidate \(k\in\{10,\ldots,20\}\), we estimate \(\alpha_0\) from the commutator criterion across 20 independent sets of probe directions, using five mean-orthogonal commutator partner directions for each run. Candidates with boundary or failed concentration estimates are discarded. Among the remaining candidates, we choose the value of \(k\) with the smallest standard deviation of \(\log \hat\alpha_0\). This first-split commutator-stability diagnostic selects \(k=15\).

After fixing this rank and the first-split compression, we run the estimator on the independent second split and inspect the second-split commutator profile, boundary behavior, compressed scale recovery, and control-adjustment conditioning as diagnostics for the selected specification. In a second-split probe diagnostic using five admissible probe draws, \(\hat\alpha_0\) ranges from 5.09 to 7.59, and none of the corresponding commutator profiles has a boundary minimum. The second-split estimator uses year fixed effects, with year 2000 as the omitted category. The coefficient estimates and standard errors are computed in compressed coordinates after the moment-based topic-scale recovery used in the observed-control adjustment; lifted word-space directions are used only to label topics. The selected specification has \[ \hat\alpha_0=5.09, \qquad n_{\rm est}=4868, \qquad \operatorname{cond}(\hat\Gamma_q)=6.56, \] where \(\hat\Gamma_q\) is the matrix in the observed-control adjustment. The minimum absolute recovered scale factor is 1.78, and the value of \(\operatorname{cond}(\hat\Gamma_q)\) is moderate in this application, so the year-control adjustment is not close to singular by this diagnostic.

Table (ref) reports the topic labels and centered coefficient estimates in the selected ordering. The labels are assigned from the leading positive words of the lifted word-space directions. The centered estimates subtract the average topic coefficient and are shown only to orient the reader; the invariant inferential objects are the joint equality test and pairwise contrasts reported below.

The topics align with recognizable fields and methods in economics. The recovered topic labeled “Econometrics / identification / asymptotics”, for example, has leading words models, estimator, estimators, estimation, asymptotic, and identification. Other topics correspond to macroeconomic shocks, monetary policy, labor search, public finance, firms and productivity, auctions, mechanism design, and private information. This is a qualitative check that the spectral directions are not merely recovering common article boilerplate.

The natural omnibus null in this setting is not \(\beta=0\), but equality of all topic coefficients, \[ H_0:\beta_1=\cdots=\beta_{15}. \] This null says that, after controlling for publication year, latent topic composition has no association with citation impact. Using the full influence-function covariance matrix, the Wald statistic for this 14-dimensional restriction is 56.26, with \(p=5.25\times10^{-7}\). Thus, within the maintained latent-topic model, abstract topic composition is associated with citation impact.

Table (ref) reports the pairwise contrasts that remain significant after Holm adjustment over all 105 pairwise comparisons. The largest contrast is between the macro/monetary shocks topic and the econometrics/identification topic. A ten-percentage-point shift of latent topic mass from the latter to the former corresponds to 0.187 log-points in the citation outcome, holding publication year fixed. The firms/trade/productivity topic also has a larger citation association than the econometrics/identification topic and the auctions/prices/products topic. These are descriptive contrasts in latent topic coordinates, not causal comparisons across research fields.

The direction of these contrasts is consistent with existing evidence on citation heterogeneity across fields in economics. CardDellaVigna2013 study the same five general-interest journals and document substantial differences in citation patterns across fields, including relatively high citation performance for more recent Development and International Economics papers and lower relative citation performance for more recent Econometrics and Theory papers. AnauatiGalianiGalvez2016 classify top-five articles into applied, applied theory, econometric methods, and theory, and find that citation life cycles differ markedly across these categories: applied and applied-theory papers have more favorable citation profiles, while theory and econometric-method papers have lower typical citation profiles, with econometric methods displaying a more heterogeneous upper tail. Our topic-level estimates should not be read as a replication of those field classifications, but the broad pattern is similar: topic content is associated with citation outcomes, and the more applied macro and firm/trade/productivity directions have higher citation associations than some econometric and theory-oriented directions.

The application is intended as a numerical illustration of the moment-based inferential workflow. It uses a single first-split selection rule, fixed before estimating the second-split coefficients, and all reported uncertainty is computed on the independent second split, conditional on the first-split compression and selected working dimension. Appendix (ref) reports second-split diagnostics, adjacent-rank checks, and a preprocessing robustness check. The main text is kept to one specification to avoid turning the illustration into a model-selection exercise.

Discussion

This paper develops direct inference for regression with latent Dirichlet topic covariates. The main distinction from the usual plug-in workflow is that document-level topic shares are not estimated and then treated as observed regressors. At fixed document length, those document-level quantities remain noisy even when the topic matrix is known. The proposed response-weighted corrected moments instead identify the downstream coefficient \(\beta\) directly from observable word--response moments under explicit low-order response-token orthogonality conditions. The same corrected moment structure also yields an estimating equation for the unknown total concentration parameter \(\alpha_0\), through the commutativity of corrected word-moment operators.

The results are model-based and should be interpreted within the stated asymptotic regime. The formal theory keeps the vocabulary dimension, number of topics and document length fixed while the number of documents increases. The commutativity identification of \(\alpha_0\) is formulated for \(k\ge3\) and relies on rank, separation and generic probe conditions. These genericity conditions should not be read as a substitute for numerical diagnostics in applied use. Boundary solutions, flat commutator profiles, small ordering eigengaps, unstable rank choices or substantial sensitivity to admissible probes are substantive warnings. They may reflect weak finite-sample identification, ill-conditioned topic separation, an unstable compression, or misspecification of the maintained LDA moment structure. The split-compression results in the Supplement identify \(\beta\) conditional on an admissible projection; compressed directions are not, by themselves, simplex-normalized topics.

Integrated latent-variable likelihood or Bayesian joint-model procedures could also be used to target \(\beta\) under fuller specifications by integrating over the latent document mixtures. The contribution here is a different tradeoff: low-order moment-based frequentist inference for \(\beta\), without document-level topic-share regression and without a full likelihood or prior specification for all unknown population parameters.

Several extensions remain. A theory allowing the vocabulary dimension or the compression dimension to grow with the number of documents would better match large-scale text applications. It would also be useful to study robustness under departures from the Dirichlet topic-mixture assumption, including non-Dirichlet latent mixtures, correlated topics and misspecified topic dimension. Finally, the same moment logic may be useful beyond topic models, in other latent simplex models where the inferential target is a downstream regression coefficient rather than recovery of each unit's latent coordinate.