2.1 Definitions of convergence of mean scores

Here, we explain why mean scores approximate true scores under Assumption 1 and consider this convergence under weaker assumptions than Assumption 1. We assume that $\theta $ is a random vector of finite dimension, which allows us to deal with conditional probabilities as though they were standard probability measures (Shiryaev, Reference Shiryaev1996, Theorems 4 and 5, Chapter II.7), as considered in more detail below as needed. We start by recalling the technical definition of the conditional independence made in Assumption 1(1)(b).

We have the following useful representation for conditionally independent variables. While the representation is a consequence of abstract representation results found, for example, in Kallenberg (Reference Kallenberg2021, Chapter 8), an elementary proof is provided for completeness. As in standard stochastic representations, we assume that the probability space is sufficiently rich to support a sequence $U_1,U_2,\ldots $ of independent and identically distributed uniform random variables on $[0,1]$ that are independent of $\theta $ . This can be achieved by possibly expanding the probability space (see Kallenberg, Reference Kallenberg2021, p. 176) and is therefore assumed essentially without loss of generality. For a cumulative distribution function F, let $F^{-1}(x) = \inf \{ t: F(t) \geq x \}$ be its generalized inverse.

Lemma 1. Suppose $X_1, X_2, \ldots $ are conditionally independent given $\theta $ . Then the distribution of $(X_1, X_2, \ldots )'$ is the same as if

(3) $$ \begin{align} X_j = F_{j|\theta}^{-1}(U_j), \qquad j \geq 1, \end{align} $$

where $U_1,U_2, \ldots $ are each uniform on $[0,1]$ and are independent of each other and of $\theta $ .

By slight abuse of notation, we will sometimes write $X_j = F_{j|\theta }^{-1}(U_j)$ , as all results in this article concern properties of distributions. For conditionally independent items, $J^{-1} \sum _{j=1}^{J} X_j = J^{-1} \sum _{j=1}^{J} F_{j|\theta }^{-1}(U_j)$ which for fixed $\theta $ will be close to the average of $E [F_{j|\theta }^{-1}(U_j) ]_U$ by the law of large numbers for independent variables, where $E[\cdot ]_U$ means integration with respect to $(U)$ only. The independence between $(U)$ and $\theta $ implies $E [F_{j|\theta }^{-1}(U_j) ]_U = E [ X_j | \theta ]$ , and we therefore expect

$$\begin{align*}D_{J} = J^{-1} \sum_{j=1}^{J} (X_j - E [X_j | \theta]) = \bar X_{J} - T_{J}(\theta) \end{align*}$$

to converge to zero. This convergence, in the sense of the upcoming Definition 2(2), defines consistent estimability, similar to Stout (Reference Stout1990) and Junker (Reference Junker1991).

Recall that an essential supremum with respect to a distribution say $F_\theta $ , is one that is taken up to $F_\theta $ -null sets, where $F_\theta $ is here the distribution of $\theta $ . Fixed values of $\theta $ are denoted by $\vartheta $ . When pointwise notation such as $P_\vartheta $ , $E_\vartheta $ , $\text {Var}\,_\vartheta $ , and $\mathrm {Cov}_\vartheta $ is used, it refers to a fixed regular conditional probability kernel $\vartheta \mapsto P_\vartheta :=P(\,\cdot \,|\theta =\vartheta )$ . Thus, for a measurable function g, $\mathrm {ess\,sup}_{\vartheta } g(\vartheta )$ denotes the smallest $c\in [-\infty ,\infty ]$ such that $g(\vartheta )\le c$ for $F_\theta $ -almost every $\vartheta $ . As far as we can tell, (3)–(5) are new at least in the present context, although they are direct adaptations of well-known definitions to conditional probability measures.

As a direct consequence of the definitions and well-known results, the convergence concepts in Definition 2 are progressively stronger, and (4) implies that X is strongly unconditionally consistently estimable.

Adapting the terminology of Stout (Reference Stout1990) (again omitting his sub-test perspective, see Remark 2), we call a sequence of items that fulfill

(4) $$ \begin{align} 2 J^{-2} \sum_{1 \leq j < k \leq J} \mathrm{Cov}(X_j, X_k |\theta) \xrightarrow{} 0 \quad \text{a.s.} \end{align} $$

essentially independent. By the same Chebyshev argument as in Stout (Reference Stout1990), binary items that are essentially independent fulfill the criteria of consistent estimability. The proof also transfers to bounded items in general.

For unconditional and uniform conditional consistency, similar simple criteria exist, which we now derive. We restrict attention to bounded items, that is, each item in X is contained within a fixed interval. For unbounded items, further moment assumptions are required.

2.2 Martingale factor models

We now propose a model class we call martingale factor models in a latent variable $\theta $ . A main conclusion of this section is that martingale factor models (and hence conditionally independent items) with bounded items are consistently estimable with exponential concentration.

As mentioned in the introduction, martingale factor models replace conditional independence with the weaker requirement that the corresponding conditional means agree. This allows conditional dependence in other aspects of the conditional distribution, such as having the conditional variation of an item depending not only on $\theta $ but also on the candidate’s previous responses. At the end of this section (Proposition 5), we provide a class of examples which are martingale factor models but may have conditional dependence. Some concrete illustrations are given in Appendix D of the Supplementary Material.

Definition 3 (Martingale factor model)

Suppose $X_1, X_2, \ldots $ are random variables, and $\theta $ is a random vector. Suppose that for all $j \geq 1,$ we have that $X_j$ is integrable. Define $\mathscr {F}_j := \sigma (\theta , X_1, \ldots , X_j)$ for $j \geq 1$ and $\mathscr {F}_0 = \sigma (\theta )$ , where $\sigma (\cdot )$ denotes the sigma-algebra generated by the indicated variables. If

(5) $$ \begin{align} E ( X_j | \mathscr{F}_{j-1}) = E [X_j | \theta], \quad \text{for all } j \geq 1, \end{align} $$

we say that $X_1, X_2, \ldots $ follow a martingale factor model in $\theta $ .

Conditional independence does not take into account the direction of time, see Definition 1 and Lemma 1, though it can be described sequentially, see eq. (2). In contrast, Definition 3 depends on the order of the items, and cannot be described without reference to the order of time. Since items are answered sequentially, this fits psychometric practice. The mathematics leaves the interpretation of this directionality open, and one natural interpretation is that the respondent’s mental state is influenced by having answered previous questions, see the discussion after Proposition 5.

