2.1 Overview of existing HMM research in education and psychology

Consider a general HMM model for an observation at time points $t = 1, \ldots , T$ . Let $\boldsymbol {X} = (\boldsymbol {X}_1, \dots , \boldsymbol {X}_T)^\top \notag $ denote the response data of an individual over time. Let $[q]$ for $q\in \mathbb N$ denote the set $\{1, 2, \dots , q\}$ . Here, $\boldsymbol {X}_t = (X_{1t}, \dots , X_{Jt})$ is the response pattern of an individual for J items at time t and $\boldsymbol {X}_t \in ([M_1] \times [M_2] \times \dots \times [M_J])$ . In other words, $M_j$ is the number of response options for an item j. Let $\boldsymbol {Z} = (Z_1, \dots , Z_T)^\top \notag $ be the latent states of an individual over time and $Z_t \in [r]$ for $r \in \mathbb {N}$ . The response patterns $\boldsymbol {X}_1, \ldots , \boldsymbol {X}_T$ are assumed to be conditionally independent given $Z_1, \ldots , Z_T$ . In other words, each $\boldsymbol {X}_t$ is directly impacted by the corresponding $Z_t$ and the conditional distribution is

$$ \begin{align*} P(\boldsymbol{X}_t \mid Z_1, \dots, Z_{t}) = P(\boldsymbol{X}_t \mid Z_t). \end{align*} $$

Furthermore, within a first-order HMM, the distribution of $Z_t$ depends only on the previous state of the immediate past $Z_{t - 1}$ . That is,

$$ \begin{align*} P(Z_t \mid Z_1, \dots, Z_{t - 1}) = P(Z_t \mid Z_{t - 1}). \end{align*} $$

HMMs consist of two components. The first element is a transition matrix that dictates the likelihood of changing between latent states over two adjacent time points. Let $\mathbf {A}$ be an $r \times r$ transition matrix with each element $(z_{t - 1}, z_{t})$ represents

$$\begin{align*}a_{z_{t - 1}, z_{t}} = P(Z_{t} = z_{t} \mid Z_{t - 1} = z_{t - 1}),\end{align*}$$

which is the probability of transition to state $Z_t = z_{t}$ given the membership during the previous time is $Z_{t - 1} = z_{t - 1}$ . Observe that each element in $\mathbf {A}$ is non-negative and each row sums to one.

Additionally, there is an emission matrix that describes the distribution of observed responses given contemporaneous latent states. Let $\mathbf {B}$ be the $q \times r$ emission matrix with each element $(\boldsymbol {x}_t, z_t)$ denoted by

$$\begin{align*}b_{\boldsymbol{x}_t, z_t} = P(\boldsymbol{X}_t = \boldsymbol{x}_t \mid Z_t = z_t).\end{align*}$$

That is, each element in $\mathbf {B}$ represents the conditional probability of a response pattern given contemporaneous hidden state $Z_t = z_t$ . The elements of $\mathbf {B}$ are also non-negative and the columns sum to one.

We define $\boldsymbol {\pi }_t = (\pi _{t, 1}, \ldots , \pi _{t, r}),$ where $\pi _{t, z_t}$ is the probability of being in state $Z_t = z_t$ at time t. Then, $\boldsymbol {\pi }_1$ is the initial distribution of the latent states. Observe that $\sum _{z_t = 1}^r \pi _{t, z_t} = 1$ for any $t \in [T]$ . Using the initial distribution, we can iterate forward and find the distribution for time $t = 2$ , $\boldsymbol {\pi }_2$ , to have elements

$$ \begin{align*} P(Z_2 = z_2) = \sum_{z_1 = 1}^r P(Z_1 = z_1) P(Z_2 = z_2 \mid Z_1 = z_1). \end{align*} $$

We can write $\boldsymbol {\pi }_2 = \boldsymbol {\pi }_1 \mathbf {A}$ . Consequently, iterating forward to time n implies that the distribution of latent states is $\boldsymbol {\pi }_n = \boldsymbol {\pi }_1 \mathbf {A}^{n - 1}$ . The stationary (or long-run) distribution, which characterizes the marginal probability of residing in state $z_t$ , is defined as $\boldsymbol {\pi } = \lim _{n \rightarrow \infty } \boldsymbol {\pi }_1 \mathbf {A}^n = (\pi _{1}, \ldots , \pi _{r})$ .

For HMMs and RHMMs, a model is identifiable when different values of the parameter set ( $\mathbf {A}, \mathbf {B}, \boldsymbol {\pi }_1$ ) correspond to distinct marginal distributions of the observed variables. Model identifiability is essential for ensuring that the underlying parameters can be consistently estimated, thereby enabling the model to accurately assess the effects of interventions. Allman et al. (Reference Allman, Matias and Rhodes2009) provide a generic identifiability condition for HMMs using results in Kruskal (Reference Kruskal1977). More recently, Bonhomme et al. (Reference Bonhomme, Jochmans and Robin2016) establish sufficient conditions for identifying HMMs, which are more feasible to verify in practice.

