2.2.1 A conditional model for RTs given item responses

To capture the effect of the ordinal response Y _ij on the RTs, we adopt a dummy coding approach. Specifically, we define a dummy variable vector as $\mathbf {Z}_{Y_{ij}}$ ,

(1) $$ \begin{align} \mathbf{Z}_{Y_{ij}}=\left(Z_{ijk}\right)_{1\leq k \leq K_{j}-1}, \quad \text{where} \quad Z_{ijk}=I_{(Y_{ij}=k)}. \end{align} $$

Let $\mathbf {z}_{y_{ij}}=(z_{ijk})_{1\leq k \leq K_{j}-1}$ denote an observation of $\mathbf {Z}_{Y_{ij}}$ . Under this encoding scheme, the relationship between y _ij and $\mathbf {z}_{y_{ij}}$ is given by

(2) $$ \begin{align} \begin{aligned} y_{ij} = 0 & \Leftrightarrow \mathbf{z}_{y_{ij}} = (0, 0, \dots, 0) \\ y_{ij} = 1 & \Leftrightarrow \mathbf{z}_{y_{ij}} = (1, 0, \dots, 0) \\ &\vdots \\ y_{ij} = K_{j} - 1 & \Leftrightarrow \mathbf{z}_{y_{ij}} = (0, 0, \dots, 1). \end{aligned} \end{align} $$

As noted above, the log-normal model is widely used for modeling RTs because it fits various data well and allows for reliable inference and future extensions. Therefore, we adopt the log-normal model in this study. To incorporate the influence of response categories on RTs, we extend the log-normal model by introducing category-specific parameters. The resulting conditional model for RTs given the observed response $Y_{ij}=y_{ij}$ is

(3) $$ \begin{align} \ln T_{ij}\vert y_{ij} =\beta_{j}-\varphi_{j}v_{i}+\mathbf{z}_{y_{ij}}{\boldsymbol{\delta}}_{j}+\epsilon_{ij}, \qquad \epsilon_{ij}{\sim} N(0, \sigma_j^2 ). \end{align} $$

Here, $v_{i}$ , $\beta _{j}$ , $\varphi _{j}$ , and $\sigma _{j}^2$ are introduced above; ${\boldsymbol {\delta }}_{j}=(\delta _{j1},\delta _{j2},\ldots ,\delta _{j(K_{j}-1)})^\top $ represents the item response-specific effect parameter for item j. The parameter ${\boldsymbol {\delta }}_{j}$ quantifies the systematic differences in RTs associated with each response category.

2.2.2 Item response model

In this study, we employ the widely-used GRM proposed by Samejima (Reference Samejima1969) to model the ordered categorical response variable $Y_{ij}$ . Under the GRM, the probability that test taker i selects response category $k~(k=0,\ldots ,K_j-1)$ for item j is given by

$$ \begin{align*}\nonumber p_{ijk}=P(Y_{ij}=k|\theta_i,a_j,b_j,{\pmb{\tau}}_j) = P(Y_{ij} \geq k|\theta_i,a_j,b_j,{\pmb{\tau}}_j) - P(Y_{ij} \geq k+1|\theta_i,a_j,b_j,{\pmb{\tau}}_j), \end{align*} $$

where

$$ \begin{align*} P(Y_{ij} \geq 0 \mid \theta_i, a_j,b_j,{\pmb{\tau}}_j) &= 1, \nonumber \\\ P(Y_{ij} \geq k \mid \theta_i, a_j,b_j,{\pmb{\tau}}_j) &= F(a_j \theta_i - b_j - \tau_{jk}), \quad k = 1, 2, \dots, K_j - 1, \nonumber \\ \ P(Y_{ij} \geq K_j \mid \theta_i, a_j,b_j,{\pmb{\tau}}_j) &= 0. \nonumber\end{align*} $$

Here, $\theta _i$ denotes the latent ability of test taker i, $a_j$ is the discrimination parameter, $b_j$ is the difficulty parameter, and ${\pmb {\tau }}_j=(\tau _{j1},\ldots ,\tau _{j(K_j-1)})$ is the thresholds for item j, which correspond to the cut-off points between the response categories. The threshold parameters must satisfy the monotonicity condition $\tau _{j1}<\tau _{j2}<\cdots <\tau _{j(K_{j}-1)}$ . In addition, to ensure identifiability between $b_j$ and $\boldsymbol {\tau }$ , we impose the constraint: $\sum _{k=1}^{K_j-1}\tau _{jk}=0 \quad \text {for } j=1,\dots ,M$ . The function $F(\cdot )$ refers to the cumulative distribution function (CDF) of either the logistic or the standard normal distribution. In this study, we use the standard normal CDF for convenience in developing the SAEM algorithm for model estimation.

