2.1 Latent class analysis

We consider N respondents answering J items. Let $Y_{ij}$ denote the response from respondent i to item j. To simplify the discussion, we assume that the responses are binary, that is, $Y_{ij} \in \{0, 1\}$ . As discussed in Section 5, the proposed method can be extended to data with polytomous items.

Suppose the data follow a latent class model with K classes. Let $\xi _i \in \{1, 2, \ldots , K\}$ indicate the latent class membership of respondent i. The model assumes that $\xi _i$ , $i = 1, \ldots , N$ , are independent and identically distributed (i.i.d.) random variables, following a categorical distribution, with $P(\xi _i = k) = \nu _k$ , where $\nu _1, \ldots , \nu _K$ are unknown parameters satisfying $\nu _k \geq 0$ and $\nu _1 + \cdots + \nu _K = 1$ . Moreover, the model assumes that an individual’s response probabilities depend only on their latent class membership. That is,

$$\begin{align*}P(Y_{ij} = 1\vert \xi_i = k) = \beta_{jk},\end{align*}$$

where $\beta _{jk}$ , $j = 1, \ldots , J$ and $k=1, \ldots , K$ , are unknown parameters. Finally, the local independence assumption is adopted, which says that for each respondent i, $Y_{i1}$ , $\ldots $ , $Y_{iJ}$ are conditionally independent given $\xi _i$ .

For the parsimony of notation, we use $\boldsymbol {\nu } = (\nu _1, \ldots , \nu _K)^\top $ , $\mathbf B = (\beta _{jk})_{J\times K}$ , and $\boldsymbol {\beta }_j = (\beta _{j1}, \ldots , \beta _{jK})^\top $ to denote the parameter vectors or matrix, and $\boldsymbol {\nu }^{\ast} = (\nu _1^{\ast}, \ldots , \nu _K^{\ast})^\top $ , $\mathbf B^{\ast} = (\beta _{jk}^{\ast})_{J\times K}$ , and $\boldsymbol {\beta }_j^{\ast} = (\beta _{j1}^{\ast}, \ldots , \beta _{jK}^{\ast})^\top $ the corresponding true values. Under the assumptions of the latent class model, the marginal log-likelihood function takes the form

$$ \begin{align*} l(\boldsymbol{\nu}, \mathbf B) = \sum_{i=1}^N\log\Big(\sum_{k=1}^K \nu_k\prod_{j=1}^J \beta_{jk}^{Y_{ij}}(1-\beta_{jk})^{1-Y_{ij}}\Big). \end{align*} $$

Classical LCA estimates the unknown model parameters by the maximum likelihood estimator (MLE), that is,

(1) $$ \begin{align} (\hat{\boldsymbol{\nu}}, \hat{\mathbf B}) = \arg\max_{(\boldsymbol{\nu}, \mathbf B)} ~l(\boldsymbol{\nu}, \mathbf B). \end{align} $$

This optimization is typically solved using an EM algorithm. Under correct model specification and mild regularity conditions, the MLE is root-N consistent, in the sense that

$$ \begin{align*} \Vert \hat{\boldsymbol{\nu}} - {\boldsymbol{\nu}}^{\ast}\Vert = O_p(1/\sqrt{N}) \text{ and } \Vert \hat{\mathbf B} - {\mathbf B}^{\ast} \Vert_F = O_p(1/\sqrt{N}), \end{align*} $$

where $\Vert \cdot \Vert $ denotes the Euclidean norm and $\Vert \cdot \Vert _F$ denotes the matrix Frobenius norm.Footnote ¹

When the number of latent classes K is unknown, the optimal number of latent classes is chosen by fitting models with increasing numbers of classes and selecting the one that minimizes the BIC value. We refer readers to Chapter 4 of Collins and Lanza (Reference Collins and Lanza2010) for a review of the inferential foundations of the classical LCA. For ease of exposition, we treat the number of classes K as known when motivating and describing the proposed refinement method in Sections 2.2 and 2.3 below.

2.2 Sparse LCA

Interpreting the latent classes is an essential step for LCA. This step is based on $\hat {\mathbf B}$ , which characterizes the item-response behaviors of all the learned latent classes. However, when J or K is moderately large, interpreting the latent classes can be a challenge, due to the large number of parameters in $\hat B$ . To tackle this challenge, Chen et al. (Reference Chen, Li, Liu and Ying2017) proposed a sparse latent class model, which assumes that the number of distinct item-response probability levels for each item can be smaller than the total number of classes.

More specifically, let $m_j^{\ast} =\text {card}(\boldsymbol {\beta }^{\ast}_{j})$ denote the true number of distinct item-response probability levels for item j. The sparse latent class model assumes that $m_j^{\ast} < K$ for many items. However, even when the true model is sparse, due to the randomness in data, each $\hat {\boldsymbol {\beta }}_{j}$ has K distinct values with probability tending to 1 as the sample size goes to infinity, and thus, cannot recover the sparse pattern in the true model. To tackle this issue, Chen et al. (Reference Chen, Li, Liu and Ying2017) proposed a regularized estimation method that simultaneously learns the number of distinct item-response probability levels, denoted by $\tilde m_j$ , and estimates the item-response probabilities, denoted by $\bar {\mathbf B} = (\bar \beta _{jk})_{J\times K}$ , satisfying that the number of distinct levels in $(\bar \beta _{j1}, \ldots , \bar \beta _{jK})^\top $ is $\tilde m_j$ , where $\tilde m_j$ can be strictly smaller than K. Consequently, $\bar {\mathbf B}$ tends to have substantially fewer parameters than those in $\hat {\mathbf B}$ , the MLE for classical LCA, and thus, tends to be much easier to interpret. Consistency for both the model selection and parameter estimation has been established for this regularized estimator.

Despite the interpretational advantage and theoretical guarantees, this regularized estimator requires maximizing a regularized log-likelihood function, which involves complex EM cycles in which each M-step involves solving multiple non-smooth optimization problems. Such a heavy computational burden renders the method impractical for many real-world applications.

Example 1. Table 1 gives a toy example of the true item-response probability matrix $\mathbf B$ for a sparse latent class model, with $J= 6$ items and $K=4$ latent classes. In this example, the elements of $\boldsymbol {\beta }_j$ have two distinct levels for items 1–5 and three distinct levels for item 6. If such a sparse pattern can be recovered from the data, one can infer that the first latent class has the lowest item-response probabilities across all items, while the fourth latent class has the highest item-response probabilities. In addition, each of the second and third classes has a combination of the lowest and highest item-response probabilities for items 1–5 and a median level of item-response probability for item 6. In a classical latent class model, however, each item would typically be assigned four distinct response-probability levels, even when the true underlying structure is sparse, making the resulting classes harder to interpret.