Here, $\mathbf {D}_{\boldsymbol {\pi }} = \text {diag}(\pi _1, \dots , \pi _r)$ . The proof of Theorem 1 involves rewriting HMMs into mixture models represented as three-way arrays, where the conditional distribution of the responses can be expressed in terms of $\mathbf {B}$ , $\mathbf {A}$ , and $\boldsymbol {\pi }$ . Then, Bonhomme et al. (Reference Bonhomme, Jochmans and Robin2016) derived conditions for finite HMMs to be identifiable using results from the simultaneous diagonalization problem (De Leeuw & Pruzansky, Reference De Leeuw and Pruzansky1978). It is worth noting that Theorem 1 applies to general HMMs and requires $\mathbf {B}$ and $\mathbf {A}$ to be full column rank. Other restrictions may be imposed to identify the parameters of an HMM. Liu et al. (Reference Liu, Culpepper and Chen2023) and Chen et al. (Reference Chen, Culpepper and Liang2020) extend the results to prove identifiability conditions for HMMs and RHMMs by relaxing some assumptions in Bonhomme et al. (Reference Bonhomme, Jochmans and Robin2016). Note that Liu et al. (Reference Liu, Culpepper and Chen2023) and Chen et al. (Reference Chen, Culpepper and Liang2020) focus on HMMs and RHMMs with binary latent attributes. Particularly, let K be the number of binary latent attributes, then the number of latent states is $r = 2^K$ corresponding to all possible attribute profiles.

Note rank $_K$ in part (b) of Theorem 2 denotes the Kruskal rank. The Kruskal rank of a matrix with r columns equals r if and only if the matrix has full column rank (see Definition A.1 in Appendix A for more details). Moreover, observe that certain assumptions are shared between Theorems 1 and 2. Namely, $\mathbf {B}$ as well as $\mathbf {A}$ need to be constant over time, and that every element of the stationary distribution is positive. We now examine what these assumptions entail in real-world settings as well as within the model.

2.2 Challenges with current HMMs

Bonhomme et al. (Reference Bonhomme, Jochmans and Robin2016), Chen et al. (Reference Chen, Culpepper and Liang2020), and Liu et al. (Reference Liu, Culpepper and Chen2023) assume $\mathbf {B}$ and $\mathbf {A}$ to be time invariant along with every element in $\boldsymbol {\pi }$ to be positive. It is essential to understand the implications of these assumptions. Having constant $\mathbf {B}$ implies that the conditional distribution $P(\boldsymbol {X}_t \mid Z_t)$ does not change with time. In other words, parallel measurements are observed across different time points. This assumption could be reasonable when a common set of items is used over time, such as diagnosing mental health conditions with a standardized questionnaire. However, the constant $\mathbf {B}$ assumption is not always feasible, such as when there may be item parameter drift or when different items are administered over time.

The assumption of constant $\mathbf {A}$ suggests that $P(Z_t \mid Z_{t - 1})$ remains the same for all t. That is, the process being modeled is assumed to be stable during the observation period. This assumption is possible when data are collected over a short time span and conditions are expected to remain consistent. However, the constant $\mathbf {A}$ assumption may be violated in studies involving extended periods or heterogeneous populations, where individual response to interventions may vary. For example, in educational settings, learning interventions may only be effective for certain subgroups, resulting in different transition dynamics between the groups.

The last assumption requires $\boldsymbol {\pi }> \boldsymbol {0}_r$ with $\boldsymbol {0}_r$ represents an r-vector of zeros and “ $>$ ” is an element-wise inequality. This pertains to the irreducibility of $\mathbf {A}$ , which is violated if $\mathbf {A}$ contains absorbing states. A state $z_{t - 1}$ is an absorbing state if it has a $100\%$ of transitioning to itself and $0\%$ to any other states. Intuitively, an absorbing state is a state such that once entered, cannot be left. The irreducibility assumption can be violated in educational research aimed at learning advancement. In such cases, interventions are designed to enhance skill development toward a mastery state, which may function as an absorbing state if skill forgetting is considered unlikely.

For a single-group pretest/posttest design, the assumptions of constant $\mathbf {B}$ and $\mathbf {A}$ are often reasonable as it is common to use similar assessments before and after an intervention. Consider a study in mathematics education examining the effects of the Enhanced Anchored Instruction program (Bottge et al., Reference Bottge, Heinrichs, Chan, Mehta and Watson2003) on mathematics achievement of middle school students (Li et al., Reference Li, Cohen, Bottge and Templin2016). The study investigates the impact of the EAI intervention by observing student’ performance on a test called fraction of the cost (FOC). The same FOC test is administered before and after the intervention. Hence, we can see that the constant $\mathbf {B}$ and $\mathbf {A}$ assumptions are satisfied here. Furthermore, the study consists of only two time points, making the likelihood of having a complete mastery state to be unlikely. This means that students can be expected to obtain or retain knowledge as well as forget some knowledge over the short periods. Therefore, the irreducibility assumption is also feasible, allowing Theorem 2 to be applied to this design. However, these assumptions may not hold for situations where item parameters change over time and there are multiple groups with heterogeneous transition parameters. The next sections discuss the specific setups and identifiability proofs of two designs that we recommend using in cases where item parameters differ over time or when there are known groups that might differ in the probability of achieving mastery.

2.3 Counterbalancing for longitudinal diagnostic design

There are several advantages for one to consider implementing a counterbalancing design. This design balances out the item positions and controls for potential order effects. Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018) proposed using a counterbalanced booklet design to address confounding factors that could obscure learning effects. Additionally, Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018) noted that counterbalancing can help avoid empirical identifiability issues associated with using RLCMs, specifically the DINA model. Particularly, they provided the example of estimated guessing parameters being inaccurate for items in later time points. This can happen when there is no guessing because all individuals within a sample have mastered the required attributes by the time they get to the items. Nevertheless, the primary aim of Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018) was not to understand the identifiability of HMMs. The current article fills this gap by formally deriving identifiability conditions for RHMMs under a counterbalancing design.

