Maximum-Likelihood Estimation of Noncompensatory IRT Models With the MH-RM Algorithm

Abstract

In “compensatory” multidimensional item response theory (IRT) models, latent ability scores are typically assumed to be independent and combine additively to influence the probability of responding to an item correctly. However, testing situations arise where modeling an additive relationship between latent abilities is not appropriate or desired. In these situations, “noncompensatory” models may be better suited to handle this phenomenon. Unfortunately, relatively few estimation studies have been conducted using these types of models and effective estimation of the parameters by maximum-likelihood has not been well established. In this article, the authors demonstrate how noncompensatory models may be estimated with a Metropolis–Hastings Robbins–Monro hybrid (MH-RM) algorithm and perform a computer simulation study to determine how effective this algorithm is at recovering population parameters. Results suggest that although the parameters are not recovered accurately in general, the empirical fit was consistently better than a competing product-constructed IRT model and latent ability scores were also more accurately recovered.

Keywords

multidimensional item response theory noncompensatory models stochastic estimation

Item response theory (IRT) is a family of statistical models that relates unobservable psychological traits (known as ability in educational testing literature) to observable test behaviors such as responding to a question on an aptitude test. IRT attempts to model how an individual’s observed response pattern can be explained by how the participant’s underlying latent trait level ( $θ$ ) interacts with various item-level characteristics such as “difficulty” or the probability of randomly guessing the correct answer. The majority of the IRT models in popular use are unidimensional, in that they focus on only one latent trait parameter in the test. However, multidimensional IRT models are becoming more popular due to increased computational power and software support on personal computers, as well as their added flexibility for capturing more complex testing phenomenon (Edwards, 2010; Reckase, 2009).

IRT models are best known for their application to dichotomous item response data (e.g., correct versus incorrect) drawn from ability and aptitude tests (Embretson & Reise, 2000). Historically, these models were estimated using a normal ogive response curve function, but software developers have often preferred to use a nearly identical logistic response curve with a scaling modification of 1.702 instead because of their simplicity to manipulate and program (Thissen & Steinberg, 2009). The classic three-parameter logistic model (3PL) introduced by Birnbaum (1968) for person $i$ responding to a dichotomous item $j$ is

P (X_{i j} = 1 | θ_{i}, α_{j}, β_{j}, g_{j}) = g_{j} + \frac{(1 - g_{j})}{1 + \exp (- α_{j} (θ_{i} - β_{j}))},

where $X_{ij} = 1$ represents a positive endorsement to the item (otherwise $X_{ij} = 0$ ), $θ_{i}$ is a single “ability” score sampled from the population of interest, $α_{j}$ is the discrimination parameter, $β_{j}$ is the item difficulty, and $g_{j}$ is the pseudo-guessing parameter. It is often convenient to re-express Birnbaum’s (1968) 3PL formulation into a manner that emphasizes the implicit “slope and intercept” relationship between the parameters; hence, an equivalent expression is

P (X_{i j} = 1 | θ_{i}, a_{j}, d_{j}, g_{j}) = g_{j} + \frac{(1 - g_{j})}{1 + \exp (- (a_{j} θ_{i} + d_{j}))},

where the slopes $a_{j} = α_{j}$ and the intercept $d_{j} = - α_{j} β_{j}$ . Equation 2 demonstrates the intimate connection between linear factor analysis and IRT more clearly. The linear factor model can be seen in the denominator of Equation 2; however, the outcome of the function is now constrained to fall between $g$ and 1 rather than from $- \infty$ to $\infty$ (McDonald, 1999). Finally, this formulation demonstrates the relationship between unidimensional and multidimensional item response theory (MIRT) models more clearly (Reckase, 2009), as explained later. For notational brevity, assume that $θ = θ_{i}$ for the remaining sections when $θ$ contains no subscripts, $θ = (θ_{(1) i}, θ_{(2) i}, \dots, θ_{(K) i})$ represents a $K$ -dimensional vector of ability parameters for individual $i$ , and item-level parameters $a = a_{j}$ , $d = d_{j}$ , and $g = g_{j}$ refer to item $j$ .

Compensatory and Noncompensatory Models

In general, there are two classes of MIRT models: compensatory (e.g., Bock & Aitkin, 1981), which are typically additive models that contain no latent variable interactions, and noncompensatory (e.g., Whitely, 1980) for incorporating products of individual item response probability curves. The former has received an overwhelming amount of attention in IRT literature due to its intimate connection to nonlinear factor analysis (see Lord, 1980). Compensatory models may be most appropriate when test items follow a disjunctive component processes, meaning that the abilities or factors can be combined additively to influence an item response probability, whereas noncompensatory models may be more appropriate for items that have a conjunctive component processes (Maris, 1999). A compensatory model for Equation 2 with $K$ -latent factors can be expressed as

