Extending an Identified Four-Parameter IRT Model: The Confirmatory Set-4PNO Model

Abstract

Given the frequent presence of slipping and guessing in item responses, models for the inclusion of their effects are highly important. Unfortunately, the most common model for their inclusion, the four-parameter item response theory model, potentially has severe deficiencies related to its possible unidentifiability. With this issue in mind, the dyad four-parameter normal ogive (Dyad-4PNO) model was developed. This model allows for slipping and guessing effects by including binary augmented variables—each indicated by two items whose probabilities are determined by slipping and guessing parameters—which are subsequently related to a continuous latent trait through a two-parameter model. Furthermore, the Dyad-4PNO assumes uncertainty as to which items are paired on each augmented variable. In this way, the model is inherently exploratory. In the current article, the new model, called the Set-4PNO model, is an extension of the Dyad-4PNO in two ways. First, the new model allows for more than two items per augmented variable. Second, these item sets are assumed to be fixed, that is, the model is confirmatory. This article discusses this extension and introduces a Gibbs sampling algorithm to estimate the model. A Monte Carlo simulation study shows the efficacy of the algorithm at estimating the model parameters. A real data example shows that this extension may be viable in practice, with the data fitting a more general Set-4PNO model (i.e., more than two items per augmented variable) better than the Dyad-4PNO, 2PNO, 3PNO, and 4PNO models.

Keywords

four-parameter model Bayesian statistics slipping guessing identification higher-order DINA model Dyad-4PNO Set-4PNO

1. Introduction

Despite having been first introduced four decades ago, the four-parameter model (4PM) has seen a relative resurgence in interest within item response theory. The 4PM is a generalization of the three-parameter model (3PM) that allows for both lower and upper asymptotes on the probability of a positive item response on a dichotomous item. Its item response function (IRF) is as follows:

P (Y_{j} = 1 | θ, α_{j}, β_{j}, g_{j}, s_{j}) = g_{j} + (1 - s_{j} - g_{j}) Ψ (α_{j} θ - β_{j}),

where $Ψ (\cdot)$ denotes either a logit link function, resulting in the four-parameter logistic (4PL) model (Barton & Lord, 1981), or a probit link function, resulting in the four-parameter normal ogive model (4PNO; Culpepper, 2016). In this formulation, $θ$ is the examinee latent trait level; $β_{j}$ and $α_{j}$ are the jth item’s threshold and slope parameters, respectively; and g_j and $1 - s_{j}$ are the jth item’s lower and upper asymptote parameters, respectively. The parameters g_j and s_j are often referred to as the guessing and slipping parameters.

The renewed interest has largely been fueled by emerging evidence that upper asymptotes may exist in a variety of contexts, such as low-stakes assessments (Culpepper, 2017), high-stakes assessments (Rulison & Loken, 2009), psychological and personality assessments (Waller & Reise, 2010), and even in the assessment and monitoring of food insecurity (Gregory, 2020). This research shows that the reasons for the existence of upper asymptotes can range from simply making a mistake during test-taking to careless responding to intentionally misresponses.

Further interest in the 4PM centers on the problems underlying its estimation. Estimation approaches that have been introduced in recent years include estimation using WinBUGS (Loken & Rulison, 2010), marginal maximum likelihood estimation (Feuerstahler & Waller, 2014), using a Gibbs sampling method to estimate the 4PNO (Culpepper, 2016), Bayesian model estimation (Waller & Feuerstahler, 2017), marginal maximum a posterior estimation using a mixture model formulation (Meng et al., 2020), using a Gibbs-slice sampler (Zhang et al., 2020), and using data augmentation in a Gibbs sampler to estimate the multidimensional 4PL (Fu et al., 2021).

Much of the interest in estimation of the 4PM stems from the serious concerns underlying its identification that are extended to it due to its generalization from the 3PM (Maris & Bechger, 2009; San Martín et al., 2015; Thissen, 2009). When a model has identification issues, several severe consequences arise. The most salient comes directly from the definition of an identified model: having an identified model means that optimal parameters are unique, so that when a model is unidentified, several parameter sets may be equally optimal. In this situation, how does one choose between them? Other problems include the consistency of parameter estimates (Gabrielsen, 1978) and the interpretability of the parameters (San Martín et al., 2009). These issues should cause any applied researcher relying on a potentially unidentified model to take pause. On the other hand, evidence suggests that not taking into account guessing effects (Aitkin & Aitkin, 2006) and slipping effects (Loken & Rulison, 2010) may lead to biased trait estimates in the extreme ends of the trait distribution.

Given the previous paragraph, it is clear that including these effects in such a way that the model is identified is a highly valuable prospect for the psychometric community. To that end, one solution proposed recently showed conditions under which a 4PM model can be identified (Kern & Culpepper, 2020). Using results from the cognitive diagnostic modeling literature (Gu & Xu, 2021), it was shown that by grouping items into pairs that measure a common attribute—which they termed as dyads—and restricting the slope and threshold parameters within that pair to be equal, that the model can be generically identified. This is as follows:

P_{j} (θ_{i}) = g_{j} + (1 - s_{j} - g_{j}) Φ (a_{d} θ_{i} - b_{d}),

where d denotes the dth dyad, and $Φ$ refers to the standard normal cumulative distribution function. Moreover, they showed that by grouping sets of at least three items that measure a common binary attribute together and restricting slopes and threshold parameter to be equal within each set that the model can be strictly identified.

In their paper, they proposed a Markov chain Monte Carlo (MCMC) algorithm to explore the space of all possible dyads to determine a likely set of dyads and estimate person, dyad, and item parameters when the dyads are not known a priori. In their approach, they used a prior distribution of potential groupings that assumed that each potential dyad was equally likely. Hence, their approach yielded an exploratory version of the proposed restricted 4PM. In the current project, we will modify the procedure to instead use a prior distribution of potential groupings that assumes one specific grouping, yielding a confirmatory version of the proposed restricted 4PM. Furthermore, the new model considered here will extend the Dyad-4PNO model further by allowing for groups of items, which we term item sets, with two or more items, rather than just the two-item dyads previously considered. These extensions will make the approach more robust, allowing for testable hypotheses with respect to item sets and improving the estimation accuracy of model parameters. Importantly, the resulting model can be thought of as a special case of the higher-order DINA (HO-DINA) model (de la Torre & Douglas, 2004) with simple structure. By exploiting this simple structure, the developed algorithm is able to accommodate many more attributes than other existing HO-DINA estimation algorithms.

The rest of the current article is structured into five sections. In the first section, the Dyad-4PNO model is introduced. Following this brief introduction is a discussion about the extension of the Dyad-4PNO model, which we term the set four-parameter normal ogive (Set-4PNO) model, to cases where a test may consist of item sets with two or more items. The second section describes a Bayesian formulation of the new model, as well as an outlined Gibbs sampling algorithm for approximating the posterior distribution. Also given in this section is a brief discussion of how the current formulation connects with the exploratory Dyad-4PNO. The third section describes and reports the results from a Monte Carlo simulation study designed to evaluate parameter recovery of the proposed Gibbs sampler under a variety of conditions. In the fourth section, we give an example of fitting several competing models, including the Set-4PNO, to an experimental IQ data set. Finally, we provide a discussion on the contributions of this article, potential future directions of this work, and limitations of the current study.

2. A Restricted 4PM

2.1. The Dyad-4PNO Model

In this section, we describe the development of the Dyad-4PNO model. As originally described, this model can be thought of as a HO-DINA model (de la Torre & Douglas, 2004) with uncertainty in which items measure a binary latent attribute, and a restriction that attributes are measured by only two items, the minimum number of items necessary for generic identifiability to hold (Kern & Culpepper, 2020). These sets of paired items are called dyads. To denote a dyad, we let $τ = (τ_{1}, \dots, τ_{J})^{⊤}$ indicate the assignment of items to pairs, such that items j and $j^{'}$ are grouped into a dyad d if $τ_{j} = τ_{j^{'}} = d$ , where $d = 1, \dots, D$ . We call $τ$ the vector of dyad indices. Each dyad d can be thought of as a group of items that measure a common binary attribute. Thus, we see that $τ$ conveys the same information as a cognitive diagnostic model’s (CDM) $Q$ matrix, where every item has simple structure. That is, every item measures one, and only one, binary attribute $W_{i d}$ .

Suppose that person i does not possess attribute d. Then, under the DINA model $P (Y_{i j} = 1 | W_{i d} = 0) = g_{j}$ , so that g_j can be interpreted as the probability of answering an item correctly, given that the examinee does not have the required attribute. Similarly, under the DINA model, $P (Y_{i j} = 0 | W_{i d} = 1) = s_{j}$ . That is, s_j is the probability of answering an item incorrectly given the examinee does have the required attribute. In this way, g_j and s_j are commonly called guessing the slipping parameters, respectively. Thus, we find that

Y_{i j} | W_{i, τ_{j}} = w_{i d}, s_{j}, g_{j} \sim Bernoulli ({(1 - s_{j})}^{w_{i d}} g_{j}^{1 - w_{i d}}) .