Finally, based on the dummy variable vector $\mathbf {Z}_{Y_{ij}}$ in Equations (1) and (2), the probability mass function of the item response variable $Y_{ij}$ can be expressed as

(4) $$ \begin{align} P(Y_{ij}=y_{ij}|\theta_i,a_j,b_j,{\pmb{\tau}}_j)=p_{ij0}^{1-\sum_{k=1}^{K_j-1}z_{ijk}}\prod_{k=1}^{K_j-1} p_{ijk}^{z_{ijk}}, \quad \text{where } y_{ij}=0,\dots,K_j-1. \end{align} $$

2.2.3 A joint model for graded responses and RTs

Let $\boldsymbol {\zeta }_i=(\theta _i,v_i)^\top $ denote the latent trait vector of examinee i. For item j, let the GRM parameters be $\boldsymbol {\eta }_j = (a_j, b_j, \boldsymbol {\tau }_j)^\top $ , and let the RT model parameters be $\boldsymbol {\alpha }_j = (\beta _j, \varphi _j, \boldsymbol {\delta }_j^\top , \sigma _j)^\top $ . The complete parameter vector for item j is given by $\boldsymbol {\xi }_j = (\boldsymbol {\eta }_j^\top , \boldsymbol {\alpha }_j^\top )^\top $ .

Based on the conditional RT model in Equation (3) and the GRM in Equation (4), the conditional joint model of $Y_{ij}$ and $T_{ij}$ given ${\boldsymbol {\zeta }}_i$ can be written as

(5) $$ \begin{align} f(t_{ij},y_{ij}| {\boldsymbol{\zeta}}_i, {\boldsymbol{\xi}}_{j})=\frac{1}{t_{ij}}\phi\left(\frac{\ln{t_{ij}}-\beta_j+\varphi_jv_i-\mathbf{z}_{y_{ij}}{\boldsymbol{\delta}}_{j}}{\sigma_j}\right) P(Y_{ij}=y_{ij}|\theta_i,a_j,b_j,{\pmb{\tau}}_j), \end{align} $$

where $\phi (\cdot )$ is the probability density function of the standard normal distribution.

Following the hierarchical modeling framework proposed by van der Linden (Reference van der Linden2007), the population-level distribution of the latent variables ${\boldsymbol {\zeta }}_i$ is assumed to follow a bivariate normal distribution

(6) $$ \begin{align} {\boldsymbol{\zeta}}_i = (\theta_i, v_i)^\top \sim N_2(\boldsymbol{\mu}_p,\boldsymbol{\Sigma}_{p}), \end{align} $$

where $\boldsymbol {\mu }_p$ and $\boldsymbol {\Sigma }_{p}$ denote the mean vector and covariance matrix, respectively. For identification purposes, these parameters are typically specified as

$$ \begin{align*}{\boldsymbol{\mu}}_{p} = (0,0)^\top, \quad \boldsymbol{\Sigma}_p=\begin{pmatrix} 1 & \sigma_{\theta v} \\ \sigma_{\theta v}& 1 \end{pmatrix}.\end{align*} $$

Finally, based on Equations (5) and (6), a two-level hierarchical joint model for $Y_{ij}$ and $T_{ij}$ is given by

