2.1 Sequential selection

In a sequential selection process, test takers select the first favorable statement from the $K$ (>1) statements in an item, the second favorable from the remaining $K-1$ statements, and so on, until all statements are ranked. A probability model models this process. For example, in a triple item ( $K$ = 3) containing three statements (Statements 1–3), the ranking pattern ( $1,2,3$ ) refers to Statement 1 as the most favorable choice and Statement 3 as the least favorable. Denote ${P}_{\left(1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)$ as the probability of the ranking pattern ( $1,2,3$ ) conditional on the column vector of the latent trait score(s) $\boldsymbol{\unicode{x3b8}}$ that the item measures. Then, based on the assumption of the sequential selection,

(1) $$\begin{align}{P}_{\left(1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)={P}_{\left(1|1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right){P}_{\left(2|2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right),\end{align}$$

where ${P}_{\left(1|1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)$ is the conditional probability of selecting Statement 1 from the three statements and ${P}_{\left(2|2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)$ is the conditional probability of selecting Statement 2 from the remaining Statements 2 and 3, both given that only one statement is selected from the available three or two statements.

Denote the selection of a statement as 1 and 0 otherwise. If multiple selections are allowed, rather than forcing the choice of one selection, eight selection patterns would exist on the three statements in the first step. However, only one selection is allowed in each step of the sequential selection. Thus, only three selection patterns, (1,0,0), (0,1,0), and (0,0,1), are legitimate. Let ${P}_{\left(1,2,3\right)}\left(1,0,0|\boldsymbol{\unicode{x3b8}} \right)$ be the conditional joint probability of three decisions made in the first step: selecting Statement 1 and declining Statements 2 and 3. Likewise, define ${P}_{\left(1,2,3\right)}\left(0,1,0|\boldsymbol{\unicode{x3b8}} \right)$ and ${P}_{\left(1,2,3\right)}\left(0,0,1|\boldsymbol{\unicode{x3b8}} \right)$ as the conditional joint probabilities of selecting Statements 2 and 3, respectively, in the first step while declining the others. Obviously, ${P}_{\left(1|1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)$ is proportional to ${P}_{\left(1,2,3\right)}\left(1,0,0|\boldsymbol{\unicode{x3b8}} \right)$ . Assuming the ratios among the three joint probabilities do not change after test takers are limited to the three selection patterns, we can normalize ${P}_{\left(1,2,3\right)}\left(1,0,0|\boldsymbol{\unicode{x3b8}} \right)$ and denote it as ${P}_{\left(1|1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)$ ,

(2) $$\begin{align}{P}_{\left(1|1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)=\frac{P_{\left(1,2,3\right)}\left(1,0,0|\boldsymbol{\unicode{x3b8}} \right)}{P_{\left(1,2,3\right)}\left(1,0,0|\boldsymbol{\unicode{x3b8}} \right)+{P}_{\left(1,2,3\right)}\left(0,1,0|\boldsymbol{\unicode{x3b8}} \right)+{P}_{\left(1,2,3\right)}\left(0,0,1|\boldsymbol{\unicode{x3b8}} \right)}.\end{align}$$

Similarly,

(3) $$\begin{align}{P}_{\left(2|2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)=\frac{P_{\left(2,3\right)}\left(1,0|\boldsymbol{\unicode{x3b8}} \right)}{P_{\left(2,3\right)}\left(1,0|\boldsymbol{\unicode{x3b8}} \right)+{P}_{\left(2,3\right)}\left(0,1|\boldsymbol{\unicode{x3b8}} \right)}.\end{align}$$

Similarly, for a pair item containing Statements 1 and 2, it is a one-step selection, and the normalized probability of selecting the first statement is

(4) $$\begin{align}{P}_{\left(1,2\right)}\left(\boldsymbol{\unicode{x3b8}} \right)=\frac{P_{\left(1,2\right)}\left(1,0|\boldsymbol{\unicode{x3b8}} \right)}{P_{\left(1,2\right)}\left(1,0|\boldsymbol{\unicode{x3b8}} \right)+{P}_{\left(1,2\right)}\left(0,1|\boldsymbol{\unicode{x3b8}} \right)},\end{align}$$

which is Andrich’s basic model for FC items (Andrich, Reference Andrich1995).

By assuming independent judgment of individual statements, the joint distributions can be simplified as the product of the probabilities of selecting or not selecting individual statements. For example, ${P}_{\left(1,2,3\right)}\left(1,0,0|\boldsymbol{\unicode{x3b8}} \right)$ is simplified as

(5) $$\begin{align}{P}_{\left(1,2,3\right)}\left(1,0,0|\boldsymbol{\unicode{x3b8}} \right)={P}_{(1)}\left(1|{\theta}_1\right){P}_{(2)}\left(0|{\theta}_2\right){P}_{(3)}\left(0|{\theta}_3\right),\end{align}$$

where ${\theta}_1$ , ${\theta}_2$ , and ${\theta}_3$ refer to the latent trait scores measured by Statements 1–3, respectively. The probability of selecting/not selecting a statement, ${P}_{(k)}\left(1|{\theta}_k\right)$ or ${P}_{(k)}\left(0|{\theta}_k\right)$ , is a function of the statement’s utility ${t}_k$ ( $k=$ 1, 2, or 3), and ${P}_{(k)}\left(0|{\theta}_k\right)=1-{P}_{(k)}\left(1|{\theta}_k\right)$ . Two popular processes are assumed to underlie the selection of individual statements: the dominance and ideal-point processes. The dominance process assumes a monotonic relationship between the degree of agreement on a statement and the underlying trait it measures. In contrast, the ideal-point process assumes that a statement has an optimal point on the trait continuum, which most likely endorses the statement, and the endorsement likelihood decreases as the trait level departs from this optimal point above or below. Different probability models have been proposed to model the two processes of judging statements. The notable model for the ideal-point process is the generalized graded unfolding model (GGUM; Roberts et al., Reference Roberts, Donoghue and Laughlin2000), whereas that for the dominance process is the two-parameter logistic IRT model (2PLM; Birnbaum, Reference Birnbaum, Lord and Novick1968). The response function of the 2PLM is

(6) $$\begin{align}{P}_{(k)}\left(1|{\theta}_k\right)=\frac{\exp \left({a}_k{\theta}_k+{b}_k\right)}{1+\exp \left({a}_k{\theta}_k+{b}_k\right)},\end{align}$$

where ${a}_k$ is the statement $k$ ’s slope/discrimination parameter and ${b}_k$ is the intercept. In Equation 6, the statement’s utility is represented by

(7) $$\begin{align}{t}_k={b}_k+{a}_k{\theta}_k+{\varepsilon}_k,\end{align}$$

where the errors ${\varepsilon}_k$ follow an independent and identical logistic distribution with a mean of 0 and a common variance across statements. Equation 6 is the cumulative distribution function of ${t}_k$ .

The IRT models under the sequential selection process were collectively named “Rank models” (de la Torre et al., Reference de la Torre, Ponsoda, Leenen and Hontangas2012; Fu et al., Reference Fu, Tan and Kyllonen2025). The specific models under the Rank models can be distinguished by the block size (i.e., the number of statements in an item) and the probability function of selecting individual statements. Notable models for the ideal-point process include Stark et al.’s (Reference Stark, Chernyshenko and Drasgow2005) multi-unidimensional pairwise preference (MUPP) model, combined with the dichotomous version of GGUM (MUPP-GGUM) for pairs and the Triplet-GGUM for triplets (Joo et al., Reference Joo, Lee and Stark2018; Lee et al., Reference Lee, Joo, Stark and Chernyshenko2019). For the dominance process, the MUPP-2PLM (Morillo et al., Reference Morillo, Leenen, Abad, Hontangas, De la Torre and Ponsoda2016) for pairs is a special case of the commonly known compensatory multidimensional two-parameter logistic IRT model (M2PLM; von Davier, Reference von Davier2008), and its response function is

(8) $$\begin{align}{P}_{\left(1,2\right)}\left(\boldsymbol{\unicode{x3b8}} \right)={\left[1+\exp \left({a}_2{\theta}_2+{b}_2-{a}_1{\theta}_1-{b}_1\right)\right]}^{-1}.\end{align}$$

The response function of the Triplet-2PLM (Fu et al., Reference Fu, Tan and Kyllonen2025; Zheng et al., Reference Zheng, Liu, Li, Xu, Zhang, Wei, Zhang, Liu and Huang2024) for triplets is

(9) $$\begin{align} &{P}_{\left(1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)\nonumber\\& \quad ={\left[1+\exp \left({a}_2{\theta}_2+{b}_2-{a}_1{\theta}_1-{b}_1\right)+\exp \left({a}_3{\theta}_3+{b}_3-{a}_1{\theta}_1-{b}_1\right)\right]}^{-1}{\left[1+\exp \left({a}_3{\theta}_3+{b}_3-{a}_2{\theta}_2-{b}_2\right)\right]}^{-1}.\end{align}$$

Equations 8 and 9 can be derived based on Equations 1–6. Note that in Equations 8 and 9, only one and two, respectively, of the intercepts ( ${b}_k$ ) can be identified. Brown (Reference Brown2016) has pointed out that the Rank models with logistic response functions for individual statements are consistent with Luce’s choice axioms (Luce, Reference Luce1959) and the Bradley–Terry model (Bradley, Reference Bradley1953; McFadden, Reference McFadden and Zarembka1973).

The estimations of the Rank models have been developed under the full-information approach. The MUPP-2PLM is a special case of traditional multidimensional IRT models; thus, the MML-EM and MCMC estimation methods developed for traditional multidimensional IRT models can be applied directly to estimate the MUPP-2PLM. Fu et al. (Reference Fu, Tan and Kyllonen2025) and Zheng et al. (Reference Zheng, Liu, Li, Xu, Zhang, Wei, Zhang, Liu and Huang2024) developed the MML-EM method for the Triplet-2PLM. Lee et al. (Reference Lee, Joo, Stark and Chernyshenko2019) estimated the MUPP-GGUM and Triplet-GGUM using MCMC, and Fu (Reference Fu2025) calibrated the MUPP-GGUM using the MML-EM method.

2.2 Thurstone’s law of pairwise comparison

According to Thurstone’s (Reference Thurstone1927) law of comparative judgment, for a pair of statements, statement selection is determined by the difference between the two statements’ utilities:

(10) $$\begin{align}{y}_{\left(1,2\right)}^{\ast }={t}_1-{t}_2,\end{align}$$

where ${t}_1$ and ${t}_2$ are independent conditional on latent traits that the two statements measure, based on the assumption of independent judgments of individual statements. If ${y}_{\left(1,2\right)}^{\ast}\ge 0$ , Statement 1 is selected; otherwise, Statement 2 is selected. Statement utilities can represent an ideal-point or dominance process. Zhang, Tu, et al.’s (Reference Zhang, Tu, Angrave, Zhang, Sun, Tay and Li2024) generalized that Thurstonian unfolding model (GTUM) is an example of an ideal-point model under Thurstone’s approach. The most popular model, the Thurstonian IRT (TIRT) model (Brown & Maydeu-Olivares, Reference Brown and Maydeu-Olivares2011), is a dominance model. In the TIRT model, each statement’s utility is represented by a one-factor model,

(11) $$\begin{align}{t}_k={u}_k+{\lambda}_k{\theta}_k+{\pi}_k,\end{align}$$

where ${u}_k$ is the intercept, ${\lambda}_k$ is the factor loading, ${\theta}_k$ is the latent trait score that the statement measures, and ${\pi}_k$ is the independent and identically normally distributed error with a mean of 0 and a variance of ${\psi}_k^2$ . Note that Equation 11 has the same basic form as Equation 7; however, we use different notations to facilitate comparison between the models under the two different modeling approaches. The difference between the two statement utilities is

(12) $$\begin{align}{y}_{\left(1,2\right)}^{\ast }={t}_1-{t}_2={\lambda}_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21}+{\pi}_1-{\pi}_2,\end{align}$$

where ${\gamma}_{21}=-\left({u}_1-{u}_2\right)$ is the threshold. The probability of selecting Statement 1 is the probability of ${y}_{\left(1,2\right)}^{\ast }>0$ and modeled by a normal ogive function

(13) $$\begin{align}{P}_{\left(1,2\right)}^{\ast}\left(\boldsymbol{\unicode{x3b8}} \right)=\Phi \left(\frac{\lambda_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21}}{\sqrt{\psi_1^2+{\psi}_2^2}}\right)\approx \frac{\exp \left[\displaystyle\frac{1.702\left({\lambda}_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21}\right)}{\sqrt{\psi_1^2+{\psi}_2^2}}\right]}{1+\exp \left[\displaystyle\frac{1.702\left({\lambda}_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21}\right)}{\sqrt{\psi_1^2+{\psi}_2^2}}\right]},\end{align}$$

where $\Phi$ is the standard normal cumulative function. For pairs, the model can be identified by setting the error variance of ${y}_{\left(1,2\right)}^{\ast }$ to a constant, for example, ${\psi}_1^2+{\psi}_2^2=1$ . The first part of Equation 13 is the TIRT model’s item response function, a special case of the compensatory two-parameter normal ogive multidimensional IRT model (McDonald, Reference McDonald, van der Linden and Hambleton1997). The approximation between the normal ogive function and the logistic function in Equation 13 is well known (Baker & Kim, Reference Baker and Kim2004). The second part of Equation 13 is equivalent to the MUPP-2PLM’s item response function (Equation 8) by parameter transformations. Furthermore, if we assume ${\pi}_1$ and ${\pi}_2$ follow a Weibull distribution so that ${\pi}_1-{\pi}_2$ follow a logistic distribution, the probability of ${y}_{\left(1,2\right)}^{\ast}\ge 0$ can be directly modeled by the MUPP-2PLM.

For items with more than two statements (i.e., $K>2$ ), the TIRT model recodes a ranking pattern into $\left(\begin{array}{c}K\\ {}2\end{array}\right)$ ranking patterns of pairs. For a triplet, the response pattern of the three statements (Statements 1–3) can be separated into three statement pairs: Pair 1 (1, 2), Pair 2 (1, 3), and Pair 3 (2, 3). The response functions (marginal probabilities of ranking patterns) for the three pairs are

(14) $$\begin{align}{P}_{1,2}^{\ast }=\Phi \left(\frac{\lambda_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21}}{\sqrt{\psi_1^2+{\psi}_2^2}}\right),\end{align}$$

(15) $$\begin{align}{P}_{1,3}^{\ast }=\Phi \left(\frac{\lambda_1{\theta}_1-{\lambda}_3{\theta}_3-{\gamma}_{31}}{\sqrt{\psi_1^2+{\psi}_3^2}}\right),\end{align}$$

(16) $$\begin{align}{P}_{2,3}^{\ast }=\Phi \left(\frac{\lambda_2{\theta}_2-{\lambda}_3{\theta}_3-{\gamma}_{32}}{\sqrt{\psi_2^2+{\psi}_3^2}}\right),\end{align}$$

where ${\gamma}_{31}={\gamma}_{21}+{\gamma}_{32}$ .Footnote ¹  Equations 14–16 can be identified by fixing the error variance of one statement, for example, setting the error variance of the first statement ( ${\psi}_1^2$ ) to 1. The TIRT model considers the local dependence among the three pairs stemming from a triplet by estimating the covariances among the three pairs’ ${y}^{\ast }$ s: $\operatorname{cov}\left({y}_{\left(1,2\right)}^{\ast },{y}_{\left(1,3\right)}^{\ast}\right)={\psi}_1^2$ , $\operatorname{cov}\left({y}_{\left(1,2\right)}^{\ast },{y}_{\left(2,3\right)}^{\ast}\right)=-{\psi}_2^2$ , and $\operatorname{cov}\left({y}_{\left(1,3\right)}^{\ast },{y}_{\left(2,3\right)}^{\ast}\right)={\psi}_3^2$ . The TIRT model is set up with mean and covariance structures under structural equation modeling (SEM) and estimated using a limited information method, unweighted least squares with mean and variance-corrected Satorra–Bentler goodness-of-fit tests (ULSMV). However, local dependence is considered only in item parameter estimation; for estimating latent trait scores, local independence is currently assumed in Equations 14–16.

The MUPP-2PLM model assumes local independence of item responses to pairs and is estimated using a full-information method. Thus, this model is unsuitable for estimating an item with more than two statements when it is separated into pairs. Some researchers have proposed adding additional random parameters to the logistic item response functions to account for dependence among pairs, as seen in Bürkner et al. (Reference Bürkner, Schulte and Holling2019, Equation 4) and Zhang, Tu, et al. (Reference Zhang, Tu, Angrave, Zhang, Sun, Tay and Li2024, Equation 14). However, the additions make the models more complicated and harder to estimate (especially for maximum likelihood estimation).

For items with more than two statements ( $K>2$ ), the difference in the modeling approaches between the Rank and Thurstonian models becomes apparent. For the Thurstonian models, the probability of a ranking pattern ( ${s}_1,{s}_2,\dots, {s}_K$ ), where ${s}_k$ denotes the statement ID with rank $k$ , is the joint distribution of selecting the first statements in the $K-1$ pairs, ( ${s}_1,{s}_2$ ), ( ${s}_2,{s}_3$ ), …, ( ${s}_{K-1},{s}_K$ ). For example, for a triplet, the probability of a ranking pattern (1, 2, 3) under the TIRT model is the joint probability of selecting the first statements in the two pairs, (1,2) and (2,3), and given by (Maydeu-Olivares, Reference Maydeu-Olivares1999)

(17) $$\begin{align}{P}_{\left(1,2,3\right)}^{\ast}\left(\boldsymbol{\unicode{x3b8}} \right)={P}^{\ast}\left({y}_{\left(1,2\right)}^{\ast }>0\cap {y}_{\left(2,3\right)}^{\ast }>0|\boldsymbol{\unicode{x3b8}} \right)=\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}\phi \left({\mathbf{y}}^{\ast }:{\mathbf{C}\boldsymbol{\unicode{x3bc}}}_{\mathbf{t}},{\mathbf{C}\boldsymbol{\Sigma}}_{\mathbf{t}}{\mathbf{C}}^{\prime}\right)d{\mathbf{y}}^{\ast },\end{align}$$

where $\mathbf{C}$ is the 2-by-3 contrast matrix $\mathbf{C}=\left(\begin{array}{l}1\kern0.5em -1\kern0.5em 0\\ {}\begin{array}{ccc}0& 1& -1\end{array}\end{array}\right)$ ; $\mathbf{t}=\left({t}_1,{t}_2,{t}_3\right)$ ; ${\boldsymbol{\unicode{x3bc}}}_{\mathbf{t}}$ and ${\boldsymbol{\Sigma}}_{\mathbf{t}}$ are the mean vector and covariance matrix of $\mathbf{t}$ ; ${\mathbf{y}}^{\ast }=\left({y}_{\left(1,2\right)}^{\ast },{y}_{\left(2,3\right)}^{\ast}\right)$ follows a bivariate normal distribution with a mean vector ${\mathbf{C}\boldsymbol{\unicode{x3bc}}}_{\mathbf{t}}$ and covariance matrix ${\mathbf{C}\boldsymbol{\Sigma}}_{\mathbf{t}}{\mathbf{C}}^{\prime }$ ; $\phi$ is the bivariate normal density function. It is obvious that Equation 17 and Equation 1 (the Rank models’ item response function for a triplet) are different: For the Rank models, Equation 1 can be rewritten using ${y}^{\ast }$ s as

(18) $$\begin{align}{P}_{\left(1,2,3\right)}\left(\boldsymbol{\unicode{x3b8}} \right)=P\left({y}_{\left(1,2\right)}^{\ast }>0\cap {y}_{\left(1,3\right)}^{\ast }>0|\boldsymbol{\unicode{x3b8}} \right)P\left({y}_{\left(2,3\right)}^{\ast }>0|\boldsymbol{\unicode{x3b8}} \right).\end{align}$$

Equation 17 represents a full joint distribution of a ranking pattern (1, 2, 3), while Equation 18 simplifies this joint distribution by separating it into two independent distributions. Specifically, in Equation 18, the comparison between Statements 2 and 3 is independent of Statement 1, which indicates that the Rank models hold a stronger assumption and, thus, are simpler than the Thurstonian models. Compared Equation 9 to Equations 14–16, the TIRT model for a triplet has two more parameters (i.e., error variances of two statements) than the Triplet-2PLM.

The FC item format discussed so far is dichotomous; test takers must select one statement from two or more statements. Sometimes, this format may be challenging because a test taker may like two statements equally. In addition, this format does not capture the degree to which a test taker prefers one statement over another. In contrast, Likert scales use graded responses (e.g., “strongly disagree,” “disagree,” “neutral,” “agree,” “disagree”) to gain more information on a statement and, thus, may generate more reliable trait scores than the dichotomous FC item format (Bürkner, Reference Bürkner2022). To address this issue, Brown and Maydeu-Olivares (Reference Brown and Maydeu-Olivares2018) extended the TIRT model to accommodate the graded FC item format (e.g., “I strongly like Statement 1,” “I like Statement 1,” “I like both statements equally,” “I like Statement 2,” “I strongly like Statement 2”). The graded TIRT model follows the approach used by the graded response model (GRM; Samejima, Reference Samejima1969) and employs $M-1$ thresholds, ${\gamma}_{21m}$ ( $m=1,\dots, M-1$ ), to divide the continuous difference of item utilities, ${y}_{\left(1,2\right)}^{\ast }$ , into $M$ ordinal response categories. Denote a test taker’s response as ${y}_{\left(1,2\right)}$ . Then, the item response function of ${y}_{\left(1,2\right)}$ is,

(19) $$\begin{align}\mathrm{for}\;m=1,\ {P}_{\left(1,2\right)}^{\ast}\left({y}_{\left(1,2\right)}=1|\boldsymbol{\unicode{x3b8}} \right)=\Phi \left(\frac{\lambda_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{211}}{\sqrt{\psi_1^2+{\psi}_2^2}}\right);\end{align}$$

(20) $$\begin{align}\mathrm{for}\kern0.24em m=2,\dots,\ M-1,\ {P}_{\left(1,2\right)}^{\ast}\left({y}_{\left(1,2\right)}=m|\boldsymbol{\unicode{x3b8}} \right)=\Phi \left(\frac{\lambda_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21m}}{\sqrt{\psi_1^2+{\psi}_2^2}}\right)-\Phi \left(\frac{\lambda_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21\left(m-1\right)}}{\sqrt{\psi_1^2+{\psi}_2^2}}\right);\end{align}$$

(21) $$\begin{align}\mathrm{for}\;m=M,\ {P}_{\left(1,2\right)}^{\ast}\left({y}_{\left(1,2\right)}=M|\boldsymbol{\unicode{x3b8}} \right)=1-\Phi \left(\frac{\lambda_1{\theta}_1-{\lambda}_2{\theta}_2-{\gamma}_{21M}}{\sqrt{\psi_1^2+{\psi}_2^2}}\right).\end{align}$$

Zhang, Tu, et al. (Reference Zhang, Tu, Angrave, Zhang, Sun, Tay and Li2024) used the same approach to handle the graded FC item format in their GTUM. Although theoretically, the graded Thurstonian models for FC items can calibrate items with any block sizes, the graded response format seems only practically feasible for pair items. Brown and Maydeu-Olivares (Reference Brown and Maydeu-Olivares2018) and Zhang, Luo, et al. (Reference Zhang, Luo and Li2023) reported improved reliabilities for using the graded format in FC items compared to the dichotomous format.

4.1 Applications to other item formats

Besides the rank item format (RANK), the MGPPM can be easily applied to other item formats or scoring schemes. Below, we discuss the applications of the MGPPM to four alternative formats: PICK (i.e., pick the statement that is most descriptive of a test taker), MOLE (i.e., choose the statements that are the most and least descriptive of a test taker) (de la Torre et al., Reference de la Torre, Ponsoda, Leenen and Hontangas2012), graded response, and a newly proposed hybrid forced-choice/Likert (HFCL; Williams et al., Reference Williams, Figueiras and Roohr2025).

4.1.1 PICK and MOLE

The PICK and MOLE formats are alternatives to the RANK format, designed to reduce test takers’ cognitive burden for FC items with larger block sizes. See Hontangas et al. (Reference Hontangas, de la Torre, Ponsoda, Leenen, Morillo and Abad2015) for examples of these two formats in real questionnaires. The response data from the PICK and MOLE formats can be considered equivalent to that from the RANK format with planned missing data. Thus, for an IRT model, the item response functions for the PICK and MOLE formats can be derived from that for the RANK format by accounting for planned missing data (Yousfi, Reference Yousfi2025). That is, the response function of an item response ${c}^{\ast }$ from the PICK or MOLE format can be written as

(36) $$\begin{align}P\left({c}^{\ast }|\boldsymbol{\unicode{x3b8}} \right)=\sum \limits_{c\in \boldsymbol{\unicode{x3c8}} \left({c}^{\ast}\right)}P\left(c|\boldsymbol{\unicode{x3b8}} \right),\end{align}$$

where $\boldsymbol{\unicode{x3c8}} \left({c}^{\ast}\right)$ is the set of full responses from the RANK format that are not in contradiction with item response ${c}^{\ast }$ from the PICK or MOLE format. For example, for a tetrad with Statements 1–4, the probability of choosing Statement 1 from the PICK format is the sum of the probabilities of the six ranking patterns with Statement 1 as rank 1 from the RANK format, and the probability of selecting Statements 1 and 2 as the most and least descriptive, respectively, from the MOLE format is the sum of the probabilities of the four ranking patterns with Statements 1 and 2 as ranks 1 and 4, respectively, from the RANK format.

4.1.2 Graded response

The MGPPM shows great flexibility in handling the graded FC item format. For example, for a pair item with five graded response categories, we can define the scoring function for both statements as shown in Table 2. Then, by treating each response category as a ranking pattern, the MGPPM for a graded pair item is defined by Equations 23–25. If we use the primary option (Equation 23), this model has the same number of model parameters as the MUPP-2PLM. If we use the alternative option (Equation 35), this model has the same number of model parameters as the graded TIRT model in Brown and Maydeu-Olivares (Reference Brown and Maydeu-Olivares2018), and their relationship is similar to that between the traditional GPCM and GRM.

Table 2 Scoring function for a pair item with five graded response categories

4.2 Relation to Thurstone’s models and rank models

The basic analysis unit of the ranking pattern approach is the utility of a statement at a rank level. Thus, the utility is within the context of comparing with other statements. In contrast, the basic analysis unit of the Rank and Thurstone’s models is the utility of a statement that is assumed to be independent of other statements. The utility of a statement is arguably constant no matter which test format (e.g., Likert vs. forced choice) or context (e.g., a statement in different FC items) this statement is administered with. This difference has significant implications for test assemblies, which we will explore further in the discussion section.

For all models, the probability of a ranking pattern is based on the joint distribution of the ranking pattern. However, the Rank models and primary MGPPM make assumptions to simplify the joint distribution. The Rank models decompose the joint probability into the product of $K-1$ independent probabilities, where statements are ranked sequentially. The primary MGPPM assumes that the differences in the probabilities of a statement’s two consecutive rankings are constant for a test taker, while the alternative MGPPM does not make this assumption. However, these assumptions only affect items with more than two statements. For pair items, the item response functions are the same or equivalent for all models, given that their formulations of statements’ utilities are the same or equivalent. For items with block sizes larger than two, the Rank-2PLM and the primary MGPPM have the same number of model parameters, while the TIRT model and the alternative MGPPM have the same number of model parameters and are more complex than the previous two models.

All three approaches provide flexible formulations to accommodate block sizes and cognitive processes underlying statement judgments (i.e., ideal point and dominance). Although not discussed here, the ranking pattern approach can handle the ideal-point process by using an ideal-point formulation of the utility of a statement at a rank. By modeling the utility at the ranking pattern level, the ranking pattern approach is more flexible in incorporating different assumptions into models than the other two approaches. Although we only discuss the two options for representing a ranked statement’s utility, many other options are possible. Using the ranking pattern approach, researchers can easily test different assumptions about a ranking process by adjusting the scoring function of a statement. For instance, one may assume that the advancement from rank 2 to 1 represents twice the utility increase from rank 3 to 2 (rather than once, as in the primary MGPPM). Additionally, the ranking pattern approach can be applied directly to the graded response format and rating score pattern, as previously demonstrated. Although Thurstone’s approach can estimate the graded response format, it cannot handle rating score patterns, whereas the sequential selection approach cannot handle both due to its modeling assumption.

The TIRT model is traditionally estimated using a limited information method under SEM, while the Rank models and MGPPM are estimated using a full information method, such as the maximum marginal likelihood estimation (MMLE). The TIRT can be estimated using the MMLE with Equation 17 as the response function. Frick (Reference Frick2023) employed Equation 17 and maximum likelihood estimation to estimate person parameters with known item parameters, thereby avoiding the local independence assumption in the person parameter estimation of the TIRT model. However, we have not found a study on the MMLE of the TIRT model’s item parameters. Under some conditions, a limited information method may be inferior to a full information method (Forero & Maydeu-Olivares, Reference Forero and Maydeu-Olivares2009). Nevertheless, all models can be easily set up for Bayesian estimations using the MCMC sampling method. The MGPPM is more elegant than the TIRT and Rank-2PLM as its conditional logit functions (e.g., Equations 32 and 33) are more concise and straightforward.

Like the traditional IRT models for cognitive tests, the models under the ranking pattern approach are statistical models to estimate the probabilities of item responses, rather than modeling a cognitive decision process or strategy that test takers are assumed to use to rank the statements in an item, because it is too cognitively demanding for test takers to consider all possible ranking patterns for FC items with a block size larger than two. Test takers likely use a heuristic strategy to simplify the ranking process. For triplets, Sass et al. (Reference Sass, Frick, Reips and Wetzel2020) showed evidence that test takers used pairwise comparisons to rank statements in a think-aloud task. However, this study also reported that 24% of the participants had difficulty remembering the information of the three statements in a triplet when making pairwise comparisons. It is conceivable that more test takers may find the pairwise comparison approach cognitively challenging when ranking four or more statements. In this respect, the sequential selection approach is more likely to be a strategy most test takers apply when responding to FC items with large block sizes. Further research is needed on the cognitive decision processes underlying FC items.

5.1 MML-EM estimation

Suppose in an FCQ, a total of $I$ test takers respond to $J$ items. Let ${x}_{ij}$ denotes test taker $i$ ’s nominal score (i.e., ranking pattern ID $c$ ) on item $j$ (Table 1); $\mathbf{X}$ is the $I\times J$ score matrix representing the total observed data. $\boldsymbol{\unicode{x3b8}}$ represents a random sample of latent trait score vectors from the test-taker population with a standard multivariate normal distribution with a $D\times D$ correlation matrix ${\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}$ . Fixing the population means and standard deviations of $\boldsymbol{\unicode{x3b8}}$ is a way to identify the model. ${\mathbf{a}}^{\#}$ and ${\mathbf{b}}^{\#}$ are the column vectors of slopes and intercepts, respectively. ${\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}$ is the intertrait correlation matrix. Assume that ${x}_{ij}$ is locally independent conditional on item parameters and $\boldsymbol{\unicode{x3b8}}$ . Then, the observed data likelihood is a marginal likelihood

(37) $$\begin{align}L\left({\mathbf{a}}^{\#},{\mathbf{b}}^{\#},{\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}|\mathbf{X}\right)=\prod \limits_{i=1}^I\left[{\int}_{\mathbf{R}}\prod \limits_{j=1}^J\prod \limits_{c=1}^C{P_{ij c}}^{\mathrm{I}\left({x}_{ij}=c\right)}\phi \left(\boldsymbol{\unicode{x3b8}} |{\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}\right)d\boldsymbol{\unicode{x3b8}} \right],\end{align}$$

where $\mathbf{R}$ is the D-dimensional area of integration, ${P}_{ijc}$ is the probability of a coded score (ranking pattern) $c$ ( $c=$ 1,…, $C$ ) (Equation 25), $\mathrm{I}\left({x}_{ij}=c\right)$ is the indicator function with a value of 1 if ${x}_{ij}=c$ and 0 otherwise, and $\phi$ is the standard multivariate normal density function with a $D\times D$ correlation matrix ${\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}$ . Because directly maximizing the observed likelihood is difficult, a common solution is the EM algorithm (Bock & Aitkin, Reference Bock and Aitkin1981) that works on the complete data log likelihood instead,

(38) $$\begin{align}l\left({\mathbf{a}}^{\#},{\mathbf{b}}^{\#},{\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}|\mathbf{X},\boldsymbol{\Theta} \right)=\sum \limits_{i=1}^I\left\{\left[\sum \limits_{j=1}^J\sum \limits_{c=1}^C\log {P_{ij c}}^{\mathrm{I}\left({x}_{ij}=c\right)}\right]+\log \phi \left({\boldsymbol{\unicode{x3b8}}}_i|{\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}\right)\right\},\end{align}$$

where $\boldsymbol{\Theta}$ is the $I\times D$ latent trait score matrix for all test takers, ${\boldsymbol{\unicode{x3b8}}}_i$ is test taker $i$ ’s latent trait score vector including $D$ traits, and $\left(\mathbf{X},\boldsymbol{\Theta} \right)$ are the complete data with observed data $\mathbf{X}$ and missing data $\boldsymbol{\Theta}$ . The EM algorithm involves an iterative process that repeatedly executes two steps: an E step and an M step.

In the E step, the expectation of the complete data log likelihood with respect to the posterior distribution of missing data (i.e., trait scores) is estimated and serves as the EM map. The expectation process involves approximating the posterior distribution of trait scores. For a small number of traits (e.g., fewer than four), the posterior distribution can be approximated by quadrature points, such as the Gauss–Hermite quadrature points (Davis & Rabinowitz, Reference Davis and Rabinowitz1984), the adaptive variants (Schilling & Bock, Reference Schilling and Bock2005), or the equidistant points within a range (e.g., −6 to 6). However, a typical FCQ includes more than three traits. The quadrature point method is not practical for high-dimensional traits due to the enormous computing memory it demands. The alternative is the stochastic approach, for example, the Metropolis-Hastings Robbins-Monro (MHRM) algorithm (Cai, Reference Cai2010), which uses MCMC to draw random samples from the posterior distribution of trait scores.

In the M step, the EM map is maximized with respect to item parameters and distributional parameters of traits using the Newton–Raphson (NR) method (Atkinson, Reference Atkinson1989) or a quasi-Newton method, for example, the BFGS algorithm (Broyden, Reference Broyden1970). The MHRM algorithm employs a method similar to the NR method to search for model parameters that maximize the EM map using the drawn samples at the E step.

The MML-EM method is a standard estimation method widely applied to IRT and many other models. The R package mirt (Chalmers, Reference Chalmers2012) implemented a general-purpose MML-EM method for IRT models, which can be used to estimate existing IRT models and users’ newly proposed IRT models. For conciseness, we do not go into the details of the MML-EM method here. Interested readers are referred to Fu et al. (Reference Fu, Tan and Kyllonen2025), which provides the technical details of the MML-EM method for estimating the Triplet-2PLM.

5.2 MCMC estimation

The MCMC is a sampling-based Bayesian estimation method that draws model parameters from their posterior distribution(s). The joint posterior distribution of all model parameters conditional on item responses in the MGPPM under the Bayesian framework is the product of the prior distributions of all parameters and the likelihood function,

(39) $$\begin{align}P\left({\mathbf{a}}^{\#},{\mathbf{b}}^{\#},\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}|\mathbf{X}\right)\propto \prod \limits_{i=1}^I\prod \limits_{j=1}^J\prod \limits_{c=1}^C{P_{ij c}}^{\mathrm{I}\left({x}_{ij}=c\right)}\phi \left({\boldsymbol{\unicode{x3b8}}}_i|{\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}\right)P\left({a}_{jk}^{\#}\right)P\left({b}_{jk}^{\#}\right)P\left({\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}\right),\end{align}$$

where ${a}_{jk}^{\#}$ is the slope of statement $k$ in item $j$ ; ${b}_{jk}^{\#}$ is the intercept of statement $k$ in item $j$ ; $P\left({a}_{jk}^{\#}\right)$ , $P\left({b}_{jk}^{\#}\right)$ , and $P\left({\boldsymbol{\Sigma}}_{\boldsymbol{\unicode{x3b8}}}\right)$ are the prior distributions of the slope, intercept, and intertrait correlations, respectively. A sampler, such as the No-U-Turn sampler (Hoffman & Gelman, Reference Hoffman and Gelman2014) implemented in the R package rstan (Stan Development Team, 2024), can be used to draw samples from the joint posterior distribution (Equation 39).

6.1 Simulation design

We conducted a simulation study to inspect the parameter recovery of the MGPPM on triplets and tetrad. We also compared the MGPPM’s model fits and trait score estimations to those of the TIRT and Triplet-2PLM on the triplet data simulated based on the MGPPM. For each simulated dataset, we had 1000 test takers and five traits, with each trait measured by 12 statements. This resulted in 20 triplets and 15 tetrads. We used two repeated combinations of three out of five to define the trait measured by each statement in the 20 triplets; for the 15 tetrads, we used three repeated combinations of four out of five traits. Trait scores were drawn from a multivariate normal distribution with a mean vector of 0 and an intertrait correlation matrix that was an excerpt of five traits (Traits 1, 6, 8, 10, and 12) from the real data, which included 14 traits (see Table 3). The intercepts were drawn from a normal distribution with a mean of −0.05 and a standard deviation (SD) of 0.54. The distributional parameters were taken from our real data. The intercept of the first statement in each triplet is set to 0 for model identification.

Table 3 Intertrait correlation matrix from real data with 14 traits

In this simulation study, we controlled two factors. The first one was the mean of the slopes. The slopes were drawn from a lognormal distribution with an SD of 0.29 and three levels of means on the normal scale: low (0.64), medium (1.28), and high (2). The SD and low mean were also taken from our real data. Previous studies (e.g., Bürkner et al., Reference Bürkner, Schulte and Holling2019; Morillo et al., Reference Morillo, Leenen, Abad, Hontangas, De la Torre and Ponsoda2016; Lee et al., Reference Lee, Joo, Stark and Chernyshenko2019) have demonstrated that slopes significantly impact the estimations of IRT models for FC items, with high slopes associated with more accurate parameter estimates. The second one was the keyed direction of statements within an item. We considered two levels. In the “same” level, all statements within a triplet/tetrad had the same direction: the first ten triplets and eight tetrads had positive statements, and the second ten triplets and seven tetrads had negative statements by multiplying −1 with the generated slopes. In the “mixed” level, for triplets, one statement was changed to a different keyed direction from the other two in half of the 20 items in the same level. The change of keyed direction on statements was balanced in terms of the five traits and positive/negative items so that there were five positive and five negative items having the direction of one statement in each item changed, and for each trait, only one statement’s direction was changed within each of the five-item groups. For tetrads, the last two statements in 7 of the 15 items at the same level (items 1–3 and 12–15) were changed to a different keyed direction from the first two statements. This design of keyed direction changes also considered balancing statement directions in each trait. Previous studies (e.g., Brown & Maydeu-Olivares, Reference Brown and Maydeu-Olivares2011; Bürkner et al., Reference Bürkner, Schulte and Holling2019; Bürkner, Reference Bürkner2022) have demonstrated that item direction is a crucial factor influencing item information, and including mixed-direction items in an FCQ enhances test reliability. In this design, we considered the proportion of mixed-direction items in a test at 0.5 and 7/15 levels for triplets and tetrads, respectively, which are probably too high in real testing. However, we intended to use a relatively high proportion level to facilitate checking the general impact of mixed-direction items on model estimations.

In sum, the simulation study has six crossed conditions (3 slope levels $\times$ 2 keyed directions) for each item format. We generated 35 replicate datasets in each condition. The item and person parameters are redrawn for each replicate dataset. Each triplet dataset was estimated using MML-EM and MCMC methods, while each tetrad dataset was estimated only using the MML-EM method. Each Triplet dataset was also estimated by TIRT with the ULSMV method and by Triplet-2PLM with the MML-EM method, and their model fits and trait score estimations were compared to those of MGPPM with the MML-EM method.

6.2 Estimation programs

The majority of the analyses was conducted in the R program (R Core Team, 2023). For the MML-EM estimations of MGPPM and Triplet-2PLM, we used the mirt package in R. Users can use the createItem function in mirt to define new item response functions. Then, the estimation and analysis modules in mirt can be used to estimate and further analyze the functions. The first and second derivatives of the MGPPM’s (Equations 25–31) and Triplet-2PLM’s (Equation 9) response functions for triplets were obtained using the symbolic derivative function, Deriv, in the R package Deriv (Clausen & Sokol, Reference Clausen and Sokol2020). We used the MHRM algorithm implemented in mirt for the MML-EM estimations. For the MHRM algorithm, mirt estimates the MHRM standard errors (SE) of model parameter estimates; however, it carries out the SE estimation separately after the parameter estimation is done, and thus, it takes much longer to obtain SE estimates. To save running time, we modified the functions in mirt so that parameters and their SE estimations were performed simultaneously. The starting values of the discrimination parameters were set to 1 for positive statements and −1 for negative statements. All the starting values of the intercepts and intertrait correlations were set to 0 and 0.25, respectively. All trait score estimates were Maximum A Posterior (MAP) estimates. The other default settings in mirt, including the convergence criterion, were used for all estimations.

The TIRT model was estimated using the 64-bit Mplus 8.4 software (Reference Muthén and MuthénMuthén & Muthén, 1998–2017) under structural equation modeling with the ULSMV method and default settings. We fixed the error variance of the first statement ( ${\psi}_1^2$ ) in a triplet to 1 to identify the model. We modified the functions in the R package thurstonianIRT (Bürkner, Reference Bürkner2024) and used them to set up and run the TIRT models in Mplus from R. Once the item parameter estimates were obtained from Mplus, we set up an M2PLM in mirt and fixed the item parameters to be those estimates from Mplus (after parameter transformation). Then, the many functions in mirt could be used for further analyses. Note that in the estimation of the TIRT model in Mplus, the covariances among the three pairs stemming from a triplet were considered only in the item parameter estimation and not in the trait score estimation (i.e., the local independence assumed in Equations 14–16). Thus, the trait score estimates from the above setup of mirt are theoretically the same as those from Mplus because their item response functions and response data are the same. Frick (Reference Frick2023) and Maydeu-Olivares and Brown (Reference Maydeu-Olivares and Brown2010) demonstrated that simplifying scoring had a minimal impact on the accuracy of latent score estimates.

The MCMC estimations of triplets were conducted using the rstan package (Stan Development Team, 2024). We applied weakly informative priors to the slopes and intercepts. For the slope prior, we used the lognormal distribution with a mean of 0 and an SD of 2 on the log scale. We assumed that the signs of slopes were known and fixed them during the estimations. For the intercept prior, we used a normal distribution with a mean of 0 and an SD of 3. For the intertrait correlation prior, we used the Lewandowski–Kurowicka–Joe (LKJ) distribution with a shape parameter of 4. Each estimation ran four chains, each with 8,000 iterations, with the first 4,000 iterations used as the burn-in period. For each parameter, the mean of the 16,000 final samples (4,000 samples per group) served as the estimate, and the SE of the estimate could also be estimated from the final samples. An estimation was deemed converged if the $\widehat{R}$ statistic (Gelman, Reference Gelman, Gilks, Richardson and Spiegelhalter1996, p. 137) was smaller than 1.01, and the bulk effective sample size was larger than 400 for every estimated parameter.

The supplemental online materials (SOM) provide the R code to simulate triplet and tetrad item responses based on the MGPPM and estimate the models using the MML-EM method for both item formats and the MCMC method for triplets. In addition, the following R code is also provided in SOM: (a) estimating the Triplet-2PLM using the MML-EM method, (b) setting up and running the TIRT model in Mplus from R, and (c) setting up the TIRT model and estimating trait scores in mirt with fixed item parameters from Mplus’s estimates (after transformation).

6.3 Results: MGPPM’s parameter recovery

We checked the recovery of slopes, intercepts, intertrait correlations, and trait score estimates in the MGPPM. The following criteria were used: mean deviation (MD), mean absolute deviation (MAD), root mean square error (RMSE), mean relative deviation (MRD), and correlation (Cor) with true values.

6.3.1 Triplets

Tables 4 and 5 show the parameter recovery results for the MML-EM and MCMC estimations of triplets, respectively. Only converged runs were included in the parameter recovery results. The converged rates were included in Tables 4 and 5. For the MML-EM method, only two runs in the same direction and one run in the mixed direction under the low-slope condition were not converged due to poor intertrait correlation estimates that caused the estimations to crash. There were more nonconvergent runs in the MCMC estimations (probably due to the different convergence rule from the MML-EM method): four runs in the same direction in the low- and high-slope conditions, nine runs in the same-direction-medium-slope condition, and two runs in the mixed-direction-high-slope condition. For the MCMC runs, the same-direction runs were more difficult to converge than the mixed-direction runs.

Table 4 Triplet MML-EM: Parameter recovery results

Note: MAD = mean absolute deviation; Cor = correlation; MD = mean deviation; MRD = mean relative deviation; RMSE = root mean square error.

Table 5 Triplet MCMC: Parameter recovery results

Note: MAD = mean absolute deviation; Cor = correlation; MD = mean deviation; MRD = mean relative deviation; RMSE = root mean square error.

Tables 4 and 5 show that the parameter estimates from the MML-EM and MCMC methods were very similar. Some general patterns in the simulation results from both estimation methods are observed as follows.

1. 1. Across all conditions, the slopes and intercepts recovered well, as the correlations with the true values were larger than 0.97, MDs and MRDs were close to 0, and MAD and RMSE were small. The slopes in the same-direction-high-slope condition performed slightly worse, as their MADs and RMSEs were the largest (in the range from 0.19 to 0.28).

2. 2. Across all conditions, the intertrait correlation estimates showed correlations with true values of greater than 0.94. However, the other four criteria in the same-direction conditions were worse than those in the mixed-direction conditions. Thus, the mixed-direction conditions had better intertrait correlation estimates than the same-direction conditions. The MCMC estimations show that the mean biases of the intertrait correlation estimates were all negative in all conditions; however, in the same-direction conditions, the mean biases were around −0.1, while in the mixed-direction conditions, they were around −0.01, indicating that trait score estimates in the mixed-direction conditions were better than the same-direction conditions. The mean biases of the intertrait correlation estimates from the MML-EM estimations were similar in magnitude to those of the MCMC; however, they had a positive direction, except in the same-direction conditions with medium or high slopes. In the low-slope-same-direction condition, the MML-EM method tended to overestimate intertrait correlations, while the MCMC method tended to underestimate them.

3. 3. For the trait score estimates, the mixed-direction conditions showed moderate improvement over the same-direction conditions, with the percentage of increase in average correlations with true values ranging from 12% to 18%. With the mixed direction, the medium-slope condition increased the average correlation in the low-slope condition from 0.83 (or 0.84 for the MCMC method) to 0.92, while the increase from the medium-slope condition to the high-slope condition was only 0.02.

The MCMC estimations took much longer than the MML-EM ones. Because our simulation runs were conducted on different computers with varying hardware configurations and operating systems, we were unable to provide an overall summary of running times for both methods. However, we ran both MML-EM and MCMC estimations on 10 datasets in the medium-slope condition on a single computer with an Intel i7-4770U CPU, 16 GB of RAM, and Windows 10 OS. Based on the 10 datasets, the MML-EM estimations took an average of 37 and 38 min for the same- and mixed-direction conditions, respectively. In contrast, the average times for the MCMC estimations were 584 and 258 min, respectively.Footnote ²

6.3.2 Tetrads

Table 6 lists the parameter recovery results for the MML-EM estimations of tetrads, including the convergence rate and running time in each condition. Like triplets, only three runs in the same direction and one in the mixed direction under the low-slope condition were not converged. The parameter recovery for tetrads was similar or slightly better than that for triplets and showed similar patterns across conditions; moderate improvements were observed in latent trait score estimates in the low-slope conditions or the same-direction conditions, where the correlations with true values for tetrads increased, on average, from 0.04 to 0.07 compared to those for triplets. The average running time on a server with Intel Xeon E5-2640 v3 CPU, 256 GB of RAM, and Microsoft Windows Server 2016 Standard OS based on 21–23 converged runsFootnote ³ was around 180 min, much longer than the average running times for triplets reported previously.

Table 6 Tetrad MML-EM: Parameter recovery results

Note: MAD = mean absolute deviation; Cor = correlation; MD = mean deviation; MRD = mean relative deviation; RMSE = root mean square error.

^a Mean running time is based on 23 runs except for the low-slope-same-direction (21 runs) and low-slope-mixed-direction (22 runs) conditions on a server with Intel Xeon E5-2640 v3 CPU, 256 GB RAM, and Microsoft Windows Server 2016 Standard OS.

6.4 Results: Model comparison

The simulated Triplet data based on MGPPM were also estimated by the TIRT using the ULSMV method and by Triplet-2PLM using the MML-EM method. Their model fits and trait score estimations were compared to those of the MGPPM using the MML-EM method. All model and trait score estimations were converged in the TIRT and Triplet-2PLM.

Table 7 presents the summary statistics of the differences in Akaike information criterion (AIC) between Triplet-2PLM and MGPPM across converged replicates in both models for each of the six simulated conditions. Note that AIC differences are equal to the differences of Bayesian information criteria (BIC), and thus, Table 7 only shows the AIC differences. For AIC and BIC, smaller is better. A difference of at least 2 for AIC or 6 for BIC is considered strong evidence in favor of a better model (Bauldry, Reference Bauldry and James2015). Table 7 shows that the MGPPM’s AIC and BIC were smaller than the Triplet-2PLM by at least 92 across all replicates, indicating that AIC and BIC prefer MGPPM over Triplet-2PLM. Note that Table 7 does not include the TIRT because the TIRT used a different dataset (i.e., recoded pair scores from triplet scores) from the other models (i.e., triplet scores), so their model comparison statistics, such as AIC and BIC, were incomparable.

Table 7 Comparisons of Triplet-2PLM and MGPPM on simulated triplet data: AIC difference (Triplet-2PLM–MGPPM)

Note: Triplet-2PLM and MGPPM were estimated by MML-EM. AIC difference is the same as BIC difference.

Table 8 presents the averages of the M2* Chi-square statistics (Cai & Hansen, Reference Cai and Hansen2013) and the related M2*-based RMSEA (root mean squared error of approximation) and SRMSR (standardized root mean squared residual) statistics for the three models across converged replicates under each condition. A p-value of M2* larger than a significant level (e.g., 0.01) indicates an overall good model fit. While M2* provides a significant test for omnibus model fit, RMSEA can be used to assess the goodness of approximation between the target model and the true model, and SRMSR indicates the magnitude of the misfit (the effect size). Maydeu-Olivares and Joe (Reference Maydeu-Olivares and Joe2014; Table 2) suggested an RMSEA and SRMSR smaller than 0.05/ $\left(m-1\right)$ and 0.027/ $\left(m-1\right)$ , respectively, for an excellent fit, smaller than 0.05 and 0.027, respectively, for a close fit, and smaller than 0.0089 and 0.05, respectively, for an adequate fit, where $m$ is the number of an item’s response categories (i.e., six for triplets). Table 8 shows that in all conditions, TIRT’s average M2* values were much larger than Triplet-2PLM’s and MGPPM’s, and TIRT’s mean p-values of the M2* tests were significant at 0.00, while the other two models’ were nonsignificant. Furthermore, TIRT’s mean RMSEA and SRMSR values were much larger than Triplet-2PLM’s and MGPPM’s, while Triplet-2PLM and MGPPM had similar values. TIRT’s mean RMSEA and SRMSR values decreased with higher item slopes, while the other two models’ values changed little across conditions. In the high-slope condition, TIRT’s mean SRMSR was just slightly higher than the other two models’. In general, the RMSEA values indicated close fits for TIRT and excellent fits for Triplet-2PLM and MGPPM; the SRMSR values indicated adequate fits for TIRT and close fits for Triplet-2PLM and MGPPM.

Table 8 Comparison of TIRT, Triplet-2PLM, and MGPPM on simulated triplet data: M2*, RMSEA, and SRMSR

Note: TIRT was estimated by ULSMV, and Triplet-2PLM was estimated by MML-EM. RMSEA = root mean squared error of approximation; SRMSR = standardized root mean squared residual.

Table 9 Comparisons of TIRT, Triplet-2PLM, and MGPPM on simulated triplet data: Reliability and correlation of latent score estimates

Note: TIRT was estimated by ULSMV, and Triplet-2PLM and MGPPM were estimated by MML-EM. Cor_MGPPM = Correlation of latent score estimates with MGPPM’s; Rel_ture = Squared correlation between true and estimated latent scores (true reliability); Rel_est = Estimated reliability of latent score estimates (based on Equation 40).

The recovery of trait score estimations for TIRT and Triplet-2PLM was very similar to MGPPM’s. The evaluation results for TIRT and Triplet-2PLM, based on the five criteria in Table 4, are presented in Table A1 in the SOM. Table 9 shows the average true reliabilities, estimated reliabilities, and correlations with MGPPM’s trait score estimates for all three estimation methods and five traits across the converged runs in each condition. True reliability is defined as the squared correlation between true and estimated trait scores. The empirical reliability of a trait’s MAP score estimates is estimated as follows (Kim, Reference Kim2012):

(40) $$\begin{align}{\widehat{\rho}}^2=\frac{\operatorname{var}\left(\widehat{\theta}\right)}{\operatorname{var}\left(\widehat{\theta}\right)+\sum \limits_{i=1}^I\operatorname{var}\left({\widehat{\theta}}_i|{\theta}_i\right)/I},\end{align}$$

where $\operatorname{var}\left(\widehat{\theta}\right)$ is the sample variance of trait score estimates, $\operatorname{var}\left({\widehat{\theta}}_i|{\theta}_i\right)$ is the error variance of subject $i$ ’s estimated trait score ${\widehat{\theta}}_i$ , and $I$ is the sample size. In Table 9, the average correlations of TIRT’s and Triplet-2PLM’s trait score estimates with MGPPM’s were at least 0.98, and their trait score reliabilities were very similar in each condition. One exception is that in the low-slope condition, TIRT’s estimated reliabilities were larger than those of the other two models and inflated compared to the true reliabilities. The inflation of estimated reliabilities was also observed in the medium-slope-same-direction condition for TIRT and Triplet-2PLM, as well as in the high-slope-same-direction condition for all three methods. Theoretically, for TIRT, due to the local independence assumption in trait score estimation, redundant information is added to the test information, resulting in an overestimated reliability (Frick, Reference Frick2023). However, our result showed that reliability overestimation occurred for TIRT only when the true reliability was low (e.g., smaller than 0.8), and other factors also impacted reliability estimation. In addition, under certain conditions, Triplet-2PLM and MGPPM also overestimate score reliabilities. Future research is needed to further investigate this issue.

In terms of estimation times, a TIRT estimation in Mplus, without estimating standard errors, completed in 1–2 s on average on a laptop with an Intel Core i5-1345U CPU and 16GB RAM, while a typical Triplet-2PLM or MGPPM run in R took from about 39 to 53 min on a Linux server. In general, an MGPPM run was slightly faster than a Triplet-2PLM one; however, MGPPM’s running times were more divergent than Triplet-2PLM’s. Table A2 in the SOM summarizes the running items of the three models in each condition.

In summary, the simulation results showed that when estimating simulated data generated from MGPPM, although the fit statistics might prefer the generating model over the TIRT and Triplet-2PLM, the trait score estimates from the three models were practically indifferent.

7.1 Results

Table 10 shows the model comparison result between the Triplet-2PLM and MGPPM based on the log likelihood, Akaike information criterion (AIC), and Bayesian information criteria (BIC). In Table 10, the MGPPM’s AIC and BIC were 253 smaller than the Triplet-2PLM’s, indicating a better fit.

Table 10 Comparison of log likelihood, AIC, and BIC

^a Model above – Model below.

Table 11 presents the M2*, RMSEA, and SRMSR statistics for the three models. The M2* statistics were significant for all models, and the SRMSR statistics were nearly identical across all models, exceeding the threshold for an adequate fit (0.5). In contrast, MGPPM’s and Triplet-2PLM’s RMSEA statistics were slightly smaller than TIRT’s, and they all were smaller than 0.05, indicating all models had close but perfect fits to the data.

Table 11 Comparison of M2* and RMSEA

Note: RMSEA = root mean squared error of approximation; SRMSR = standardized root mean squared residual.

We use the ${z}_h$ statistic (Drasgow et al., Reference Drasgow, Levine and Williams1985) to assess the fit of an individual item. Note that items in the TIRT model refer to the pairs stemming from the triplets. The asymptotic distribution of the ${z}_h$ statistic is standard normal, assuming true model parameters. The hypothesis test is conducted structurally as a lower one-sided test. Because the estimated model parameters from these IRT models were used in estimating ${z}_h$ , the asymptotic distributions of the ${z}_h$ statistics were unknown. However, we used ${z}_h<-1.64$ for the Triplet-2PLM and MGPPM and ${z}_h<-2.13$ for the TIRT model (due to the Bonferroni adjustment) as the rough criteria to flag potential misfit items corresponding to a one-sided significance level of 0.05 based on the standard normal distribution. Table 12 shows only one misfitted pair (stemming from a triplet) in the TIRT model and none in the other two models.

Table 12 Comparison of number (N) and percentage (%) of misfit items

^a Flagging criterion ${z}_h$ < −1.64 for the Triple-2PLM and MGPPM and ${z}_h$ < −2.13 for the TIRT model.

^b Flagging criterion Cramer’s V > 0.3.

An item pair fit test has been used to check the local independence assumption in IRT models. We used Chen and Thissen’s (Reference Chen and Thissen1997) Pearson ${X}^2$ local dependence ( $\mathrm{LD}\hbox{-} {\mathrm{X}}^2$ ) index to check the fits of item pairs. Because the asymptotic distribution of $\mathrm{LD}\hbox{-} {\mathrm{X}}^2$ is unknown, we transformed the $\mathrm{LD}\hbox{-} {\mathrm{X}}^2$ index to Cramer’s V index (Cramer, Reference Cramer1946, p. 282), a normalized version of the ${X}^2$ test statistic, and used the criterion, Cramer’s V larger than 0.3, as a rough indicator of strong local dependence among items. The justification for applying the $\mathrm{LD}\hbox{-} {\mathrm{X}}^2$ test to the TIRT models lies in the local independence assumption among the three pairs that stem from a triplet in trait score estimation. In Table 12, a tiny (0.36%) proportion of item pairs (i.e., 58 out of 16,110 pairs of items) in the TIRT model showed strong local dependence, while none did for the other two models.

Table 13 shows that the Triplet-2PLM and MGPPM had almost identical reliability estimates across the 14 traits, and both had the mean and SD of the reliability estimates of 0.68 and 0.07, respectively. In contrast, the TIRT model had slightly and consistently higher reliability estimates than the other two models across the 14 traits; the mean and SD of the reliability estimates were 0.74 and 0.06, respectively, for the TIRT model. The reliability estimates in the TIRT model were likely inflated due to the local independence assumption among the pairs stemming from the triplets in trait score estimations. This result is consistent with that in the low-slope condition of the previous simulation study, as the slope’s distribution parameters in the low-slope condition were derived from real data.

Table 13 Comparison of trait score reliability estimates

Table 14 lists the correlations of trait score estimates between the MGPPM and the other two models. The Triplet-2PLM showed nearly perfect correlations with the MGPPM, with a mean correlation of 0.99 and a standard deviation of 0.00. The TIRT’s correlations were slightly lower, with a mean and standard deviation of 0.98 and 0.01, respectively. This result is also consistent with that in the simulation study.

Table 14 Correlations of trait score estimates with MGPPM

In summary, all three models provided a good fit to the real data, as indicated by the model and item fit statistics. Their trait score estimates were almost perfectly correlated. The Triplet-2PLM and MGPPM had almost identical reliability estimates on trait scores, while the TIRT model’s reliability estimates were slightly inflated. Based on the model comparison statistics, the MGPPM fitted the data better than the Triplet-2PLM.