In what we refer to as a counterbalancing design, a fixed number of questions are arranged into different orderings. These orderings correspond to the various versions of the test, and examinees are then categorized into groups depending on the version they receive. That is, counterbalancing here refers to an experimental design in which the same set of items is administered in different orders across examinees to control for potentially confounding factors, such as order and time effects. This design is relatively popular in both classroom and large-scale assessments (Johnson, Reference Johnson1992; Shute et al., Reference Shute, Hansen and Almond2008). While counterbalancing can be implemented in many ways, we focus only on those configurations that are relevant for ensuring model identifiability, which will be specified later in this article. An example is the experiment conducted by Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018), where participants responded to 50 items organized into five 10-item blocks, with a learning intervention administered between blocks. Here, the five equal-size partitions of the 50 items are referred to as booklets 1–5. Counterbalancing was implemented by rotating each of the five booklets to be administered first in the sequence (Table 1).

Table 1

Balanced block design in Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018)

Table 1 Long description

The table consists of six columns and six rows. The first row contains the header Time points, spanning columns two through six, labeled 1, 2, 3, 4, and 5. The first column lists the Groups from 1 to 5.

Row 2 (Group 1): Time 1 is Booklet 1, Time 2 is Booklet 2, Time 3 is Booklet 3, Time 4 is Booklet 4, Time 5 is Booklet 5.

Row 3 (Group 2): Time 1 is Booklet 2, Time 2 is Booklet 3, Time 3 is Booklet 4, Time 4 is Booklet 5, Time 5 is Booklet 1.

Row 4 (Group 3): Time 1 is Booklet 3, Time 2 is Booklet 4, Time 3 is Booklet 5, Time 4 is Booklet 1, Time 5 is Booklet 2.

Row 5 (Group 4): Time 1 is Booklet 4, Time 2 is Booklet 5, Time 3 is Booklet 1, Time 4 is Booklet 2, Time 5 is Booklet 3.

Row 6 (Group 5): Time 1 is Booklet 5, Time 2 is Booklet 1, Time 3 is Booklet 2, Time 4 is Booklet 3, Time 5 is Booklet 4.

The design follows a Latin Square pattern where each booklet appears exactly once in each group and once at each time point.

Following the experimental setup of Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018), this article focuses on the case with binary response data, but the proof can be extended to the case where items are polytomous. The remainder of this section introduces the RHMM framework and presents a new identifiability theorem for RHMMs specific to the counterbalancing design.

2.3.1 General model

Let there be G total groups where each group corresponds to a unique ordering of the booklets. Observe that, for a given time point, two groups may have different booklets. That is, examinees from separate groups may answer multiple different items at varying time points, resulting in different responses. In the HMM framework, this implies that the emission matrix $\mathbf {B}$ may change over time and we can define the multi-item emission matrix.

Definition 1 (Multi-item emission matrix).

Suppose there are J total items with $M_j$ representing the number of response options for item j. Let $\mathbb {J}_t \subset [J]$ to be the set of items administered at time t. The emission matrix of an item is defined as

$$ \begin{align*} \underset{M_j \times r}{\mathbf{B}_j} = \begin{bmatrix} b_{j,01} & \cdots& b_{j,0r} \\ \vdots & \cdots & \vdots \\ b_{j,M_{j-1},1} & \cdots & b_{j,M_{j-1},r} \end{bmatrix}, \end{align*} $$

and the emission matrix at time t is

$$ \begin{align*} \underset{(\prod_{j \in \mathbb{J}_t} M_j) \times r}{\mathbf{B}_t} = \mathop{{{{\circledast}}}}_{j \in \mathbb J_t} \mathbf{B}_j, \end{align*} $$

where $\mathop {{{{\circledast }}}}_{j \in \mathbb J_t}$ denotes the Khatri–Rao product (i.e., a column-wise Kronecker product) of matrices $\mathbf {B}_j$ for all $j \in \mathbb J_t$ (see Definition A.3 in Appendix A).

Let K be the total number of binary latent attributes and J be the total number of items. The conditional distribution of the response $\boldsymbol {X}$ in the emission matrix is written in terms of $\boldsymbol {\alpha } \in \{0, 1\}^K$ , as opposed to $\boldsymbol {Z}$ . That is, for the remainder of this section on counterbalancing, the number of latent states is $r = 2^K$ corresponding to all possible attribute profiles. Here, $\boldsymbol {\alpha } = (\boldsymbol {\alpha }_1, \dots , \boldsymbol {\alpha }_T)^\top \notag $ denotes the learning trajectory of individual i, where $\boldsymbol {\alpha }_t = (\alpha _{1t}, \dots , \alpha _{Kt})$ represents the binary attribute profile at time t. The element $\alpha _{kt} = 1$ indicates that the respondent possesses attribute k for $k \in [K]$ , and $\alpha _{kt} = 0$ otherwise. When the items are dichotomous, the emission probability conditional on the attribute profile $\boldsymbol {\alpha }_t$ is

(1) $$ \begin{align} P(\boldsymbol{X}_t = \boldsymbol{x}_t \mid \boldsymbol{\alpha}_t, \boldsymbol{\Omega}) = \prod_{j = 1}^{\mathbb{J}_t} b_{j, x_t, \boldsymbol{\alpha}_t}^{y_{jt}} (1 - b_{j, x_t, \boldsymbol{\alpha}_t})^{1 - y_{jt}}, \end{align} $$