(7) $$ \begin{align} \left\{ \begin{aligned} f(t_{ij}, y_{ij} | {\boldsymbol{\zeta}}_i, {\boldsymbol{\xi}}_j) &= \frac{1}{t_{ij}} \phi \left( \frac{\ln{t_{ij}}-\beta_j + \varphi_j v_i - \mathbf{z}_{y_{ij}} {\boldsymbol{\delta}}_j}{\sigma_j} \right) P(Y_{ij}=y_{ij}|\theta_i,a_j,b_j,{\pmb{\tau}}_j), \\ f({\boldsymbol{\zeta}}_i|\sigma_{\theta v})&= \frac{1}{2\pi \sqrt{1-\sigma_{\theta v}^{2}}} \exp\left\{ -\frac{\theta_{i}^{2} + v_{i}^2 - 2\sigma_{\theta v}\theta_{i}v_{i}}{2(1-\sigma_{\theta v}^{2})} \right\}. \end{aligned} \right. \end{align} $$

In this article, the proposed joint model is referred to as the GR–RT model, as it accounts for both graded responses and response times. Notably, when ${\boldsymbol {\delta }}_{j}=\mathbf {0}$ , the GR–RT model reduces to the basic HM framework proposed by Klein Entink et al. (Reference Klein Entink, Kuhn, Hornke and Fox2009) for graded scoring items, where responses and RTs are conditionally independent. Moreover, since binary scoring is a special case of graded scoring, the GR–RT model provides a unified framework for analyzing test data that include both dichotomous and polytomous items.

The practical applicability of a measurement model depends on the feasibility of accurately estimating its parameters. The complex structure of the GR–RT model presents significant challenges for parameter estimation. Therefore, developing an effective computational method to estimate the parameters of the GR–RT model is a central focus of this study. In the following section, we present an SAEM algorithm designed to compute the MMLE of the GR–RT model.

3.2.1 EM algorithm

With the introduction of ${\mathbf {X}}_{ij}$ , the E-step of EM algorithm simplifies to working with sufficient statistics, and the M-step directly updates the estimate of $\boldsymbol {\Omega }$ using the updated values of these sufficient statistics. Below, we present the general procedure of the EM algorithm, with a detailed implementation provided in Appendix A of the Supplementary Material.

Let $\boldsymbol {\Omega }^{0}=(\boldsymbol {\xi }_{1}^{0},\ldots , \boldsymbol {\xi }_{M}^{0}, \sigma _{\theta v}^{0})$ be the initial values of $\boldsymbol {\Omega }$ , and $\boldsymbol {\Omega }^{r-1}=(\boldsymbol {\xi }_{1}^{r-1},\ldots , \boldsymbol {\xi }_{M}^{r-1}, \sigma _{\theta v}^{r-1})$ denote the parameter values updated at the end of the $(r-1)$ th iteration. The rth EM iteration consists of the following steps:

1. 1. E-step: Compute the conditional expectations of the sufficient statistics $S_{\boldsymbol {\eta }_{j}}$ , $S_{\boldsymbol {\alpha }_{j}}$ , and $ S_{\sigma _{\theta v}}$ given $\mathbf {y}$ , $\mathbf {t}$ , and $\boldsymbol {\Omega }^{r-1}$ :

   $$ \begin{align*} {s}_{\boldsymbol{\eta}_j}^{r}=E_{(\mathbf{x}_j, \boldsymbol{\theta})}(S_{\boldsymbol{\eta}_{j}}|\mathbf{y}, \mathbf{t},\boldsymbol{\Omega}^{r-1}), s_{\boldsymbol{\alpha}_j}^{r}=E_{\mathbf{v}}(S_{\boldsymbol{\alpha}_{j}}|\mathbf{y}, \mathbf{t},\boldsymbol{\Omega}^{r-1}), \quad \text{for } j=1,\ldots,M, \nonumber \end{align*} $$
   and
   $$ \begin{align*} s_{\sigma_{\theta v}}^{r}=E_{\boldsymbol{\zeta}}(S_{\sigma_{\theta v}}|\mathbf{y}, \mathbf{t},\boldsymbol{\Omega}^{r-1}).\nonumber \end{align*} $$

2. 2. M-step: Update the parameter estimate of $\boldsymbol {\Omega }$ by

   $$ \begin{align*}\boldsymbol{\eta}^r_j=\hat{\boldsymbol{\eta}}_j(s_{\boldsymbol{\eta}_j}^{r}), \boldsymbol{\alpha}^r_j=\hat{\boldsymbol{\alpha}}_j(s_{\boldsymbol{\alpha}_j}^{r}),\text{for } j=1,\ldots,M,\end{align*} $$
   and
   $$ \begin{align*}\sigma_{\theta v}^r=\hat{\sigma}_{\theta v}(s_{\sigma_{\theta v}}^{r}),\end{align*} $$

   where the expressions of $\hat {\boldsymbol {\eta }}_j(s_{\boldsymbol {\eta }_j}^{r})$ , $\hat {\boldsymbol {\alpha }}_j(s_{\boldsymbol {\alpha }_j}^{r})$ , and $\hat {\sigma }_{\theta v}(s_{\sigma _{\theta v}}^{r})$ are provided in Equations (A10)–(A13) in the Supplementary Material.

While the M-step in the above EM iteration is straightforward and avoids the need for tedious optimization procedures, the E-step involves complex multidimensional integrals. To address this challenge, we propose an SAEM algorithm, inspired by existing studies (Camilli & Geis, Reference Camilli and Geis2019; Lavielle & Mbogning, Reference Lavielle and Mbogning2014; Meng & Xu, Reference Meng and Xu2023), which replaces the E-step computation with stochastic approximation, thereby reducing the computational burden.

3.2.2 A dimensionality reduction computation of $E_{(\mathbf {x}_j, \boldsymbol {\theta })}(S_{\boldsymbol {\eta }_{j}} \mid \mathbf {y}, \mathbf {t}, \boldsymbol {\Omega }^{r-1})$

Using the law of total expectation (Kolmogorov, Reference Kolmogorov2018),

$$ \begin{align*} E_{(\mathbf{x}_j, \boldsymbol{\theta})}(S_{\boldsymbol{\eta}_{j}} \mid \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1}) = E_{\boldsymbol{\theta}} \left( E_{\mathbf{x}_{j}}(S_{\boldsymbol{\eta}_{j}} \mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1}) \mid \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1} \right).\nonumber \end{align*} $$