The conditional expectations in Definition 3 are (nonparametric) least squares best predictors (e.g., Klenke, Reference Klenke2020, Corollary 8.17), so a martingale factor model says that once $\theta $ is known, no additional predictive power for item j is gained by knowing all earlier answers to the items. In contrast, conditional independence in $\theta $ means that once $\theta $ is known, additionally knowing all earlier items is completely uninformative for the current item, not merely at the level of its conditional mean, but for its entire conditional distribution.

Recall (e.g., Davidson, Reference Davidson1994, Chapter 15.1) that a martingale $(X_j)$ with respect to a filtration (i.e., an increasing sequence of sigma-algebras) $(\mathscr {F}_j)$ fulfills $E(X_j|\mathscr {F}_{j-1}) = X_{j-1}$ for all j and is integrable (meaning $E |X_j| < \infty $ for all j). Definition 3 is therefore a martingale-like assumption.

Recall (e.g., Davidson, Reference Davidson1994, Chapter 15.2) that $(Z_j)$ is a martingale difference sequence under an expectation operator $\mathbf {E}$ with respect to a filtration $(\mathscr {A}_j)$ if (i) $Z_j$ is $\mathscr {A}_{j}$ -measurable, (ii) $\mathbf {E} |Z_j| < \infty $ for all j, and (iii) $\mathbf {E}[Z_j | \mathscr {A}_{j-1}] = 0$ .

Due to Lemma 3, we may apply martingale theory to study martingale factor models. We now give two such applications, for simplicity restricted to bounded items, as this is of main interest in psychometric applications. The Azuma–Hoeffding bound (e.g., Davidson, Reference Davidson1994, Theorem 15.20, p. 245) gives the following result.

We may also apply martingale central limit theorems (see Hall & Heyde, Reference Hall and Heyde1980, Chapter 3). For central limit theorems, some restriction on second-order moments must be made. Here, we apply a central limit theorem from Chapter 24.2 in Davidson (Reference Davidson1994), using a general filtration $(\mathscr G_j)$ that is required for a later illustration. Lemma 3 provides a canonical choice for a filtration for martingale factor models.

Consider items following Assumption 1(1)(b). Then the assumptions of Proposition 4 hold with $\alpha _j = 1$ for all j when $\sum _{j=1}^J \text {Var}\,_\theta (X_j) \xrightarrow {} \infty $ as $J \xrightarrow {} \infty $ , yielding a conditional and unconditional central limit theorem for conditionally independent items.

Proposition 3 can be refined into multivariate and covariance bounds using, for example, results from Tropp (Reference Tropp2011). Proposition 4 can be extended in various directions (e.g., multivariate versions and functional central limit theorems) using the available martingale literature. We consider these extensions outside the scope of the present article.

We now show that a certain class of combinations of martingale factor models and exogenous randomness results in a new martingale factor model under certain assumptions. Since conditionally independent items are martingale factor models that can be represented as in Lemma 1, this provides a way to create a large space of martingale factor models. The construction is similar to that of discrete stochastic integral theory (e.g., Klenke, Reference Klenke2020, Chapter 9.3).

Proposition 5 allows $(A_{j,m})$ to depend on a randomization device $(V_j)$ , such as a Markov chain, which allows the formulation of complex models for item generation switching mechanisms. This may be used to create rich martingale factor models, as illustrated in Appendix D of the Supplementary Material, where we consider a hidden Markov model which switches between locally independent items in a way that may be interpreted as having the respondent’s levels of stress vary dynamically. In this situation, the martingale factor model has time-dependent dynamics due to test-taking stress, though the items themselves are locally independent items, and therefore time-invariant. Appendix D of the Supplementary Material also briefly sketches that it is possible to empirically detect deviations from the martingale factor model, and how Proposition 5 can be used to deal with missing values in martingale factor models.

Remark 5. Martingale factor models have applications also outside the strict study of mean scores. As discussed also in Junker (Reference Junker1991, Sections 4 and 5), likelihood-like estimators are often connected to weighted means. For example, suppose X is a martingale factor model in $\theta $ and let the conditional marginal item distributions follow the generalized partial credit model of Muraki (Reference Muraki1992) (see also, e.g., Van der Linden, Reference Van der Linden2016, Chapter 8). That is, let

(6) $$ \begin{align} P_\theta(X_j=k) = p_{j,\theta}(k) := \frac{\exp\{\gamma_{jk}+w_j\,k\,\theta\}} {\sum_{h=0}^{K}\exp\{\gamma_{jh}+w_j\,h\,\theta\}}, \qquad k=0,\ldots,K, \end{align} $$

where, for identifiability, we assume $\gamma _{j,0} = 0$ , and gather the remaining parameters in $\beta _j = (\gamma _{j,1}, \ldots , \gamma _{j,K}, w_j)$ for $j \geq 1$ . For a given respondent, we may consider the conditional item distribution and form the resulting independence-based likelihood, treating $\theta $ as a parameter. When solving the score equation for $\theta $ , the resulting estimator is a non-linear function of a weighted mean score, as calculated in Appendix C of the Supplementary Material. This non-linear relation to the mean score implies, with some additional arguments and assumptions, concentration inequalities and asymptotic normality for the resulting estimator.

As a final result, we characterize conditional independence through the martingale factor model property, showing that conditional independence is the transformation-invariant closure of the martingale factor model restriction. Closely related characterizations of conditional independence, phrased in terms of conditional expectations of bounded measurable transforms, appear in the abstract conditional independence literature (see, e.g., Dawid, Reference Dawid1980, Section 5, especially Definition 5.2). We provide a simple self-contained proof for completeness.

The martingale factor model is not invariant under non-linear transformations and therefore allows features such as conditional variances that depend on previous items. Proposition 6 makes precise how this model differs from conditional independence: conditional independence requires the martingale factor model assumption to hold for every bounded item transformation, whereas the martingale factor model itself requires this only for the identity transformation.

3.1 Sufficient conditions for convergent estimability