Table 1

An example of the item-response probabilities for a sparse latent class model

Table 1 Long description

The table consists of five columns. The first column is labeled Item, followed by four columns under the header Latent class, numbered 1, 2, 3, and 4.

* Item 1: Latent class 1 is 0.10, class 2 is 0.10, class 3 is 0.10, and class 4 is 0.90.

* Item 2: Latent class 1 is 0.20, class 2 is 0.20, class 3 is 0.80, and class 4 is 0.80.

* Item 3: Latent class 1 is 0.15, class 2 is 0.75, class 3 is 0.15, and class 4 is 0.75.

* Item 4: Latent class 1 is 0.30, class 2 is 0.30, class 3 is 0.30, and class 4 is 0.85.

* Item 5: Latent class 1 is 0.25, class 2 is 0.25, class 3 is 0.90, and class 4 is 0.90.

* Item 6: Latent class 1 is 0.20, class 2 is 0.55, class 3 is 0.55, and class 4 is 0.95.

2.3 Proposed refinement procedure

We propose a refinement procedure to estimate a sparse latent class model, which remains statistically consistent while avoiding the high computational burden of the regularized estimator. Like the regularization approach in Chen et al. (Reference Chen, Li, Liu and Ying2017), our goal is to find an estimate of the latent class model, denoted by $(\bar {\boldsymbol \nu }, \bar {\mathbf B})$ , such that $\bar {\mathbf B}$ is sparse, in the sense that each row of $\bar {\mathbf B}$ has a small number of distinct values, that is, each item has a small number of unique item-response probability levels.

The proposed estimator has two building blocks, including (1) a stepwise search algorithm for exploring sparse models and (2) an EBIC for comparing models with different sparsity levels, both of which rely on an item-specific pseudo-likelihood function constructed based on the initial model estimate from the first stage. In what follows, we first define the pseudo-likelihood, then describe the stepwise search algorithm, and finally the EBIC.

2.3.1 Item-specific pseudo-likelihood

Let

$$ \begin{align*} \hat \gamma_{ik} = \frac{\hat\nu_k\prod_{j=1}^J \hat\beta_{jk}^{Y_{ij}}(1-\hat\beta_{jk})^{1-Y_{ij}}}{\sum_{h=1}^K\hat\nu_h\prod_{j=1}^J \hat\beta_{jh}^{Y_{ij}}(1-\hat\beta_{jh})^{1-Y_{ij}}} \end{align*} $$

be the posterior probability of respondent i belonging to the kth latent class evaluated under the MLE for the classical LCA, and we denote $\hat {\boldsymbol {\Gamma }} = \{\hat {\gamma }_{ik}\}_{N\times K}$ . The pseudo-likelihood function for item j is defined as

(2) $$ \begin{align} Q_j(\boldsymbol{\beta}_j; \hat{\boldsymbol{\Gamma}}) = \sum_{i=1}^N \sum_{k=1}^K \hat \gamma_{ik} \left( Y_{ij} \log(\beta_{jk}) + (1 - Y_{ij}) \log(1 - \beta_{jk}) \right). \end{align} $$

This is a weighted Bernoulli log-likelihood with the latent class memberships replaced by the corresponding posterior probabilities $\hat {\gamma }_{ik}$ ; for brevity, we refer to it as the pseudo-likelihood throughout. We should also note that this function is the objective function for $\boldsymbol {\beta }_j$ in the M-step of the EM algorithm for solving (1).

This item-specific pseudo-likelihood enables independent model selection item by item, substantially reducing the complexity of the model selection problem.

2.3.2 Stepwise search algorithm

For each item j and each possible number of distinct item-response probability levels $m_j \in \{1, \dots , K\}$ , ideally, we hope to solve the following maximization problem:

(3) $$ \begin{align} \tilde{\boldsymbol{\beta}}_{j}^{(m_j)} = \arg\max_{\boldsymbol{\beta}_{j}} ~Q_j(\boldsymbol{\beta}_{j}; \hat{\boldsymbol{\Gamma}}), \quad \text{s.t. } \text{card}(\boldsymbol{\beta}_{j}) = m_j. \end{align} $$

Recall that $\text {card}(\boldsymbol {\beta }_{j})$ denotes the number of distinct values within the j-th item-response probability vector $\boldsymbol {\beta }_{j}$ . In other words, $\tilde {\boldsymbol {\beta }}_{j}^{(m_j)}$ re-estimates the response probabilities for item j, under the constraint that the item has $m_j$ distinct item-response probability levels. The solutions $\tilde {\boldsymbol {\beta }}_{j}^{(1)}$ , …, $\tilde {\boldsymbol {\beta }}_{j}^{(K)}$ are then compared using the EBIC introduced below. Unfortunately, (3) becomes computationally infeasible when K is sufficiently large, due to the combinatorial nature of the discrete constraint $\text {card}(\boldsymbol {\beta }_{j}) = m_j$ . For example, when $K = 10$ , solving (3) for $m_j = 1, \ldots , 10$ requires comparing a total of 115,975 models, which is computationally demanding.

To tackle this computational challenge, we instead use a stepwise search algorithm that leverages the ordering of the initial estimate $\hat {\boldsymbol \beta }_j$ . The algorithm iteratively identifies partitions of the entries of $\boldsymbol {\beta }_j$ such that the parameters in the same subset of the partition share the same value. More specifically, let $\mathcal C = \{C_1, \ldots , C_{m_j}\}$ be a partition of $\{1, \ldots , K\}$ , satisfying that (1) each subset $C_l$ is nonempty, (2) $C_l$ s are mutually exclusive, and (3) the union of $C_l$ is $\{1, \ldots , K\}$ . A partition defines a constrained parameter space for $\boldsymbol {\beta }_j$ :