From Equation (14), the inner expectation $E_{\mathbf {x}_{j}}(S_{\boldsymbol {\eta }_{j}} \mid \boldsymbol {\theta }, \mathbf {y}, \mathbf {t}, \boldsymbol {\Omega }^{r-1})$ is given by

(17) $$ \begin{align} E_{\mathbf{x}_{j}}(S_{\boldsymbol{\eta}_{j}} \mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r - 1})=\left\{ \sum_{i = 1}^{N}\sum_{k = 1}^{K_{j}-1}\theta_{i}E(X_{ijk} \mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r - 1}), \sum_{i = 1}^{N}E({\mathbf{X}}_{ij}\mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r - 1}), \sum_{i = 1}^{N}\theta_{i}, \boldsymbol{\theta}^\top\boldsymbol{\theta}\right\}, \end{align} $$

where

(18) $$ \begin{align} \sum_{i=1}^{N}E({\mathbf{X}}_{ij}\mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1})=\left(\sum_{i=1}^N E(X_{ij1}\mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1}),\dots,\sum_{i=1}^N E(X_{ij(K_j-1)}\mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1})\right)^\top. \end{align} $$

According to the relationship between $X_{ijk}$ and $Y_{ij}$ defined in Equation (9), the conditional distribution of $X_{ijk}$ given $Y_{ij}=y_{ij}$ is given by

$$ \begin{align*}\nonumber X_{ijk}|y_{ij}\sim \left\{ \begin{aligned} & N(a_j\theta_i-b_j-\tau_{jk},1)I_{x_{ijk}<0}, \quad &\text{if } y_{ij} < k, \\ &N(a_j\theta_i-b_j-\tau_{jk},1)I_{x_{ijk}\geq0}, \quad &\text{if }y_{ij} \geq k. \end{aligned} \right. \end{align*} $$

Further, it can be derived that the conditional expectation of $X_{ijk}$ given $\boldsymbol {\theta }$ , $\mathbf {y}$ , $\mathbf {t}$ , and $\boldsymbol {\Omega }^{r-1}$ is

(19) $$ \begin{align} E(X_{ijk} \mid \boldsymbol{\theta}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1}) = \begin{cases} a_j^{r-1} \theta_{i} - b_j^{r-1} - \tau_{jk}^{r-1} - \frac{\phi(-a_j^{r-1} \theta_{i} + b_j^{r-1} + \tau_{jk}^{r-1})}{\Phi(-a_j^{r-1} \theta_{i} + b_j^{r-1} + \tau_{jk}^{r-1})}, & \text{if } y_{ij} < k, \\ a_j^{r-1} \theta_{i} - b_j^{r-1} - \tau_{jk}^{r-1} + \frac{\phi(-a_j^{r-1} \theta_{i} + b_j^{r-1} + \tau_{jk}^{r-1})}{1 - \Phi(-a_j^{r-1} \theta_{i} + b_j^{r-1} + \tau_{jk}^{r-1})}, & \text{if }y_{ij} \geq k. \end{cases} \end{align} $$