In the current model, we assume that if items j and $j^{'}$ are a dyad pair d, then responses to these items are locally independent given $W_{i d}$ . This implies that

P (Y_{i d} = y_{i d} | τ, W_{i d}, s_{d}, g_{d}) = \prod_{k \in {j, j^{'}}} {({(1 - s_{k})}^{w_{i d}} g_{k}^{1 - w_{i d}})}^{y_{i k}} {(s_{k}^{w_{i d}} {(1 - g_{k})}^{1 - w_{i d}})}^{1 - y_{i k}},

where $Y_{i d}$ denotes the binary item responses for items $(Y_{i j}, Y_{i j^{'}})$ in the dth dyad with realization $y_{i d} = (y_{i j}, y_{i j^{'}})$ , and $s_{d}$ and $g_{d}$ contain the slipping and guessing parameters for the corresponding items in dyad d.

Treating $W_{i d}$ s as latent categorical responses, we might assume that they are controlled by a single latent trait $θ_{i}$ . A simple model for this is the two-parameter normal ogive (2PNO) model

Φ_{i d} \equiv P (W_{i d} = 1 | θ_{i}, a_{d}, b_{d}) = Φ (a_{d} θ_{i} - b_{d}),

where a_d is a slope parameter describing the strength of the relationship between $θ_{i}$ and $W_{i d}$ , and b_d is a threshold parameter. A discrimination-difficulty parameterization can be obtained by setting $a_{d}^{*} = a_{d}$ and $b_{d}^{*} = b_{d} / a_{d}$ so that $a_{d} θ_{i} - b_{d} = a_{d}^{*} (θ_{i} - b_{d}^{*})$ . Here, $a_{d}^{*}$ is a discrimination parameter describing how much the attribute—and, thus, the dyad—distinguishes between persons with low and high levels of $θ_{i}$ . Accordingly, $b_{d}^{*}$ is a difficulty parameter describing the location of the binary attribute $W_{i d}$ on the $θ$ -scale, wherein a person is equally likely to possess as not possess attribute d; those with higher values of $θ_{i}$ are more likely to possess $W_{i d}$ .

Here, we can define a so-called dyad response function (DRF). The DRF is the joint distribution of the item responses to items in dyad d. The DRF is

\begin{array}{l} P (Y_{i d} = y_{i d} | τ, a_{d}, & b_{d}, s_{d}, g_{d}, θ_{i}) = \sum_{w = 0}^{1} P (Y_{i d} = y_{i d} | τ, W_{i d} = w, s_{d}, g_{d}) P (W_{i d} = w | θ_{i}, a_{d}, b_{d}) \\ = [\prod_{k \in {j, j^{'}}} g_{k}^{y_{i k}} {(1 - g_{k})}^{1 - y_{i k}}] (1 - Φ_{i d}) + [\prod_{k \in {j, j^{'}}} {(1 - s_{k})}^{y_{i k}} s_{k}^{1 - y_{i k}}] Φ_{i d} \\ = g_{d} (y_{i d}) + [s_{d} (y_{i d}) - g_{d} (y_{i d})] Φ_{i d}, \end{array}

where

g_{d} (y_{i d}) = \prod_{k \in {j, j^{'}}} g_{k}^{y_{i k}} {(1 - g_{k})}^{1 - y_{i k}},

and

s_{d} (y_{i d}) = \prod_{k \in {j, j^{'}}} {(1 - s_{k})}^{y_{i k}} s_{k}^{1 - y_{i k}},

are the asymptotes as $θ \to - \infty$ and $θ \to + \infty$ , respectively, of the probability of a given response pattern $y_{i d}$ on dyad d.

The marginal IRF for a single item can be found by summing over the response categories of the other item in a dyad. Doing this shows that the resulting marginal item characteristic curve (ICC) is

P_{j} (θ_{i}) = g_{j} + (1 - s_{j} - g_{j}) Φ (a_{d} θ_{i} - b_{d}) .

Notice that this is similar to the ICC of the usual 4PNO model in Equation 1 with the difference that the item slope and intercept parameters are restricted to be equal to each other within each dyad pair d. The resulting model is generically identified.

Assuming local independence of the responses between dyads, the likelihood over all J item responses for the entire set of n examinees is as follows:

P (Y | τ, θ, s, g, a, b) = \prod_{i = 1}^{n} \prod_{d = 1}^{D} P (Y_{i d} = y_{i d} | τ, a_{d}, b_{d}, s_{d}, g_{d}, θ_{i}) .

Note that in general, this likelihood is not equivalent to the one often given in the context of the 4PM. This is because item responses within an item dyad are not assumed independent; that is, the dyad makes an explicit dependence assumption between items within the dyad. However, this dyadic relationship is necessary for the identifiability conditions to minimally hold.

2.2. The Confirmatory Set-4PNO Model

Treating the dyad index vector $τ$ in the Dyad-4PNO model as an unknown parameter to be estimated gives us the exploratory version of the model. This is how the model was originally described. Specifically, the vector $τ$ was explored over the space of all potential sets of item pairings. This approach was taken to satisfy the minimal requirements for generic identifiability, that is, that all items are grouped into pairs of item, or dyads, such that $τ_{j} = τ_{j^{'}} = d$ . This is a particularly powerful approach when the structure of the test is unknown.

However, here we observe two things. First, while a minimum of two items forming a dyad loading onto a common attribute $W_{i d}$ is sufficient for identifiability (Corollary 2, Kern & Culpepper, 2020), this does not preclude the possibility that item sets may contain three of more items. In fact, if each item set contains at least three items, then it was shown that the model is strictly identified (Corollary 1, Kern & Culpepper, 2020). Second, it is quite possible for the test structure to be known or theorized, either through expert evaluation or otherwise. This is important for being able to make testable hypotheses regarding the test structure.

Given these observations, the Dyad-4PNO model can be generalized. First, instead of only allowing for item sets of two items, we can instead allow for sets of size two or larger. To do this, the definition of an item set in the vector $τ$ can be expanded, so that $m \geq 2$ items $j_{1}, j_{2}, \dots, j_{m}$ are grouped together in item set d if $τ_{j_{1}} = τ_{j_{2}} = \dots = τ_{j_{m}} = d$ , where $d = 1, \dots, D$ . The items grouped together in item set d can be denoted by $D_{d}$ . The m items in item set d all load onto a common latent binary indicator $W_{i d}$ . Assuming the items $j_{1}, j_{2}, \dots, j_{m}$ are locally independent given $W_{i d}$ , then

P (Y_{i d} = y_{i d} | τ, W_{i d}, s_{d}, g_{d}) = \prod_{k \in D_{d}} {({(1 - s_{k})}^{w_{i d}} g_{k}^{1 - w_{i d}})}^{y_{i k}} {(s_{k}^{w_{i d}} {(1 - g_{k})}^{1 - w_{i d}})}^{1 - y_{i k}},

where $Y_{i d} = (Y_{i j_{1}}, \dots, Y_{i j_{m}})$ , $s_{d} = (s_{j_{1}}, \dots, s_{j_{m}})$ , and $g_{d} = (g_{j_{1}}, \dots, g_{j_{m}})$ . Furthermore, the corresponding DRF is expanded to be

P (Y_{i d} = y_{i d} | τ, a_{d}, b_{d}, s_{d}, g_{d}, θ_{i}) = g_{d} (y_{i d}) + [s_{d} (y_{i d}) - g_{d} (y_{i d})] Φ_{i d},

where

g_{d} (y_{i d}) = \prod_{k \in D_{d}} g_{k}^{y_{i k}} {(1 - g_{k})}^{1 - y_{i k}},

and

s_{d} (y_{i d}) = \prod_{k \in D_{d}} {(1 - s_{k})}^{y_{i k}} s_{k}^{1 - y_{i k}},

with interpretation as given previously. With this, we can find the expected sum-score S_d for an item set d to be

E (S_{i d} | θ_{i}) = \sum_{y_{i d} \in Y_{d}} (1^{⊤} y_{i d}) P (Y_{i d} = y_{i d} | τ, a_{d}, b_{d}, s_{d}, g_{d}, θ_{i}),

where $Y_{d}$ is the set of all possible response patterns $y_{d}$ for item set d.

Approaching the second observation that test structure may be known or theorized, we can simply treat the dyad index vector $τ$ —which is now better described as the item set index vector—as known and fixed. The only restriction that should be placed on $τ$ is that each item set d should have at least two items, so the model is identified. By fixing $τ$ according to the known test structure, we are no longer looking to explore the space of potential groupings but rather are seeking to examine the hypothesis that the test follows the given structure. It is in this sense that we call this the confirmatory Set-4PNO model. Figure 1 gives an example of the confirmatory Set-4PNO model implied by setting the index vector to $τ = (1,1,2,2,2,3,3)$ . If there are competing test structures $τ_{1}$ and $τ_{2}$ , we can simply fit the models implied by $τ_{1}$ and $τ_{2}$ and compare their fits.

Figure 1.

Path diagram of the confirmatory set four-parameter normal ogive (Set-4PNO) model implied by the item set index vector $τ = (1,1,2,2,2,3,3)$ .

2.3. Singlet Items

An important point missing from the above discussion is that ungrouped, or singlet, items may be included in the test as well. However, to guarantee test model identifiability, these items must have their guessing and slipping parameters set to zero, that is, $s_{j} = g_{j} = 0$ . Assuming that the singlet items are locally independent of each other and the item sets, their inclusion into the test likelihood is simple. This is done by multiplying the implied singlet item likelihood by the likelihoods of the item sets. For example, suppose there are $J - N$ items grouped into D item sets and N singlet items for a total of J items. Assuming a 2PNO model for the singlet items, the full likelihood for the test with n persons is

P (Y | τ, θ, s, g, γ, ξ) = \prod_{i = 1}^{n} [\prod_{d = 1}^{D} P (Y_{i d} = y_{i d} | τ, γ_{d}, s_{d}, g_{d}, θ_{i}) \prod_{j = 1}^{N} P (Y_{i j} = y_{i j} | ξ_{j}, θ_{i})],

where $γ_{d} = (a_{d}, b_{d})$ contains the slope and intercept parameters for item set d, $y_{i j}$ is the person i response to singlet item j, and $ξ_{j} = (a_{j}, b_{j})$ contains the slope and intercept parameters for singlet item j.

3. Bayesian Estimation of the Confirmatory Set-4PNO Model

In this section, a Bayesian formulation for the confirmatory Set-4PNO model is presented. First, a presentation of the model and a description of Bayesian priors are given. Second, the full conditional distribution for implementing an MCMC Gibbs sampler is provided.

3.1. Bayesian Formulation

A directed acyclic graph representing the confirmatory Set-4PNO model is given in Figure 2. The following Bayesian formulation is proposed for the confirmatory Set-4PNO model described above:

Y_{i j} | W_{i, τ_{j}} = w_{i d}, s_{j}, g_{j} \sim Bernoulli ({(1 - s_{j})}^{w_{i d}} g_{j}^{1 - w_{i d}}),

W_{i d} = I (W_{i d}^{*} > 0),

W_{i d}^{*} | θ_{i}, τ_{j} = d, a_{d}, b_{d} \sim N (a_{d} θ_{i} - b_{d},1),

θ \sim N (μ_{θ}, σ_{θ}^{2}),

ξ_{j} = (s_{j}, g_{j}) \sim Beta (a_{s}, b_{s}) Beta (a_{g}, b_{g}) I (0 \leq s_{j} + g_{j} < 1),

γ_{d} = (a_{d}, b_{d}) \sim N_{2} (μ_{γ}, Σ_{γ}) I (a_{d} > 0) .

Figure 2.

Direct acyclic graph of the set four-parameter normal ogive (Set-4PNO) model. Note. Squares and circles distinguish observed random variables and random parameters, respectively, rectangles are plates that denote variables that share common indices, and solid and dashed lines distinguish between stochastic and deterministic relationships.

The choice of priors is both a reflection of uncertainty in the parameters and of a desired simplicity in sampling. The chosen priors lend themselves well to a Gibbs sampling approach.

Equation 6 is the Bernoulli IRF for the observed response $Y_{i j}$ of individual i on item j, given the latent response $W_{i d}$ on dyad d and item slipping parameters, s_j , and guessing parameters, g_j . This is simply the IRF for the DINA model. Equations 7 and 8 constitute the 2PNO model portion of the model. In Equation 7, it is assumed that there is a continuous latent response variable $W_{i d}^{*}$ underlying $W_{i d}$ , where $W_{i d}$ takes on a value of one if $W_{i d}^{*} > 0$ and a value of zero otherwise. Equation 8, in turn, shows $W_{i d}^{*}$ following a normal distribution given the latent trait $θ_{i}$ , dyad slope parameter a_d , and dyad threshold parameter b_d . Together, Equations 6 –8 constitute the Set-4PNO model.

The rest of the above equations complete the Bayesian formulation of the Set-4PNO model by giving prior distributions for person, item, and item set parameters. Equation 9 shows that the latent trait $θ$ has a normal prior. This is a common prior for latent person parameters in many contexts. Because there is no natural unit for a latent variable, it is necessary to set a scale for $θ$ . Here, we adopt the common practice of standardizing $θ$ by setting $μ_{θ} = 0$ and $σ_{θ}^{2} = 1$ . This allows for parameter interpretations in terms of standardized units.

For the items parameters s_j and g_j , we give a joint beta prior distribution, as shown in Equation 10. If we have a simple monotonicity assumption for the model—that is, that $P (Y_{i j} = 1 | W_{i d} = 0) < P (Y_{i j} = 1 | W_{i d} = 1)$ —then we must include the restriction on the space of s_j and g_j , where $0 \leq s_{j} + g_{j} < 1$ (Junker & Sijtsma, 2001). This restriction induces a dependence between s_j and g_j and is the reason that independent beta prior distributions are not appropriate. Hence, the above monotonicity restriction is reflected in Equation 10. The hyperparameters a_s , b_s , a_g , and b_g are shape and scale parameters for the joint beta distribution. These can be set to reflect prior information as desired. For example, an uninformative uniform prior over the region ${(s_{j}, g_{j}) : 0 \leq s_{j}, g_{j} \leq 1; 0 \leq s_{j} + g_{j} < 1}$ is obtained when $a_{s} = b_{s} = a_{g} = b_{g} = 1$ .

Finally, for the item set parameters a_d and b_d , a truncated bivariate normal distribution prior is given. As shown in Equation 11, the distribution has mean $μ_{γ}$ and covariance $Σ_{γ}$ as well as the restriction that the item set slopes must be positive. The restriction given is commonly found in educational measurement, where it is assumed that the item score and the latent trait $θ$ have a positive relationship and could be relaxed in other measurement contexts. If prior information about a_d and b_d is known, then $μ_{γ}$ and $Σ_{γ}$ can be set to account for this. Conversely, if little is known about a_d and b_d , then this can be accounted for by setting $Σ_{γ}$ to have large elements on the diagonal (i.e., setting large variances) and zero-valued off-diagonal elements (i.e., setting the covariance between a_d and b_d to be 0).

3.2. Connection With Exploratory Dyad-4PNO

The difference between the current estimation approach and the one given by Kern and Culpepper (2020) is that they proposed using a uniform prior over the space of potential item dyads. That is, they proposed using $p (τ) \propto 1$ as a prior. In the current approach, it is possible to view the Bayesian formulation given as a special case of the previous method, where for a given item set index vector $τ_{1}$ , $P (τ = τ_{1}) = 1$ . In this way, it is possible to treat the index vector as fixed. One possible relaxation, when there are multiple competing item structures, say $τ_{1}, \dots, τ_{M}$ , but uncertainty about which is better, is to specify prior probabilities $P (τ = τ_{m})$ for each, such that $\sum_{m = 1}^{M} P (τ = τ_{m}) = 1$ . In the case where each is equally likely, then it would be natural to set $P (τ = τ_{m}) = \frac{1}{M}$ for all M competing models.

3.3. Posterior Approximation

Using the formulation given, we can construct a simple Gibbs sampling algorithm to approximate the full posterior distribution of the model parameters given the observed responses. The posterior distribution, in turn, can be used to determine point estimates, variances, and other summary statistics for each of the parameters in the model.

Gibbs sampling is an approach to finding the full posterior distribution of a set of parameters by iteratively sampling from each parameter’s conditional distribution, given the remaining parameters in the set and the observed data. Following is the series of sampling steps implied by the Bayesian formulation given above. Importantly, the choices for prior distributions generally allow for sampling from well-known distributions, which can greatly increase the efficiency of the MCMC algorithm. For each iteration of the Gibbs sampling algorithm, each step is taken in order. This is done sequentially a number of times until the Markov chain converges in distribution to the joint posterior distribution; this is the so-called burn-in sample. Once the chain has converged, each additional sample is said to be representative of the posterior distribution. This post-burn-in sample can then be used for describing the posterior distributions for each parameter.

1. Update for $W_{i d}$ , $i = 1, \dots, n$ , $d = 1, \dots, D$ :

W_{i d} | Y_{i d}, τ_{j} = d, s_{d}, g_{d}, a_{d}, b_{d} \sim Bernoulli ({\tilde{p}}_{i d}),

where

{\tilde{p}}_{i d} = \frac{P (Y_{i d} = y_{i d} | W_{i d} = 1, s_{d}, g_{d}) P (W_{i d} = 1 | θ_{i}, τ_{j} = d, a_{d}, b_{d})}{\sum_{w = 0}^{1} P (Y_{i d} = y_{i d} | W_{i d} = w, s_{d}, g_{d}) P (W_{i d} = w | θ_{i}, τ_{j} = d, a_{d}, b_{d})} .

2. Update for augmented data $W_{i d}^{*}$ , $i = 1, \dots, n$ , $d = 1, \dots, D$ :

W_{i d}^{*} | W_{i d} = w_{i d}, θ_{i}, τ_{j} = d, a_{d}, b_{d} \sim N (a_{d} θ_{i} - b_{d},1) I {(W_{i d}^{*} > 0)}^{w_{i d}} I {(W_{i d}^{*} < 0)}^{1 - w_{i d}} .

3. Update for $θ_{i}$ , $i = 1, \dots, n$ :

θ_{i} | W_{i}^{*}, τ, a, b \sim N (\frac{\sum_{d = 1}^{D} a_{d} (W_{i d}^{*} + b_{d}) + μ_{θ} σ_{θ}^{- 2}}{\sum_{d = 1}^{D} a_{d}^{2} + σ_{θ}^{- 2}}, {(\sum_{d = 1}^{D} a_{d}^{2} + σ_{θ}^{- 2})}^{- 1}) .

4. Update for $γ_{d} = (a_{d}, b_{d})$ , $d = 1, \dots, D$ :

γ_{d} | W_{d}^{*}, τ, θ \sim N_{2} (μ_{γ}^{*}, Σ_{γ}^{*}) I (a_{d} > 0),

where

X_{θ} = (θ, - 1_{n}),

Σ_{γ}^{*} = (X_{θ}^{'} X_{θ} + Σ_{γ}^{- 1})^{- 1},

μ_{γ}^{*} = Σ_{ξ}^{*} (X_{θ}^{'} W_{d}^{*} + Σ_{γ}^{- 1} μ_{γ}) .

Here, note that the added constraint that $a_{d} > 0$ makes the joint distribution above a truncated bivariate normal distribution. To sample $γ_{d}$ , first sample $a_{d} | W_{d}^{*}, θ$ from a univariate truncated normal distribution whose mean is the first element of $μ_{γ}^{*}$ and whose variance is the $(1, 1)$ element of $Σ_{γ}^{*}$ . Next, sample $b_{d} | a_{d}, W_{d}^{*}, θ$ from a normal distribution conditioned on a_d . The mean and variance of this distribution is simply

μ_{b_{d} | a_{d}}^{*} = μ_{b_{d}}^{*} + \frac{σ_{a_{d} b_{d}}^{*}}{σ_{a_{d}}^{2 *}} (a_{d} - μ_{a_{d}}^{*}),

and

σ_{b_{d} | a_{d}}^{2 *} = σ_{b_{d}}^{2 *} - \frac{σ_{a_{d} b_{d}}^{2 *}}{σ_{a_{d}}^{2 *}},

respectively. Using this approach yields a sample value of $γ_{2 d}$ from its truncated bivariate normal distribution.

5. Update $ξ_{j} = (s_{j}, g_{j})$ , $j = 1, \dots, J$ :

ξ_{j} | τ_{j}, W_{τ_{j}}, Y_{j} \propto g_{j}^{{\tilde{a}}_{g j} - 1} {(1 - g_{j})}^{{\tilde{b}}_{g j} - 1} s_{j}^{{\tilde{a}}_{s j} - 1} {(1 - s_{j})}^{{\tilde{a}}_{s j} - 1} I (0 \leq g_{j} < 1 - s_{j}) I (0 \leq s_{j} < 1),

where

{\tilde{a}}_{g j} = a_{g} + \sum_{i = 1}^{n} y_{i j} (1 - w_{i, τ_{j}}),

{\tilde{b}}_{g j} = b_{g} + \sum_{i = 1}^{n} (1 - y_{i j}) (1 - w_{i, τ_{j}}),

{\tilde{a}}_{s j} = a_{s} + \sum_{i = 1}^{n} (1 - y_{i j}) w_{i, τ_{j}},

{\tilde{b}}_{s j} = b_{s} + \sum_{i = 1}^{n} y_{i j} w_{i, τ_{j}} .

In this instance, as in Step 4, the constraint that $0 \leq g_{j} < 1 - s_{j}$ and $0 < s_{j} < 1$ has the effect of making it somewhat more difficult to sample from the resulting analytical joint distribution of $ξ_{j}$ . To do so, one may take the approach given in either Culpepper (2015) or Culpepper (2016).

4. Monte Carlo Simulation Studies

In this section, the proposed model is investigated using two Monte Carlo simulation studies. First, the efficacy of the above MCMC is examined. Specifically, this is an examination into how well the algorithm recovers model parameters when the structure of estimated model is properly specified as the one given in the data generating model. In the second simulation study, the ability of the deviance information criterion (DIC)—both the standard version and the marginal version—to choose the correct model is investigated.

4.1. Simulation 1

4.1.1. Method

In this simulation study, we seek to determine how well the algorithm recovers model parameters in a variety of contexts when the model structure is properly specified in estimation. Five factors were manipulated, largely following the first study in Kern and Culpepper (2020). First, guessing parameters were either generated from $g_{j} \sim Unif (.15, .25)$ or from $g_{j} \sim Unif (.35, .45)$ . These two conditions are labeled the low guessing condition and the high guessing condition, respectively, with values centered at .20 and .40. Second, the slipping parameters were generated based on a low slipping condition—that is, $s_{j} \sim Unif (.05, .15)$ —or a high slipping condition—that is, $s_{j} \sim Unif (.15, .25)$ . Together, the guessing and slipping conditions increase the overall level of uncertainty in the responses and are expected to be related to decreases in model parameter recovery. Third, the number of item sets was manipulated to be either $D = 10$ or $20$ . Fourth, the number of items per item set was $G = 2$ , 3, 4, or 5. Together, the number of item sets and number of items per item set determine the overall number of items in the test, ranging from $J = 20$ (10 sets of two items) to $J = 100$ (20 sets of five items). Finally, the sample size was $n = 500$ , 1,000, or 2,000. All five factors were fully crossed for a total of $2 (guessing) \times 2 (slipping) \times 2 (item sets) \times 4 (items per item set) \times 3 (sample size) = 96$ conditions.

Latent trait values $θ$ were sampled as $θ \sim N (0, 1)$ , and the item set parameters a_d and b_d were sampled as $a_{d} \sim Unif (1.0, 3.5)$ and $b_{d} \sim N (0, 1)$ . These distributions were chosen to approximate estimated values seen in the investigations of this modeling approach with real data. Each of the 96 conditions was replicated 100 times. Preliminary investigations showed that a chain length of 10,000 with a burn-in of 5,000 was sufficient to achieve parameter convergence using the Gelman–Rubin $\hat{R}$ statistic (Gelman & Rubin, 1992), so this setting is used throughout the simulation. Finally, the hyperparameters for the priors are as follows: $μ_{θ} = 0$ , $σ_{θ}^{2} = 1$ , $μ_{γ} = (2, 0)$ , $Σ_{γ} = I_{2}$ , and $a_{g} = b_{g} = a_{s} = b_{s} = 1$ . Because the prior for a_d is a truncated normal, the choice of hyperparameters $μ_{a}$ and $σ_{a}^{2}$ must be chosen carefully, so that the expected value and standard deviation of the distribution of a_d are as desired. Here, with $μ_{a} = 2$ and $σ_{a}^{2} = 1$ , it can be shown that $E (a_{d}) \approx 2.06$ and $SD (a_{d}) \approx 0.94$ .

4.1.2. Parameter recovery

The efficacy of model parameter recovery was determined using root mean squared error (RMSE), bias, and correlation with the true parameter. Each of these measures is averaged over all the $R = 100$ replications. Specifically, for a given parameter $ζ$ , we calculated these as follows:

RMSE (\hat{ζ}) = \frac{1}{M} \sum_{m = 1}^{M} \sqrt{\frac{1}{R} \sum_{r = 1}^{R} {({\hat{ζ}}_{m}^{(r)} - ζ_{m}^{(r)})}^{2}},

Bias (\hat{ζ}) = \frac{1}{M} \sum_{m = 1}^{M} \frac{1}{R} \sum_{r = 1}^{R} ({\hat{ζ}}_{m}^{(r)} - ζ^{(r)}),

Cor (\hat{ζ}) = \frac{1}{R} \sum_{r = 1}^{R} \frac{1}{M - 1} \sum_{m = 1}^{M} \frac{({\hat{ζ}}_{m}^{(r)} - {\bar{\hat{ζ}}}^{(r)}) (ζ^{(r)} - {\bar{ζ}}^{(r)})}{s_{\hat{ζ}}^{(r)} s_{ζ}^{(r)}},

where M is either the number of items, item sets, or examinees (dependent upon the given parameter $ζ$ ); ${\hat{ζ}}_{m}^{(r)}$ is the mean of the posterior samples for the mth parameter from replication r, $ζ_{m}^{(r)}$ ; the values ${\bar{\hat{ζ}}}^{(r)}$ and ${\bar{ζ}}^{(r)}$ are

{\bar{\hat{ζ}}}^{(r)} = \frac{1}{M} \sum_{m = 1}^{M} {\hat{ζ}}_{m}^{(r)},

{\bar{ζ}}^{(r)} = \frac{1}{M} \sum_{m = 1}^{M} ζ_{m}^{(r)};

and $s_{ζ^}^{(r)}$ and $s_{ζ}^{(r)}$ are within replication standard deviations for the estimated and true parameters, respectively.

4.1.3. Results

As overall result patterns are as expected for sample size, to not clutter this space, only results for the sample size $N = 2, 000$ are shown; the interested reader may contact the author for this information.

Simulation results for the item set parameters a and b are given in Table 1. First, while not shown, it was found that with increasing sample size, item set parameters show increasingly better recovery. This is especially the case with the b-parameters, which show near perfect recovery for $N = 2, 000$ . Second, as the number of item sets D and items per item set J increase, accurate recovery of the item set parameters increases, though this effect is stronger for increasing J than increasing D. This makes sense, because increasing J essentially increases the “sample size” of each item set parameter—that is, there are more observations associated with the given item set parameter. On the other hand, increasing D does not have this same effect, instead increasing the number of item set parameters being estimated. Finally, with increasing levels of slipping and guessing, item set parameter recover decreases. This is due to the increasing uncertainty in the responses that results for higher levels of slipping and guessing.

Table 1.

Simulation 1 Recovery of Item Set Parameters a and b for $N = 2,000$

D	G	Guessing	Slipping	Bias(a)	RMSE(a)	Cor(a)	Bias(b)	RMSE(b)	Cor(b)
10	2	Low	Low	.02	.25	.94	−.01	.13	0.99
10	2	Low	High	.03	.32	.89	.00	.18	0.99
10	2	High	Low	.03	.40	.84	−.02	.22	0.98
10	2	High	High	.02	.55	.67	.01	.36	0.93
10	3	Low	Low	.00	.19	.96	.00	.11	1.00
10	3	Low	High	.01	.25	.95	.01	.13	0.99
10	3	High	Low	.02	.30	.91	.00	.16	0.99
10	3	High	High	.01	.41	.80	.00	.24	0.97
10	4	Low	Low	.01	.18	.97	.00	.09	1.00
10	4	Low	High	.03	.21	.96	.00	.10	1.00
10	4	High	Low	.02	.24	.94	.01	.12	0.99
10	4	High	High	.00	.36	.86	.01	.19	0.98
10	5	Low	Low	.00	.16	.97	.01	.09	1.00
10	5	Low	High	.01	.19	.97	.00	.10	1.00
10	5	High	Low	.02	.22	.96	−.01	.12	0.99
10	5	High	High	.00	.29	.91	.00	.16	0.99
20	2	Low	Low	.00	.21	.96	.00	.12	0.99
20	2	Low	High	.02	.26	.94	−.01	.15	0.99
20	2	High	Low	.01	.33	.90	−.02	.20	0.98
20	2	High	High	.02	.47	.78	−.01	.28	0.96
20	3	Low	Low	.00	.17	.97	.00	.10	1.00
20	3	Low	High	.00	.21	.96	.00	.12	0.99
20	3	High	Low	.02	.25	.94	.01	.14	0.99
20	3	High	High	.00	.37	.87	−.01	.20	0.98
20	4	Low	Low	.00	.15	.98	.01	.09	1.00
20	4	Low	High	.01	.19	.97	.01	.10	1.00
20	4	High	Low	.00	.22	.95	.00	.13	0.99
20	4	High	High	.01	.29	.92	−.01	.17	0.99
20	5	Low	Low	.00	.15	.98	.00	.09	1.00
20	5	Low	High	.01	.17	.97	.00	.10	1.00
20	5	High	Low	.00	.19	.97	−.01	.11	1.00
20	5	High	High	.01	.26	.94	−.01	.14	0.99

Results for the item parameters s and g are given in Table 2. First, we note that bias and RMSE are especially small for s and g, which is notable in that these parameters are notoriously difficult to estimate in IRT models. This result shows some of the usefulness of this model in accounting for guessing and slipping effects. Furthermore, as with the item set parameters, it was found that increasing sample size improves item parameter recovery. Other factors that affect item parameter recovery include the number of items per item set J and the levels of slipping and guessing. Specifically, as J increases or as slipping and guessing levels decrease, item parameters show better recovery. Interestingly, it seems that increasing the number of item sets D does not affect item parameter recovery very much.

Table 2.

Simulation 1 Recovery of Item Parameters g and s for $N = 2,000$

D	G	Guessing	Slipping	Bias(g)	RMSE(g)	Cor(g)	Bias(s)	RMSE(s)	Cor(s)
10	2	Low	Low	.00	.02	.83	.00	.02	.86
10	2	Low	High	.00	.02	.81	.00	.02	.79
10	2	High	Low	.00	.03	.74	.00	.02	.82
10	2	High	High	.00	.03	.66	.00	.03	.71
10	3	Low	Low	.00	.02	.88	.00	.01	.92
10	3	Low	High	.00	.02	.86	.00	.02	.86
10	3	High	Low	.00	.02	.82	.00	.02	.88
10	3	High	High	.00	.02	.76	.00	.02	.81
10	4	Low	Low	.00	.02	.89	.00	.01	.92
10	4	Low	High	.00	.02	.87	.00	.02	.87
10	4	High	Low	.00	.02	.83	.00	.01	.91
10	4	High	High	.00	.02	.78	.00	.02	.83
10	5	Low	Low	.00	.01	.90	.00	.01	.93
10	5	Low	High	.00	.02	.88	.00	.01	.89
10	5	High	Low	.00	.02	.85	.00	.01	.91
10	5	High	High	.00	.02	.82	.00	.02	.86
20	2	Low	Low	.00	.02	.85	.00	.02	.88
20	2	Low	High	.00	.02	.83	.00	.02	.81
20	2	High	Low	.00	.02	.77	.00	.02	.83
20	2	High	High	.00	.03	.70	.00	.03	.74
20	3	Low	Low	.00	.02	.88	.00	.01	.92
20	3	Low	High	.00	.02	.86	.00	.02	.86
20	3	High	Low	.00	.02	.82	.00	.01	.90
20	3	High	High	.00	.02	.79	.00	.02	.81
20	4	Low	Low	.00	.01	.89	.00	.01	.93
20	4	Low	High	.00	.02	.88	.00	.02	.88
20	4	High	Low	.00	.02	.84	.00	.01	.91
20	4	High	High	.00	.02	.81	.00	.02	.84
20	5	Low	Low	.00	.01	.90	.00	.01	.94
20	5	Low	High	.00	.01	.89	.00	.02	.89
20	5	High	Low	.00	.02	.84	.00	.01	.92
20	5	High	High	.00	.02	.82	.00	.02	.86

Lastly, results for the latent trait $θ$ are in Table 3. As with the dyad and item parameters, bias is very small—nearly zero—in all conditions. Furthermore, we found that sample size has a negligible effect on the recovery of the latent trait. On the other hand, it was found that increasing the number of item sets and the number of items per item set had a positive impact on the estimation of $θ$ . This is expected, because increasing either of these factors amounts to an increase in test length. Lastly, increasing levels of guessing and slipping has detrimental effects on latent trait estimation. Together, the evidence is as expected and shows that the algorithm can effectively recover model parameters.

Table 3.

Simulation 1 Recovery of Person Parameters $θ$ for $N = 2,000$

D	G	Guessing	Slipping	Bias( $θ$ )	RMSE( $θ$ )	Cor( $θ$ )
10	2	Low	Low	.00	.45	.89
10	2	Low	High	.00	.50	.86
10	2	High	Low	.00	.54	.84
10	2	High	High	.00	.64	.77
10	3	Low	Low	.00	.42	.91
10	3	Low	High	.00	.46	.89
10	3	High	Low	.00	.49	.87
10	3	High	High	.00	.57	.82
10	4	Low	Low	.00	.40	.92
10	4	Low	High	.00	.43	.90
10	4	High	Low	.00	.46	.89
10	4	High	High	.00	.53	.85
10	5	Low	Low	.01	.39	.92
10	5	Low	High	.00	.41	.91
10	5	High	Low	.00	.44	.90
10	5	High	High	.00	.50	.87
20	2	Low	Low	.00	.37	.93
20	2	Low	High	.00	.41	.91
20	2	High	Low	.00	.44	.90
20	2	High	High	.00	.52	.86
20	3	Low	Low	.00	.34	.94
20	3	Low	High	.00	.37	.93
20	3	High	Low	.01	.40	.92
20	3	High	High	.00	.46	.89
20	4	Low	Low	.01	.32	.95
20	4	Low	High	.00	.35	.94
20	4	High	Low	.00	.37	.93
20	4	High	High	.00	.43	.90
20	5	Low	Low	.00	.32	.95
20	5	Low	High	.00	.34	.94
20	5	High	Low	.00	.35	.94
20	5	High	High	.00	.41	.91

4.2. Simulation 2

4.2.1. Method

In this simulation study, the ability of the DIC to select the correct model in several contexts is examined. In this simulation, DIC and its marginal counterpart (i.e., marginal DIC) are studied. The DIC was chosen for this simulation study because of its ubiquity in Bayesian model selection. The marginal DIC was chosen because it has been argued that it is the more appropriate index in latent variable modeling contexts, that is, that an index selecting between latent variable models should have the latent variable integrated out (Merkle et al., 2019).

This study was conducted in two parts, with each part consisting of a different generating model. In the first part, the generating model was the Set-4PNO model, whereas in the second part, the generating model was the 2PNO model. This approach was used as an approximation to model complexity, because the Set-4PNO model has more parameters than the 2PNO model. Ideally, the criterion indices would do as they are designed to do and choose the generating model more often, regardless of model complexity. In both generating model situations, five competing models are fit: the Set-4PNO, exploratory Dyad-4PNO, 4PNO, 3PNO, and 2PNO models.

Using the Set-4PNO model as the generating model, four factors were manipulated. The first factor is the number of item sets $D = 10$ or $20$ , and the second factor is the number of items per item set G (i.e., 3 or 4). Together, these two factors control the number of items in the overall test as well the structure for the estimated Set-4PNO model. The third and fourth factors control the guessing and slipping levels (either low or high for each factor) in the generating model. These were as in Simulation 1. That is, guessing was either low [i.e., $g_{j} \sim U n i f (.15, .25)$ ] or high [i.e., $g_{j} \sim U n i f (.35, .45)$ ], and slipping was either low [i.e., $s_{j} \sim U n i f (.05, .15)$ ] or high [i.e., $s_{j} \sim U n i f (.15, .25)$ ]. Together, these control the overall response uncertainty. The four factors were fully crossed yielding a total of $2 \times 2 \times 2 \times 2 = 16$ conditions.

Using the 2PNO model as the generating model, two factors were manipulated. As in the previous paragraph, these were the number of item sets $D = 10$ or $20$ , and the number of items per item set $G = 3$ or 4. These controlled the total number of items in the overall test and also the structure for estimating the Set-4PNO model. With the 2PNO model as the generating model, the fitted Set-4PNO model was randomly chosen among the possible Set-4PNO models. This was done given that each one is expected to be poorly fit compared to the 2PNO, though there will be some distribution of model fits across all possible Set-4PNO models.

In all cases, latent trait values were sampled as $θ \sim N (0, 1)$ . Furthermore, the a- and b-parameters were sampled as $a \sim Unif (1.0, 3.5)$ and $b \sim N (0, 1)$ . As in Simulation 1, a chain length of 10,000 with a burn-in of 5,000 was used throughout the simulation. Also like in Simulation 1, the hyperparameters for the priors used in estimation are as follows: $μ_{θ} = 0$ , $σ_{θ}^{2} = 1$ , $μ_{γ} = (2, 0)$ , $Σ_{γ} = I_{2}$ , and $a_{g} = b_{g} = a_{s} = b_{s} = 1$ . Finally, as the computation of the marginal DIC is relatively time-consuming—it requires the integration of $θ$ out of the likelihood at every post-burn-in MCMC sample—each condition was only replicated 50 times.

4.2.2. Results

The main results in this section concern the proportion of times the correct model was chosen under the various conditions. These results are given in Tables 4 and 5. Some important trends can be found. First, when the 2PNO model was the generating model, both DIC and marginal DIC correctly selected the 2PNO model in every instance of the simulation. Thus, when no slipping or guessing effects are present, this seems to be easily detected via both indices.

Table 4.

Proportion of Models Chosen by Deviance Information Criterion by Condition

D	G	Guessing	Slipping	Set-4PNO	Dyad-4PNO	4PNO	3PNO	2PNO
Generating model: Set-4PNO
10	3	Low	Low	1.00	.00	.00	.00	0.00
10	3	Low	High	1.00	.00	.00	.00	0.00
10	3	High	Low	0.94	.06	.00	.00	0.00
10	3	High	High	0.14	.50	.36	.00	0.00
10	4	Low	Low	1.00	.00	.00	.00	0.00
10	4	Low	High	1.00	.00	.00	.00	0.00
10	4	High	Low	1.00	.00	.00	.00	0.00
10	4	High	High	0.42	.48	.10	.00	0.00
20	3	Low	Low	1.00	.00	.00	.00	0.00
20	3	Low	High	1.00	.00	.00	.00	0.00
20	3	High	Low	1.00	.00	.00	.00	0.00
20	3	High	High	0.80	.18	.02	.00	0.00
20	4	Low	Low	1.00	.00	.00	.00	0.00
20	4	Low	High	1.00	.00	.00	.00	0.00
20	4	High	Low	1.00	.00	.00	.00	0.00
20	4	High	High	1.00	.00	.00	.00	0.00
Generating model: 2PNO
10	3	—	—	0.00	.00	.00	.02	0.98
10	4	—	—	0.00	.00	.00	.00	1.00
20	3	—	—	0.00	.00	.00	.00	1.00
20	4	—	—	0.00	.00	.00	.00	1.00

Note. 4PNO = four-parameter normal ogive; 3PNO = three-parameter normal ogive; 2PNO = two-parameter normal ogive.

Table 5.

Proportion of Models Chosen by Marginal Deviance Information Criterion by Condition

D	G	Guessing	Slipping	Set-4PNO	Dyad-4PNO	4PNO	3PNO	2PNO
Generating model: Set-4PNO
10	3	Low	Low	1.00	.00	.00	.00	0.00
10	3	Low	High	1.00	.00	.00	.00	0.00
10	3	High	Low	1.00	.00	.00	.00	0.00
10	3	High	High	1.00	.00	.00	.00	0.00
10	4	Low	Low	1.00	.00	.00	.00	0.00
10	4	Low	High	1.00	.00	.00	.00	0.00
10	4	High	Low	1.00	.00	.00	.00	0.00
10	4	High	High	1.00	.00	.00	.00	0.00
20	3	Low	Low	1.00	.00	.00	.00	0.00
20	3	Low	High	1.00	.00	.00	.00	0.00
20	3	High	Low	1.00	.00	.00	.00	0.00
20	3	High	High	1.00	.00	.00	.00	0.00
20	4	Low	Low	1.00	.00	.00	.00	0.00
20	4	Low	High	1.00	.00	.00	.00	0.00
20	4	High	Low	1.00	.00	.00	.00	0.00
20	4	High	High	1.00	.00	.00	.00	0.00
Generating model: 2PNO
10	3	—	—	0.00	.00	.00	.02	0.98
10	4	—	—	0.00	.00	.00	.00	1.00
20	3	—	—	0.00	.00	.00	.00	1.00
20	4	—	—	0.00	.00	.00	.00	1.00

Note. 4PNO = four-parameter normal ogive; 3PNO = three-parameter normal ogive; 2PNO = two-parameter normal ogive.

The same results do not hold in general when the Set-4PNO model is the generating model. Consider the case using the DIC for model selection. In this case, we find two things. It seems that increasing the effect sizes for guessing and slipping, that is, increasing response uncertainty, increases the chances for the index to incorrectly choose the wrong model. This is especially true with decreasing numbers of items and smaller item set sizes. In particular, with high guessing and slipping, the exploratory Dyad-4PNO is incorrectly chosen much of the time. For lower values of D and G, the 4PNO model is incorrectly chosen often as well, though less often that the Dyad-4PNO. In the simulation, neither the 3PNO nor the 2PNO were incorrectly chosen.

Using marginal DIC for model selection given the Set-4PNO is the generating model looks to be a generally better approach than using the DIC. First, increasing guessing and slipping effects does not affect model selection uncertainty with the marginal DIC nearly as much as with the DIC. In all conditions, the marginal DIC never chose the 2PNO, 3PNO, or 4PNO models. In particular, when D is low, and guessing and slipping are high, the Dyad-4PNO was never incorrectly chosen a very small amount of the time. This is clearly much less often than with the standard DIC.

Taken together, these results show that the Set-4PNO model is distinguishable from other potential competing models. In particular, the marginal DIC is a good metric for making a model selection among the competing models. Furthermore, this shows that if a nondyadic structure for the items is present, then the set-4PNO model will generally fit better than the exploratory Dyad-4PNO model does, as one would expect.

5. Application

In this section, we provide a real data example to demonstrate the use of this model. The data consist of an experimental IQ data set found on Open Psychometrics website. The data include the responses of 400 examinees to 25 different items. For each item, examinees are shown a patterned $3 \times 3$ grid of tiles, in which the lower left tile is missing. The examinee is asked to select the tile from eight options that best fits in the missing spot in the grid. In accordance with the use of these data in prior literature (Chen et al., 2020; Kern & Culpepper, 2020), a 20-item subset comprised of the 20 items with a correct response rate of at least 25% will be used. This is a voluntary Internet sample made up of 71% males and 29% females with an average age of 33.23. There was an overall missingness rate of 3.25%. For purposes of this example, missing responses were coded as incorrect responses. This will be revisited later.

To demonstrate that allowing for item sets with more than two items can be beneficial, two confirmatory models were fit. For the first model, the dyads chosen by the exploratory Dyad-4PNO in Kern and Culpepper (2020) were fit to the data. The model was chosen as it is presumably the best possible set of two-item sets for these data. However, as shown in Table 6, while several dyad groupings are quite stable in terms of percentage of time items are matched in the exploratory algorithm, four dyads—namely, Dyads 3, 6, 8, and 10—exhibit much less stability. It can be hypothesized that the reason for this instability might be that some residual dependencies between dyads may be present. In this way, two item sets of size four, each consisting of the items from two dyads in the exploratory solution, can be constructed. Implicitly, choosing to merge a pair of dyads imposes an equality constraint on the item set parameters, so to choose which pairs of dyads to merge, we chose the pairs with the smallest differences between the dyad parameters. This leads to an eight-set Set-4PNO model with six item sets of two items and two item sets of four items. Additionally, the data were also fit to the more standard 2PNO, 3PNO, and 4PNO models. This is done as a check to determine whether the Set-4PNO is a better fit to the data than other standard IRT models. For all five models, the data were fit using Gibbs sampling.

Table 6.

Exploratory Dyad-4PNO Solution

Dyad	Items	% Matched	a_d	b_d
1	${1, 2}$	100.00	1.74	−2.22
2	${3, 8}$	99.77	2.48	−2.47
3	${4, 6}$	82.37	2.26	1.20
4	${5, 7}$	99.83	1.70	−0.32
5	${9, 10}$	99.96	1.76	−0.26
6	${11, 12}$	81.98	2.23	1.36
7	${13, 16}$	97.02	3.20	−0.11
8	${14, 17}$	81.60	3.07	−0.41
9	${15, 20}$	99.83	2.56	−1.35
10	${18, 19}$	84.25	2.60	0.01

Note. 4PNO = four-parameter normal ogive.

To implement the MCMC algorithms, it is necessary to set the hyperparameters for the prior distributions. For the person parameters $θ_{i}$ , we used a standard normal prior. For the item parameters s_j and g_j , the joint beta prior as described above was set with $a_{s} = b_{s} = a_{g} = b_{g} = 1$ , yielding a uniform distribution over the lower left triangular region of the unit square. For the item set parameters a_d and b_d , a prior mean of $μ_{γ_{2}} = (2, 0)$ and prior covariance of $Σ_{γ_{2}} = I_{2}$ are used; this choice was made recognizing the slopes to often be somewhat larger than 0, and because there is no information on a particular correlative relationship between the slopes and intercepts prior to estimation. Furthermore, these were the priors used in Kern and Culpepper (2020). Similar decisions for prior distributions are made for each model. The sensitivity of the model choice to the prior is also investigated in the next section.

For all models, parameter convergence was assessed using the potential scale reduction factor $\hat{R}$ (Gelman & Rubin, 1992). They suggested that if $\hat{R} < 1.1$ for all model parameters, then the chain has successfully converged to the posterior distribution. In this study, we used five independent chains with random start values to compute $\hat{R}$ . To make model comparisons as fair as possible, we determined the number of MCMC iterations as the approximate number of iterations needed, so that the MCMC chains of the parameters for each of the five competing models (i.e., the Dyad-4PNO, the Set-4PNO, 2PNO, 3PNO, and 4PNO) simultaneously converged to their respective posteriors. It was found that the criterion for convergence was met with approximately 40,000 iterations with a burn-in sample of 20,000. Model parameter estimates are taken to be the simple averages of the samples from the posterior distribution, that is, those MCMC samples taken post-burn-in. Model parameter estimates for the Dyad-4PNO and Set-4PNO models are given in Table 7, and model parameter estimates for the 2PNO, 3PNO, and 4PNO models are given in Table 8.

Table 7.

Item Set and Item Parameter Estimates for the Dyad-4PNO and the Set-4PNO Models

Dyad-4PNO Model							Set-4PNO Model
Dyad	Item	a_d	b_d	$b_{d}^{*}$	g_j	s_j	Set	Item	a_d	b_d	$b_{d}^{*}$	g_j	s_j
1	1	1.70	−2.20	−1.29	.18	.02	1	1	1.76	−2.23	−1.26	.18	.02
	2				.35	.02		2				.35	.02
2	3	2.45	−2.47	−1.00	.25	.14	2	3	2.53	−2.47	−0.97	.24	.14
	8				.47	.11		8				.48	.11
3	4	2.31	1.21	0.52	.22	.22	3	4	1.84	0.94	0.52	.22	.24
	6				.46	.03		6				.47	.06
6	11	2.28	1.38	0.61	.15	.19		11				.16	.29
	12				.17	.43		12				.16	.48
4	5	1.63	−0.33	−0.21	.51	.16	4	5	1.76	−0.36	−0.17	.51	.16
	7				.18	.21		7				.18	.21
5	9	1.75	−0.26	−0.16	.16	.19	5	9	1.86	−0.22	−0.09	.17	.19
	10				.36	.15		10				.35	.14
7	13	3.25	−0.13	−0.04	.15	.21	6	13	3.35	−0.09	−0.02	.16	.21
	16				.21	.25		16				.20	.24
8	14	3.02	−0.48	−0.14	.21	.04	7	14	3.12	−0.32	−0.09	.26	.05
	17				.10	.61		17				.10	.60
10	18	2.63	−0.02	−0.01	.15	.38		18				.15	.42
	19				.12	.33		19				.14	.38
9	15	2.56	−1.36	−0.53	.17	.04	8	15	2.51	−1.25	−0.50	.17	.04
	20				.29	.23		20				.30	.23

Note. 4PNO = four-parameter normal ogive.

Table 8.

Item Parameter Estimates for the 4PNO, 3PNO, and 2PNO Models

	4PNO					3PNO				2PNO
Item	a_j	b_j	$b_{j}^{*}$	g_j	s_j	a_j	b_j	$b_{j}^{*}$	g_j	a_j	b_j	$b_{j}^{*}$
1	1.88	−2.35	−1.25	.15	.02	1.23	−1.60	−1.30	0.17	1.14	−1.70	−1.49
2	2.01	−2.86	−1.43	.14	.02	1.05	−1.65	−1.57	0.18	0.97	−1.73	−1.79
3	1.68	−1.66	−0.99	.16	.11	0.71	−0.65	−0.92	0.14	0.65	−0.79	−1.22
4	2.31	1.54	0.67	.23	.15	1.52	1.35	0.88	0.22	0.51	0.29	0.57
5	1.49	−0.18	−0.12	.42	.11	0.76	−0.04	−0.05	0.33	0.48	−0.58	−1.21
6	1.12	−0.19	−0.17	.18	.06	1.03	−0.06	−0.06	0.18	0.75	−0.39	−0.52
7	1.35	−0.07	−0.05	.16	.14	1.00	0.23	0.23	0.15	0.75	−0.08	−0.11
8	1.94	−2.50	−1.28	.21	.10	0.51	−0.79	−1.55	0.22	0.49	−1.00	−2.05
9	1.91	−0.37	−0.20	.12	.20	1.05	0.14	0.14	0.09	0.89	−0.05	−0.05
10	1.95	−0.58	−0.30	.28	.15	0.99	−0.01	−0.01	0.24	0.69	−0.41	−0.59
11	2.20	1.40	0.64	.14	.15	1.39	1.14	0.82	0.13	0.71	0.49	0.68
12	1.90	1.43	0.75	.15	.31	0.99	1.34	1.35	0.14	0.48	0.63	1.33
13	2.82	0.31	0.11	.15	.13	1.72	0.54	0.31	0.15	0.97	0.04	0.05
14	2.84	−0.52	−0.18	.19	.04	1.90	−0.26	−0.14	0.16	1.27	−0.56	−0.44
15	2.37	−1.24	−0.52	.14	.04	1.67	−0.85	−0.50	0.11	1.48	−0.99	−0.67
16	2.82	−0.46	−0.16	.14	.25	1.08	0.28	0.26	0.11	0.83	0.03	0.04
17	2.70	−0.03	−0.01	.09	.57	0.68	0.96	1.41	0.06	0.53	0.72	1.35
18	2.03	−0.19	−0.09	.10	.35	0.87	0.57	0.66	0.08	0.67	0.34	0.51
19	2.54	0.31	0.12	.13	.28	1.12	0.65	0.58	0.10	0.79	0.32	0.40
20	2.09	−0.99	−0.48	.23	.19	0.73	−0.16	−0.22	0.15	0.62	−0.36	−0.59

Note. 4PNO = four-parameter normal ogive; 3PNO = three-parameter normal ogive; 2PNO = two-parameter normal ogive.

As a first diagnostic, we investigated the latent trait estimates given across the alternative models. Specifically, a scatter plot matrix containing plots of the $θ$ estimates for each of the five models against each other is given in Figure 3. If slipping and guessing effects are not present, there should not be any major systematic differences among the estimates. In particular, the presence of guessing effects should bias those on the low end of the latent trait scale and the presence of slipping effects should bias those on the high end of the latent trait scale in models that do not account for these effects. These patterns would manifest in plots that look nonlinear with flattening at the high and low ends of the scale, similar to an s-shaped ogive curve. Indeed, both of these patterns are observed in the scatter plots, giving support for the 4PNO models. Moreover, interestingly, when comparing the restricted 4PNO models (i.e., Dyad-4PNO and Set-4PNO) to the standard unrestricted 4PNO model, there also seems to be evidence of potential unaccounted-for low-end and high-end bias. On the other hand, the Dyad-4PNO and Set-4PNO have similar trait estimates throughout the scale.

Figure 3.

Scatter plot matrix of latent trait estimates for each of the five competing models.

Models were also compared using a marginal DIC that was computed for each. DIC is the average deviance of the posterior distribution with a penalty for the effective number of parameters. As previously described, the marginal DIC is like the standard DIC, but the deviance is computed with $θ$ integrated out. This approach was taken following arguments that it is better able to determine the best-fitting model than the standard unmarginalized DIC in latent variable modeling contexts (Merkle et al., 2019) and simulation results that give evidence that this approach is better for the current context. Specifically, the model with the smallest marginal DIC is taken to the one that best fits the data taking into account a trade-off between model fit and model complexity. Table 9 shows the model fit results, including marginal DIC and the effective number of parameters p_d , for each of the five competing models. Here, the eight-set Set-4PNO model ( $m a r g i n a l D I C = 8651.73$ and $p_{d} = 50.86$ ) fits the IQ data better than the other four models.

Table 9.

Model Fit for the Dyad-4PNO, Set-4PNO, 4PNO, 3PNO, and 2PNO Models

Index	Dyad-4PNO	Set-4PNO	4PNO	3PNO	2PNO
Marginal DIC	8,658.02	8,651.73	8,693.23	8,757.94	8,750.98
p_d	51.95	50.86	53.72	47.20	39.51

Note. 4PNO = four-parameter normal ogive; 3PNO = three-parameter normal ogive; 2PNO = two-parameter normal ogive; DIC = deviance information criterion.

To show that it was the particular proposed Set-4PNO model that fit the data well and not just any eight sets of items, two cases of competing models were compared with the proposed model. In the first case of models, an item set index vector $τ_{m}$ was chosen (uniform) randomly from the set of all possible item set index vectors for a set of 20 items partitioned into six sets of two items and two sets of four items. Here, 1,000 random models were chosen. For the second competing model case, competing models were chosen to be as similar as possible to the proposed Set-4PNO model. Specifically, these models are identical to the proposed model except for two items with swapped item set indices. In this case, there were 172 unique close competing models. For each competing model in both cases, the Set-4PNO model implied by the competing model’s item set index vector $τ_{m}$ was estimated using the previously presented Gibbs sampling approach with 40,000 MCMC iterations and a burn-in sample of 20,000. Finally, the marginal DIC was computed.

A competing model is said to fit better than the proposed Set-4PNO if it has a smaller marginal DIC. In 1,000 randomly chosen models, it was found that none fit better than the proposed model. Similarly, in 172 unique close models, none fit better than the proposed model, though in general, the difference in marginal DIC was smaller than for the randomly chosen models. However, the smallest difference among the close models to the proposed model was 9.98. These results are shown in Figure 4. This shows (1) that there is a great deal of heterogeneity in fit among the possible models with a similar structure, (2) that the chosen model is most likely better than one that we could have arrived upon by chance alone, and (3) that there are no extremely close models that are better. Note that this does not rule out the possibility that some model is better. More generally, this approach shows how one could assess the goodness of a confirmatory model against those of a similar structure, where the smaller the proportion of random models with a smaller DIC than the proposed model is better.

Figure 4.

Histograms of marginal deviance information criteria (DICs) for models with similar item set structure to the proposed set four-parameter normal ogive (Set-4PNO) model. Plot A is the histogram for the randomly chosen item set. Plot B is the histogram for the close item sets. The black vertical line shows the marginal DIC of the proposed Set-4PNO model.

One may question the reason for increased a-parameters in the 4PMs. Experience shows that a-parameters tend to be most affected by the inclusion of guessing and slipping effects into the response model. Note that model fitting in the IRT context is equivalent to finding the set of ICCs under the constraints of a given model that best describe the patterns in the data. Then, when fitting two separate IRT models to the same data, the differences in the models’ implied ICCs would typically be close subject to model constraints. When minimizing the differences between the 2PNO model ICC to the 4PM ICC, the 2PNO slope must be smaller than the 4PM slope. In this way, if guessing and slipping effects are truly present, then four-parameter-type models’ slopes larger than a 2PNO model fit on the same data.

5.1. Alternative Prior Investigation

It is possible that parameter estimates and outcomes of model selection could be sensitive to the choice of a prior distribution. This has been cited before as a reason to be cautious in using the 4PM (Asparouhov & Muthén, 2020). To investigate this, the five models (Dyad-4PNO, Set-4PNO, 4PNO, 3PNO, and 2PNO) were also estimated using several less informative priors than in the previous section. The priors were constructed by crossing two factors: (1) a-parameter prior ( $μ_{a} = 0$ and $σ_{a}^{2} = 6$ vs. $μ_{a} = - 4$ and $σ_{a}^{2} = 15$ ) and (2) b-parameter prior ( $σ_{b}^{2} = 2$ vs. $σ_{b}^{2} = 3$ ). Recall that the a-parameter prior is a normal distribution truncated below at zero. The choice of a-parameter prior levels was made, so that $E (a) \approx 2$ and $Var (a) \approx 2$ or 3, respectively. The joint prior for the slipping and guessing parameters was already uninformative, so it was not adjusted. Also, for model identifiability purposes, the prior mean and variance of the latent trait distribution was maintained at $μ_{θ} = 0$ and $σ_{θ}^{2} = 1$ .

Two measures are used to investigate sensitivity. The first is to compute marginal DIC for each of the models being compared and see which model corresponds with the smallest marginal DIC in each case. Model selection is then not sensitive to the choice of prior to the extent that the same model is chosen in each case. The second measure of interest is the root mean squared difference between of the expected item set proportion correct S_d given the prior in the previous section (i.e., $μ_{a} = 2$ , $σ_{a}^{2} = 1$ , and $σ_{b}^{2} = 1$ ) and an alternative prior. This is done as a simple approach to evaluating the closeness of the model predictions made for each item set under alternative priors. Proportion correct $p_{c} (θ) = E (S_{d} | θ) / n_{d}$ , where n_d is the number of items in item set d, is used, so that the measure is comparable across item sets of varying sizes. To get a total measure across the scale of $θ$ , this difference is integrated over $θ$ . Thus, the root mean squared difference in expected proportion correct (RMSDPC) is

RMSDPC = \sqrt{\int {(p_{c}^{(1)} (θ) - p_{c}^{(2)} (θ))}^{2} f (θ) d θ},

where $p_{c}^{(m)} (θ)$ , $m = 1, 2$ , is the expected item set proportion correct for the two models being compared and $f (θ)$ is the standard normal density.

First, as shown in Table 10, the choice of prior does not seem to change the decision of which of the five models to select. However, decreasing the prior information generally seems to decrease the value of the marginal DIC for each model. This might be expected as the search space for the parameter values would be widened with a larger prior variance. However, the decrease in each case is generally small, though the decrease for the 4PNO is a bit larger (i.e., a decrease of less than 4 for all models except 4PNO).

Table 10.

Marginal Deviance Information Criterion for the Dyad-4PNO, Set-4PNO, 4PNO, 3PNO, and 2PNO Models Under Several Alternative Priors

	$μ_{a} = 0$ and $σ_{a}^{2} = 6$		$μ_{a} = - 4$ and $σ_{a}^{2} = 15$
Model	$σ_{b}^{2} = 2$	$σ_{b}^{2} = 3$	$σ_{b}^{2} = 2$	$σ_{b}^{2} = 3$
Dyad-4PNO	8,658.22 (54.40)	8,655.37 (53.79)	8,655.82 (53.92)	8,654.89 (53.73)
Set-4PNO	8,648.90 (51.51)	8,648.96 (51.69)	8,649.16 (52.06)	8,647.66 (51.62)
PNO	8,683.50 (53.56)	8,680.11 (52.95)	8,682.77 (52.66)	8,681.37 (54.79)
PNO	8,754.40 (48.97)	8,754.47 (48.72)	8,754.06 (49.19)	8,754.48 (49.34)
PNO	8,749.53 (39.93)	8,749.60 (39.96)	8,749.41 (39.96)	8,749.48 (39.99)

Note. Bold-face indicates the model chosen by the given prior. 4PNO = four-parameter normal ogive; 3PNO = three-parameter normal ogive; 2PNO = two-parameter normal ogive.

Next, Table 11 shows the RMSDPC of each item set for each of the alternative priors compared to the prior given in the previous section. In every case, this is nearly zero, indicating a close fit in the expected scores within each item set under the Set-4PNO model regardless of prior. This can further be seen in Figure 5, which shows the expected item set scores for each prior. These results, along with the results of the marginal DIC values, show that the choice of prior does not have much effect on the outcome of results, given that the prior is suitably wide.

Table 11.

Squared Deviation of the Expected Proportion Correct for Each Item Set for Alternative Priors From the Prior $μ_{a} = 2$ , $σ_{a}^{2} = 1$ , and $σ_{b}^{2} = 1$ Under the Set-4PNO Model

Set	$μ_{a} = 0$ and $σ_{a}^{2} = 6$		$μ_{a} = - 4$ and $σ_{a}^{2} = 15$
	$σ_{b}^{2} = 2$	$σ_{b}^{2} = 3$	$σ_{b}^{2} = 2$	$σ_{b}^{2} = 3$
1	.01	.01	.00	.01
2	.01	.02	.01	.02
3	.00	.01	.00	.00
4	.01	.01	.01	.01
5	.00	.01	.01	.01
6	.02	.03	.03	.03
7	.01	.01	.01	.01
8	.01	.01	.01	.01

Note. 4PNO = four-parameter normal ogive.

Figure 5.

Plots of expected proportion correct for each of the eight item sets and each of the priors.

5.2. Investigating Missingness

One thing noted earlier is that there was an overall missingness rate of 3.25% in the sample. To handle missingness above, it was decided to code missing responses as incorrect responses. However, upon further investigation, a couple observations were made. First, there was a fair amount of heterogeneity in the level of missingness between items, ranging from 0.5% to 13.75% missing, with a median of 2.25% missing. Second, the level of missingness was correlated with the item slipping parameter estimates ${\hat{s}}_{j}$ at a level of $r = .76$ and with the item guessing parameter estimates ${\hat{g}}_{j}$ at a level of $r = - .52$ . This suggests that missingness and the parameter estimates are somewhat confounded.

To investigate this further, it was decided to consider the case when the missing data values were imputed. To do this, the mice package in R was used, using the predictive mean matching imputation method. After imputation, the data were again used to fit the five competing models. As shown in Table 12, the Set-4PNO model is again chosen as the best model. Further analysis of the Set-4PNO model solution shows that estimated slipping parameters decreased by 11.67% on average and estimated guessing parameters grew by 9.46% on average. This suggests that missingness did play a small factor in potential over- and underestimation of slipping and guessing parameters, respectively. This seems to be due to general missingness effects seen in IRT modeling, in this case induced by the data being collected in a voluntary low-stakes testing context. However, future investigations should look at this issue further.

Table 12.

Model Fit for the Dyad-4PNO, Set-4PNO, 4PNO, 3PNO, and 2PNO Models Fit With Imputed Data

Index	Dyad-4PNO	Set-4PNO	4PNO	3PNO	2PNO
Marginal DIC	8,609.02	8,588.65	8,628.98	8,700.12	8692.44
p_d	50.28	43.66	50.60	58.58	48.33

Note. DIC = deviance information criterion; 4PNO = four-parameter normal ogive; 3PNO = three-parameter normal ogive; 2PNO = two-parameter parameter ogive.

6. Discussion

6.1. Contributions

In this article, we offer several contributions to the literature. First, we show that the algorithm presented is able to handle large numbers of items and item sets. Notably, each item set is equivalent to an attribute in the HO-DINA model. Initially, the intention was to compare the newly presented Gibbs sampling algorithm to existing HO-DINA estimation software. Unfortunately, it was quickly found out that the existing software was unable to handle a large number of attributes. While theoretically this number of attributes could be found in real assessments, practical constraints presented by the software can hinder investigation of them. The current algorithm, however, can easily estimate a model with 20 or more item sets or attributes. The only constraint on this is that test structure should be in simple structure. While a seemingly restrictive constraint, we consider this to be a positive in that only knowledge of similarity amongst sets of items is necessary for this to hold approximately.

Second, the algorithm presented is relatively fast. Specifically, for 20 item sets of five items (100 items overall) and 2,000 persons, estimation took approximately 30 minutes on a PC with an Intel i7-6700HQ processor with 16 GB of RAM. In comparison, an attempt to use the GDINA package in R to estimate the model with the same data took 12 hours for one iteration before stalling on the same PC. Two reasons why this would be the case are posited here. First, because of simple structure, the model presented has a simplified representation of the likelihood. Specifically, as shown in Equation 5, it is possible to form the joint response function for the items in item set d marginalized over the attribute $W_{i j}$ in a very simplified form. This allows for quick computation of the full likelihood over all responses for all persons. Second, the algorithm uses an efficient Bayesian formulation rather than using a marginal maximum likelihood estimation (MMLE) approach. It is well-known that the performance of MMLE can suffer as dimensionality increases. Thus, as the number of attributes greatly increases, MMLE may be inefficient compared to the Gibbs sampling approach for this model. With this in mind, while beyond the scope of the current article, it is possible that a stochastic approximation method such as Metropolis-Hastings Robbins-Monro (Cai, 2010; Gu & Kong, 1998) may be faster with that same estimation task. Future research could investigate this possibility.

Third, the current model extends the Dyad-4PNO by allowing for confirmatory settings. In a confirmatory setting, a researcher would be able to specify a priori an item structure as a set of items and fit the implied model to the data. This grouping could be obtained through prior exploration or expert knowledge. Furthermore, if several competing structures have been obtained a priori, then each of the implied models can be fit to the data, which can then be compared via fit indices such as the marginal DIC. This is particularly powerful when moving from an exploratory stage to a theory testing stage.

Fourth, the current model further extends the Dyad-4PNO by allowing item sets that are larger than two. The Dyad-4PNO model was developed with minimally sufficient conditions for identification in mind. While this is incredibly theoretically useful, it may be quite limited in many situations. Identification conditions, however, allow for the possibility that items can be grouped into sets of more than two items. For instance, while two items per item set (i.e., a dyad) is sufficient for generic identifiability to hold, three items per item set are sufficient for strict identifiability to hold (Kern & Culpepper, 2020). If there are any more than three items per item set, then strict identifiability will hold. Noting that an item set is equivalent to items matching on a binary attribute, it seems reasonable that item sets of more than two items may occur in practice. The Set-4PNO model can accommodate these situations, and as shown in the simulation study, evidence shows that increasing numbers of items in an item set can improve parameter estimation.

Fifth, it was shown that the model can be applied to real data settings. To do so, a subset of an experimental IQ data set was used. A comparison of five competing models (i.e., 2PNO, 3PNO, 4PNO, Dyad-4PNO, and Set-4PNO) showed that the models with a grouped item structure (Dyad-4PNO and Set-4PNO) fit better than the standard IRT models (2PNO, 3PNO, and 4PNO). This is unsurprising given previous findings using these data (Kern & Culpepper, 2020). However, it was also shown that allowing an item set structure with some sets of items with more than two items rather than simply forcing all sets to have two items, like in the Dyad-4PNO model, can increase model fit. Moreover, this model fit increase is not simply a byproduct of the kind of grouping allowed but rather of the specific sets of items chosen to group together. In other words, knowledge of the sets of items to be grouped together can go a long way in terms of maximizing model fit. Further investigations with these data showed that the conclusion that the Set-4PNO model fit best among the competing models was not due to choice of prior, nor was it due to a decision to rescore missing data as incorrect. One thing that was shown, however, was that the estimation of guessing and slipping parameters may be affected by missingness in the data, though the affect of missingness was ultimately small.

6.2. Future Directions and Limitations

While the Set-4PNO model is an important extension onto the Dyad-4PNO, several future directions can be identified. First, something not explored is the estimation of attribute mastery patterns as in a traditional CDM. This has not been pursued due to the focus on the continuous latent trait, but given the growing interest in the psychometric literature on formative assessments and attribute mastery, future study should investigate the potential of estimating attribute patterns. Importantly, it seems that given the capability of the model to efficiently accommodate large numbers of attributes, this may open up possibilities for more fine-grained assessments.

Second, the model and its accompanying algorithm should be expanded in several ways. One important inclusion would be to consider the possibility of a multidimensional structure of the continuous latent trait. With the model as currently posed with simple structure on the submodel for the attributes, this should be fairly straightforward, though the efficacy of latent trait and parameter estimation would have to be examined, and conditions for identification would need to be established. While noting that criticism of the simple structure requirement of this model, another consideration is determining the extent to which attribute simple structure can be relaxed while still maintaining estimation accuracy and algorithmic speed.

One possible way to expand as briefly discussed is to allow partial information on the item-attribute structure. This would result in a partially exploratory Set-4PNO model. Conceptually, this partially exploratory approach could be used as a way of implementing a Bayesian approach to comparing models. For example, suppose there are two competing models $τ_{1}$ and $τ_{2}$ and a prior is specified, such that $P (τ = τ_{2}) = 1 - P (τ = τ_{1})$ . Then, after convergence of a corresponding MCMC algorithm, subsequent samples from the full posterior distribution would inform on the value of $P (τ | \cdot)$ , the posterior distribution of $τ$ . If $P (τ = τ_{1} | \cdot) > P (τ = τ_{2} | \cdot)$ , then the model specified by $τ_{1}$ is preferred.

Another possible way to expand is to allow for larger item sets in an exploratory approach. Currently, the original exploratory Dyad-4PNO allows for only simple two-item item sets. As shown, this is the minimum necessary for guaranteed identification of slipping and guessing parameters. However, as the current article shows in the example, it is possible to find better fitting solutions with item sets larger than two. The approach taken to finding this solution was to use previous findings from the exploratory Dyad-4PNO and an observation of less certainty regarding two of the dyads. While this approach, using data-based observation worked well, it seems possible that other groupings among the items could fit better yet, which a more general exploratory Set-4PNO approach may be able to discover.

Finally, as suggested by the example, it would be fruitful to further investigate the effects of missingness on the estimation of guessing and slipping parameters. It does seem that level of missingness may have a small overall effect on these parameter estimates. Further work could be done to modify the estimation routine to allow for missingness similarly to how it is currently handled in standard IRT model estimation. In this instance, special attention must be given to estimating the item set parameters.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Justin L. Kern

References

Aitkin

(2006). Investigation of the identifiability of the 3PL model in the NAEP 1986 math survey (Technical Report). National Center for Educational Statistics.

Asparouhov

Muthén

(2020). IRT in Mplus (Technical Report). http://www.statmodel.com/download/MplusIRT. pdf (Version 4).

Barton

M. A.

Lord

F. M.

(1981). An upper asymptote for the three-parameter logistic item-response model (Technical Report No. 80-20). Educational Testing Service.

Cai

(2010). High-dimensional exploratory item factor analysis by a metropolis-hastings Robbins-Monro algorithm. Psychometrika, 75(1), 33–57.

Chen

Culpepper

Liang

(2020). A sparse latent class model for cognitive diagnosis. Psychometrika, 85, 1–33.

Culpepper

S. A.

(2015). Bayesian estimation of the DINA model with Gibbs sampling. Journal of Educational and Behavioral Statistics, 40(5), 454–476.

Culpepper

S. A.

(2016). Revisiting the 4-parameter item response model: Bayesian estimation and application. Psychometrika, 81(4), 1142–1163.

Culpepper

S. A.

(2017). The prevalence and implications of slipping on low-stakes large-scale assessments. Journal of Educational and Behavioral Statistics, 42, 706–725.

de la Torre

Douglas

J. A.

(2004). Higher-order latent trait models for cognitive diagnosis. Psychometrika, 69(3), 333–353.

10.

Feuerstahler

L. M.

Waller

N. G.

(2014). Abstract: Estimation of the 4-parameter model with marginal maximum likelihood. Multivariate Behavioral Research, 49(3), 285.

11.

Zhang

Y.-H.

Shi

Tao

(2021). A Gibbs sampler for the multidimensional four-parameter logistic response model via a data augmentation scheme. British Journal of Mathematical and Statistical Psychology, 74(3), 427–464.

12.

Gabrielsen

(1978). Consistency and identifiability. Journal of Econometrics, 8(2), 261–263.

13.

Gelman

Rubin

D. B.

(1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472.

14.

Gregory

C. A.

(2020). Are we underestimating food insecurity? Partial identification with a Bayesian 4-parameter IRT model. Journal of Classification, 37(3), 632–655.

15.

M. G.

Kong

F. H.

(1998). A stochastic approximation algorithm with Markov chain Monte-Carlo method for incomplete data estimation problems. Proceedings of the National Academy of Sciences, 95(13), 7270–7274.

16.

(2021). Sufficient and necessary conditions for the identifiability of the Q-matrix. Statistica Sinica, 31, 449–472. https://doi.org/10.5705/ss.202018.0420

17.

Junker

B. W.

Sijtsma

(2001). Cognitive assessment models with few assumptions, and connections with nonparametric item response theory. Applied Psychological Measurement, 25(3), 258–272.

18.

Kern

J. L.

Culpepper

S. A.

(2020). A restricted four-parameter IRT model: The dyad four-parameter normal ogive (Dyad-4PNO) model. Psychometrika, 85(3), 575–599.

19.

Loken

Rulison

K. L.

(2010). Estimation of a four-parameter item response theory model. British Journal of Mathematical and Statistical Psychology, 63, 509–525.

20.

Maris

Bechger

(2009). On interpreting the model parameters for the three parameter logistic model. Measurement, 7(2), 75–88.

21.

Meng

Zhang

Tao

(2020). Marginalized maximum a posteriori estimation for the four-parameter logistic model under a mixture modelling framework. British Journal of Mathematical and Statistical Psychology, 73(S1), 51–82. https://doi.org/10.1111/bmsp.12185.

22.

Merkle

E. C.

Furr

Rabe-Hesketh

(2019). Bayesian comparison of latent variable models: Conditional versus marginal likelihoods. Psychometrika, 84(3), 802–829.

23.

Rulison

K. L.

Loken

(2009). I’ve fallen and I can’t get up: Can high-ability students recover from early mistakes in CAT? Applied Psychological Measurement, 33, 83–101.

24.

San Martín

González

Tuerlinckx

(2009). Identified parameters, parameters of interest and their relationships. Measurement, 7(2), 97–105.

25.

San Martín

González

Tuerlinckx

(2015). On the unidentifiability of the fixed-effects 3PL model. Psychometrika, 80(2), 450–467.

26.

Thissen

(2009). On interpreting the parameters for any item response model. Measurement, 7(2), 106–110.

27.

Waller

N. G.

Feuerstahler

L. M.

(2017). Bayesian modal estimation of the four-parameter item response model in real, realistic, and idealized data sets. Multivariate Behavioral Research, 52, 350–370.

28.

Waller

N. G.

Reise

S. P.

(2010). Measuring psychopathology with non-standard item response theory models: Fitting the four-parameter model to the Minnesota Multiphasic Personality Inventory. In Embretson

S. E.

(Ed.), Measuring psychological constructs: Advances in model-based approaches (pp. 147–173). American Psychological Association.

29.

Zhang

(2020). Gibbs-slice sampling algorithm for estimating the four-parameter logistic model. Frontiers in Psychology, 11, 1–18. https://doi.org/10.3389/fpsyg.2020.02121