where $\boldsymbol {\Omega }$ represents the model parameters, $\mathbb {J}_t$ denotes the set of items administered at time t with elements denoted by $\mathbb {j}_1,\dots , \mathbb {j}_{|\mathbb J_t|}$ , and $b_{j, x_t, \boldsymbol {\alpha }_t}$ is an element of an item-specific matrix, which denotes the probability of a response $x_t$ to item j at time t for a subject with attribute profile $\boldsymbol {\alpha }_t$ .

Equation (1) is a general form of the emission probability that contains the response probabilities for all $2^K$ latent classes for each item. The current model is unrestricted in the sense that there are no constraints on the emission probabilities. This lack of constraint, however, makes the general model susceptible to the label switching problem, a well-known issue in latent class modeling (Jasra et al., Reference Jasra, Holmes and Stephens2005; Redner & Walker, Reference Redner and Walker1984). Specifically, across different groups and time points, the latent classes associated with the items can vary in order or even differ entirely. However, this problem is mitigated in RHMMs, where some structure is introduced by an expert-specified $\boldsymbol {Q}$ matrix (see below) and imposed on the emission matrix.

The restrictions can be connected to alternative representations where CDMs are written as a mixture of generalized linear models (de la Torre, Reference de la Torre2011; Henson et al., Reference Henson, Templin and Willse2009; von Davier, Reference von Davier2008). The response probability for an item j is now modeled as

(2) $$ \begin{align} b_{j, 1, \boldsymbol{\alpha}_t} = P(X_{jt} = 1 \mid \boldsymbol{\alpha}_{t}, \boldsymbol{\beta}_{j}) = \Psi(\boldsymbol{\alpha}_{t}^{*\top} \boldsymbol{\beta}_{j}), \end{align} $$

where $\Psi (\cdot )$ is an arbitrary cumulative distribution function,

(3) $$ \begin{align} \boldsymbol{\alpha}_{t}^{*\top} = \left( 1, \alpha_{1t}, \dots, \alpha_{Kt}, \alpha_{1t} \alpha_{2t}, \dots, \alpha_{(K - 1)t} \alpha_{Kt}, \dots, \prod_{k = 1}^K \alpha_{kt} \right) \end{align} $$

is an alternative form of the binary vector $\boldsymbol {\alpha }_{t}$ , and

(4) $$ \begin{align} \boldsymbol{\alpha}_{t}^{*\top} \boldsymbol{\beta}_{j} = \beta_{j, 0} + \sum_{k = 1}^K \beta_{j, k} \alpha_{kt} + \sum_{k> k'} \beta_{j, kk'} \alpha_{kt} \alpha_{k't} + \dots + \beta_{j, 12 \dots K} \prod_{k = 1}^K \alpha_{kt}. \end{align} $$

The regression coefficients create a sparse vector $\boldsymbol {\beta }_{j}$ , where the nonzero elements represent the effects of latent skills or combination of skills on the response of item j. It can be shown that CDMs can be reparameterized to Equation (2) with different sparsity patterns of the vector $\boldsymbol {\beta }_{j}$ . Following notations from Chen et al. (Reference Chen, Culpepper and Liang2020), let $\boldsymbol {\delta }_j = (\delta _{j0},\dots ,\delta _{j,2^K-1})^\top \notag \in \{0, 1\}^{2K}$ be the structure vector which indicates the relevant main effects and/or interactions of the K attributes for item j. The elements $\delta _{jc} = 0$ correspond to $\beta _{jc} = 0$ and are inactive, whereas $\delta _{jc} = 1$ implies that $\beta _{jc}$ is active and non-zero. The matrix of structure vectors for a given time is denoted by $\Delta _t = (\boldsymbol {\delta }_{\mathbb {j}_1} \;\dots \; \boldsymbol {\delta }_{\mathbb {j}_{|\mathbb J_t|}} )^\top $ . The structure matrices impose restrictions on the emission probabilities by providing information about the relationship between items and attributes. Using this general model to describe CDMs, we have the following new identifiability theorem for RHMMs in the context of a counterbalancing design.

Theorem 3 (Identifiability for RHMMs in a counterbalancing design).

Suppose there are J items, K attributes, T time points, G groups, and L booklets. Let $\mathbb {J}_t \subset [J]$ be the set of items administered at time t and $\boldsymbol {\pi }_1 = (\pi _{1, 1}, \ldots , \pi _{1, 2^K})$ is the initial distribution. For an HMM where time-varying emission matrices $\mathbf {B}_t$ contain probabilities formed by the general model in Equation (2) and $T \geq 3$ , the parameters are identifiable if:

1. (a) $rank(\mathbf {A}) = 2^K$ ;

2. (b) for all $t \in [T]$ , $\boldsymbol {\Delta }_{t}$ has the form $\boldsymbol {\Delta }_{t} = \begin {pmatrix} \boldsymbol {D} \\ \boldsymbol {\Delta }' \end {pmatrix}$ after row swapping, where $\boldsymbol {\Delta }^{\prime }_{t}$ is any $|\mathbb {J}_{t}| \times 2^K$ binary matrix and $\boldsymbol {D}$ is a $K \times 2^K$ binary matrix with the structure

   $$\begin{align*}\boldsymbol{D} = \begin{pmatrix} 1 & 1 & 0 & \dots & 0 & \dots & 0 \\ 1 & 0 & 1 & \dots & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & & \vdots \\ 1 & 0 & 0 & \dots & 1 & \dots & 0 \end{pmatrix},\end{align*}$$

   which includes a diagonal pattern of 1’s for the main-effect coefficients.

3. (c) $G \geq L$ and there is a $t \in (1, T)$ such that all booklets show up once across G groups;

4. (d) $\pi _{1, l}> 0$ for all $l \in [2^K]$ .

Proof. The proof is found in Appendix B.

To clarify Part (c), observe, for instance, time point $3$ in Table 1. The design satisfies Part (c) because all five booklets appear once across the different groups. On the other hand, consider the design in Table 2. This design does have counterbalancing but it does not satisfy condition (c) in Theorem 3 because none of the time points in Table 2 contains all four booklets across the groups.

Table 2

Unidentifiable design for $T = 4$

Table 2 Long description

The table consists of five columns and five rows. The first column is labeled Groups, and the subsequent four columns are labeled Time points 1, 2, 3, and 4.

Row 1: Group 1 receives Booklet 2 at Time 1, Booklet 1 at Time 2, Booklet 3 at Time 3, and Booklet 4 at Time 4.

Row 2: Group 2 receives Booklet 4 at Time 1, Booklet 2 at Time 2, Booklet 3 at Time 3, and Booklet 1 at Time 4.

Row 3: Group 3 receives Booklet 4 at Time 1, Booklet 3 at Time 2, Booklet 2 at Time 3, and Booklet 1 at Time 4.

Row 4: Group 4 receives Booklet 3 at Time 1, Booklet 1 at Time 2, Booklet 4 at Time 3, and Booklet 2 at Time 4.

Part (d) of Theorem 3 requires the initial distribution to be strictly positive. That is, the samples are representative and include observations from all attribute profiles. Parts (a) and (b) assume that the emission and transition matrices are full rank. These are analogous to the requirement of a full-rank design matrix in linear regression. In practice, $\mathbf {A}$ and $\mathbf {B}$ may fail to be full rank when certain latent states are indistinguishable in their transition probabilities and response processes. An example for why this may happen is that the quality of items is bad such that they are unable to differentiate between the attribute profiles. Moreover, it should be noted that Part (b) is a sufficient condition for identifying model parameters and weaker conditions can be imposed for more general models (see, e.g., Culpepper, Reference Culpepper2023). Hence, different CDMs and RLCMs may require different assumptions for identifiability of the parameters. We will consider the case with the DINA as an example given that the DINA was used in Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018).