Finally, by substituting the expression for $E(X_{ijk} \mid \boldsymbol {\theta }, \mathbf {y}, \mathbf {t}, \boldsymbol {\Omega }^{r-1})$ in Equation (19) into Equations (17) and (18), the closed form of $E_{\mathbf {x}_{j}}(S_{\boldsymbol {\eta }_{j}} \mid \boldsymbol {\theta }, \mathbf {y}, \mathbf {t}, \boldsymbol {\Omega }^{r - 1})$ is obtained. Let $h_j^r({\boldsymbol {\theta }})=E_{\mathbf {x}_{j}}(S_{\boldsymbol {\eta }_{j}} \mid \boldsymbol {\theta }, \mathbf {y}, \mathbf {t}, \boldsymbol {\Omega }^{r-1})$ , then we have

(20) $$ \begin{align} E_{(\mathbf{x}_j, \boldsymbol{\theta})}(S_{\boldsymbol{\eta}_{j}} \mid \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1}) = E_{\boldsymbol{\theta}} \left(h_j^r({\boldsymbol{\theta}}) \mid \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1} \right). \end{align} $$

This reduces the expectation of $S_{\boldsymbol {\eta }_{j}}$ with respect to $\boldsymbol {\theta }$ and ${\mathbf {X}}$ to the expectation of $h_j^r({\boldsymbol {\theta }})$ with respect to $\boldsymbol {\theta }$ , thereby achieving dimensionality reduction.

3.2.3 General framework of the SAEM algorithm

Based on the E-step in Equation (20), and following Delyon et al. (Reference Delyon, Lavielle and Moulines1999), the rth SAEM iteration consists of the following three steps:

1. 1. S-step: Randomly sample $\boldsymbol {\zeta }_{i}^{l}=(\theta _i^l, v_i^l) (l=1,\ldots ,m_r)$ from $f(\boldsymbol {\zeta }_{i}|\mathbf {y},\mathbf {t},\boldsymbol {\Omega }^{r-1})$ , for $i=1,\ldots ,N$ . The detailed sampling procedure is given in Appendix B of the Supplementary Material.

2. 2. SA-step: Let $\boldsymbol {\theta }^l=(\theta _1^l,\ldots ,\theta _N^l)^\top $ and $\mathbf {v}^l=(v_1^l,\ldots ,v_N^l)^\top $ . Using the generated $\boldsymbol {\theta }^l$ and $\mathbf {v}^l$ , compute $S_{\boldsymbol {\alpha }_{j}}^{l}$ , $S_{\sigma _{\theta v}}^{l}$ , and $h_j^{r}(\boldsymbol {\theta }^{l})$ according to Equations 15–17, for $l=1,\ldots ,m_r$ . The sufficient statistics are then updated as follows:

   $$ \begin{align*} s_{\boldsymbol{\eta}_j}^{r}&=s_{\boldsymbol{\eta}_j}^{r-1}+\gamma_r\left(\frac{\sum_{l=1}^{m_r}h_j^{r}(\boldsymbol{\theta}^{l})}{m_r}-s_{\boldsymbol{\eta}_j}^{r-1}\right),\nonumber\\ s_{\boldsymbol{\alpha}_j}^{r}&= s_{\boldsymbol{\alpha}_j}^{r-1}+\gamma_r\left(\frac{\sum_{l=1}^{m_r}S_{\boldsymbol{\alpha}_j}^{l}}{m_r}-s_{\boldsymbol{\alpha}_j}^{r-1}\right), \text{ for } j=1,\ldots,M, \nonumber \end{align*} $$
   and
   $$ \begin{align*} s_{\sigma_{\theta v}}^{r}=s_{\sigma_{\theta v}}^{r-1}+\gamma_r\left(\frac{\sum_{l=1}^{m_r}S_{\sigma_{\theta v}}^{l}}{m_r}-s_{\sigma_{\theta v}}^{r-1}\right), \nonumber \end{align*} $$
   where $ \gamma _r>0$ , $\sum _{r=1}^{+\infty }\gamma _r=+\infty $ and $\sum _{r=1}^{+\infty }\gamma _r^2<+\infty $ .

3. 3. M-Step: Compute

   $$ \begin{align*}\boldsymbol{\eta}^r_j=\hat{\boldsymbol{\eta}}_j(s_{\boldsymbol{\eta}_j}^{r}),\boldsymbol{\alpha}^r_j=\hat{\boldsymbol{\alpha}}_j(s_{\boldsymbol{\alpha}_j}^{r}), \text{for } j=1,\ldots,M\end{align*} $$
   and
   $$ \begin{align*}\sigma_{\theta v}^r=\hat{\sigma}_{\theta v}(s_{\sigma_{\theta v}}^{r}),\end{align*} $$
   according to Equations (A10)–(A13) in the Supplementary Material.