$$\begin{align*}\Theta_j(\mathcal C) = \{\boldsymbol{\beta}_j:\beta_{jk}=\beta_{jk'}, \forall\, k,k'\in C_l,\ l=1, \ldots, m_j, \text{ and } \beta_{jk} < \beta_{jk'}, \text{ if } k\in C_l, k' \in C_{l'}, l < l'\}, \end{align*}$$

satisfying $\text {card}(\boldsymbol {\beta }_{j}) = m_j$ for all $\boldsymbol {\beta }_{j} \in \Theta _j(\mathcal C)$ . As $\hat {\boldsymbol \beta }_j$ is a statistically consistent estimate of ${\boldsymbol \beta }_j^{\ast}$ , the ordering of the estimated parameters in $\hat {\boldsymbol \beta }_j$ is useful when we search for the partition based on the true parameters $\boldsymbol {\beta }_j^{\ast}$ .

For ease of exposition, we assume that $\hat {\boldsymbol \beta }_j$ has K distinct entries, which is typically the case in practice and asymptotically guaranteed. Let $\pi _j$ be a permutation of $\{1,\ldots ,K\}$ that defines the order statistic for $\hat \beta _{j1}$ , …, $\hat \beta _{jK}$ , satisfying that $\hat \beta _{j,\pi _j(1)} < \cdots < \hat \beta _{j,\pi _j(K)}$ . The algorithm starts with $\tilde {\boldsymbol {\beta }}_{j}^{(K)} \ = \ \hat {\boldsymbol \beta }_j$ . This solution defines a partition $\mathcal C_j^{(K)}= \big \{\{\pi _j(1)\},\ldots ,\{\pi _j(K)\}\big \}$ , where each subset contains exactly one latent class and the subsets are ordered based on the initial estimate $\hat {\boldsymbol {\beta }}_j$ . It is easy to see that

$$\begin{align*}\tilde{\boldsymbol{\beta}}_{j}^{(K)} = \arg\max_{\boldsymbol{\beta}_j} Q_j(\boldsymbol{\beta}_{j};\hat{\boldsymbol{\Gamma}}) \quad \text{s.t. } \boldsymbol{\beta}_{j} \in \Theta_j(\mathcal C_j^{(K)}).\end{align*}$$

The algorithm then proceeds to find a partition $\mathcal C_j^{(K-1)}$ by merging two subsets of $\mathcal C_j^{(K)}$ and then learns

$$\begin{align*}\tilde{\boldsymbol{\beta}}_{j}^{(K-1)} = \arg\max_{\boldsymbol{\beta}_j} Q_j(\boldsymbol{\beta}_{j};\hat{\boldsymbol{\Gamma}}) \quad \text{s.t. } \boldsymbol{\beta}_{j} \in \Theta_j(\mathcal C_j^{(K-1)}).\end{align*}$$

More specifically, we consider merging two adjacent subsets, with the intuition that $\beta _{jk}^{\ast}$ s that share the same value should cluster together in $\hat {\boldsymbol \beta }_j$ . We consider $K-1$ candidate partitions

$$ \begin{align*} \begin{aligned} \mathcal C_{j1}^{(K-1)} &= \left\{\{\pi_j(1), \pi_j(2)\}, \{\pi_j(3)\}, \ldots, \{\pi_{j}(K)\}\right\},\\ \mathcal C_{j2}^{(K-1)} &= \left\{\{\pi_j(1)\}, \{\pi_j(2), \pi_j(3)\}, \ldots, \{\pi_{j}(K)\}\right\},\\ &\dots\\ \mathcal C_{j,K-1}^{(K-1)} &= \left\{\{\pi_j(1)\}, \{\pi_j(2)\}, \ldots, \{\pi_{j}(K-1), \pi_{j}(K)\}\right\}, \end{aligned} \end{align*} $$

each of which contains $K-1$ subsets and is obtained by merging two adjacent subsets in $\mathcal C_j^{(K)}$ . For each candidate partition, $\mathcal C_{jb}^{(K-1)}$ , $b = 1, \ldots , K-1$ , we compute a constrained solution

$$\begin{align*}\tilde{\boldsymbol{\beta}}_{jb}^{(K-1)} = \arg\max_{\boldsymbol{\beta}_j} Q_j(\boldsymbol{\beta}_{j};\hat{\boldsymbol{\Gamma}}) \quad \text{s.t. } \boldsymbol{\beta}_j\in\Theta_j(\mathcal C_{jb}^{(K-1)}).\end{align*}$$

By construction, $\text {card}(\tilde {\boldsymbol {\beta }}_{jb}^{(K-1)}) = K-1$ for all b. The best $(K-1)$ -level solution is given by the one that yields the largest pseudo-likelihood. That is, $\tilde {\boldsymbol {\beta }}_{j}^{(K-1)} \ = \ \tilde {\boldsymbol {\beta }}_{jb_{j}^{(K-1)}}^{(K-1)},$ with associated partition $\mathcal C_j^{(K-1)} = \mathcal C_{j, b_{j}^{(K-1)}}^{(K-1)}$ , where

$$\begin{align*}b_{j}^{(K-1)} = \arg\max_{1\le b\le K-1} Q_j\!\left(\tilde{\boldsymbol{\beta}}_{jb}^{(K-1)};\hat{\boldsymbol{\Gamma}}\right).\end{align*}$$

Iteratively, suppose that the best s-level solution $\tilde {\boldsymbol {\beta }}_j^{(s)}$ has been obtained, with associated partition $\mathcal C_j^{(s)} = \{C_{j1}^{(s)},\ldots ,C_{js}^{(s)}\}$ , for $s = 2, 3, \ldots , K$ . The derivation of the best $(s-1)$ -level solution $\tilde {\boldsymbol {\beta }}_j^{(s-1)}$ is similar to that of $\tilde {\boldsymbol {\beta }}_j^{(K-1)}$ . That is, the candidate partitions for the $(s-1)$ -level solution are given by

$$ \begin{align*} \begin{aligned} \mathcal C_{j1}^{(s-1)} &= \left\{C_{j1}^{(s)}\cup C_{j2}^{(s)},C_{j3}^{(s)},\ldots,C_{js}^{(s)}\right\},\\ \mathcal C_{j2}^{(s-1)} &= \left\{C_{j1}^{(s)}, C_{j2}^{(s)}\cup C_{j3}^{(s)},\ldots,C_{js}^{(s)}\right\},\\ &\dots\\ \mathcal C_{j,s-1}^{(s-1)} &= \left\{C_{j1}^{(s)}, {C_{j2}^{(s)}}, \ldots, C_{j,s-1}^{(s)}\cup C_{js}^{(s)}\right\}. \end{aligned} \end{align*} $$

For each candidate partition, we compute a constrained solution

$$\begin{align*}\tilde{\boldsymbol{\beta}}_{jb}^{(s-1)} = \arg\max_{\boldsymbol{\beta}_j} Q_j(\boldsymbol{\beta}_{j};\hat{\boldsymbol{\Gamma}}) \quad \text{s.t. } \boldsymbol{\beta}_j\in\Theta_j(\mathcal C_{jb}^{(s-1)}), ~b = 1, \ldots, s-1.\end{align*}$$

Then, $\tilde {\boldsymbol {\beta }}_{j}^{(s-1)} \ = \ \tilde {\boldsymbol {\beta }}_{jb_{j}^{(s-1)}}^{(s-1)},$ with associated partition $\mathcal C_j^{(s-1)} = \mathcal C_{j b_{j}^{(s-1)}}^{(s-1)}$ , where

$$\begin{align*}b_{j}^{(s-1)} = \arg\max_{1\le b\le s-1} Q_j\!\left(\tilde{\boldsymbol{\beta}}_{jb}^{(s-1)};\hat{\boldsymbol{\Gamma}}\right).\end{align*}$$

By running the above calculation for $s = K, K-1, \ldots , 2$ , we obtain $\tilde {\boldsymbol {\beta }}_{j}^{(1)}$ , …, $\tilde {\boldsymbol {\beta }}_{j}^{(K)}$ as the candidate 1- to K-level models. This process is summarized in Algorithm 1. Compared with an exhaustive search over all possible sparse models, the stepwise algorithm above is significantly more computationally efficient. That is, Algorithm 1 compares only $K(K-1)/2$ models, whereas the exhaustive search compares a number of models that is superexponential in K. When $K=10$ , the stepwise method reduces the number of models to compare from 115,975 in the exhaustive search to 45.

Although the algorithm does not exhaustively search over all possible sparse models, the theoretical results in Section 2.4 show that, as the sample size N goes to infinity, the true sparse structure for item j is captured by $\tilde {\boldsymbol {\beta }}_{j}^{(m_j^{\ast})}$ with probability tending to 1. Consequently, together with the EBIC introduced below, consistent model selection can be achieved.

Algorithm 1

Stepwise search for item-level sparsity.

Algorithm 1 Long description

The algorithm is presented in a structured list format.

1. Input. Item j, number of latent classes K, unrestricted item-response estimate beta-hat sub j, and pseudo-likelihood function Q sub j.

2. Let pi sub j be a permutation of 1 to K such that beta-hat sub j, pi sub j (1) is less than beta-hat sub j, pi sub j (K).

3. Set beta-tilde sub j super (K) as beta-hat sub j.

4. Initialize the ordered partition C sub j super (K) as the set of pi sub j (1) through pi sub j (K).

5. A for loop starts for s equals K, K minus 1, down to 2.

6. Write the current ordered partition as C sub j super (s) equals the set of C sub j 1 super (s) through C sub j s super (s).

7. A nested for loop starts for b equals 1 to s minus 1.

8. Form C sub j b super (s minus 1) by merging adjacent blocks C sub j b super (s) and C sub j, b plus 1 super (s) in C sub j super (s).

9. Compute beta-bar sub j b super (s minus 1) equals arg max over beta sub j of Q sub j (beta sub j; Gamma-hat) such that beta sub j is an element of Theta sub j (C sub j b super (s minus 1)).

10. End of nested for loop.

11. Select b sub j super (s minus 1) equals arg max for 1 less than or equal to b less than or equal to s minus 1 of Q sub j (beta-bar sub j b super (s minus 1); Gamma-hat).

12. Set beta-tilde sub j super (s minus 1) as beta-bar sub j, b sub j super (s minus 1).

13. Set C sub j super (s minus 1) as C sub j, b sub j super (s minus 1) and keep the inherited block order.

14. End of main for loop.

15. Output. Candidate solutions beta-tilde sub j super (m) for m equals 1 to K.

Example 1 (continued). To further illustrate Algorithm 1, we use a simulated dataset generated from the toy example in Table 1, with sample size $N=1,000$ . We consider Item 1, for which the true item-response probability vector $\boldsymbol {\beta }_1^{\ast}$ is $(0.10,0.10,0.10,0.90)$ . The initial estimate gives $\tilde {\boldsymbol {\beta }}_1^{(4)} = \hat {\boldsymbol {\beta }}_1 = (0.141, 0.103, 0.083, 0.883)$ , where the estimates of the first three parameters are close to each other but not exactly the same. Sorting the entries of $\hat {\boldsymbol {\beta }}_1$ in increasing order gives $\pi _1(1)=3$ , $\pi _1(2)=2$ , $\pi _1(3)=1$ , and $\pi _1(4)=4$ , as well as the ordered partition $ \mathcal C_1^{(4)}=\big \{\{3\},\{2\},\{1\},\{4\}\big \}. $

To find $\tilde {\boldsymbol {\beta }}_1^{(3)},$ there are three candidate partitions obtained by merging adjacent blocks:

$$ \begin{align*} \mathcal C_{11}^{(3)}=\big\{\{2,3\},\{1\},\{4\}\big\},\quad \mathcal C_{12}^{(3)}=\big\{\{3\},\{1,2\},\{4\}\big\},\quad \mathcal C_{13}^{(3)}=\big\{\{3\},\{2\},\{1,4\}\big\}. \end{align*} $$

These partitions yield the constrained estimates $ \tilde {\boldsymbol {\beta }}_{11}^{(3)} =(0.141,\,0.094,\,0.094,\,0.883), $ $ \tilde {\boldsymbol {\beta }}_{12}^{(3)} =(0.124, 0.124,\,0.083,\,0.883), $ and $ \tilde {\boldsymbol {\beta }}_{13}^{(3)} =(0.490,\,0.103,\,0.083,\,0.490), $ with pseudo-likelihood values $-353.713$ , $-354.332$ , and $-532.043$ , respectively. Since $\tilde {\boldsymbol {\beta }}_{11}^{(3)}$ gives the largest pseudo-likelihood, we obtain $\tilde {\boldsymbol {\beta }}_1^{(3)} = \tilde {\boldsymbol {\beta }}_{11}^{(3)}=(0.141,\,0.094,\,0.094,\,0.883)$ , with associated partition $ \mathcal C_1^{(3)}=\mathcal C_{11}^{(3)}=\big \{\{2,3\},\{1\},\{4\}\big \}. $

By further executing Algorithm 1, we obtain $ \tilde {\boldsymbol {\beta }}_1^{(2)} =(0.114,\,0.114,\,0.114,\,0.883) $ and $ \tilde {\boldsymbol {\beta }}_1^{(1)} =(0.323,\,0.323,\,0.323,\,0.323). $ We note that when the number of levels is correctly set at 2, the estimate $\tilde {\boldsymbol {\beta }}_1^{(2)}$ recovers the same sparsity pattern as the true parameter in Table 1, which is $\mathcal C_{1}^{(2)} = \{\{1,2,3\},\{4\}\}$ .

2.3.3 Extended BIC

To determine the optimal sparsity level for each item, we construct an EBIC as

(4) $$ \begin{align} \text{EBIC}_{j}(m_j) = -2Q_j(\boldsymbol{\tilde{\beta}}_{j}^{(m_j)}; \hat{\boldsymbol{\Gamma}}) + m_j\log(N) + 2m_j \log(\rho). \end{align} $$

This EBIC is constructed following the idea proposed in the EBIC for high-dimensional regression models (Chen & Chen, Reference Chen and Chen2008, Reference Chen and Chen2012), where $\rho $ is a constant that characterizes the user’s prior belief on model sparsity. More specifically, the first two terms in (4) mimic the standard BIC (Schwarz, Reference Schwarz1978), but replace a real likelihood function with the pseudo-likelihood. These terms arise from the Laplace approximation to the marginal likelihood in a Bayesian formulation in which a prior distribution is imposed on all possible models and their model parameters. The last term arises from the prior distribution over models with different sparsity levels, with $\rho \geq 1$ being the ratio of the prior probability of a randomly chosen m-level model to that of a randomly chosen $(m+1)$ -level model. We set $\rho \geq 1$ to ensure that sparse models are preferred. When $\rho $ is set to be 1, the third term becomes zero, so that the EBIC becomes a pseudo-BIC. The role of $\rho $ will be clarified in Section 3.3.

Based on the EBIC, we select the number of distinct levels for item j as

$$ \begin{align*} \tilde m_j = \arg\min_{m_j} \text{EBIC}_{j}(m_j). \end{align*} $$

Given $\tilde m_j$ , we infer the sparse item-response probability pattern based on the refined estimate $\boldsymbol {\tilde {\beta }}_{j}^{(\tilde {m}_{j})}$ .

Next, we comment on the EBIC. As we set $\rho $ to be a constant, the third term in (4) is dominated by the second term $m_j\log (N)$ as the sample size N goes to infinity. Consequently, it will not affect the asymptotic property of the EBIC. However, this term is useful when the number of classes is relatively large compared with the sample size. In that case, the model space is large, and the effective sample size is relatively small. As a result, the standard BIC tends to over-select, as in high-dimensional regression settings (Chen & Chen, Reference Chen and Chen2008, Reference Chen and Chen2012). Our simulation results in Section 3 show that this additional penalty term substantially improves model selection accuracy. In particular, $\rho =20$ provides a practical finite-sample compromise: it reduces over-selection relative to weaker penalties while avoiding the pronounced under-selection induced by larger penalties. Moreover, as shown in Section 2.4, the EBIC consistently learns the true number of levels $m_j^{\ast}$ for each item j, that is, as the sample size N goes to infinity, with probability tending to 1, $\tilde m_j = m_j^{\ast}$ . Combined with the model selection consistency of $\tilde {\boldsymbol {\beta }}_{j}^{(m_j^{\ast})}$ from the stepwise search algorithm, the refined estimator $\boldsymbol {\tilde {\beta }}_{j}^{(\tilde {m}_{j})}$ consistently recovers the sparse pattern of $\boldsymbol {\beta }^{\ast}_j$ .

Example 1 (continued). We illustrate the use of EBIC by continuing the analysis of Example 1. From the stepwise search, we obtained the four candidate models with pseudo-likelihood values

$$\begin{align*}Q_1\!(\tilde{\boldsymbol{\beta}}_1^{(4)}\!;\hat{\boldsymbol{\Gamma}})=-353.47,\ Q_1\!(\tilde{\boldsymbol{\beta}}_1^{(3)}\!;\hat{\boldsymbol{\Gamma}})=-353.71,\ Q_1\!(\tilde{\boldsymbol{\beta}}_1^{(2)}\!;\hat{\boldsymbol{\Gamma}})=-355.60,\ Q_1\!(\tilde{\boldsymbol{\beta}}_1^{(1)}\!;\hat{\boldsymbol{\Gamma}})=-629.11. \end{align*}$$

With $N=1,000$ and $\rho =20$ , the EBIC values for the four candidate models are

$$\begin{align*}\mathrm{EBIC}_1(4)=758.54,\quad \mathrm{EBIC}_1(3)=746.12,\quad \mathrm{EBIC}_1(2)=736.99,\quad \mathrm{EBIC}_1(1)=1,271.12. \end{align*}$$

Hence,

$$\begin{align*}\tilde m_1=\arg\min_{m\in\{1,2,3,4\}} \mathrm{EBIC}_1(m)=2, \end{align*}$$

which correctly recovers the true number of distinct levels for Item 1. And the corresponding estimate $ \tilde {\boldsymbol {\beta }}_1^{(2)}=(0.114,\,0.114,\,0.114,\,0.883) $ recovers the same sparsity pattern as the true parameter vector in Table 1.

2.3.4 Final estimator

We obtain the final estimate of the latent class model via constrained marginal maximum likelihood estimation, using the item-level clustering pattern determined by EBIC. That is,

(5) $$ \begin{align} (\bar{\boldsymbol\nu}, \bar{\mathbf B}) = \arg\max_{(\boldsymbol{\nu}, \mathbf B)} ~l(\boldsymbol{\nu}, \mathbf B), s.t. \ \beta_{jg} = \beta_{jh} \text{ if } \tilde \beta_{jg}^{(\tilde{m}_{j})} = \tilde \beta_{jh}^{(\tilde{m}_{j})}, j=1, \ldots, J, \text{ and } h,g = 1, \ldots, K. \end{align} $$

The standard error for each $\bar \beta _{jk}$ and $\bar \nu _k$ is also obtained based on the Fisher information matrix of this marginal likelihood. We summarize the full proposed procedure in Algorithm 2.

Algorithm 2

Proposed refinement procedure

Algorithm 2 Long description

The pseudocode is structured in 11 numbered steps.

Step 1: Input includes binary response matrix Y, number of latent classes K, initial L C A estimate (nu-hat, B-hat), posterior probability matrix Gamma-hat, and E B I C tuning parameter rho.

Step 2: A ‘for’ loop begins for j = 1 to J.

Step 3: Construct item-specific pseudo-likelihood Q sub j (beta sub j; Gamma-hat).

Step 4: Apply an algorithm to obtain candidate estimates beta-tilde sub j super (m) for m = 1 to K.

Step 5: A nested ‘for’ loop begins for m = 1 to K.

Step 6: Compute E B I C sub j (m) = -2 Q sub j (beta-tilde sub j super (m); Gamma-hat) + m log N + 2 m log rho.

Step 7: End nested loop.

Step 8: Select m-bar sub j = arg min over 1 <= m <= K of E B I C sub j (m), and let C-tilde sub j be the partition induced by beta-tilde sub j super (m-bar sub j).

Step 9: End outer loop.

Step 10: Estimate (nu-tilde, B-bar) by maximizing the marginal likelihood subject to beta sub j in Theta sub j (C-tilde sub j) for all j = 1 to J.

Step 11: Output includes selected item-level sparsity pattern {m-bar sub j, C-tilde sub j : j = 1 to J} and final sparse L C A estimate (nu-tilde, B-bar).

2.4 Theoretical results

In what follows, we show that the sparse item-response probability pattern implied by $\boldsymbol {\tilde {\beta }}_{j}^{(\tilde {m}_{j})}$ for each item j is statistically consistent, for all $j = 1, \ldots , J$ . This result further implies that the statistical inference based on the final estimate $(\bar {\mathbf B}, \bar {\boldsymbol \nu })$ is asymptotically valid. We make the following assumptions.

Assumption 1. The MLE in (2.1) satisfies

$$ \begin{align*} \Vert \hat{\boldsymbol{\nu}} - {\boldsymbol{\nu}}^{\ast}\Vert = O_p(1/\sqrt{N}) \text{ and } \Vert \hat{\mathbf B} - {\mathbf B}^{\ast} \Vert_F = O_p(1/\sqrt{N}), \end{align*} $$

where ${\boldsymbol {\nu }}^{\ast}$ and ${\mathbf B}^{\ast}$ are the corresponding true model parameters.

Since J and K are fixed, Assumption 2 implies that there exists a constant $\varepsilon>0$ such that $\varepsilon \le \beta _{jk}^{\ast}\le 1-\varepsilon $ and $\nu _k^{\ast}\ge \varepsilon $ for all $j=1,\dots ,J$ and $k=1,\dots ,K$ .

The proof of the theorem is given in Section A of the Supplementary Material. Theorem 1 gives two results about the statistical consistency of the proposed method. The first result, $P(\tilde m_j=m_j^{\ast})$ converges to 1, states that the EBIC selection step consistently identifies the true number of distinct response-probability levels for each item. The second result, which concerns $P(\tilde {\beta }_{jg} = \tilde {\beta }_{jh})$ and $P(\tilde {\beta }_{jg} \neq \tilde {\beta }_{jh})$ , further establishes the consistency in recovering the sparse pattern for each item: latent classes with the same true item-response probability are asymptotically grouped together, whereas latent classes with different true item-response probabilities are asymptotically kept in separate groups. Thus, the refinement procedure consistently recovers both the item-specific model complexity and the equality pattern among the K latent classes.

3.1 Simulation settings

We consider two simulation settings. Setting I has $K=4$ latent classes and $J=32$ binary items, and Setting II has $K=8$ latent classes and $J=64$ binary items. For each setting, we generate one true class proportion vector $\boldsymbol {\nu }$ and one true item-response probability matrix $\mathbf B$ and then keep these true parameters fixed across all sample sizes and replications. The true class proportions are $\boldsymbol {\nu }=(0.241,0.259,0.202,0.298)$ in Setting I and $\boldsymbol {\nu }=(0.106,0.137,0.130,0.101,0.124,0.117,0.150,0.136)$ in Setting II. For each setting, we consider sample sizes $N\in \{500,750,1,000,1,500,2,000\}$ and run 100 independent replications for each sample size. The two settings also differ in the level configurations of the true item-response probabilities. In Setting I, all 32 items have two distinct response-probability levels, where a level denotes a distinct value among the K class-specific response probabilities for a given item rather than an observed response category. In Setting II, the 64 items consist of 48 two-level items and 16 three-level items, yielding a richer and more complex item-level structure. The full true class proportions and item-response probabilities are reported in Section B of the Supplementary Material.

In the refinement stage, we consider multiple candidate values of the EBIC tuning parameter $\rho $ , namely, 1, 5, 10, 20, 40, 80, 160, and 320. Recall that the additional EBIC penalty is $2m_j\log \rho $ , so this term is zero when $\rho =1$ and increases with $\rho $ , placing a stronger penalty on denser models.

In the main simulation analysis, the unrestricted first-stage model for $(\hat {\mathbf B},\hat {\boldsymbol {\nu }})$ is fitted using the true value of K, so that the reported results focus on the proposed item-level refinement.Footnote ² Furthermore, we also conducted additional simulations to compare the proposed procedure with the regularized method of Chen et al. (Reference Chen, Li, Liu and Ying2017); the results are reported in Section E of the Supplementary Material.

3.2 Evaluation criteria

Multiple criteria are used to evaluate the performance of the proposed method. First, we evaluate the selection of the item-specific numbers of probability levels. For each simulated dataset, we record the numbers of items that are under-selected, correctly selected, and over-selected relative to the true values $\{m_j^{\ast}\}_{j=1}^J$ . We also record the total number of incorrectly selected items per replication, where an item is counted as incorrect whenever $\tilde m_j \neq m_j^{\ast}$ .

Second, we assess recovery of the sparse item-response probability pattern itself. This is done by comparing the partition of the latent classes induced by $\boldsymbol {\beta }_j^{\ast}$ , denoted $\mathcal C_j^{\ast}$ , with that induced by $\tilde {\boldsymbol {\beta }}_j$ , denoted $\tilde {\mathcal {C}}_j$ , using the adjusted Rand index (ARI; Hubert & Arabie, Reference Hubert and Arabie1985). For each item j, we compute

$$\begin{align*}\mathrm{ARI}(\mathcal C_j^{\ast},\tilde{\mathcal C}_j) = \frac{\sum_{g,h}\binom{n_{gh}}{2}-\frac{\sum_g\binom{n_{g\cdot}}{2}\sum_h\binom{n_{\cdot h}}{2}}{\binom{K}{2}}} {\frac12\left[\sum_g\binom{n_{g\cdot}}{2}+\sum_h\binom{n_{\cdot h}}{2}\right]-\frac{\sum_g\binom{n_{g\cdot}}{2}\sum_h\binom{n_{\cdot h}}{2}}{\binom{K}{2}}}. \end{align*}$$

Here, $n_{gh}$ denotes the number of latent classes that fall into group g under $\mathcal C_j^\star $ and into group h under $\tilde {\mathcal C}_j$ , and $n_{g\cdot }=\sum _h n_{gh}$ and $n_{\cdot h}=\sum _g n_{gh}$ are the corresponding marginal counts. For each simulated dataset and each value of $\rho $ , we compute the average item-level ARI, $(1/J)\sum _{j=1}^J \mathrm {ARI}(\mathcal C_j^\star ,\tilde {\mathcal C}_j)$ .

Finally, we assess parameter estimation accuracy for both the item-response probabilities and the class proportions. Specifically, for each simulated dataset and each value of $\rho $ , we compute the mean squared error of the final refined estimator:

$$\begin{align*}\mathrm{MSE}(\bar{\mathbf B}) = \frac{1}{JK}\min_{\pi}\sum_{j=1}^J\sum_{k=1}^K\left(\bar{\beta}_{j,\pi(k)}-\beta_{jk}^{\ast}\right)^2,\end{align*}$$

where the minimum is taken over all permutations $\pi $ on $\{1, \ldots , K\}$ to account for label indeterminacy. Let $\hat {\pi }$ denote the permutation that attains this minimum. We then compute the corresponding error for the class proportions using the same alignment:

$$\begin{align*}\mathrm{MSE}(\bar{\boldsymbol{\nu}}) = \frac{1}{K}\sum_{k=1}^K\left(\bar{\nu}_{\hat{\pi}(k)}-\nu_k^{\ast}\right)^2. \end{align*}$$

To evaluate the effect of the refinement step, we compute the same quantities for the unrestricted first-stage estimator $(\hat {\mathbf B},\hat {\boldsymbol {\nu }})$ and compare them with those of the final refined estimator $(\bar {\mathbf B},\bar {\boldsymbol {\nu }})$ .

3.3 Simulation results

We first discuss the choice of $\rho $ , which controls a clear trade-off between over-selection and under-selection in the refinement step. As $\rho $ increases, denser models are penalized more heavily, so over-selection tends to decrease. At the same time, a larger $\rho $ makes the procedure more prone to under-selection. Table 2 shows the average numbers of under-selected and over-selected items under different values of $\rho $ . In Setting I, this trade-off is mild. Because this setting is relatively simple, the overall performance is strong across all sample sizes. Under-selection is essentially absent across all sample sizes and all candidate values of $\rho $ , while over-selection decreases steadily as $\rho $ increases and also becomes smaller as N increases. In contrast, the trade-off is much more pronounced in Setting II. For each fixed $\rho $ , both types of selection error generally decrease as N increases. For each fixed sample size, increasing $\rho $ sharply reduces over-selection, but this improvement is accompanied by more under-selection, especially when N is small. For example, at $N=500$ , over-selection drops from 8.94 items at $\rho =1$ to 0.40 at $\rho =20$ , while under-selection increases from 0.27 to 2.33.

Table 2

Item-level selection errors under different EBIC tuning values $\rho $

Table 2 Long description

The table is divided into two main sections: Setting I and Setting II. Each section contains columns for rho and two main categories: Under-selected items and Over-selected items. Both categories are further subdivided by sample sizes: 500, 750, 1,000, 1,500, and 2,000.

In Setting I:

* For all values of rho (1, 5, 10, 20, 40, 80, 160, 320), the number of under-selected items is consistently 0.00 across all sample sizes.

* Over-selected items decrease as sample size and rho increase. For rho equals 1, values drop from 1.45 at n equals 500 to 0.69 at n equals 2,000. For rho equals 320, values drop from 0.03 at n equals 500 to 0.00 at n equals 2,000.

In Setting II:

* Under-selected items increase as rho increases but decrease as sample size increases. For rho equals 1, values drop from 0.27 at n equals 500 to 0.00 at n equals 1,000 and above. For rho equals 320, values drop from 6.24 at n equals 500 to 0.00 at n equals 1,500 and above.

* Over-selected items decrease as both rho and sample size increase. For rho equals 1, values drop from 8.94 at n equals 500 to 4.25 at n equals 2,000. For rho equals 320, values drop from 0.03 at n equals 500 to 0.01 at n equals 2,000.

The row for rho equals 20 is highlighted in bold in both settings as a reference point for further simulations.

Note: For each setting and sample size, the table reports the average numbers of under-selected and over-selected items per replication. The

row for $\rho =20$ is highlighted in bold because this value is used in the remaining simulation summaries.

For the remaining analyses, we use $\rho =20$ . This choice is motivated by the fact that the real-data application in Section 4 has sample and item sizes similar to those of the $N=750$ case under Setting II. For this case, $\rho =20$ provides a reasonably balanced compromise between over-selection and under-selection. In particular, at $N=750$ in Setting II, increasing $\rho $ from 1 to 20 reduces the average number of over-selected items from 7.13 to 0.41, while the average number of under-selected items is still moderate at 0.32; larger values of $\rho $ reduce over-selection only slightly further but induce more under-selection. This asymmetry is important because under-selection can introduce substantial bias in the refined estimates of $\bar {\mathbf B}$ , whereas moderate over-selection is less harmful for estimation accuracy. We therefore use $\rho =20$ in the remaining simulation summaries and in the real-data analysis. Results for the other candidate values of $\rho $ are reported in Section D of the Supplementary Material.

Figure 1 summarizes, for each replication, the proportion of items for which the number of levels, $m_j$ , is correctly selected. In Setting I, this proportion is concentrated very close to one across all sample sizes, indicating that EBIC almost always recovers the true item-level complexity. Even in the smallest samples, the worst replications incorrectly select at most two items out of 32, and the distribution becomes increasingly concentrated at perfect selection as N increases. Setting II is more challenging, but the boxplots still shift steadily upward with N. At $N=500$ , the median proportion of correctly selected items is already above 0.95, meaning that the typical replication incorrectly selects no more than three items out of 64. Once N reaches 750 or above, the median is essentially at perfect selection, and by $N=2,000$ , most replications achieve perfect selection. Overall, these results show that the EBIC step accurately identifies the true number of item-level probability levels, and that its performance improves with sample size, in line with theoretical results.

Figure 1

Proportion of items with correctly selected $m_j$ in each replication at $\rho =20$ . Columns correspond to Settings I and II.

Figure 1 Long description

The y-axis represents Correctly selected item proportion ranging from 88% to 100% in 4% increments. The x-axis represents Sample size N with discrete values at 500, 750, 1000, 1500, and 2000.

* Setting I (Left Panel): Shows near-perfect performance. For all sample sizes, the median is at 100%. There are a few outlier dots at approximately 97% and 94% for sample sizes 500, 750, and 1000, and a single outlier at 97% for sample sizes 1500 and 2000.

* Setting II (Right Panel): Shows more variability that decreases as sample size increases.

- At N = 500, the box spans from 94% to 97% with a median around 95.5% and whiskers extending from 89% to 100%.

- At N = 750, 1000, and 1500, the median and upper quartile are at 100%, with the lower quartile at approximately 98.5% and whiskers extending down to 97%. Outliers are present between 92% and 96%.

- At N = 2000, the box collapses to a single line at 100% with three outliers between 95% and 98%.

Figure 2 provides a complementary summary based on the mean item-level ARI, which measures how well the estimated class partition matches the true sparse partition once $m_j$ has been selected. In Setting I, the ARI remains high, and the median is essentially at perfect partition recovery across all sample sizes; as N increases, the remaining outliers become less frequent. In Setting II, recovery is weaker at $N=500$ , with a mean item-level ARI of around 0.95, but it improves quickly with sample size and exceeds 0.98 once N reaches 750. By $N=2,000$ , most replications recover the true sparse pattern almost perfectly, apart from a small number of outliers. Taken together, these results show that the selected item-level sparsity pattern converges rapidly to the truth, with most remaining recovery errors concentrated in the smallest-sample Setting II scenario.

Figure 2

Mean item-level ARI in each replication at $\rho =20$ . Columns correspond to Settings I and II.

Figure 2 Long description

A two-panel box plot visualization. The Y-axis is labeled Mean item-level A R I and ranges from 0.925 to 1.000. The X-axis is labeled Sample size N with discrete values at 500, 750, 1000, 1500, and 2000.

* Setting I (Left Panel): Shows near-perfect performance across all sample sizes. The median lines are consistently at 1.000. A few outlier dots appear below the main line, primarily between 0.960 and 0.990.

* Setting II (Right Panel): Shows an asymptotic increase in performance as sample size increases.

- At N = 500, the median is approximately 0.953 with a wide interquartile range extending from 0.942 to 0.962.

- At N = 750, the median rises to approximately 0.985.

- At N = 1000, the median is approximately 0.993.

- At N = 1500, the median is nearly 1.000 with a very narrow box.

- At N = 2000, the median is 1.000 with minimal variance and a few outliers just below the line.

Figure 3 compares the MSE of the unrestricted first-stage estimate with that of the final refined estimate, separately for $\mathbf B$ and $\boldsymbol {\nu }$ . For both parameters, the MSE decreases with increasing sample size, consistent with theoretical results. For $\mathbf B$ , the refined estimator uniformly improves on the unrestricted estimator across both simulation settings and all sample sizes. In Setting I, the mean MSE decreases from $1.73\times 10^{-3}$ to $8.32\times 10^{-4}$ at $N=500$ , and from $4.27\times 10^{-4}$ to $2.04\times 10^{-4}$ at $N=2,000$ . In Setting II, the corresponding decrease is from $3.27\times 10^{-3}$ to $1.69\times 10^{-3}$ at $N=500$ , and from $7.90\times 10^{-4}$ to $2.20\times 10^{-4}$ at $N=2,000$ . These gains indicate that, once the item-level sparse structure is recovered with sufficient accuracy, imposing the selected equality constraints substantially improves estimation of the item-response probability matrix.

Figure 3

Mean squared error of the item-response probability matrix $\mathbf B$ and class proportions $\boldsymbol {\nu }$ across replications at $\rho =20$ . The top panel reports the MSE for $\mathbf B$ , and the bottom panel reports the MSE for $\boldsymbol {\nu }$ . Within each setting and sample size, the two boxplots compare the unrestricted first-stage estimate with the final refined estimate.

Figure 3 Long description

The figure consists of two main horizontal rows.

Top Row: M S E of B.

This row contains two panels labeled Setting I and Setting II. The y-axis represents M S E of B on a logarithmic scale from 0.0001 to 0.0030. The x-axis shows Sample size N with values 500, 750, 1000, 1500, and 2000.

- In Setting I, M S E decreases as N increases. For every N, the Refined B (dark blue) shows significantly lower M S E and less variance than the Unrestricted B hat (light blue).

- In Setting II, the trend is similar but the gap between Unrestricted and Refined estimates is wider, with Refined B consistently lower.

Bottom Row: M S E of nu.

This row also contains Setting I and Setting II panels. The y-axis represents M S E of nu on a logarithmic scale from 0.00001 to 0.00100.

- In Setting I, both Unrestricted nu hat and Refined nu bar show a gradual decrease in M S E as N increases. The two boxplots are nearly identical in height and median for all sample sizes, though Refined nu bar has slightly fewer extreme outliers.

- In Setting II, the M S E values are lower overall than in Setting I. The Unrestricted and Refined estimates perform almost identically across all sample sizes, with a steady downward trend as N increases from 500 to 2000.

For the class proportions $\boldsymbol {\nu }$ , the unrestricted and refined estimates are much closer, as expected, because the refinement step is designed primarily to improve estimation of $\mathbf B$ rather than the class proportions. The MSE of $\boldsymbol {\nu }$ nevertheless decreases with sample size in both settings. In Setting I, the mean refined MSE decreases from $4.53\times 10^{-4}$ at $N=500$ to $1.18\times 10^{-4}$ at $N=2,000$ ; in Setting II, it decreases from $2.14\times 10^{-4}$ to $6.23\times 10^{-5}$ over the same range. Thus, the main estimation gain from refinement occurs for $\mathbf B$ , while estimation of $\boldsymbol {\nu }$ is largely preserved.

4.1 Data

We applied the proposed refinement procedure to the Social Role Performance (SRPPER) item bank from the PROMIS Wave 1 calibration study (Cella et al., Reference Cella, Riley, Stone, Rothrock, Reeve, Yount, Amtmann, Bode, Buysse, Choi, Cook, DeVellis, DeWalt, Fries, Gershon, Hahn, Lai, Pilkonis, Revicki and Hays2010). The analysis included respondents with complete responses to all SRPPER items, yielding a response matrix with $N=737$ respondents and $J=56$ items, listed in Table G1 in the Supplementary Material.

The original SRPPER items in the PROMIS dataset were measured on a five-point ordinal scale with response categories Never, Rarely, Sometimes, Often, and Always. Because the items included both positively and negatively worded statements, the numeric scores from 1 to 5 were assigned to align all items in the same direction, with higher values indicating more severe limitations in social role performance. In the analysis, we dichotomized each item as $Y_{ij}=1$ for responses coded 4 or 5 and $Y_{ij}=0$ otherwise. Therefore, the item-response probabilities reported below can be interpreted as the probabilities of endorsing a more severe response under this binary coding.

4.2 Analysis workflow

The real-data analysis followed the strategy developed in Section 2.3. We first fitted unrestricted latent class models with candidate numbers of classes ranging from 1 to 8 and selected $K=8$ latent classes based on BIC. We then applied the proposed item-level refinement procedure to this selected unrestricted solution. The refinement step used $\rho =20$ , as suggested by the simulation results (see Section 3.3 and Setting II in Table 2), to select item-specific equality groups among the eight latent classes. Finally, the sparse latent class model was re-estimated under these selected equality constraints, producing the refined item-response probability matrix $\bar {\mathbf B}$ and the estimated class proportions.