2.3.2 DINA model

We now introduce another binary matrix, which is more commonly used with the DINA, that also defines the underlying structure between the latent skills and items.

Definition 2 ( $\boldsymbol {Q}$ matrix).

Suppose there are J total items and K attributes. A $\boldsymbol {Q}$ matrix is a $J \times K$ binary matrix that can be written as $\boldsymbol {Q} = (\boldsymbol {q}_1, \dots , \boldsymbol {q}_J)^\top \notag \in \{0, 1\}^{J \times K}$ . The $\boldsymbol {q}_j^\top \notag = (q_{j1}, \dots , q_{jK})$ represents the j-th row of the Q-matrix. Element $q_{jk} = 1$ indicates that answering item j requires the mastery of skill k while $q_{jk} = 0$ means that skill k is not required for item j.

It is possible to derive the $\boldsymbol {Q}$ matrix when the structure matrix is known. Different from $\boldsymbol {\Delta }_{t}$ , the $\boldsymbol {Q}$ matrix only contains partial information about the item–attributes relationship (Fang et al., Reference Fang, Liu and Ying2019). The $\boldsymbol {Q}$ matrix does not provide information on how the attributes may interact with each other without specifying a particular CDM. Despite this limitation, the $\boldsymbol {Q}$ matrix is widely used in CDM research, and its identifiability has been extensively explored (Chen et al., Reference Chen, Liu, Xu and Ying2015; Gu & Xu, Reference Gu and Xu2019; Xu, Reference Xu2017; Xu & Shang, Reference Xu and Shang2018). We consider the DINA model along with the assumption of a known $\boldsymbol {Q}$ matrix.

For a DINA, the probability of a correct response is

(5) $$ \begin{align} b_{j, 1, \boldsymbol{\alpha}_t} = P(X_{jt} = 1 \mid \boldsymbol{\alpha}_t, \boldsymbol{q}_t, s_j, g_j) = (1 - s_j)^{\eta_{jt}} g_j^{1 - \eta_{jt}}, \end{align} $$

where $\eta _{jt} = 1$ if an individual mastered the required attributes for item j at time t and $\eta _{jt} = 0$ if at least one required attribute is not mastered. Additionally, we have $s_j = P(X_{jt} = 0 \mid \eta _{jt} = 1)$ and $g_j = P(X_{jt} = 1 \mid \eta _{jt} = 0)$ correspond to the slipping and guessing parameter of an item j. Now, Theorem 3 can be rewritten for the current setting with the DINA as

Corollary 1 (Identifiability for RHMMs-DINA in a counterbalancing design).

Suppose there are J items, K attributes, T time points, G groups, and L booklets. For an HMM where time-varying emission matrices $\mathbf {B}_t$ contain probabilities formed by the DINA model shown in Equation (5) and $T \geq 3$ , the parameters are identifiable if:

1. (a) $rank(\mathbf {A}) = 2^K$ ;