Remark 2. A crucial element of our SAEM algorithm is the use of $h_j^r({\boldsymbol {\theta }})=E_{\mathbf {x}_{j}}(S_{\boldsymbol {\eta }_{j}} \mid \boldsymbol {\theta }, \mathbf {y}, \mathbf {t}, \boldsymbol {\Omega }^{r-1})$ , which eliminates the need to sample $\mathbf {x}_j$ during the S-step. This enhances computational efficiency and, more importantly, ensures that the estimates of $\tau _{jk}^r (k=1,\ldots ,K_j-1)$ obtained in the M-step reliably preserve the desired ascending order: $\tau ^r_{j1}<\tau ^r_{j2}<\cdots <\tau ^r_{j(K_j-1)}$ . In contrast, if $\mathbf {x}_j$ were randomly sampled during S-step, this order could be violated, as highlighted by Camilli and Geis (Reference Camilli and Geis2019). To illustrate this advantage, we first define $\tilde {\tau }_{jk,l}^r$ according to Equation (A11) in the Supplementary Material, as follows:

$$ \begin{align*} \tilde{\tau}_{jk,l}^r=\frac{- \sum_{i=1}^N E(X_{ijk}\mid \boldsymbol{\theta}^l, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1})}{N}+\frac{\left(\sum_{i=1}^{N}E({\mathbf{X}}_{ij}\mid\boldsymbol{\theta}^l, \mathbf{y},\mathbf{t},\boldsymbol{\Omega}^{r-1})\right)^\top\mathbf{1}_{K_j-1}}{N(K_{j}-1)}, \nonumber \end{align*} $$

where $\mathbf {1}_{K_j-1}$ denotes a $(K_j-1)\times 1$ column vector with all entries equal to $1$ . The term $\left (\sum _{i=1}^{N}E({\mathbf {X}}_{ij}\mid \boldsymbol {\theta }^l, \mathbf {y},\mathbf {t},\boldsymbol {\Omega }^{r-1})\right )^\top \mathbf {1}_{K_j-1}$ represents the summation of the elements within the vector $\sum _{i=1}^{N}E({\mathbf {X}}_{ij}\mid \boldsymbol {\theta }^l, \mathbf {y},\mathbf {t},\boldsymbol {\Omega }^{r-1})$ . Using $\tilde {\tau }_{jk,l}^r$ , the estimate $\tau _{jk}^r$ in M-step of our SAEM can be equivalently computed as

$$ \begin{align*}\tau_{jk}^r=\tau_{jk}^{r-1}+\gamma_r \left(\frac{\sum_{l=1}^{m_r}\tilde{\tau}_{jk,l}^r}{m_r}-\tau_{jk}^{r-1}\right).\end{align*} $$

Moreover, it can be easily shown that if the initial values of $\boldsymbol {\tau }_j$ satisfy $\tau _{j1}^{0}<\cdots <\tau _{j(K_{j}-1)}^{0}$ , the following order holds:

$$ \begin{align*}E(X_{ij1}\mid \boldsymbol{\theta}^{l}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1})>E(X_{ij2}\mid \boldsymbol{\theta}^{l}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1})>\dots>E(X_{ij(K_{j}-1)}\mid \boldsymbol{\theta}^{l}, \mathbf{y}, \mathbf{t}, \boldsymbol{\Omega}^{r-1}).\end{align*} $$

Thus, the order $\tilde {\tau }_{j1,l}<\dots <\tilde {\tau }_{j(K_j-1),l}$ is preserved, ensuring that $\tau _{j1}^r<\dots <\tau _{j(K_j-1)}^r$ holds.

It is important to note that the SAEM procedure does not require the number of response categories $K_j$ to be the same across items, which makes it well-suited for estimating the GR–RT model with mixed-scoring items. In addition, because the HM is nested within the GR–RT model, the SAEM algorithm developed here can be directly applied to estimate the HM.

Large-scale assessments often employ matrix-sampling designs, leading to sparse response matrices with non-random missing data. This missingness is not a flaw in data collection but rather a natural consequence of the test booklet assignment process. To facilitate the flexible application of the GR–RT model in this context, we extended the SAEM algorithm to accommodate non-random missing data by introducing an indicator matrix $\mathbf {D} = \{d_{ij}\}_{N \times M}$ , where $d_{ij} = 1$ denotes that item j is administered to respondent i, and $d_{ij} = 0$ otherwise. This extension explicitly accounts for the non-ignorable missingness mechanism, allowing for valid parameter estimation under complex large-scale assessment designs. The detailed procedures of the extended SAEM algorithm are provided in Appendix C of the Supplementary Material.