Here, we use standard convergence results (see, e.g., Billingsley, Reference Billingsley1995, Section 16) to provide two sufficient conditions for convergent estimability that are based on an item sequence being consistently estimable. These sufficient conditions are based on assuming that $\vartheta \mapsto E_\vartheta [\bar X_J]$ converges to a function pointwise. Conditions that imply this pointwise convergence are given in Section 3.2.

Under a monotonicity assumption, we provide the following sufficient condition for convergent estimability. Recall that if X is an unconditionally consistently estimable item sequence and $T_J(\theta )=T(\theta )+o_P(1)$ , then $\bar X_{J} = T(\theta ) + o_P(1)$ .

3.2 Items with conditional expectations parameterized by sequences which fulfill a law of large numbers result

We now present a general technique which can show that $T_J$ converges pointwise. Choose, for each $j\ge 1$ , a measurable version $m_j:\mathbb {R}^{d_\theta }\to \mathbb {R}$ , where $d_\theta $ is the dimension of $\theta $ , such that almost surely, $ m_j(\theta )=E[X_j|\theta ]$ . Such versions exist by Kallenberg (Reference Kallenberg2021, Theorem 8.1).

We may always represent $m_j(x)=H(x,\beta _j)$ for a function H and a sequence $(\beta _j)$ ; for example, one may take $d_\beta =1$ , $\beta _j=j$ , and $H(x,j)=m_j(x)$ . This representation is, however, non-unique, and for parametric models we can often express H in terms of standard parameter sequences.

We may often approximate the average $T_J(\vartheta ) = J^{-1} \sum _{j=1}^{J} H(\vartheta , \beta _j)$ by an integral, implying the pointwise convergence assumed in Lemmas 5 and 6.

Technically, Definition 6 says that the empirical probability measure of $(\beta _j)_{j=1}^J$ converges weakly (Klenke, Reference Klenke2020, Definition 13.12) to some probability measure. There is a large class of item-parameter sequences for which this property is natural. This includes sequences with a finite number of possible outcomes whose relative proportions stabilize as the sequence grows ( $F_\beta $ is then a discrete distribution), many low-discrepancy or quasi-Monte Carlo constructions (Niederreiter, Reference Niederreiter1992, Chapters 2 and 3), and the random-item setting considered in Appendix E of the Supplementary Material after conditioning on the realized item-selection sequence. Showing the convergence of such sums to an integral from first principles may be highly non-trivial. Variants of Definition 6 are a standard tool for dealing with limits of non-identically distributed random variables. See, e.g., the asymptotic framework of Section 2.2 in Claeskens and Hjort (Reference Claeskens and Hjort2008) and e.g. Ploberger and Krämer (Reference Ploberger and Krämer1996, Assumption A1), and many similar papers in econometrics).

Under parametric assumptions, $\beta _1, \ldots , \beta _{J}$ may be estimated empirically. Any observed stability among such estimates must be interpreted with caution, as Definition 6 is asymptotic, and the stability is required to continue indefinitely for the law of large numbers property to hold. Substantive knowledge of the item sequence, such as it being generated by drawing at random from an infinite item pool, is required to justify the law of large numbers property.

The law of large numbers property can be recast as a between-form stability property similarly to Lemma 4. We assume that $(\beta _j)$ is bounded, which is the case in practice.

Lemma 7 allows the psychometrician to think about the law of large numbers property as a between-form stability condition: as the number of items increases, averages of bounded continuous functions of the item parameters are asymptotically unchanged by removing a vanishing fraction of the items.

Since Definition 6 is asymptotic, it cannot be used to derive finite J bounds, such as those we considered earlier in the article. This may be improved with further assumptions. For example, when drawing items from an infinite item pool, one can derive finite-item guarantees based on the probability measure of the randomization device. We do not consider such refinements in the present article.

As a direct consequence of Definition 6 and Lemmas 5 and 6, we get the following conclusion. In its first statement, we recall that $\bar X_{J} = T_J(\theta ) + o_{\text {a.s.}}(1)$ if X is consistently estimable with exponential concentration (Proposition 1), which holds for bounded martingale factor models by Proposition 3 with $\alpha _j\equiv 1$ —and therefore also for bounded conditionally independent items by Proposition 2.

We now consider the conditionally independent case. For conditionally independent items, a representation of the form assumed in the upcoming result is available at the level of distributions by Lemma 1. The substantial assumption in the corollary is therefore that the parameter sequence fulfills the law of large numbers property.

In Appendix E of the Supplementary Material, we show that, after conditioning on the realized item-selection sequence, randomly selecting from a large pool of conditionally independent items yields an item sequence that is conditionally independent and parameterized by a sequence that fulfills the law of large numbers property, thereby providing a more substantial illustration of Corollary 1.

4.1 Unconditional and person-specific confidence intervals for the true score

We now consider confidence intervals for a respondent’s true score, computed based only on that respondent’s answers.

We only consider coverage rates under the unconditional probability measure, thereby taking into account that respondents are sampled at random, which is often the relevant perspective (see Chang & Stout, Reference Chang and Stout1993, Section 4, p. 50 and Holland, Reference Holland1990, p. 581 onward).

For a given $0 < \alpha < 1$ , let $z_\alpha =\Phi ^{-1}(1-\alpha /2)$ , where $\Phi $ is the cumulative distribution function of the standard normal distribution. Under the conditions of Proposition 4, we have that $ \lim _{J \xrightarrow {} \infty } P\!\left ( T_J(\theta ) \in \bar X_{J} \pm z_\alpha \hat V_{J}^{1/2}/\sqrt {J} \right ) = 1 - \alpha , $ where $\bar V_J(\theta ) = J^{-1} \sum _{j=1}^J \text {Var}\,_\theta X_j$ , and $\hat V_{J}$ is assumed to fulfill $\hat V_{J} = \bar V_{J}(\theta ) + o_P(\bar V_J(\theta ))$ . Finding such a consistent variance estimator may be possible under weak assumptions, for example, using methods such as those in Douglas (Reference Douglas1997) and Ramsay (Reference Ramsay1991). While this would be highly advantageous, it requires large J and a large number of respondents n. We here restrict attention to the less studied case where n is small and may even be $1$ , and use a variance estimator that is inconsistent though with a known positive bias.