2. (b) $\boldsymbol {Q}_t$ takes the form (after a row permutation):

   $$\begin{align*}\boldsymbol{Q}_t = \begin{pmatrix} \boldsymbol{I}_K \\ \boldsymbol{Q}^*_t \end{pmatrix}_{J \times K}\end{align*}$$

   with $\boldsymbol {I}_K$ as a $K \times K$ identity matrix and $\boldsymbol {Q}^*_t$ can be any form except zero rows;

3. (c) $G \geq L$ and there exists $t \in (1, T)$ such that all booklets show up once across G groups;

4. (d) $\pi _{1, l}> 0$ for all $l \in [2^K]$ .

Proof. The proof is found in Appendix C.

Remark 1. Part (b) in Corollary 1 is specific to the DINA case and corresponds to Part (b) in Thereom 3. It is worth noting that this condition is related to the reliance of the proof on Kruskal’s theorem, as it imposes some requirements on $\boldsymbol {Q}_t$ when the DINA model is used (see Remark C.1 in Appendix C). Moreover, Part (b) in Corollary 1 also indicates that the $\boldsymbol {Q}_t$ matrix at each time point is complete (Chiu et al., Reference Chiu, Douglas and Li2009), which is a common condition for identifying the DINA model when there is only one time point (Gu & Xu, Reference Gu and Xu2019). The completeness condition guarantees that the $\boldsymbol {Q}_t$ matrix can identify all potential attribute patterns.

The identifiability theorems introduced thus far are limited to RHMMs. There are certain advantages to using RHMMs over conventional HMMs. Having either the $\boldsymbol {Q}$ or the structure matrix, some information can be obtained regarding the underlying structure of the latent states to help reduce the number of estimated parameters. It is worth noting that the identifiability conditions in Thereom 3 and Corollary 1 do not require the $\boldsymbol {Q}$ matrix or the structure matrix to be known but rather to take certain form. Whether or not this assumption is feasible is beyond the scope of this article; however, there has been substantial research on estimating the $\boldsymbol {Q}$ matrix more accurately (Balamuta & Culpepper, Reference Balamuta and Culpepper2022; Chen et al., Reference Chen, Culpepper, Wang and Douglas2018, Reference Chen, Culpepper and Liang2020; Köhn et al., Reference Köhn, Chiu, Oluwalana, Kim and Wang2025).

Another benefit of using RHMMs, as mentioned briefly before, is that using RHMMs with known $\boldsymbol {Q}$ matrices avoids the issue of label switching by mapping the responses $\boldsymbol {X}$ onto a more parsimonious collection of binary attributes. In conventional HMMs, the time-varying emission matrices are susceptible to the label switching problem. Specifically, for different groups, $\mathbf {B}_t$ can represent different orderings of the latent classes, or even entirely different underlying latent states across groups. Methods that enforce unique labeling across groups are thus necessary to ensure meaningful interpretations and inference.

2.3.3 Empirical investigation of counterbalanced design identifiability

Theorem 3 establishes sufficient conditions for identifying the parameters of RHMMs within a counterbalancing research design. Part (c) is particularly relevant to the research setting for which we want to provide additional insights. Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018) offer intuition on how the counterbalancing helps with identifying the DINA model parameters by exposing the items to individuals with varying mastery attribute profiles. This corresponds to the assumption of nonzero population proportions (Part (d) of Theorem 3). Moreover, counterbalancing addresses potential confounding factors such as test fatigue and systematic changes in item difficulty over time that are unrelated to changes in individuals’ mastery. For instance, one might observe the items administered at later time points appear more difficult (e.g., lower guessing parameters in the DINA model). Without counterbalancing, it would be unclear whether this phenomenon reflects test fatigue, genuine increases in item difficulty, proficiency loss, or some combinations of these factors. Consequently, the lack of counterbalancing is likely to introduce potential issues with parameter recovery of RHMMs.

To demonstrate the potential consequences of failing to meet the identifiability conditions, we conducted a simple simulation. The simulation design followed Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018) and included $350$ individuals, $50$ items, and $5$ time points. We employed the same Q-matrix from Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018) and used the DINA first-order HMM (DINA-FOHM; Chen et al., Reference Chen, Culpepper, Wang and Douglas2018) to generate the data. The DINA-FOHM is conceptually similar to the model described in the previous section which integrates DINA with RHMM. Data were generated using the HMCDM package (Zhang et al., Reference Zhang, Wang, Chen and Kwon2024). Two sets of item parameters were generated, which we call the normal and easier conditions. In the normal condition, the slipping and guessing parameters for all items were drawn from a $U(0.1, 0.2)$ distribution. In the easier condition, the guessing parameters for the items increased over time by generating from uniform distributions with progressively higher bounds. This condition accounts for a potential confound where improvements in performance are not solely a result of changes in students’ mastery levels. Additionally, we manipulated the research design such that there exists counterbalancing in one condition, similar to Wang et al. (Reference Wang, Yang, Culpepper and Douglas2018) and Table 1. The other situation does not have counterbalancing, so that all individuals receive same booklet of items at each time point and the item parameters change over time. Only one replication was conducted for each condition.

We estimated the item parameters using the generated data and computed the absolute difference between the estimated and true item parameters. Figure 1 presents the results for the normal condition. Observe that with counterbalancing, the item parameters were accurately recovered. On the other hand, when no counterbalancing is present, parameter recovery was poorer, particularly for some items where the absolute differences in the slipping parameters are around $0.6$ . This result is even more pronounced for the easier condition in Figure 2. Without counterbalancing, several items exhibited large absolute differences for both slipping and guessing parameters. With counterbalancing, the item parameters remained well estimated/recovered. From this small simulation, we argue that it is crucial to consider model’s assumptions to ensure identifiability and to enable meaningful interpretation of the model’s results.