P (X_{i j} = 1 | θ, a, d, g) = g + \frac{(1 - g)}{1 + \exp [- (a^{'} θ + d)]} .

Here, we see that $a = a_{1}, \dots, a_{K}$ and $θ$ are now vectors that linearly combine to form a scalar value. The inspiration for naming this a compensatory model is that a linear combination of $a' θ$ can yield the same scalar value as other $a' θ$ combinations, creating identical expected probability values on the response surface. For example, if $a_{1} = a_{2} = 1$ , and $d = 0$ , then a subject with ability scores $θ_{1} = θ_{2} = 0$ will have exactly the same probability of endorsement as a subject with scores $θ_{1} = - 3$ and $θ_{2} = 3$ . Note that an identical relationship holds for linear factor analysis (Mulaik, 2010) and is intimately related to the well-known problem of rotational indeterminacy.

The compensatory nature of Equation 3 raises theoretical issues for testing situations that require simultaneous contribution from multiple abilities, but where the probability of correct endorsement should be limited by the lowest ability. For compensatory models, if a subject has one very low factor score, she can still have a very high probability of correct response if the remaining factor scores are sufficiently high. Sympson (1977) argued that this hypothesized compensation is not realistic for some types of test items and gives an example of a mathematics test item that requires both arithmetic computation and reading skills. Sympson believed that a subject with poor reading skills would not be able to determine the problem that needed to be solved, and, conversely, that a subject with poor mathematical skills would not be able to solve the problem despite adequately understanding the question. To correct this inherent limitation, Sympson proposed that the probability of correct endorsement for an item be modeled as the product of the individual response curves from separate unidimensional models. The product of these models results in a bounded nonlinear response surface where the probability of correct response cannot exceed the lowest of the unidimensional models, and is in this sense noncompensatory with respect to increasing values of $θ_{k}$ .

More formally, the noncompensatory version of Equation 2 can be expressed as

P (X_{i j} = 1 | θ, a, d, g) = g + (1 - g) \prod_{k = 1}^{K} \frac{1}{1 + \exp [- (a_{k} θ_{k} + d_{k})]},

where $k$ represents the $k$ th element of each parameter vector. Again, if $g = 0$ , the equation reduces to a noncompensatory version of the 2PL model, and if $a = (a_{1} = a_{2} = \dots = a_{K}) = 1$ the equation further reduces the response curve to a noncompensatory version of 1PL model. Henceforth, multidimensional versions of the 1-, 2-, and 3PL models that follow a compensatory structure will be prefixed with a “C,” while noncompensatory models will be prefixed with an “N.”

To determine the observed data likelihood function for compensatory and noncompensatory models, let $P_{ij} = P (X_{ij} = 1 | Ψ_{j}, θ_{i})$ and $Q_{ij} = 1 - P_{ij}$ , where $Ψ_{j}$ is the collection of the item parameters applicable to item $j$ , and let $Ψ$ be the collection of all $Ψ_{j}$ sets and group-level parameters. We also require a fundamental assumption that is necessary for all IRT models known as local independence (Lazarsfeld, 1950). Local independence occurs when $P (X_{ij} = 1)$ is completely specified by the item and ability parameters, and, after controlling for the $θ$ parameters, $P (X_{ij} = 1 | Ψ_{j}, θ_{i})$ is strictly independent of the remaining response vector. When local independence is assumed, the likelihood of the $i$ th ( $J \times 1$ ) response vector, $x_{i}$ , is

L_{ℓ} (X_{i} | Ψ, θ_{i}) = \prod_{j = 1}^{J} P_{i j}^{χ_{i j}} Q_{i j}^{(1 - χ_{i j})},

Here, $χ_{ij}$ is equal to $1$ when the $j$ th item is positively endorsed by individual $i$ , else it is $0$ . Assuming a population distributional form $g (θ)$ —most often a multivariate normal—the marginal distribution becomes

L_{ℓ} (X_{i} | Ψ) = \int_{- \infty}^{\infty} \dots \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} L_{ℓ} (X_{i} | Ψ, θ) g (θ) d θ,

where there are $K$ -fold integrals to be evaluated. Finally, this brings us to the observed data likelihood function. Letting $X$ represent the complete $N \times J$ data matrix, the observed likelihood equation is

L (X | Ψ) = \prod_{i = 1}^{N} [\int_{- \infty}^{\infty} \dots \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} L_{ℓ} (X_{i} | Ψ, θ) g (θ) d θ] .

The IRT parameter estimation for exploratory compensatory models has progressed over that past 60 years, moving from heuristic estimation techniques to more computationally intensive Bayesian Markov chain Monte Carlo (MCMC) methods (Baker & Kim, 2004). The early focus was on estimating the item-specific parameters for unidimensional models, and until Bock and Aitkin (1981) introduced an Expectation-Maximization estimation solution, IRT applications were largely limited to small testing situations (Baker & Kim, 2004). The Expectation-Maximization (EM) algorithm appeared to be a reasonable solution for lower dimensional models without forfeiting numerical accuracy. Unfortunately, this technique quickly becomes inefficient as the number of dimensions increases because the number of quadrature points required for estimating the “E-step” increases exponentially and must be accommodated by decreasing the number of quadrature nodes. A solution for a moderate number of dimensions was described by Schilling and Bock (2005), who demonstrated that adaptive quadrature could be used for better accuracy when a smaller number of quadratures per dimension is used. However, the problem of high-dimensional solutions still remained, and without supplementary EM techniques parameter standard error computations from the parameter information matrix were not possible.

Although interesting in their own right, noncompensatory models have received relatively little attention in theoretical and applied research mainly because they present a greater estimation challenge, especially in exploratory applications. Knol and Berger (1991) state that “the disadvantage of noncompensatory models is that no efficient algorithms for estimation of the item parameters are available,” and to date this statement largely remains true. Since noncompensatory models often include factor-specific difficulty parameters, their estimation requires sufficient variability in the relative difficulties of factors across items to identify the dimensions (Bolt & Lall, 2003). When Sympson (1977) first introduced noncompensatory models, he attempted to estimate a noncompensatory two-parameter logistic (N2PL) model by way of a heuristic regression method. At the same conference proceeding, Lord (1977) suggested that a maximum-likelihood or Bayesian estimation framework would be more appropriate for obtaining model parameters. Shortly thereafter, Whitely (1980) estimated the item coefficients for a N1PL model by maximum likelihood (ML), but this particular implementation required that the factor scores were known a priori and hence error-free.

The first attempt at estimating the item parameters without prior knowledge of the factor scores was introduced by Maris (1995). Maris derived a feasible EM solution for a N1PL model termed a “conjunctive Rasch model.” Unfortunately, the standard errors for the estimated parameters were not available for this method, the factors were required to be orthogonal, and the method itself posed difficulties in its implementation. The main problem was that the EM algorithm was not used to marginalize missing factor scores, as is typical in IRT estimation (e.g., Bock and Aitkin, 1981) and was instead used to deal with the missing latent response processes (Maris, 1995). Consequently, this approach has made the integration of the compensatory and noncompensatory frameworks difficult.

Bolt and Lall (2003) were the first to demonstrate that a feasible and flexible estimation approach for noncompensatory models was to use Bayesian MCMC methods. In their study, the C2PL model was compared with the N1PL model using simulation data and an empirical data set from an English usage test. Whereas the C2PL model was estimated accurately in many scenarios, the N1PL suffered when latent correlations were greater than .6, often requiring larger sample sizes to estimate the parameters to within a reasonable tolerance (n > 3,000). Also, the standard errors for the N1PL model were much larger than the compensatory model for nearly all conditions, and the estimation times took approximately 12 hr to complete on a 1.7 GHz processor for 31 items. Babcock (2011) is the most recent to use MCMC for estimating noncompensatory models and also was the first to attempt to estimate an N2PL model. His simulation results again revealed that higher correlations between the factors led to overall less effective parameter recovery, that at least six unidimensional items per factor were needed to help stabilize the model, and that large sample sizes (e.g., n = 4,000) were necessary for adequate parameter estimation. As Babcock demonstrated, a Bayesian implementation for an N2PL model required substantial modifications to typical prior parameter distributions, and even ad hoc stochastic acceptance schemes were required to control the factor correlation parameter (determined by many prior trail-and-error runs).

Noncompensatory IRT Models With the MH-RM Algorithm

Only recently have confirmatory IRT methods arose in the literature, mainly because programming effective confirmatory methods is difficult due to the inherent “curse of dimensionality,” but also because appropriate standard errors for estimated parameters have not been readily available. Confirmatory IRT models imply that certain restrictions are placed on the model a priori and subsequently tested for their fit to data, such as restricting certain item slopes to equal zero. Estimation methods have mainly been drawn from Bayesian MCMC techniques as they are very flexible and can readily handle user-declared restrictions (see Edwards, 2010).

Recently, however, Cai (2010a) introduced a flexible framework for parameter estimation by combining the Metropolis–Hastings (MH; Hastings, 1970; Metropolis, Rosenbluth, Teller, & Teller, 1953) and Robbins–Monro (RM; Robbins & Monro, 1951) algorithms to form a joint estimation framework that circumvents many of the less attractive features of strictly Bayesian MCMC or ML approaches. The Metropolis–Hastings Robbins–Monro (MH-RM) algorithm estimates the item- and group-level parameters by using a stochastically imputed complete-data solution with an assumed population distributional form (typically, multivariate normal) to capitalize on the more manageable complete-data likelihood:

L (X | Ψ, Θ) = \prod_{i = 1}^{N} L_{ℓ} (X_{i} | Ψ, θ_{i}) g (θ_{i} | μ, Σ) .

For exploratory item factor analysis, the population mean vector $μ$ is usually assumed to be a fixed $K \times 1$ vector of zeros and $∑$ is assumed to be a fixed $K \times K$ identity matrix. However, in confirmatory item analysis, various elements in $μ$ and $∑$ may be estimated in a way analogous to confirmatory factor analysis in the structural equation modeling framework (e.g., see Bollen, 1989, Chapter 7). As can be seen in Equation 8, the complete-data log likelihood is composed of two additive components: the log likelihood for a multivariate logistic regression and the log-likelihood relationships among the factor scores.

Cai (2010a) demonstrated that when using a stochastically imputed complete-data model, and properly accounting for the error in the stochastic imputations, ML estimation and observed parameter standard errors can be calculated. The iterative algorithm works well when partitioned into three stages: perform $b$ burn-in stochastic EM cycles (Diebolt & Ip, 1996), collect $c$ cycles of $Ψ^{(c)}$ and find the average of this set, $\hat{Ψ}$ , and finally perform the MH-RM cycle until the model converges to within a predetermined tolerance. The MH-RM stage begins estimation with the initial parameter estimates $Ψ^{(t = 1)} = \hat{Ψ}$ and recursively completes the following steps:

Given $Ψ^{(t)}$ , stochastically impute an $N \times K$ matrix of missing ability scores $Θ^{(t + 1)}$ with an MH sampler after a small number of burn-ins (say, 5-10).

Using the imputed $Θ^{(t + 1)}$ , find a recursive approximation for the information matrix in place of the noise-corrupted Hessian matrix $H (Ψ | X, Θ^{(t + 1)})$ ,

Γ_{t + 1} = Γ_{t} + γ_{t} (H (Ψ | X, Θ^{(t + 1)}) - Γ_{t}),

where $γ_{t}$ is an RM gain constant (note that $Γ_{1} = 0$ and $γ_{1} = 1$ ).

Update the parameter estimates using the noise-corrected Hessian and gradient vector $\nabla_{t + 1}$ ,

Ψ^{(t + 1)} = Ψ^{(t)} + γ_{t} ({[Γ_{t + 1}]}^{- 1} \nabla_{t + 1}) .

If $| | Ψ^{(t + 1)} - Ψ^{(t)} | | < tolerance$ on $q$ consecutive iterations, then stop. Otherwise, repeat Steps 1 to 3. The $γ_{t}$ parameter is known as a gain constant and slowly decreases to zero as $t \to \infty$ . In practice, $γ_{t}$ can be understood as any reasonable decreasing function so long as it satisfies the following properties:

γ_{t} \in [0, 1], \sum_{t = 1}^{\infty} γ_{t} = \infty, and \sum_{t = 1}^{\infty} γ_{t}^{2} < \infty .

The MH-RM algorithm handles the inherent noise-corrupted complete-data imputations using the RM root-finding algorithm to stabilize both the updates and the information matrix. In this way, the inaccuracies borne from the MH sampler are properly accounted for when attempting to maximize the parameter estimates, and subsequent standard errors can be computed appropriately using Louis’s (1982) complete-data method (see Cai, 2010a). Also, during the MH-RM stage, multiple sets of $θ$ s can be drawn and averaged for constructing the gradient and Hessian matrices, but doing so is rarely required for item factor analysis.

For the noncompensatory model with a lower asymptote parameter, let

P = P (X_{i j} = 1 | θ, a, d, g) = g + (1 - g) \cdot P_{1} (θ, a, d) \cdot P_{2} (θ, a, d) \dots P_{K} (θ, a, d),

be the probability of a positive item endorsement and $Q (X_{i j} = 0 | θ, a, d, g) = 1 - P (X_{i j} = 1 | θ, a, d, g)$ , where

P_{k} = P_{k} (X_{i j} = 1 | θ, a, d) = \frac{1}{1 + \exp (- (a_{k} θ_{k} + d_{k}))},

for $k = 1, 2, \dots, K$ . Assuming that all $θ$ are known, the likelihood function in Equation 7 (or, more appropriately, the $\log_{e}$ values) can be differentiated with respect to each item parameter to obtain the gradient and Hessian functions that are needed to find ML estimates. The analytical evaluation of the noncompensatory model’s gradient and the Hessian functions with and without a lower asymptote parameter can be found in Appendix A. Following the calculations for each noncompensatory item, these values are combined to form a large gradient vector and block-diagonal Hessian matrix and are then used to update all of the item parameters simultaneously using a multivariate Newton–Raphson root finding algorithm, as well as to apply the RM gain constant. The complete-data log likelihood in Equation 8 is also a straightforward generalization when assuming that the item parameters are provisionally known and only requires assuming an appropriate multivariate density function.

Unfortunately, the numerical identification of models that contain noncompensatory items may require slightly more care compared with their compensatory counterparts. While noncompensatory models may be mathematically identified by constraining item and group parameters appropriately (say, constraining the factors to be orthogonal, two slopes to be set to 1, and one additional slope to 0) due to the stochastic nature of how the $\hat{θ}$ values are obtained, it may be beneficial to include sets of unidimensional items that load on only one factor (Babcock, 2011). In this way, the $θ$ values are less likely to switch axes during estimation, and the solution will converge in the correct orientation. The unidimensional items do not necessarily need to be dichotomous to be effective but in the following simulation study they are to keep the models and data less complicated. Strictly speaking, the items need not necessarily be unidimensional to be beneficial in controlling the $θ$ metric (i.e., cross-loadings are permissible), but unidimensional items are used in this simulation study for simplicity and comparability with Babcock’s (2011) simulation study.

The MH-RM algorithm is also suitable when the latent factors are nonlinear or functionally related by products. McDonald (1962) was the first to derive a method to estimate nonlinear factor effects for continuous variables, but modeling categorical variables with nonlinear terms was soon recognized to be problematic due to the “curse of dimensionality” that occurs when many integrals must be evaluated numerically. However, when drawing the $θ$ values instead of using numerical quadrature schemes, this problem quickly diminishes since this is precisely what the stochastic methods were designed to do well: Evaluate high-dimensional and difficult to evaluate integrals. For example, if we wished to model a simple two-factor product function using a C2PL model given the item parameters,

P (X_{i j} = 1 | θ, a, d) = \frac{\exp (a_{1} θ_{1} + a_{2} θ_{2} + a_{3} θ_{1} θ_{2} + d)}{1 + \exp (a_{1} θ_{1} + a_{2} θ_{2} + a_{3} θ_{1} θ_{2} + d)},

then this is accomplished by simply drawing two random $θ_{1}$ and $θ_{2}$ values, multiplying the values together to create the respective product value, and evaluating whether these newly sampled values should be accepted or rejected with an MH sampler when considering the product variable contribution as well. This particular model can have a probability response surface that is similar to a noncompensatory model (see Figure 1 for a visual comparison of noncompensatory and product models), and given that the C2PL model can often model noncompensatory configurations quite well (see Ackerman, 1989, for more details), this interaction-based model may offer subtle but important improvements. However, if the interaction coefficient ( $a_{3}$ ) becomes too large, then the response surface may actually increase as $θ$ values decrease, which may not be theoretically desirable for ability testing applications.

Figure 1.

The probability surface plots on the top represent two-dimensional compensatory (left) and noncompensatory (right) models with the parameters $a_{1} = . 5$ , $a_{2} = 1$ , and $g = . 2$ , where $d = 0$ for the compensatory, and $d_{1} = - 0.5$ and $d_{2} = 0.5$ for the noncompensatory. The bottom figures represent probability surfaces for the product model (Equation 13), where $a_{1} = . 5$ , $a_{2} = . 5$ , and $d = 1$ for both plots, and the interaction coefficients are $a_{3} = . 3$ (left) and $a_{3} = . 1$ (right).

Noncompensatory models offer important and desirable characteristics not inherent in common multidimensional IRT, and should be utilized in situations where the additive property in compensatory models is theoretically inappropriate. However, the estimation of these models clearly has been difficult, with mixed results regarding parameter recovery and ease of implementation. The purpose of this article is to approach the estimation of noncompensatory models by using a flexible MH-RM algorithm, which has several desirable properties when estimating mixed item types as well as multidimensional models. To evaluate the performance of the MH-RM algorithm with noncompensatory models, a simulation study was designed and is described in detail below.

Simulation Study

The simulation design was organized for investigating two- and three-dimensional noncompensatory IRT models that were identified by including several packets of unidimensional models. Noncompensatory items that included non-zero $g$ parameters were also investigated to determine the effects of freely estimating this lower asymptote parameter. The purpose of this simulation study was to determine how effective the MH-RM algorithm could recover the known parameters for compensatory and noncompensatory IRT models, and determine which characteristics were most effective in improving the overall parameter recovery accuracy. In addition, the product model described in Equation 13 was contrasted with the N2PL model, and expected a posteriori (EAP) scores were compared with determine which estimated model could recover the $θ$ values most effectively.

Following Babcock (2011), item response data were simulated for noncompensatory and compensatory items to determine how well parameters could be recovered under known conditions. For simplicity and potentially better stability, the compensatory item types were constrained to load only on one designated factor (i.e., unidimensional), and did not include the lower asymptote parameter, $g$ , so that 2PL models could be used throughout. The design described below was organized with the purpose of creating a “simple structure” factor loadings pattern consisting of compensatory IRT items, and a set of 10-noncompensatory IRT item types estimated with and without the $g$ parameters.

Design and Simulated Parameters

The simulation design had a $2 \times 2$ format, where two- and three-factor models were examined with and without lower asymptote parameters. Within each of these four conditions, there were fixed effects for sample size ( $n = 1, 000, 2, 500, 4, 000$ ), number of unidimensional indicators (5, 10, 15), and factor correlation values (0, .2, .4, .6, .8). Every fixed-factor combination was replicated 30 times for numerical stability. Each cell of the design consisted of 10 noncompensatory items. Furthermore, each fixed-factor combination was segregated by the item class (noncompensatory versus compensatory) since the recovery behavior was expected to be inconsistent based on previous research (e.g., Bolt & Lall, 2003).

The factor scores were drawn from a multivariate normal distribution with mean $0$ and covariance $∑$ ( $MVN ~ (0, ∑)$ ), where the diagonal elements in $∑$ all equaled 1.0 and all the off-diagonals were one of the five previously named fixed values for the inter-factor correlations. The unidimensional item parameters were obtained by drawing $d$ s from a unit normal distribution ( $d ~ N (0, 1)$ ), $a$ s from a log-normal distribution ( $a ~ \ln N (. 2, . 2)$ ), while the $g$ s were fixed at 0. The noncompensatory items were drawn using $d ~ uniform (- 0.5, 3.00)$ , $a ~ \ln N (. 2, . 2)$ , and $g$ s from a beta distribution ( $g ~ β (8, 50)$ ). For the three-factor model, if the $d$ parameter combinations summed to greater than 5.0, they were redrawn to accommodate for the fact that the product of unidimensional probability terms quickly shrinks to zero across all levels of $θ$ . Relatively little is known about the parameters in noncompensatory models, so this parameter set seemed to be a good place to begin since it reflected typical test data characteristics.

Estimation Details

During preliminary analyses, it was noticed that the $d$ parameters for the noncompensatory items occasionally reached unreasonable levels (e.g., $| d | > 6$ ) and had very large standard errors; thus, to better control the range, the $d$ parameters were estimated using a maximum a posteriori (MAP) scheme rather than ML by incorporating a $d ~ N (0, 3)$ prior. In addition, to control the $g$ parameters, a $β (16, 100)$ prior was included so that these parameters could also be recovered using MAP estimation. All IRT models were estimated using the MH-RM algorithm found in the mirt software package (Chalmers, 2012) in R (R Core Team, 2013). The MH-RM controlling parameters were set to have 150 burn-in iterations, 50 iterations to be collected, and averaged over to be used as the starting values for the MH-RM, and a maximum of 2,000 iterations for the MH-RM algorithm. The algorithm was stopped before 2,000 iterations were reached if all of the parameters changed less than 0.001 across three consecutive iterations. In addition, there were five MH draws of the $θ$ parameters per iteration, and the gain constant $(γ_{t})$ used in the RM damper used the function,

γ_{t} = {(\frac{p_{1}}{t - 1})}^{p_{2}},

where the constants $p_{1} = 0.15$ and $p_{2} = 0.65$ for $t = 2, 3, \dots, 2000$ were used for all conditions, with $γ_{1} = 1$ set automatically to initialize the RM damper on the first iteration.

Results

Two statistics were computed to determine how effectively the item parameters were recovered: bias and root mean-squared deviation (RMSD). Bias and RMSD were defined as

Bias = \frac{\sum_{p = 1}^{P} ({\hat{x}}_{p} - x_{p})}{P},

and

RMSD = \sqrt{\frac{\sum_{p = 1}^{P} {({\hat{x}}_{p} - x_{p})}^{2}}{P}},

where $x$ is an arbitrary population parameter, $\hat{x}$ is the respective sample estimate, and $P$ is the number of estimated parameters. To obtain more stable estimates of these statistics, each fixed-factor combination cell was formed by calculating the average of the 30 repeated cells. The marginalized average of these cells with respect to item class (compensatory vs. noncompensatory) and item parameter type can be seen in Table 1.

Table 1.

Compensatory and Noncompensatory Parameter Recovery Statistics.

		Bias			RMSD
Item class	Design	a	d	g	a	d	g
Compensatory	Two-factor, N2PL	0.005	−0.004	—	0.088	0.067	—
	Three-factor, N2PL	0.000	0.016	—	0.093	0.082	—
	Two-factor, N3PL	−0.038	−0.088	—	0.103	0.126	—
	Three-factor, N3PL	−0.007	−0.014	—	0.095	0.078	—
Noncompensatory	Two-factor, N2PL	0.043	−0.029	—	0.327	0.432	—
	Three-factor, N2PL	0.156	0.093	—	0.958	1.131	—
	Two-factor, N3PL	0.555	−0.309	0.059	3.176	2.174	0.121
	Three-factor, N3PL	0.335	−0.449	0.020	2.639	2.251	0.051

Note. N2PL = noncompensatory two-parameter logistic model; RMSD = root mean-squared deviation; N3PL = noncompensatory two-parameter logistic model.

The recovery of the compensatory item parameters on average was adequate, where the RMSD ranged from .067 to .126, and had little bias. Compensatory parameters were recovered with an average RMSD of .077 in the two-factor N2PL design, which is consistent with previous findings (e.g., Cai, 2010a, 2010b; Chalmers, 2012). Noncompensatory item parameters, however, were not recovered well for all cells of the design. The two-factor N2PL had the most acceptable parameter recovery for compensatory and noncompensatory designs, although the RMSD for the noncompensatory models were still greater than 0.3. However, the inter-factor correlation parameters were recovered well for the designs with a non-zero $g$ parameter (average RMSD = 0.045), and very well for the designs without (average RMSD = 0.025). Table 1 makes apparent the increased bias and extreme RMSD values for the noncompensatory two-parameter logistic (N3PL) model compared with the N2PL designs, indicating that the $g$ parameters were highly influential and negatively affected the recovery of the noncompensatory parameters. Given the extremely poor recovery of the N3PL models, and the poor recovery of the three-factor N2PL model, the tables in Appendix B present results demonstrating the effects of varying the number of unidimensional indicators, inter-factor correlation, and sample size for only the two-factor N2PL.

The tables in Appendix B indicate that increasing the sample size reduced the amount of bias and RMSD. For the compensatory parameters, increasing the number of indicators tended to reduce the RMSD, while increasing the inter-factor correlation did not appear to have an influence in reducing RMSD or bias. For the noncompensatory parameters, increasing the inter-factor correlation and number of indicators were highly influential. Using more indicators tended to effectively decrease the RMSD, and higher inter-factor correlations increased the RMSD. In addition, the recovery of the noncompensatory intercepts and slopes appeared to be consistently biased toward zero, possibly due to the normal prior imposed on the $d$ parameters, though this bias tended to decrease as the sample size increased. The compensatory slope parameters were not recovered as effectively as the compensatory intercepts, but for noncompensatory models the opposite trend appeared. However, as can be seen in Table 1, this trend was not consistent when including the $g$ parameter.

The N3PL designs demonstrated more difficulty converging compared with the N2PL designs, often reaching the maximum number of allocated MH-RM cycles. The two- and three-factor N2PL designs on average converged in 1,153.2 ( $SD = 256.3$ ) and 1,831.7 ( $SD = 118.6$ ) iterations, respectively, after an average of $879.3$ ( $SD = 343.0$ ) and 1,938.8 ( $SD = 883.4)$ s. The two- and three-factor N3PL designs, however, converged after an average of 1,876.8 ( $SD = 72.91$ ) and 1,995.5 ( $SD = 6.94$ ) iterations, respectively, after 1,324.1 ( $SD = 644.8$ ) and 2,231.1 ( $SD = 1, 142.7$ ) s. As can be seen from the very small standard deviation, the three-factor N3PL model almost never converged after 2,000 MH-RM cycles, indicating that the model was the most unstable of the four designs.

Noncompensatory and Product Models

The product model from Equation 13 was also estimated to evaluate how effectively this model could fit the simulated N2PL data. Estimation of the product-constructed model appeared to be more stable than the N2PL, converging in an average of 282.9 MH-RM iterations ( $SD = 124.4$ ) compared with 1,153.2 iterations ( $SD = 256.3$ ) for the N2PL model. This in turn cause the product-constructed models to converge in less than half the estimation time compared with the N2PL models.

Given that the N2PL and product model have the same degrees of freedom, only the log-likelihood and $G^{2}$ goodness-of-fit statistic were compared. For these data, the two-factor N2PL model had a higher log-likelihood (average difference of 19.5) and lower $G^{2}$ value (average difference of −39.4) for all fixed-factor combinations, indicating that the N2PL provided a better empirical fit to the simulated data than the product model across all conditions. Increasing the sample size and number of unidimensional indicators tended to increase the log-likelihood differences, while increasing the inter-factor correlation tended to decrease the difference. The same trend also occurred with $G^{2}$ values. It appeared that regardless of the number of indicators, sample size, and inter-factor correlation, the N2PL model better represented the data despite the fact that the estimates of the noncompensatory parameters were fairly discrepant from the population values.

In addition to comparing model statistics, recovery of the $θ$ parameters is also important in determining the effectiveness of the product-constructed model approximation. Therefore, EAP estimates of the $θ$ parameters were estimated for both models. EAP scores were compared with the simulated $Θ$ values, and RMSD statistics were computed for the noncompensatory and product model designs separately. The noncompensatory model recovered the ability parameters with an RMSD = 0.463, while the product model had RMSD = 0.467. A paired t test indicated that the noncompensatory ability parameters were recovered with a lower RMSD across the experimental conditions, $t_{44} = - 10.51$ , $p < . 001,$ $η^{2} = . 715$ , and of the 45 design pairings, the RMSD was lower than the product model in every condition. Therefore, despite the lack of precision in recovering the population parameters, the N2PL managed to outperform the competing product-constructed approximation in recovering the $Θ$ s. However, from a practical perspective this small improvement may not be important in applied settings, and therefore, the product model can still serve as a useful approximation to the noncompensatory model provided that the interaction coefficient is not too large in magnitude.

Discussion

This study examined the accuracy of recovering compensatory and noncompensatory IRT parameters using the MH-RM algorithm. While previous studies have demonstrated the complexities of estimating noncompensatory models using only ML and Bayesian methods, the MH-RM was believed to provide a suitable compromise between the two frameworks since it borrows strength from both disciplines. The algorithm provided a flexible framework for estimating mixed IRT item types for two- and three-dimensional 2PL and 3PL noncompensatory models, and demonstrated some promising results for noncompensatory models that did not freely estimate the lower asymptote parameter.

The results from the simulation study indicated that the N2PL model could provide a better empirical fit to data simulated from the N2PL model compared with a competing product model. In addition, the N2PL model was able to more accurately recover the population ability ( $θ$ ) parameters than was the approximate but more stable product model. These findings suggests that if items theoretically should be modeled from noncompensatory factors, then estimating a noncompensatory IRT model is justified, and, if appropriate, may even be included into computerized adaptive testing procedures. This discovery was somewhat surprising given the elevated RMSD statistics obtained in Tables B1 to B3, and appears to demonstrate the fungibility of item-level parameters when computing incidental estimates such as ability scores.

Unfortunately, there were several concerns for fitting noncompensatory models by ML with the MH-RM algorithm as well. To begin, researchers interested in recovering the population coefficients for noncompensatory models will experience difficulties when estimating all noncompensatory parameters, although the most promising is the two-factor N2PL model. The overall recovery statistics clearly indicate that the population parameters cannot be recovered accurately without substantial effort, and although increasing the sample size and number of unidimensional indicators does help with increasing precision, the sample sizes might have to be in the tens of thousands before adequate confidence in the estimates can be obtained. Also for the noncompensatory model, the parameter standard errors were very large under all conditions, which is consistent with past estimation results (e.g., Babcock, 2011; Bolt & Lall, 2003). Increasing the number of latent factors in the model also decreased the estimation recovery accuracy, suggesting that higher dimensional solutions require even larger data sets, and including lower bound parameters ( $g$ ) caused extremely inaccurate recovery of the noncompensatory slopes and intercepts. In empirical applications where lower bounds should be present it may be beneficial to treat these parameters as fixed at prespecified values and only estimate the slopes, intercepts, and inter-factor correlations. This approach has been adopted when estimating compensatory IRT models with lower bound parameters to avoid instability issues (cf. TESTFACT 4; Wood et al., 2003), and may be even more important for noncompensatory models do to their natural instability.

The size of the factor correlation also played a negative role in the parameter estimation accuracy, although this problem was not as severe as Babcock (2011) noticed when using an MH within Gibbs sampling estimation approach. In fact, the factor correlation estimates were relatively close to the simulated population values, and, unlike the MH within Gibbs sampling approach, did not require ad hoc prior distributions to be imposed to help facilitate stability. Factor correlations greater than .8, however, may create further estimation problems and should be investigated in the future. Another potential issue in noncompensatory IRT models is that different sets of item parameters could produce very similar item characteristic surfaces to compared with the population values, therefore research into checking the closeness of estimated item response surface (IRS) recovery should be investigated (Zhang, 2012).

An interesting area to investigate in future work is the possibility of combining compensatory and noncompensatory models and applying them to uniquely suited items that would not be adequately fitted by either class in isolation. One example of this type of model could be

P (X_{i j} = 1 | θ, a, d, g) = g + \frac{(1 - g)}{(1 + \exp [- (a_{1} θ_{1} + d_{1})]) (1 + \exp [- (a_{k \neq 1}^{’} θ_{k \neq 1} + d_{2})])},

where $K - 1$ factors are modeled to have a compensatory relationship but jointly hold a noncompensatory relationship with the values from $θ_{1}$ . This item type can also easily be included within the aforementioned MH-RM estimation framework, and may be more stable than the N3PL models for three dimensions since a compensatory component is included. Empirically, this model would be ideal for a situation where multiple responses that require different forms of reasoning are required, but still must be limited by another latent factor. For instance, imagine the following mathematics test item:

One of the most important theorems for understanding Euclidean space and geometry is the Pythagorean theorem, $a^{2} + b^{2} = c^{2}$ . The theorem can be generalized to include higher dimensional spaces, to objects that are not right triangles, and to objects that are in fact not triangles at all. Prove that the two-dimensional version of the theorem is true using an algebraic, geometric, or differential method. Show your work.

This question is meant to measure a student’s ability to prove an important theorem using whatever method they are most comfortable with; however, it requires some additional information before the question can be attempted. A potential limiting factor may be the student’s ability to comprehend what the question is asking, for without this important ability the probability of correct endorsement drops rapidly to zero, regardless of how talented the student may be at successfully completing proofs in other contexts.

Although examples such as these can be accommodated by noncompensatory and compensatory combinations, this article showed that accurately recovering population parameters for items that contain a noncompensatory component is very difficult. However, if researchers are more interested in specifying models that optimally fit their data both empirically and theoretically, rather than recovering population parameters accurately, then these item response models may be warranted, and the flexibility of the MH-RM algorithm for compensatory and noncompensatory response models would be an effective estimation approach for estimating these models by full-information ML.

Footnotes

Appendix A

The observed-data log-likelihood equation used to estimate model parameters is of the form

L = C + \sum_{i = 1}^{U} \sum_{j = 1}^{J} (r_{ij} \log (P_{ij}) + (f_{ij} - r_{ij}) \log (Q_{ij})),

where $C$ is a constant, $P_{ij}$ is defined by (Equation 4) with $g = 0$ and is composed of the elements in (Equation 12), $Q_{ij} = 1 - P_{ij}$ , $r_{ij}$ represents the number of subjects who correctly endorsed the item in response pattern $i$ , $f_{ij}$ is the total number of subjects who responded to item $j$ with pattern $i$ , and $U$ is the number of unique response patterns. In the MH-RM algorithm, the summation in $U$ can be taken over the total sample size $(N)$ rather than over the unique response patterns, $f_{ij}$ can be set to a constant of $1$ , and $r_{ij}$ represents whether subject $i$ responded correctly ( $1$ ) or incorrectly ( $0$ ) for item $j$ . For notation clarity, the indices for the unique response patterns (or sample size) and items will be omitted. It is also useful to note a few analytic identities before expressing the gradient and Hessian, where for noncompensatory models without a lower asymptote parameter $\partial P / \partial a_{k} = P Q_{k} θ_{k}$ , $\partial P / \partial d_{k} = P Q_{k}$ , $\partial Q_{k} / \partial a_{k} = - P_{k} Q_{k} θ_{k}$ , and $\partial Q_{k} / \partial d_{k} = - P_{k} Q_{k}$ , where $Q_{k} = 1 - P_{k}$ .

Differentiating the log likelihood with respect to $d_{k}$ and $a_{k}$ for item $j$ , where $k$ is the index for each $1, 2, \dots, K$ dimensions in $θ$ , gives the gradient values

\frac{\partial L}{\partial d_{k}} = \sum Q_{k} [r - (f - r) \frac{P}{Q}],

\frac{\partial L}{\partial a_{k}} = \sum θ_{k} Q_{k} [r - (f - r) \frac{P}{Q}],

which can be collected into the gradient vector

\nabla L' = (\begin{matrix} \frac{\partial L}{\partial d_{1}} & \frac{\partial L}{\partial a_{1}} & \dots & \frac{\partial L}{\partial d_{K}} & \frac{\partial L}{\partial a_{K}} \end{matrix}),

Further differentiating the log likelihood, and using $s$ to index off-diagonal elements across the $K$ dimensions (where $s \neq k$ ) gives the Hessian values

\frac{\partial^{2}}{\partial^{2} d_{k}} = - \sum Q_{k} [P_{k} [r - (f - r) \frac{P}{Q}] + Q_{k}^{2} (f - r) \frac{PQ + P^{2}}{Q^{2}}],

\frac{\partial^{2} L}{\partial^{2} a_{k}} = - ​ \sum θ_{k}^{2} Q_{k} [P_{k} [r - (f - r) \frac{P}{Q}] + Q_{k} (f - r) \frac{P Q + P^{2}}{Q^{2}}],

\frac{\partial^{2} L}{\partial a_{k} \partial d_{k}} = - \sum θ_{k} Q_{k} [P_{k} [r - (f - r) \frac{P}{Q}] + Q_{k} (f - r) \frac{P}{Q^{2}}],

\frac{\partial^{2} L}{\partial d_{k} \partial d_{s}} = - \sum Q_{k} Q_{s} (f - r) \frac{P}{Q^{2}},

\frac{\partial^{2} L}{\partial d_{k} \partial a_{s}} = - \sum θ_{s} Q_{k} Q_{s} (f - r) \frac{P}{Q^{2}},

\frac{\partial^{2} L}{\partial a_{k} \partial a_{s}} = - \sum θ_{k} θ_{s} Q_{k} Q_{s} (f - r) \frac{P}{Q^{2}} .

These values are then collected into the $2 K \times 2 K$ Hessian matrix:

\nabla^{2} L = (\begin{matrix} \frac{\partial^{2} L}{\partial^{2} d_{1}} & \frac{\partial^{2} L}{\partial a_{1} \partial d_{1}} & \dots & \frac{\partial^{2} L}{\partial d_{1} \partial d_{s}} & \frac{\partial^{2} L}{\partial d_{1} \partial a_{s}} \\ \frac{\partial^{2} L}{\partial a_{1} \partial d_{1}} & \frac{\partial^{2} L}{\partial^{2} a_{1}} & \dots & \frac{\partial^{2} L}{\partial a_{1} \partial d_{s}} & \frac{\partial^{2} L}{\partial a_{1} \partial a_{s}} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ \frac{\partial^{2} L}{\partial d_{1} \partial d_{s}} & \frac{\partial^{2} L}{\partial a_{1} \partial d_{s}} & \dots & \frac{\partial^{2} L}{\partial^{2} d_{K}} & \frac{\partial^{2} L}{\partial a_{K} \partial d_{K}} \\ \frac{\partial^{2} L}{\partial d_{1} \partial a_{s}} & \frac{\partial^{2} L}{\partial a_{1} \partial a_{s}} & \dots & \frac{\partial^{2} L}{\partial a_{K} \partial d_{K}} & \frac{\partial^{2} L}{\partial^{2} a_{K}} \end{matrix}) .

For a noncompensatory model with a lower asymptote parameter, let $P = g + (1 - g) P_{1} P_{2} \dots P_{K}$ , $Q = 1 - P$ , and let $P^{*}$ be identical to Equation 12. Some useful derivatives for this model are $\partial P / a_{k} = (1 - g) θ_{k} P^{*} Q_{k}$ , $\partial P / \partial d_{k} = (1 - g) P^{*} Q_{k}$ , and $\partial P / \partial g = Q^{*}$ . The gradient vector is then

\frac{\partial L}{\partial d_{k}} = \sum (1 - g) P^{*} Q_{k} [\frac{r}{P} - \frac{(f - r)}{Q}],

\frac{\partial L}{\partial a_{k}} = \sum (1 - g) θ_{k} P^{*} Q_{k} [\frac{r}{P} - \frac{(f - r)}{Q}],

\frac{\partial L}{\partial g_{k}} = \sum Q^{*} [\frac{r}{P} - \frac{(f - r)}{Q}],

\nabla L' = (\begin{matrix} \frac{\partial L}{\partial d_{1}} & \frac{\partial L}{\partial a_{1}} & \dots & \frac{\partial L}{\partial d_{K}} & \frac{\partial L}{\partial a_{K}} & \frac{\partial L}{\partial g} \end{matrix}) .

Further differentiating the log likelihood gives

\frac{\partial^{2} L}{\partial^{2} d_{k}} = \sum (1 - g) P^{*} Q_{k} [(\frac{r}{P} - \frac{(f - r)}{Q}) [(1 - g) Q_{k} - P_{k}] - P^{*} Q_{k} (1 - g) (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

\frac{\partial^{2} L}{\partial^{2} a_{k}} = \sum (1 - g) θ_{k}^{2} P^{*} Q_{k} [(\frac{r}{P} - \frac{(f - r)}{Q}) [(1 - g) Q_{k} - P_{k}] - P^{*} Q_{k} (1 - g) (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

\frac{\partial^{2} L}{\partial d_{k} \partial a_{k}} = \sum (1 - g) θ_{k} P^{*} Q_{k} [(\frac{r}{P} - \frac{(f - r)}{Q}) [(1 - g) Q_{k} - P_{k}] - P^{*} Q_{k} (1 - g) (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

\frac{\partial^{2} L}{\partial^{2} g} = - \sum Q^{* 2} [\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}}],

\frac{\partial^{2} L}{\partial d_{k} \partial g} = - \sum P^{*} Q_{k} [(\frac{r}{P} - \frac{(f - r)}{Q}) + (1 - g) Q^{*} (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

\frac{\partial^{2} L}{\partial a_{k} \partial g} = - \sum θ_{k} P^{*} Q_{k} [(\frac{r}{P} - \frac{(f - r)}{Q}) + (1 - g) Q^{*} (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

\frac{\partial^{2} L}{\partial d_{k} \partial d_{s}} = \sum (1 - g) P^{*} Q_{k} Q_{s} [(\frac{r}{P} - \frac{(f - r)}{Q}) - (1 - g) P^{*} (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

\frac{\partial^{2} L}{\partial a_{k} \partial a_{s}} = \sum (1 - g) θ_{k} θ_{s} P^{*} Q_{k} Q_{s} [(\frac{r}{P} - \frac{(f - r)}{Q}) - (1 - g) P^{*} (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

\frac{\partial^{2} L}{\partial d_{k} \partial a_{s}} = \sum (1 - g) θ_{s} P^{*} Q_{k} Q_{s} [(\frac{r}{P} - \frac{(f - r)}{Q}) - (1 - g) P^{*} (\frac{r}{P^{2}} + \frac{(f - r)}{Q^{2}})],

which is collected into the $(2 K + 1) \times (2 K + 1)$ Hessian matrix

\nabla^{2} L = (\begin{matrix} \frac{\partial^{2} L}{\partial^{2} d_{1}} & \frac{\partial^{2} L}{\partial a_{1} \partial d_{1}} & \dots & \frac{\partial^{2} L}{\partial d_{1} \partial d_{s}} & \frac{\partial^{2} L}{\partial d_{1} \partial a_{s}} & \frac{\partial^{2} L}{\partial d_{1} \partial g} \\ \frac{\partial^{2} L}{\partial a_{1} \partial d_{1}} & \frac{\partial^{2} L}{\partial^{2} a_{1}} & \dots & \frac{\partial^{2} L}{\partial a_{1} \partial d_{s}} & \frac{\partial^{2} L}{\partial a_{1} \partial a_{s}} & \frac{\partial^{2} L}{\partial a_{1} \partial g} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ \\ \frac{\partial^{2} L}{\partial d_{1} \partial d_{s}} & \frac{\partial^{2} L}{\partial a_{1} \partial d_{s}} & \dots & \frac{\partial^{2} L}{\partial^{2} d_{K}} & \frac{\partial^{2} L}{\partial a_{K} \partial d_{K}} & \frac{\partial^{2} L}{\partial d_{K} \partial g} \\ \frac{\partial^{2} L}{\partial d_{1} \partial a_{s}} & \frac{\partial^{2} L}{\partial a_{1} \partial a_{s}} & \dots & \frac{\partial^{2} L}{\partial a_{K} \partial d_{k}} & \frac{\partial^{2} L}{\partial^{2} a_{K}} & \frac{\partial^{2} L}{\partial a_{K} \partial g} \\ \frac{\partial^{2} L}{\partial d_{1} \partial g} & \frac{\partial^{2} L}{\partial a_{1} \partial g} & \dots & \frac{\partial^{2} L}{\partial d_{K} \partial g} & \frac{\partial^{2} L}{\partial a_{K} \partial g} & \frac{\partial^{2} L}{\partial^{2} g} \end{matrix}) .

Appendix B

Table B3.

N2PL for $n = 4, 000$ .

			Compensatory		Noncompensatory
Parameters	Correlation	No. of indicators	Bias	RMSD	Bias	RMSD
$a$	0.0	5	0.004	0.066	0.009	0.206
		10	0.003	0.061	0.028	0.204
		15	0.005	0.061	0.016	0.162
	0.2	5	0.007	0.068	0.030	0.226
		10	−0.004	0.063	0.026	0.195
		15	0.002	0.062	0.016	0.178
	0.4	5	0.000	0.069	0.039	0.258
		10	0.002	0.061	0.024	0.208
		15	0.001	0.060	0.010	0.187
	0.6	5	0.002	0.063	0.025	0.327
		10	0.004	0.060	0.025	0.246
		15	0.000	0.059	0.015	0.206
	0.8	5	0.002	0.061	0.032	0.389
		10	0.008	0.058	0.042	0.341
		15	0.004	0.058	0.040	0.294
$d$	0.0	5	−0.007	0.051	0.003	0.298
		10	0.002	0.047	0.024	0.306
		15	0.001	0.048	0.013	0.241
	0.2	5	−0.003	0.047	0.022	0.320
		10	0.002	0.048	0.026	0.288
		15	0.002	0.048	0.021	0.269
	0.4	5	0.001	0.046	0.026	0.367
		10	0.002	0.046	0.019	0.299
		15	−0.002	0.049	0.007	0.282
	0.6	5	−0.007	0.046	0.009	0.460
		10	−0.004	0.043	0.012	0.346
		15	0.001	0.047	0.012	0.283
	0.8	5	0.004	0.046	0.019	0.509
		10	−0.001	0.046	0.021	0.437
		15	−0.005	0.049	0.020	0.407

Note. N2PL = noncompensatory two-parameter logistic model; RMSD = root mean-squared deviation.

Acknowledgements

The authors would like to thank four anonymous reviewers for providing insightful feedback, which greatly improved the quality of this manuscript.

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.

References

Ackerman

T. A.

(1989). Unidimensional IRT calibration of compensatory and noncompensatory multidimensional items. Applied Psychological Measurement, 13, 113-127.

Babcock

(2011). Estimating a noncompensatory IRT model using Metropolis within Gibbs sampling. Applied Psychological Measurement, 35, 317-329.

Baker

F. B.

Kim

S. H.

(2004). Item response theory: Parameter estimation techniques (2nd ed.). New York, NY: Marcel Dekker.

Birnbaum

(1968). Some latent trait models and their use in inferring an examinee’s ability. In Lord

F. M.

Novick

M. R.

(Eds.), Statistical theories of mental test scores (pp. 395-479). Reading, MA: Addison-Wesley.

Bock

R. D.

Aitkin

(1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46, 443-459.

Bollen

K. A.

(1989). Structural equations with latent variables. New York, NY: Wiley.

Bolt

D. M.

Lall

V. F.

(2003). Estimation of compensatory and noncompensatory multidimensional item response models using Markov Chain Monte Carlo. Applied Psychological Measurement, 27, 395-414.

Cai

(2010a). High-dimensional exploratory item factor analysis by a Metropolis-Hastings Robbins-Monro algorithm. Psychometrika, 75, 33-57.

Cai

(2010b). Metropolis-Hastings Robbins-Monro algorithm for confirmatory item factor analysis. Journal of Educational and Behavioral Statistics, 35, 307-335.

10.

Chalmers

R. P.

(2012). mirt: A multidimensional item response theory package for the R environment. Journal of Statistical Software, 48(6), 1-29. Retrieved from http://www.jstatsoft.org/v48/i06

11.

Diebolt

E. H. S.

(1996). Stochastic EM: Method and application. In Gilks

Richardson

Spiegel-halter

(Eds.), Markov chain Monte Carlo in practice (pp. 259-273). London, England: Chapman & Hall.

12.

Edwards

M. C.

(2010). A Markov Chain Monte Carlo approach to confirmatory item factor analysis. Psychometrika, 75, 474-497.

13.

Embretson

S. E.

Reise

S. P.

(2000). Item response theory for psychologists. Mahwah, NJ: Erlbaum.

14.

Hastings

W. K.

(1970). Monte Carlo simulation methods using Markov chains and their applications. Biometrika, 57, 97-109.

15.

Knol

D. L.

Berger

M. P. F.

(1991). Empirical comparison between factor analysis and multidimensional item response models. Multivariate Behavioral Research, 26, 457-477.

16.

Lazarsfeld

P. F.

(1950). The logical and mathematical foundation of latent structure analysis. In Stouffer

S. A.

Guttman

Suchman

E. A.

Lazarsfeld

P. F.

Star

S. A.

Clausen

J. A.

(Eds.), Measurement and prediction (pp. 362-412). New York, NY: Wiley.

17.

Lord

F. M.

(1977, July). Discussion: Section 2. In Weiss

D. J.

(Ed.), Proceedings of the 1977 Computerized Adaptive Testing Conference (pp. 99-103). Minneapolis, MN.

18.

Lord

F. M.

(1980). Applications of item response theory to practical testing problems. Hillsdale, NJ: Lawrence Erlbaum.

19.

Louis

T. A.

(1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society: Series B, 44, 226-233.

20.

Maris

(1995). Psychometric latent response models. Psychometrika, 60, 523-547.

21.

Maris

(1999). Estimating multiple classification latent class models. Psychometrika, 64, 187-212.

22.

McDonald

R. P.

(1962). A general approach to nonlinear factor analysis. Psychometrika, 27, 397-415.

23.

McDonald

R. P.

(1999). Test theory: A unified treatment. Mahwah, NJ: Erlbaum.

24.

Metropolis

Rosenbluth

A. W.

Teller

A. H.

Teller

(1953). Equations of state space calculations by fast computing machines. Journal of Chemical Physics, 21, 1087-1091.

25.

Mulaik

S. A.

(2010). Foundations of factor analysis (2nd ed.). Boca Raton, FL: Chapman & Hall.

26.

R Core Team. (2013). R: A language and environment for statistical computing [Computer software]. Vienna, Austria. Available from http://www.R-project.org/

27.

Reckase

M. D.

(2009). Multidimensional item response theory. New York, NY: Springer.

28.

Robbins

Monro

(1951). A stochastic approximation method. The Annals of Mathematical Statistics, 22, 400-407.

29.

Schilling

Bock

R. D.

(2005). High-dimensional maximum marginal likelihood item factor analysis by adaptive quadrature. Psychometrika, 70, 533-555.

30.

Sympson

J. B.

(1977, July). A model for testing with multidimensional items. In Weiss

D. J.

(Ed.), Proceedings of the 1977 Computerized Adaptive Testing Conference (pp. 82-98). Minneapolis, MN.

31.

Thissen

Steinberg

(2009). Item response theory. In Millsap

R. E.

Maydeu-Olivares

(Eds.), The SAGE handbook of quantitative methods in psychology (pp. 148-177). Los Angeles, CA: Sage.

32.

Whitely

S. E.

(1980). Multicomponent latent trait models for ability tests. Psychometrika, 60, 181-198.

33.

Wood

Wilson

Gibbons

Schilling

Muraki

Bock

(2003). TESTFACT 4 for Windows: Test scoring, item statistics, and full-information item factor analysis [Computer software]. Scientific Software International. Chicago, IL: Scientific Software International.

34.

Zhang

(2012). Calibration of response data using MIRT models with simple and mixed structures. Applied Psychological Measurement, 36, 375-398.