For a fixed $0 < \alpha < 1$ , let $z_\alpha =\Phi ^{-1}(1-\alpha /2)$ and consider the confidence interval

$$\begin{align*}I_J = \left[ \bar X_J \pm z_\alpha s_J/\sqrt{J} \right], \quad \text{where} \quad s_J^2:=\frac{1}{J-1}\sum_{j=1}^J (X_j-\bar X_J)^2. \end{align*}$$

If the items were independent and identically distributed, this interval would have asymptotic nominal coverage, but under a martingale factor model (and hence for conditionally independent items) that satisfies the assumptions of Proposition 4, the interval is instead asymptotically conservative.

We revisit this result in Section 4.1.3 where we introduce assumptions that lead to a limit of $B_J(\theta )$ as J increases, yielding a concrete expression for the asymptotic coverage of $I_J$ .

4.1.1 Small J improvements: A hybrid approach using a Clopper–Pearson-type argument

While Theorem 2 is an asymptotic result, its empirical relevance concerns finite J. We now consider finite-J improvements. An important limitation of $I_J$ is that if $s_J = 0$ , the resulting singleton interval will almost surely not contain $T_J(\theta )$ . In short tests, it is highly likely that some candidates answer all questions either correctly or incorrectly, yielding no empirical variance. This indicates that $I_J$ should not be used in such cases.

Clopper and Pearson (Reference Clopper and Pearson1934) provide a finite-J conservative confidence interval for independent and identically distributed binary observations that is valid and non-singleton even when $s_J = 0$ . From Hoeffding (Reference Hoeffding1956, Section 5), Clopper–Pearson intervals are finite-J conservative as confidence intervals for the average success probability also in the non-identically distributed case. Under conditional independence, this means that the interval is a finite-J conservative confidence interval under $P_\theta $ for the true score. Taking expectations, this holds also under P. However, numerical experiments reported in Section 5.3 show that the Clopper–Pearson interval is excessively conservative unless J is large. In our simulations, the hybrid approach we now propose outperforms Clopper–Pearson for all J. Further, the hybrid approach can be used for non-binary items, and holds under martingale factor models and does not require the stricter assumption of conditional independence. We therefore do not further consider the Clopper–Pearson interval in the main text.

The hybrid interval starts by considering the case $X_1 = X_2 = \cdots = X_J = x$ . For $K=1$ , the Clopper–Pearson interval is $[0, 1 - c_{\alpha ,J}]$ when $x=0$ and $[c_{\alpha ,J}, 1]$ when $x=1$ (Thulin, Reference Thulin2014, p. 820), where $c_{\alpha ,J} = (\alpha /2)^{1/J}$ . These intervals have endpoints that are convex combinations of the common value x and the interval endpoints, and for $x = 0, 1$ are of the form $ [c_{\alpha ,J} x + (1-c_{\alpha ,J})\cdot 0, c_{\alpha ,J}x + (1-c_{\alpha ,J})\cdot 1] = [c_{\alpha ,J}x, c_{\alpha ,J}x+(1-c_{\alpha ,J})]. $ Let us therefore consider, for general $K \geq 1$ , intervals of the form $[c_{\alpha ,J} x + (1-c_{\alpha ,J})\cdot 0, c_{\alpha ,J}x + (1-c_{\alpha ,J})\cdot K] = [c_{\alpha ,J}x, x + (K-x)(1-c_{\alpha ,J})]$ . We obtain the following bound.

Proposition 8. Suppose that $X=(X_1,X_2,\ldots )$ follows a martingale factor model in $\theta $ and that $0\le X_j\le K$ almost surely for all j. Let $0<\alpha <1$ , $c_{\alpha ,J} = (\alpha /2)^{1/J}$ , and let $I_{\text {CP}} = \big [\,X_1 c_{\alpha ,J},\; X_1+(K-X_1)(1-c_{\alpha ,J})\,\big ]$ . Then $ P\!\left ( s_J=0,\; T_J(\theta )\notin I_{\text {CP}} \right )\le \alpha. $

The lengths of these intervals are $|I_J| = 2z_\alpha s_J/\sqrt {J}$ and $|I_{\text {CP}}| = K(1-c_{\alpha ,J})$ . Switching to $I_{\text {CP}}$ when $s_J=0$ creates a discontinuity in length as $s_J \to 0$ , where the length suddenly jumps to $K(1-c_{\alpha ,J})$ . This discontinuity reflects a finite-sample limitation of the normal approximation, and we therefore suggest a continuity-corrected hybrid interval