Figure 1

Absolute difference between the estimated and true item parameters in the normal case (similar items over time).

Figure 1 Long description

A two-panel multi-plot labeled a and b. Both panels share a common Y axis titled Absolute difference and an X axis titled Situation with two categories: With counterbalancing in purple and No counterbalancing in light green.

Panel a, titled Slipping parameters, shows Y axis values from 0.0 to 0.6. The With counterbalancing group has a very low, tight distribution of data points near 0.0. The No counterbalancing group also has a low median, but contains several high-value outliers clustered between 0.4 and 0.6.

Panel b, titled Guessing parameters, shows Y axis values from 0.00 to 0.10. Both groups show a wider spread of data compared to panel a. The With counterbalancing group has a median absolute difference around 0.025 with a few outliers reaching up to 0.09. The No counterbalancing group has a slightly lower median near 0.02 but features a higher density of extreme outliers reaching above 0.10.

Figure 2

Absolute difference between the estimated and true item parameters where items get easier over time.

Figure 2 Long description

A two-panel figure labeled a and b. Both panels share a common Y-axis titled Absolute difference and an X-axis titled Situation with two categories: With counterbalancing in dark purple and No counterbalancing in light green.

Panel a, titled Slipping parameters, shows the Y-axis ranging from 0.0 to 0.6. The With counterbalancing group has a very low, tight distribution of points near 0.0. The No counterbalancing group also has a low median, but features several high-value outliers reaching up to approximately 0.65.

Panel b, titled Guessing parameters, shows the Y-axis ranging from 0.0 to 0.3. The With counterbalancing group has a median near 0.05 with points spread up to 0.15. The No counterbalancing group shows a higher median and a much wider vertical spread of points, with one significant outlier reaching above 0.3.

2.4 Multiple-group longitudinal design

Another commonly employed experimental design is the multiple-group longitudinal design. This design is widely used to evaluate the effectiveness of interventions across various domains, including psychology (van Hoorn et al., Reference van Hoorn, van Dijk, Meuwese, Rieffe and Crone2016), education (Madison & Bradshaw, Reference Madison and Bradshaw2018b; Tang & Zhan, Reference Tang and Zhan2021), and clinical trials (Khori et al., Reference Khori, Shalamzari, Isanejad, Alizadeh, Alizadeh, Khodayari, Khodayari, Shahbazi, Zahedi, Sohanaki, Khaniki, Mahdian, Saffari and Fayad2015). The design serves as a tool for investigating whether the effects observed can be attributed to the intervention itself rather than to extraneous factors. Figure 3 shows the research flowchart of two consecutive time points within the multiple-group longitudinal design when there are three groups. Individuals are randomly assigned to treatment and control groups and exposed to different interventions depending on the group membership. All groups of individuals take the same assessment at time $t - 1$ . Conceptually, this design is somewhat similar to the multi-armed bandit problem in reinforcement learning where there are multiple options with different reward distributions. In contrast to adaptive bandit algorithms, the experimental design fixes group membership at the outset as individuals are randomly assigned to a single intervention arm and remain in that group for the duration of the study. The control group functions as a baseline and typically receives a placebo or no intervention which should yield minimal to no improvement compared to the treatment groups. Changes in the performance are tracked over time both within each group and between groups. The intervention(s) is considered to be effective when the changes observed in the treatment group(s) are determined to be practically and/or statistically better than that of the control group.

Figure 3

Flowchart for two consecutive points of the multiple-group longitudinal design.

Consider, for example, the experiment in Tang and Zhan (Reference Tang and Zhan2021) that assessed the effectiveness of an intervention called CDF in promoting learning. The authors compared this feedback mode to two other interventions: correct–incorrect response feedback (CIRF) and no feedback. To implement this, they used a quasi-experimental design in which different classes within a school were randomly assigned to one of the feedback types. Students took three parallel tests at different time points, with the interventions applied between testing occasions. The data were analyzed using RHMMs but the identifiability was not investigated. Using this example as a guide, we adapt a more general HMM for this design and derive an identifiability theorem for HMM in this context.

In the experiment by Tang and Zhan (Reference Tang and Zhan2021), the control group consisted of individuals who did not receive any feedback after a test. Individuals in the control group are assumed to exhibit no improvement at all relative to those in the treatment groups. In this context, the assumption is equivalent to fixing the transition matrix of the control group to be the identity matrix. In contrast, one would expect the treatment groups to exhibit changes in the latent attributes. Their transition matrices would differ not only from the identity matrix, but potentially also from each other. Unlike the counterbalancing design, individuals from different groups answer the same items at each time point. However, the items can differ between time points. Therefore, we still consider an HMM where the emission matrix $\mathbf {B}_t$ changes over time and the transition matrix $\mathbf {A}$ differs between groups.

Following the same setting in Section 2.1, suppose there are N subjects, J items, T time points, and r latent states. Let $\mathbf {A}^{(g)}$ represent the transition matrix for group g with $\mathbf {A}^{(c)}$ denote the transition matrix for the control group. We have the following theorem for HMM identifiability with respect to this multiple-group longitudinal design.

2.4.1 Empirical investigation of multigroup longitudinal design identifiability