3.2.4 Some specifications for implementing the SAEM algorithm

Initialization: Several studies, including Delyon et al. (Reference Delyon, Lavielle and Moulines1999), Lavielle and Mbogning (Reference Lavielle and Mbogning2014), and Meng and Xu (Reference Meng and Xu2023), have demonstrated that the SAEM algorithm is robust to the selection of initial values. To evaluate this robustness of our SAEM algorithm, we conducted a simulation study comparing parameter estimation performance with three different sets of initial values: a set of randomly selected values, a set of common values, and the true values of $\boldsymbol {\Omega }$ . The parameter recovery results across these different initial values were nearly identical, providing additional evidence of the SAEM algorithm’s robustness in terms of initial value selection. For the sake of conciseness, the details of this simulation study are not reported here. In both our simulation and empirical studies, we used the set of common values as the initial values.

The specifications of $\gamma _r$ : The choice of step size $\gamma _r$ plays an important role in determining the convergence of the SAEM algorithm. Drawing on existing studies (Camilli & Geis, Reference Camilli and Geis2019; Lavielle & Mbogning, Reference Lavielle and Mbogning2014; Meng & Xu, Reference Meng and Xu2023), we adopt a two-stage step size sequence selection approach. Specifically, we set $\gamma _r=1$ during the initial R iterations (the first stage) and $\gamma _r=\frac {1}{(r-R)^\kappa }$ for $r> R$ (the second stage), where $0.5 < \kappa \leq 1$ . In this study, we use $\kappa = 1$ . The step size sequence is thus expressed as $\gamma _r = I_{(r \le R)} + I_{(r> R)}\,\frac {1}{r - R}$ . In addition, through simulation studies, we verify that the convergence of the SAEM algorithm is ensured when $R = 500$ .

The specifications of $m_r$ : Delyon et al. (Reference Delyon, Lavielle and Moulines1999) have demonstrated that the convergence of the SAEM algorithm is guaranteed with $m_r = 1$ . Notably, increasing $m_r$ does not improve computational accuracy but instead adds to the computational burden. Therefore, in practical applications of the SAEM algorithm, $m_r$ is typically set to $1$ .

Convergence check: The convergence of the SAEM iteration is monitored using $\Delta ^r=\|\boldsymbol {\Omega }^{r} - \boldsymbol {\Omega }^{r-1}\|_{\max }$ . Here, $\|\cdot \|_{\max }$ denotes the maximum of the absolute values within a vector. In this study, the SAEM iteration is recommended to terminate when $\Delta ^r<10^{-4}$ or the maximum number of iterations exceeds 1,000.

4.1.1 Design

This study investigated the performance of the proposed SAEM algorithm under mixed-format testing conditions that included both dichotomous and graded-response items. Without loss of generality, the graded items were specified to have four response categories. To assess the impact of item-type composition, five mixture proportions (MPs) of binary and graded-scoring items were considered, as shown in Table 2.

Table 2

Five mixture proportion (MP) structures of four-category graded and binary items

Table 2 Long description

The table consists of six columns and three rows.

* The header row contains the column labels: M P 1, M P 2, M P 3, M P 4, and M P 5.

* The second row, labeled Four-category graded items, contains the following percentages: 100 percent, 70 percent, 50 percent, 30 percent, and 0 percent.

* The third row, labeled Binary items, contains the following percentages: 0 percent, 30 percent, 50 percent, 70 percent, and 100 percent.

The data shows an inverse relationship between the two item types across the five mixture proportions.

Since the sample size (N) and test length (M) are two important factors influencing the accuracy of the parameter estimation, this simulation considered three sample sizes ( $N = 500$ , 1,000, and 2,000) and two test lengths ( $M = 20$ and 40). The three factors (sample size, test length, and MP) were crossed, resulting in $ 3 \times 2 \times 5 = 30$ simulation conditions. For each simulation condition, the item response and RT data were generated under two true models: the GR–RT model and the HM. The true parameter values were set as follows.

Let $M_g$ denote the number of graded-scored items, where $j = 1, \dots , M_g$ indexes the graded items and $j = M_g + 1, \dots , M$ indexes the binary items. Let $U(l_1,l_2)$ denote the uniform distribution between $l_1$ and $l_2$ , and $LN(\mu ,\sigma ^2)$ denote the log-normal distribution with parameters $\mu $ and $\sigma ^2$ .