(7) $$ \begin{align} \tilde I_J = \left\{ \begin{matrix} \bar X_J \pm z_\alpha \tilde s_J/\sqrt{J}, &\quad \text{if } s_J> 0, \\ [X_1\,c_{\alpha,J},\; X_1+(K-X_1)(1-c_{\alpha,J})] & \quad \text{if } s_J = 0, \end{matrix} \right. \end{align} $$

where $\tilde s_J = \max (s_J,\eta _J^{\text {len}})$ is chosen to preserve continuity, with

$$\begin{align*}\eta_J^{\text{len}} = \inf \{ \eta \geq 0 : 2z_\alpha \eta/\sqrt{J} \geq K(1-c_{\alpha,J}) \} = \sqrt{J}|I_{\text{CP}}|/(2 z_\alpha). \end{align*}$$

Numerical experiments reported in Section 5.3 show that for the data generating processes we consider, the hybrid approach is superior to $I_J$ and to a discontinuous hybrid approach using $\tilde s_J = s_J$ , and that this holds for all J.

In the proof of Theorem 2, we show that, conditionally on each admissible value of $\theta $ , $ s_J^2 = \bar V_J(\theta ) + W_J(\theta ) + o_P(\bar V_J(\theta )). $ Moreover, the assumptions of the theorem imply $J\bar V_J(\theta ) \to \infty $ , whereas $ (\eta _J^{\text {len}})^2 = O(J^{-1}). $ It follows that $P(s_J < \eta _J^{\text {len}}) \to 0$ . Hence, for large J, the hybrid interval coincides with the asymptotic interval with probability tending to one, and therefore has the same asymptotic conservative coverage as $I_J$ . When $s_J < \eta _J^{\text {len}}$ , $I_J$ is shorter than $I_{\text {CP}}$ ; the continuity bridge defining $\tilde I_J$ corrects this short-length regime without the excessive conservatism of finite-J conservative intervals such as the Clopper–Pearson interval of Clopper and Pearson (Reference Clopper and Pearson1934).

4.1.3 Simplifications under consistent estimability

We now identify the limiting coverage probability for the confidence interval considered in Theorem 2 under consistent estimability, and a certain second moment convergence which can be justified in similar ways as convergent estimability.

Lemma 8. Suppose that $ T_J(\theta ) = T(\theta ) + o_P(1), $ that $J^{-1} \sum _{j=1}^J E_\theta [X_j^2] = T^{(2)}(\theta ) + o_P(1), $ and that $ J^{-1} \sum _{j=1}^J \big (E_\theta [X_j]\big )^2 = M^{(2)}(\theta ) + o_P(1). $ Also, assume $T^{(2)}(\theta )-T(\theta )^2>0$ with probability one. Then

(8) $$ \begin{align} B_J(\theta) = B(\theta) + o_P(1), \quad \text{where} \quad B(\theta) = \frac{T^{(2)}(\theta)-M^{(2)}(\theta)} {T^{(2)}(\theta)-T(\theta)^2}. \end{align} $$

Under eq. (8), and with the convention that $B(\theta )^{-1/2}=\infty $ on $\{B(\theta )=0\}$ , the coverage probability from Theorem 2 converges to $1-\alpha + p_B$ , where

(9) $$ \begin{align} p_B = 2 P \left( z_\alpha \leq Z \leq z_\alpha B(\theta)^{-1/2} \right). \end{align} $$

This is close to $1-\alpha $ when $B(\theta )$ is concentrated around $1$ . A simple numerical illustration is given in Section 5.2.

4.2 Convergence of mean-score correlations

Here, we consider the empirical practice of using mean scores in correlation-based methods, such as factor analysis or structural equation models (see Bollen, Reference Bollen1989, for a textbook introduction to these models). This practice may be motivated by linear factor models of the form $(X_1, \ldots , X_J)' = \Lambda \theta + u$ , where $\theta $ is a vector, $\Lambda $ is a full column rank matrix, and $\theta $ is uncorrelated with u. Then, for suitably chosen linear transformations of the form $H \cdot (X_1, \ldots , X_J)'$ , such as a left inverse of $\Lambda $ , one obtains variables that are linear transformations of $\theta $ plus an error term. Under additional assumptions making this error small, the correlation matrix of $H \cdot (X_1, \ldots , X_J)'$ is therefore close to that of $\theta $ , which motivates this empirical practice.

This argument does not hold for item response models. In contrast, non-linear factor scores are required to consistently approximate $\theta $ . Yet we may inquire into what this common empirical practice approximates from a nonparametric IRT perspective, and we will show in the numerical illustrations that using mean scores as input into covariance models may work well, and even be superior to using misspecified parametric estimation methods.

Several models, including linear regression, factor analysis, and structural equation models are, under standard assumptions, dependent only on covariances and means (Bollen, Reference Bollen1989, Chapter 8). For standardized variables, the consistency of these methods only depends on correlations, which depend on pairs of items. We therefore consider only two scales in the following result.

Proposition 9 shows that inputting sum or mean scores into correlation methods will be approximately the same as inputting functions of their underlying latent variables. Proposition 9(1) and (2) show that this holds when comparing the correlation of two mean scores, and Proposition 9(3) and (4) show that this holds when inputting a mix of, for example, continuous explanatory variables and mean scores. Under convergent estimability, Proposition 9(2) and (4) identify the limit of mean-score correlations. Importantly, Proposition 9(1) and (3) identify what mean-score correlations approximate without assuming convergent estimability, that is, without making an assumption on a theoretical infinite item sequence.

Inference in structural equation models typically relies on fourth-order moments (see, e.g., Satorra, Reference Satorra1989), and extensions of Proposition 9 lead to the convergence properties also of higher-order moments. We consider this outside the scope of the present article.

As an illustration of Proposition 9, we study the phenomenon that mean-score correlations have been reported to have an absolute value less than the true latent correlation under normality, as also seen in the upcoming numerical illustrations. For large J, this effect is explained by Proposition 9: suppose $(\theta _X, \theta _Y)$ are jointly normal. Then $|\rho (T^{(X)}(\theta _X), T^{(Y)}(\theta _Y))| \leq |\rho (\theta _X, \theta _Y)|$ , as $|\rho (\theta _X, \theta _Y)|$ is the maximal correlation coefficient under normality (see, e.g., the discussion on maximal correlation in Klaassen & Wellner, Reference Klaassen and Wellner1997, Section 2). If $(\theta _X, \theta _Y)$ has a non-normal copula, there is no simple monotonic relationship between the magnitude of $\rho (T^{(X)}(\theta _X), T^{(Y)}(\theta _Y))$ and $\rho (\theta _X, \theta _Y)$ without further assumptions.

4.3 Nonparametric factor scores

In IRT, estimating $\theta $ is of practical importance, though it may be both computationally and theoretically difficult. Computationally, the difficulty lies in the non-linear optimization problems required, while theoretically, the difficulty lies in identifying realistic assumptions resulting in robust methods which perform well in realistic sample and item number settings.

Alternatively, we may consider the mean score as an ordinal proxy for $\theta $ , and map its empirical ranks to a chosen latent scale. This idea appears, for example, as ordinal ability estimation in Douglas (Reference Douglas1997) and Ramsay (Reference Ramsay1991) and is closely related to the equipercentile equating argument (see, e.g., Lord, Reference Lord1982). In our setting, convergent estimability provides a direct justification for such rank-based scoring rules under weak assumptions, without requiring estimation of item response curves.

Suppose X is convergently estimable to $T(\theta )$ and $\vartheta \mapsto T(\vartheta )$ is a strictly increasing map, and suppose $\theta $ is univariate and has a known distribution function G. Then for each person $i \geq 1$ in a random sample,

(10) $$ \begin{align} \bar X_{i,J} \;=\; T(\theta_i) + o_P(1). \end{align} $$

We may estimate $F(t) = P(T(\theta ) \leq t)$ with the empirical distribution function $F_{n,J}(t)$ based on independently sampled mean scores. Let $P[\cdot ]_\theta $ mean a probability, where only $\theta $ is considered random. When $n,J$ increase, $F_{n,J}(\bar X_{i,J})$ is close to

$$\begin{align*}F(T(\theta_i)) = P[T(\theta) \leq T(\theta_i)]_\theta = P[\theta \leq \theta_i]_\theta = G(\theta_i), \end{align*}$$

so that the NPFS $\hat \theta _i = G^{-1}(F_{n,J}(\bar X_{i,J}))$ is close to $\theta _i$ . To avoid boundary issues in evaluating $G^{-1}$ , we re-scale the standard empirical distribution function and work with

(11) $$ \begin{align} F_{n,J}(x):=\frac{1}{n+1}\sum_{r=1}^n I \{\bar X_{r,J}\le x\}\in\Big[\frac{1}{n+1},\frac{n}{n+1}\Big] \end{align} $$

for $x \geq \min _r \bar X_{r,J}$ .

This method, which can be refined in various directions that we do not study here, is a large-J, large-n method. It does not require knowledge of the limiting function T, and it avoids item-parameter estimation and optimization once the latent distribution G has been specified.

While G is assumed known, this is also the case in most other methods, such as parametric IRT. An advantage of the analytical simplicity of $\hat \theta $ is that it provides a transparent analysis of what is estimated if G is misspecified. Indeed, suppose $\theta _i$ has distribution function H, though $\hat \theta _i$ erroneously assumes $\theta \sim G$ . Then by Proposition 10, we get $\hat \theta _i = G^{-1}(F_{n,J}(\bar X_{i,J})) = G^{-1} \circ H \circ H^{-1}(F_{n,J}(\bar X_{i,J})) = G^{-1} \circ H (\theta _i) + o_P(1) = \tilde \theta _i + o_P(1)$ , where $\circ $ is function composition and $\tilde \theta _i = G^{-1} \circ H (\theta _i)$ . Interestingly, since $P(\tilde \theta _i \leq t) = P( G^{-1} \circ H (\theta _i) \leq t ) = P( \theta _i \leq H^{-1} \circ G(t)) = H \circ H^{-1} \circ G(t) = G(t)$ , the erroneous assumption that $\theta _i \sim G$ in fact reaches a transformed latent variable where this assumption is fulfilled.

Sample covariances or correlations computed from the resulting NPFSs can be input into covariance methods and will, under conditions similar to Proposition 9, converge to the true covariance or correlation between latent variables. Section 5.4 provides a numerical illustration of this idea. In contrast to mean-score correlations, the resulting estimate is consistent.

As an illustration, recall the martingale factor model with generalized partial credit marginals considered in Remark 5 under the additional assumption that its parameter sequence fulfills the law of large numbers property (Definition 6). In contrast to standard parametric latent variable estimates, the NPFS of Proposition 10 requires neither nuisance parameter estimation nor optimization for its computation.

5.1 A simple analytical illustration for calculating a limiting mean-score correlation

Here, we consider a simple model class with parameters drawn from an infinite item pool and having a particularly simple limiting true score.

Figure 1(a) shows the limiting function T of Proposition 11 for $c = 0.9, \lambda = 0.5 = \psi $ . Its curvature induces a small bias into the mean-score correlations as an estimator of the latent correlation (see Figure 1(c)). In contrast, with $c = 2.5$ , the limiting function is almost exactly linear, see Figure 1(b), which in turn induces an almost consistent mean-score correlation, see Figure 1(d).

Figure 1

Illustrations for Proposition 11. For panels (c) and (d), $(\theta _X, \theta _Y)$ is jointly normal with specified correlation $\rho $ .

Figure 1 Long description

The image consists of four panels arranged in a two by two grid.

Top-left panel A. The x axis is labeled with the Greek letter theta and ranges from negative 3 to 3. The y axis ranges from 0.0 to 1.0. A black curve shows a sigmoidal S-shaped increase, starting near 0 at x equals negative 3 and approaching 1.0 at x equals 3.

Top-right panel B. The x axis is labeled with the Greek letter theta and ranges from negative 3 to 3. The y axis ranges from 0.2 to 0.8. A black line shows a perfectly linear increase from the bottom-left to the top-right of the plot area.

Bottom-left panel C. The x axis is labeled rho open parenthesis xi comma eta close parenthesis and ranges from negative 1.0 to 1.0. The y axis is labeled rho open parenthesis X bar comma Y bar close parenthesis and ranges from negative 1.0 to 1.0. A red dashed line labeled Truth represents a perfect diagonal from negative 1 to 1. A solid black line labeled Estimate follows a similar path but has a slightly shallower slope, crossing the red line at the origin.

Bottom-right panel D. This panel uses the same axes and labels as panel C. The red dashed Truth line and the solid black Estimate line are nearly identical, showing a much tighter fit and a steeper slope than in panel C, almost perfectly overlapping along the diagonal.

5.2 The limiting bias of the naive confidence interval for the true score

Consider items of the form $X_j = I \{ \varepsilon _j \leq \xi \}$ , where $\xi \sim N(0,1)$ and $\varepsilon _j \sim N(\mu _j, \sigma _j^2)$ , and where $\xi , \varepsilon _1, \ldots $ are mutually independent. Before any items are drawn, we sample $\mu _1, \ldots , \mu _{J} \sim U[-1,0], \sigma _1, \ldots , \sigma _{J} \sim U[0.2, 10/6]$ independently, which is as if we sample from a random item pool. While written as a threshold model, the items follow a two-parameter probit IRT model by the equivalence shown in, for example, Takane and De Leeuw (Reference Takane and De Leeuw1987).

We study the coverage of the confidence interval from Theorem 2. Lemma 8 identifies the asymptotic coverage rate as $1-\alpha + p_B$ for $p_B$ as in eq. (9). We numerically approximate $p_B$ to be $0.01$ , giving an asymptotic coverage rate of $96 \%$ , compatible with the simulation illustration given next.

5.3 A numerical illustration of the person-specific confidence intervals

We extend the numerical illustration from Section 5.2 to polytomous items using a normal-threshold model, which reduces to the binary setup when $K=1$ . Let $\xi \sim N(0,1)$ , and independently draw $\sigma _j\sim U[0.2,10/6]$ and $\eta _j\sim N(0,\sigma _j^2)$ for each $j \geq 1$ , independently of $\xi $ . Given ordered thresholds $\tau _{j,1}<\cdots <\tau _{j,K}$ , we let $ X^{(K)}_j = \sum _{k=1}^{K} I\{\xi +\eta _j>\tau _{j,k}\}\in \{0,1,\ldots ,K\}. $ Thresholds are chosen to avoid nearly empty categories: For each item j, we draw a central location $m_j\sim U[-1,0]$ and set the K thresholds to be equally spaced on $[m_j-c,\;m_j+c]$ with $c=2.5$ . We consider $2,3,5,7,$ and $9$ response categories. The resulting items are conditionally independent. In Appendix D.1 of the Supplementary Material, we show that the hybrid confidence interval is well-behaved also under a martingale factor model which is not locally independent when $K \geq 2$ .

We let J vary from $2$ to $200$ . For each J and for a fixed parameter sequence $(m_j), (\sigma _j)$ , we compute coverage rates with 10,000 repetitions for the asymptotic confidence interval in Theorem 2, as well as the hybrid interval described in Section 4.1.1. Since the coverage rates depend on the parameter draws of $(m_j), (\sigma _j)$ , we repeat this whole process $500$ times (each time generating new $(m_j), (\sigma _j)$ parameters as described above). We then obtain a distribution of coverage rates for each J and each method.

The mean, and the lower and upper $2.5 \%$ quantiles, of the coverage rate distributions of the hybrid and asymptotic method, as well as the Clopper–Pearson method for $K=1$ , are shown in Figure 2, and tabulated in Appendix F.3 of the Supplementary Material.

Figure 2

Coverage rates as functions of J. The mean and the lower and upper $2.5 \%$ quantiles of the rate distributions are indicated.

Figure 2 Long description

A two-panel line graph. Both panels share a horizontal x-axis labeled J ranging from 0 to 200 and a vertical y-axis labeled Coverage rate ranging from 0.85 to 1.00. A horizontal dashed line is fixed at 0.95.

Panel A is titled Unconditional coverage probabilities of the 95 percent hybrid confidence interval, and the Clopper-Pearson interval for K = 1. The data lines start near 1.00 at low J and decrease. Most lines (K = 2, 4, 6, 8) drop sharply before stabilizing just above 0.95. The line for K = 1 (asterisks) stays slightly below 0.95. A dashed line with open circles for K = 1 (C P) remains higher, around 0.97. Shaded regions represent the lower and upper 2.5 percent quantiles.

Panel B is titled Unconditional coverage probabilities of the 95 percent asymptotic conservative interval I sub J. The data lines start at low coverage rates (below 0.85) and increase asymptotically toward 0.95 as J increases. The line for K = 1 (asterisks) is the lowest, while lines for higher K values (2, 4, 6, 8) reach the 0.95 threshold faster. Shaded regions indicate the distribution spread.

Legends in both panels identify:

* Asterisk: K = 1

* Blue circle: K = 2

* Red triangle: K = 4

* Green square: K = 6

* Purple diamond: K = 8

* Open circle (Panel A only): K = 1 (C P)

In the present simulation setup, the hybrid method performs well already from around $J \geq 13$ . In contrast, the Clopper–Pearson method and the uncorrected asymptotic interval $I_J$ perform poorly: Clopper–Pearson is excessively conservative, while $I_J$ has poor performance when $s_J$ is small, which occurs at high rates for small J (see Appendix F.3 of the Supplementary Material).

The coverage rate of the hybrid method is slightly conservative for sufficiently large J and for the values of K considered here. For J = 5,000 and $K = 1$ , the mean coverage rate was found to be $95.89 \%$ , compatible with the asymptotic bias calculated in Section 5.2.

5.4 Latent correlation estimates: Mean-score correlations, NPFS correlations, and ordinal factor analysis—A simple numerical illustration

We consider a standard threshold-based ordinal factor model with binary items, given by

(12) $$ \begin{align} X_k = I \left\{ \lambda_k^{(X)} \xi + \epsilon^{(X)}_{k}> \tau^{(X)}_{k} \right\}, \quad Y_j = I \left\{ \lambda_j^{(Y)} \eta + \epsilon^{(Y)}_{j} > \tau^{(Y)}_{j} \right\}, \end{align} $$

for $k=1,\ldots ,J$ and $j=1,\ldots ,J$ (see Jöreskog, Reference Jöreskog1994, for background). The error terms $(\epsilon ^{(X)}_{k})$ and $(\epsilon ^{(Y)}_{j})$ are assumed mutually independent and also independent of the latent variables $(\xi ,\eta )$ . Under latent normality, this ordinal factor model is equivalent to a normal-link IRT model (Takane & De Leeuw, Reference Takane and De Leeuw1987).

In applied work, one is often interested in the latent correlation $\rho = \rho (\xi ,\eta )$ , perhaps as input to estimate a structural equation model. Under joint normality of $(\xi , \eta )$ and the error terms, $\rho $ can be estimated by fitting a factor model to the empirical tetrachoric correlations (Jöreskog, Reference Jöreskog1994). We call this the ordinal factor model estimator. In this illustration, we compare three estimates of $\rho $ : the ordinal factor model estimator, mean-score correlations, and correlations between the NPFSs of Section 4.3.

The NPFS requires knowledge of the marginal distributions of $\xi $ and $\eta $ . These marginals will be kept standard normal, and the non-normal scenarios we consider change only the copula of $(\xi , \eta )$ . The mean-score correlation does not rest on explicit distributional assumptions. In contrast, the ordinal factor model estimator requires latent joint normality and is generally inconsistent outside this assumption (Foldnes & Grønneberg, Reference Foldnes and Grønneberg2020, Reference Foldnes and Grønneberg2022). Even when all empirically testable normality assumptions for the latent variables appear to hold, latent normality does not need to hold, and model-based correlation estimates may deviate substantially from their population values under non-normality (Grønneberg & Foldnes, Reference Grønneberg and Foldnes2024). Latent non-normality is empirically relevant, as illustrated in Grønneberg and Foldnes (Reference Grønneberg and Foldnes2024), who applied a test for underlying normality from Grønneberg et al. (Reference Grønneberg, Moss and Foldnes2020) to a collection of psychometric data, and found that latent normality was strongly rejected for the seven psychometric datasets they considered. Testing for latent non-normality is not commonly done in practical applications even when using estimators that rest on this assumption.

For the simulations, we fix factor loadings $\lambda _j^{(X)}=\lambda _j^{(Y)}=0.5$ and set the error-term standard deviations also to $0.5$ . Thresholds $\tau _j^{(X)}$ and $\tau _j^{(Y)}$ are drawn once, uniformly from $[-0.9,0.9]$ , and kept constant across replications. In contrast to Section 5.3, we do not consider re-drawing the thresholds and redoing all the simulations per drawn threshold configuration, as this would be computationally too demanding due to the tetrachoric correlation calculations.

Extending the threshold range to $[-2.5,2.5]$ would significantly improve the performance of the mean-score correlation approach, as seen in Section 5.1. We did not consider this setup here, as it makes some item pairs so strongly dependent that the tetrachoric correlations become numerically unstable.

We consider $\rho $ either equal to $0.8$ or $0.5$ . A latent correlation of $0.5$ is common in practice, while a strong correlation of $0.8$ emulates cases where tests of unidimensionality are relevant, that is, that it is of interest to test whether $\eta $ and $\xi $ coincide. Further, we consider non-normal latent variables as described below, and normal latent variables. We consider normally distributed latent data only for $\rho = 0.8$ , as the mean-score correlations with normal latent variables were investigated systematically in Section 5.1, and the model-based approach is well-known to perform well under normality.

For the non-normal scenarios, $(\xi , \eta )'$ has standard normal marginals and the non-normal Joe copula with parameters chosen via the R-package covsim (Grønneberg et al., Reference Grønneberg, Foldnes and Marcoulides2022) to achieve the target correlation $\rho $ ( $0.5$ or $0.8$ ). This type of copula-induced non-normality is known to detrimentally affect ordinal factor analysis (Foldnes & Grønneberg, Reference Foldnes and Grønneberg2020). Furthermore, the item error terms $\epsilon ^{(X)}_{k}$ are drawn from a $\chi ^2_2$ -distribution, and $\epsilon ^{(Y)}_{j}$ from a log-normal distribution, both subsequently scaled and shifted to have a mean of zero and a standard deviation of $0.5$ .

We generate data for $n=3{,}000$ respondents and set $J=50$ items per scale. Each distributional configuration is replicated $1{,}000$ times. For $50$ binary items per scale, computing the $4{,}950$ unique tetrachoric correlations for each replication is computationally intensive, and the resulting correlation matrix is frequently non-positive-definite (occurring in $37.4\%$ of the replications in the non-normal, $\rho =0.80$ case; for comparison, the corresponding rates are $0.4\%$ under normality with $\rho =0.8$ and $0.3\%$ in the non-normal, $\rho =0.5$ case). When a non-positive-definite matrix occurs, we apply the nearPD function from the R Matrix package (Bates et al., Reference Bates, Maechler and Jagan2024) to find the closest positive definite matrix before passing it to lavaan (Rosseel, Reference Rosseel2012). For ordinal data estimated via tetrachoric correlations, either diagonally weighted least squares (DWLS) or unweighted least squares (ULS) are recommended (Li, Reference Li2016; Forero et al., Reference Forero, Maydeu-Olivares and Gallardo-Pujol2009) when fitting factor models. We choose the ULS option in lavaan, as the DWLS option requires the additional computation of the asymptotic covariance matrix of the vectorized tetrachoric correlation matrix.

The mean absolute errors (MAEs) for each condition are shown in Table 1. Sampling distributions are shown in Appendix F.1 of the Supplementary Material; the appendix figures illustrate the bias directions discussed below.

Table 1

Mean absolute errors for three estimators of the latent correlation based on 1,000 replications with $n=3{,}000$ and $J=50$ items per scale

Table 1 Long description

The table consists of four columns and four rows including the header.

Header Row:

- Estimator

- Normal, rho = 0.8

- Non-normal, rho = 0.8

- Non-normal, rho = 0.5

Row 1 (Ordinal factor model):

- Normal, rho = 0.8: 0.006 (bold-faced)

- Non-normal, rho = 0.8: 0.076

- Non-normal, rho = 0.5: 0.034

Row 2 (Nonparametric factor score):

- Normal, rho = 0.8: 0.042

- Non-normal, rho = 0.8: 0.022

- Non-normal, rho = 0.5: 0.014 (bold-faced)

Row 3 (Mean-score):

- Normal, rho = 0.8: 0.043

- Non-normal, rho = 0.8: 0.018 (bold-faced)

- Non-normal, rho = 0.5: 0.025

A note indicates that bold-faced values represent the best performing method for each condition.

Note: The best performing method is bold-faced.

We first consider the high-correlation scenario ( $\!\! \rho =0.8$ ). Under latent normality, the ordinal factor model estimator performs as expected (MAE = 0.006) and outperforms the sum-score estimator and the NPFS, both of which exhibit downward bias (MAE = 0.042 and 0.043, respectively). However, when the normality assumption is violated, ordinal factor analysis substantially overestimates the latent correlation (MAE = 0.076). In contrast, the mean-score estimator, while remaining slightly conservative, is far more accurate under non-normality (MAE = 0.018), and slightly more accurate than the NPFS (MAE = 0.022).

In the moderate correlation scenario ( $\rho = 0.5$ ) under non-normality, ordinal factor analysis maintains a positive bias (MAE = 0.034), whereas the sum-score estimator is slightly negatively biased but considerably more precise (MAE = 0.025). The best performing method is the NPFS, which is close to unbiased (MAE = 0.014).

Although this small illustration has limited external validity, it demonstrates that even when a binary factor model structure is correctly specified, the simple mean-score correlation can outperform a standard parametric estimator when its underlying distributional assumptions are violated, and that the NPFS estimator either performs similarly or better than the mean-score correlation. We note here that NPFS is applied with the correct latent variable distribution, which is in contrast to the ordinal factor model, and that the mean score does not directly rely on such a specification. Finally, the performance of the mean-score correlation is driven by the formula of the true score, which except under strong assumptions cannot be estimated (though with strong assumptions, it can be estimated as illustrated in Appendix F.2 of the Supplementary Material).