Similar to the counterbalancing design, we now want to present additional clarification for the identifiability conditions in Theorem 4. Specifically, we focus on Part (c), its interpretation, and how it facilitates the identification of HMM parameters. Part (c) assumes that individuals in the control group do not exhibit any changes in their latent classes. In essence, the control group serves as an anchor that, similar to the use of anchor items to fix the scale of item parameters in item response theory, ensures identification of both the parameters and latent classes. The existence of a control group helps to address the potential issue of label switching that arises when employing a Bayesian framework to estimate HMM parameters (Jasra et al., Reference Jasra, Holmes and Stephens2005). Additionally, the inclusion of a control group mitigates potential confounding factors introduced when different items are administered across time points. Without a control group with constant latent states, it might be difficult to evaluate whether changes in the estimated item parameters between time points are due to different mastery levels and/or genuine differences in item parameters.

We again performed a simple simulation to illustrate what would happen when the identifiability conditions are not met. The simulation design consisted of $350$ individuals, $30$ items, and $3$ time points. Consistent with the previous simulation settings, the data were generated under the DINA-FOHM model. More specifically, a single dataset was simulated with a control group whose latent classes remained constant over time. Using this dataset, we evaluated the role of the control group by comparing two scenarios: one in which the software incorporated the control group information, and another where the software did not. In the first condition, all individuals answered the same set of $10$ items at each time point. As existing algorithms do not allow the specification of individuals’ group membership, part (c) in Theorem 4 is violated, because all individuals are implicitly assumed to have a nonzero probability of transitioning to different latent classes over time. To reflect the notion of a non-changing control group within the constraints of current algorithms, one possible approach is to require those in the control group to answer all items at the first time points, rather than administering sets of $10$ items across three time points. Hence, in the second condition, the existence of a control group is specified by having a number of individuals answer all $30$ items at the first time point. Following the setup in Section 2.3.3, only one replication was conducted for each condition.

Data were fitted using the HMCDM and TDCM (Madison et al., Reference Madison, Jeon, Cotterell, Haab and Zor2025) packages, and the absolute difference between the estimated and true item parameters was recorded. The results are illustrated in Figure 4. When there is a control group, the item parameters were accurately recovered. On the other hand, without a control group, parameter recovery was poorer for both the slipping and guessing parameters. These findings further support the necessity to ascertain that the conditions for identification of models’ parameters are satisfied.

Figure 4

Absolute difference between the estimated and true item parameters for the multiple-group design with a control group.

Figure 4 Long description

A multi-panel figure with two side-by-side box plots labeled a and b. Both plots share a common Y-axis labeled Absolute difference ranging from 0.0 to 0.4 and an X-axis labeled Situation with two categories: With control group and No control group.

Panel a, titled Slipping parameters, shows a dark purple box for the With control group situation centered near 0.02 with low variance and a few outliers reaching 0.12. The No control group situation is represented by a light green box with a higher median around 0.07, a larger interquartile range, and several high outliers, including one data point exceeding 0.4.

Panel b, titled Guessing parameters, follows a similar trend. The With control group box is dark purple and compressed near the 0.0 baseline. The No control group box is light green, showing a median near 0.05 and a much wider spread of data points, with extreme outliers reaching above 0.4. In both panels, the presence of a control group significantly reduces the absolute difference and variability of the parameters.

5.1 Settings

We next report results from a Monte Carlo experiment to evaluate the performance of the Gibbs sampling algorithm in recovering the parameters under the identification conditions of Theorem 4. We conducted the simulation study with three time points, and two groups: one treatment and one control. Other design factors that were manipulated include number of latent attributes (i.e., $K = 2, 4, 6,$ and $8$ ), sample sizes per group (i.e., $N = 100, 500,$ and $1,000$ ), and number of items per time point ( $10, 20,$ and $100$ ). For each setting, we repeated the simulation study $500$ times.

By the identifiability condition, the transition matrix for the control group was specified as the identity matrix. For the treatment group, the transition matrix was simulated such that it is full rank and each row included entries with relatively higher probabilities compared to the others. This structure captures potential real-life situations where certain latent states might be more likely to transition into one another, indicating potential correlations among these states. Moreover, all elements of the treatment group’s transition matrix were nonzero. That is, we allowed transitions between any pairs of states to demonstrate the flexibility of the model. In practical settings, researchers may evaluate whether this transition matrix is realistic and, if necessary, impose additional constraints.

The emission matrices for the items were simulated to be full rank across time points and allow for potential separation between the latent classes. Specifically, all items were dichotomous, and the probabilities of a correct response of an item are simulated from a Beta distribution. The beta parameter is fixed at $10$ and the alpha parameter is a function of the latent classes. This corresponds to different latent classes having different means or locations. That is, the latent classes, on average, have different probabilities of answering a question correctly. The emission matrices across time points were generated from the same underlying distributions to reflect the use of parallel items similar to the conditions described in Tang and Zhan (Reference Tang and Zhan2021). The initial distributions for the two groups were set to be uniform from all possible latent classes. That is, the initial distributions have strictly positive elements, which satisfies conditions (d) of Theorem 4.

Lastly, we utilized the identifiability assumption to find good initial values for the Gibbs sampling algorithm. Particularly, we performed a latent class analysis on the control group to estimate the emission matrix. Based on the estimates, we identified the most likely latent classes for individuals in the treatment group at each time point. These classifications were then used to derive the transition matrix and initial distribution. Providing the algorithm with more informed starting values may help the Markov chain converge more quickly and, therefore, reduce the length of the burn-in period and overall computational burden. We used Markov chains of length $70,000$ with the first $15,000$ discarded as burn-in for all conditions.