1. 1. GR–RT model:

   1. (a) Graded-scored items:

      • • GRM: For $j=1, \dots , M_g$ , the discrimination parameter $a_j$ was generated from a uniform distribution, $a_j \sim U(0.5,2.5)$ , allowing for moderate to high discrimination. The location parameter $b_j$ was generated from a standard normal distribution, $b_j\sim N(0,1)$ , to center the items on the ability scale. These specifications follow previous studies (e.g., Zhang & Chen, Reference Zhang and Chen2024; Wang & Xu, Reference Wang and Xu2015; Meng & Xu, Reference Meng and Xu2023). The threshold parameters $\tau _{j1}, \tau _{j2}$ , and $\tau _{j3}$ were generated as follows: $\tau _{j1} \sim U(-1,-0.3)$ , $\tau _{j2} \sim U(-0.25,0.25)$ , and $\tau _{j3}=-\tau _{j1}-\tau _{j2}$ , such as $\sum _{k=1}^{3}\tau _{jk}=0$ . These specifications based on empirical results from the data of one of the test booklets of the PISA 2022 science literacy assessment.

      • • Log-normal RT model: For $j=1, \dots , M_g$ , $\beta _j \sim U(3,5)$ , ${\boldsymbol {\delta }}_j \sim U(-1,1)$ , $\sigma _j \sim U(0.3,0.7)$ , and $\varphi _j\sim LN(0,0.5)$ . These specifications based on empirical results from the data of one of the test booklets of the PISA 2022 science literacy assessment.

   2. (b) Binary-scored items:

      • • For $j=M_g+1,\dots , M$ , the true values of $a_j$ , $b_j$ , $\sigma _j$ , $\varphi _j$ , $\beta _j$ , and ${\boldsymbol {\delta }}_j$ were generated in the same way as the parameters for the graded-scored items described in (a). Note that there is no $\boldsymbol {\tau }$ -parameter, since the GRM is reduced to the 2PNO model.

2. 2. Hierarchical model:

   • • The true values for the parameters in the GRM and log-normal RT models were generated in the same manner as for GR–RT model. Note that there is no $\boldsymbol {\delta }$ -parameter in the HM.

3. 3. To account for a diverse range of correlations between ability and speed, the correlation coefficient $\sigma _{\theta v}$ was drawn from a uniform distribution, $\sigma _{\theta v} \sim U(-1,1)$ , for each simulation condition. Given $\sigma _{\theta v}$ , the person parameters were then sampled from the following bivariate normal distribution:

   $$ \begin{align*}(\theta_i,v_i)^T\sim N_2\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix},\begin{pmatrix} 1 & \sigma_{\theta v}\\ \sigma_{\theta v} & 1 \end{pmatrix}\right),\end{align*} $$
   for $i = 1,\ldots ,N$ .

As noted by previous studies (Molenaar et al., Reference Molenaar, Tuerlinckx and van der Maas2015a, Reference Molenaar, Tuerlinckx and van der Maas2015c), the HM can be estimated using Mplus. Following the reviewers’ suggestion, we successfully implemented Mplus to estimate the GR–RT model. Thus, in this simulation, Mplus was used as a comparison tool to evaluate the performance of our proposed method. The Mplus code used for this purpose is provided in Appendix D of the Supplementary Material.

A total of 100 test datasets were generated for each of the two true models across all 30 testing scenarios. The GR–RT and HM models were fitted to their respective simulated data using both the SAEM algorithm and Mplus. Parameter estimation accuracy was assessed using the root mean square error (RMSE) and absolute bias (ABias), both computed as averages over 100 replications and all items. Taking the “a-parameter” as an example, we denote $\hat {a}_{jg}$ as the estimate of $a_j$ obtained from the gth replication. The RMSE and ABias for the a-parameter were calculated as follows:

$$ \begin{align*} \text{RMSE}_{\hat{a}} =\frac{1}{M}\sum_{j=1}^M\sqrt{ \frac{1}{100} {\sum_{g=1}^{100}\left( \hat{a}_{jg}-a_{j}\right) ^{2}} }, \quad \text{ABias}_{\hat{a}} = \frac{1}{M}\sum_{j=1}^M \left| \frac{1}{100} \sum_{g=1}^{100} \left( \hat{a}_{jg} - a_{j} \right) \right|. \end{align*} $$

The computing time of both the SAEM algorithm and Mplus was also recorded for all simulation conditions to evaluate computational efficiency.