Fast Bayesian Estimation for the Four-Parameter Logistic Model (4PLM)

Abstract

There is a rekindled interest in the four-parameter logistic item response model (4PLM) after three decades of neglect among the psychometrics community. Recent breakthroughs in item calibration include the Gibbs sampler specially made for 4PLM and the Bayes modal estimation (BME) method as implemented in the R package mirt. Unfortunately, the MCMC is often time-consuming, while the BME method suffers from instability due to the prior settings. This paper proposes an alternative BME method, the Bayesian Expectation-Maximization-Maximization-Maximization (BE3M) method, which is developed from by combining an augmented variable formulation of the 4PLM and a mixture model conceptualization of the 3PLM. The simulation shows that the BE3M can produce estimates as accurately as the Gibbs sampling method and as fast as the EM algorithm. A real data example is also provided.

The four-parameter logistic item response model (4PLM) was first mentioned as far back as McDonald (1967) and formally proposed by Barton and Lord (1981). Despite this, for many years, it received very little attention in the literature due to doubts about its utility, as well as technical difficulties related to parameter estimation. Recently, however, the psychometrics community has rekindled its interest in the 4PLM after three decades of neglect.

Increasingly, applications of the 4PLM have appeared in the literatures as supporting evidence of its usefulness. Several studies argued that the slipping parameter from the 4PLM can be estimated to account for early careless errors to improve ability estimation (Liao et al., 2012; Rulison & Loken, 2009). In psychology, it has been noted that there is value in modeling the lower and upper asymptotes of item response functions (Reise & Waller, 2003; Waller & Reise, 2010) in instances where subjects with higher levels of psychopathology may be reluctant to self-report attitudes, behaviors, and/or experiences which may indicate a social desirability bias (Rouse et al., 1999). In education, recent studies have shown the prevalence and implication of the slipping effect in low-stakes large-scale educational assessments (Culpepper, 2016, 2017).

Despite this work, one daunting challenge regarding the 4PLM pertains to the difficulty in estimating item parameters. The lack of a fast and stable implementation of the 4PLM remains a major barrier to its widespread application in practice. Few studies to date have addressed the item calibration issue for the 4PLM. All the methods fall into two approaches: MCMC methods and marginal maximum likelihood (Feuerstahler & Waller, 2014). For MCMC methods, Loken and Rulison (2010) presented Bayesian model formulations with Metropolis-Hastings (MH) sampling using the general-purpose software (e.g., “bayesmh” module in Stata). More recently, Culpepper (2016) proposed a flexible Gibbs sampling method. Both MCMC methods can provide accurate item estimates, but are inherently time-consuming. Even for Culpepper’s (2016, 2017) Gibbs sampling approach, calibration is expected to take at least 1 hour, which also motivates us to find a faster Bayesian estimation.

Another recent implementation of the 4PLM was done by Waller and Feuerstahler (2017). In their work, they extended the Bayesian modal estimation method (BME) to the 4PLM using the R package mirt (Chalmers, 2012). Compared to MCMC methods, BME is a fast calibration method, but several issues need to be explored to evaluate the utility of BME and its implementation in the mirt package. For example, Guo and Zheng (2019) found that the item parameter estimates yielded by BME were unstable for the 3PLM when changing priors of guessing parameters. Moreover, in the mirt package of version 1.30, the standard errors of all item parameters cannot be calculated as a consequence of program errors in computing the Hessian matrix when a Beta distribution prior (“expbeta” in mirt) is imposed on guessing or slipping parameters.

At present, we propose an alternative Bayesian modal estimation approach that we term Bayesian Expectation-Maximization-Maximization-Maximization (BE3M) method, which is developed by using a latent variable augmentation approach for the 4PLM in conjunction with the mixture modeling approach for the 3PLM (Zheng et al., 2017). The remainder of the manuscript is as follows: First, we present a detailed derivation of the BME and discuss potential difficulties in item calibration. Second, we present a latent variable augmentation approach to formulating the 4PLM (Béguin & Glas, 2001; Culpepper, 2016), and the resulting BE3M estimation method. Third, we present two simulation studies assessing the accuracy of parameter recovery for the four calibration methods (BE3M, mirt, fourPNO, Stata). Fourth, we apply the four methods to bullying item responses of 7,491 adolescents from the 2005 to 2006 Health Behavior in School-Aged Children (HBSC) study. Lastly, we discuss future research directions and concluding remarks.

Bayesian Modal Estimation for the 4PLM

Following Culpepper (2016), the slope-intercept form IRT formulation of the 4PLM is:

$$P\left(Y_{ij}=1|{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)\equiv P_j\left({\mathrm{\theta }}_i\right)={\mathrm{\gamma }}_j+\left(1-{\mathrm{\varsigma }}_j-{\mathrm{\gamma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right),$$

(1)

with

$$P_j^{*}\left({\mathrm{\theta }}_i\right)=\frac{1}{1+\exp \left(-D\left({\mathrm{\alpha }}_j{\mathrm{\theta }}_i-{\mathrm{\beta }}_j\right)\right)},$$

(2)

where $Y_{ij}$ is the observed score of examinee $i$ on item $j$ ; ${\mathrm{\alpha }}_j$ , ${\mathrm{\beta }}_j$ , ${\mathrm{\gamma }}_j$ , and ${\mathrm{\varsigma }}_j$ are item slope, intercept, guessing, and slipping parameter of item $j\left(j=1,2,3\dots ,J\right)$ , respectively; ${\mathrm{\theta }}_i\in \mathrm{\theta }$ is the ability parameter of examinee $i\left(i=1,2,3\dots ,N\right)$ ; $D$ is a constant, equals 1.702. Let ${\mathrm{\psi }}_j={\left({\mathrm{\alpha }}_j,{\mathrm{\beta }}_j,{\mathrm{\gamma }}_j,{\mathrm{\varsigma }}_j\right)}^T$ be the $j$ th column vector of the matrix $\mathbf{\Psi }$ collecting the item parameters for all J items. Similarly, let $y_i$ be the row vector of the matrix $\mathbf{Y}$ collecting all item responses of examinee $i$ . Then, the marginal likelihood function should be

$$\begin{array}{l}L\left(\mathbf{Y}|\mathbf{\Psi }\right)={\displaystyle \prod }_{i=1}^{\mathrm{N}}P\left(y_i|\mathbf{\Psi }\right)\\ =\prod _{i=1}^N{\int }_{{\mathrm{\theta }}_i}P\left(y_i|{\mathrm{\theta }}_i,\mathbf{\Psi }\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i\\ =\prod _{i=1}^N{\int }_{{\mathrm{\theta }}_i}\prod _{j=1}^J\left[P_j{\left({\mathrm{\theta }}_i\right)}^{Y_{ij}}{\left(1-P_j\left({\mathrm{\theta }}_i\right)\right)}^{1-Y_{ij}}\right]g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i\end{array}$$

(3)

with

$$P\left(y_i|\mathbf{\Psi }\right)={\int }_{{\mathrm{\theta }}_i}P\left(y_i|{\mathrm{\theta }}_i,\mathbf{\Psi }\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i,$$

(4)

$$P\left(y_i|{\mathrm{\theta }}_i,\mathbf{\Psi }\right)={\displaystyle \prod }_{j=1}^J{P}_j{\left({\mathrm{\theta }}_i\right)}^{Y_{ij}}{\left(1-P_j\left({\mathrm{\theta }}_i\right)\right)}^{1-Y_{ij}},$$

(5)

where $g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)$ is a density function for each examinee $i$ , and $\mathrm{\tau }$ is the vector containing the parameters of the examinee population ability distribution. Please refer to Supplemental Appendix A for a detailed derivation. It is worth pointing out that the BME is a natural extension of the one for 3PLM, but that the Fisher information matrix is slightly larger—it is size 4 × 4 rather than size 3 × 3—with small diagonal elements for the guessing and slipping parameters. Taken together, we find a great deal of numerical instability in the calculation of the Fisher information matrix, resulting in the BME requiring a much larger sample size than the 3PLM to avoid convergence problems. According to Waller and Feuerstahler (2017), the convergence rate for a sample of 1,000, which is always feasible for 3PLM, is only 9%. The new method in the current paper can address this issue by reformulating the 4PLM as a latent mixture model.

The Latent Mixture Modeling Reformulation of the 4PLM

From inspection, we can determine that the 4PLM can be regarded as a latent mixture model. First, note that

$$\begin{array}{l}{P}_j\left({\mathrm{\theta }}_i\right)={\mathrm{\gamma }}_j+\left(1-{\mathrm{\varsigma }}_j-{\mathrm{\gamma }}_j\right)\times P_j^{*}\left({\mathrm{\theta }}_i\right)\\ =P_j^{*}\left({\mathrm{\theta }}_i\right)\times \left(1-{\mathrm{\varsigma }}_j\right)+\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right)\times {\mathrm{\gamma }}_j.\end{array}$$

(6)

Next, we define a latent indicator variable $W_{ij}$ to indicate whether examinee $i$ knows the answer to item $j$ , where $W_{ij}=1$ if he knows it and 0 if he does not (Béguin & Glas, 2001; Culpepper, 2016). Furthermore, we model $W_{ij}$ probabilistically as $W_{ij}~Bernoulli\left(P_j^{*}\left({\mathrm{\theta }}_i\right)\right),$ meaning that $P\left(W_{ij}=1\right)=P_j^{*}\left({\mathrm{\theta }}_i\right)$ .

From the mixture-modeling perspective, depending on the value of $W_{ij}$ , the 4PLM can be decomposed into two distinct parts: $1-{\mathrm{\varsigma }}_j$ and ${\mathrm{\gamma }}_j$ . We can further obtain the conditional probabilities $P\left(Y_{ij}|W_{ij},{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)$ as:

$$\begin{array}{l}P\left(Y_{ij}=1|W_{ij}=1,{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)=1-{\mathrm{\varsigma }}_j,\\ P\left(Y_{ij}=0|W_{ij}=1,{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)={\mathrm{\varsigma }}_j,\\ P\left(Y_{ij}=1|W_{ij}=0,{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)={\mathrm{\gamma }}_j,\\ P\left(Y_{ij}=0|W_{ij}=0,{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)=1-{\mathrm{\gamma }}_j.\end{array}$$

(7)

Thus, we find the joint distribution of $Y_{ij}=1$ and $W_{ij}$ (conditional on item and person parameters) as follows:

$$\begin{array}{l}P\left(Y_{ij},W_{ij}|{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)\\ =\{\begin{array}{cc}\left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right),& \mathrm{for}{Y}_{ij}=1,W_{ij}=1\\ {\mathrm{\varsigma }}_j{P}_j^{*}\left({\mathrm{\theta }}_i\right),& \mathrm{for}{Y}_{ij}=0,W_{ij}=1\\ {\mathrm{\gamma }}_j\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right),& \mathrm{for}{Y}_{ij}=1,W_{ij}=0\\ \left(1-{\mathrm{\gamma }}_j\right)\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right),& \mathrm{for}{Y}_{ij}=0,W_{ij}=0\end{array}.\end{array}$$

(8)

Let $w_i$ be the latent indicator vector for examinee i. Then the joint distribution for the new augmented complete data $\left(\mathbf{Y},\mathbf{W},\mathrm{\theta }\right)$ is

$$P\left(y_i,w_i,{\mathrm{\theta }}_i|\mathbf{\Psi },\mathrm{\tau }\right)=P\left(y_i,w_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)$$

(9)

where

$$\begin{array}{l}P\left(y_i,w_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)\\ ={\displaystyle \prod }_{j=1}^J{\left[\left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right)\right]}^{Y_{ij}{W}_{ij}}{\left[{\mathrm{\gamma }}_j\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right)\right]}^{Y_{ij}\left(1-W_{ij}\right)}\\ \times {\left[{\mathrm{\varsigma }}_j{P}_j^{*}\left({\mathrm{\theta }}_i\right)\right]}^{\left(1-Y_{ij}\right)W_{ij}}{\left[\left(1-{\mathrm{\gamma }}_j\right)\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right)\right]}^{\left(1-Y_{ij}\right)\left(1-W_{ij}\right)}.\end{array}$$

(10)

Just like the MMLE/EM algorithm, the BE3M algorithm also marginalizes out the latent ability variable, so the likelihood function of the BE3M without priors is:

$$\begin{array}{l}L\left(\mathbf{Y},\mathbf{W}|\mathbf{\Psi }\right)=\prod _{i=1}^{N}P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi }\right)\\ =\prod _{i=1}^N{\int }_{{\mathrm{\theta }}_i}P\left({\mathbf{y}}_i,{\mathbf{w}}_i,{\mathrm{\theta }}_i|\mathbf{\Psi },\mathrm{\tau }\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i\\ =\prod _{i=1}^N{\int }_{{\mathrm{\theta }}_i}P\left(y_i,w_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i.\end{array}$$

(11)

This formulation follows Culpepper (2016) for 4PNO. In contrast to his work, however, we now develop a Bayesian modal method to estimate model parameters for the 4PLM. This approach follows the basic premise of an EM algorithm and proposes to split the maximization step into three parts: first maximizing over ${\mathrm{\gamma }}_j$ , then maximizing over ${\mathrm{\varsigma }}_j$ , and finally jointly maximizing over ${\mathrm{\alpha }}_j$ and ${\mathrm{\beta }}_j$ . We further use prior distributions for the item parameters in the maximization step. As such, we call this approach the Bayesian Expectation-Maximization-Maximization-Maximization (BE3M) method.

Expectation Step and Artificial Data

In the expectation step, we take the expectation of the log-likelihood over the incomplete portion of the data, the latent variables $W_{ij}$ and ${\mathrm{\theta }}_i$ . In the maximization step, we maximize the posteriors for the item parameters. This maximization process is equivalent to finding the zero of the summation of the first derivatives of the expected log-likelihood and the first derivatives of the priors. Thus, here we find the expected value of the first derivatives of the log-likelihood function. To begin, the first derivatives of the log-likelihood for item j is

$$\begin{array}{l}\frac{\partial \ln L\left(\mathbf{Y},\mathbf{W}|\mathbf{\Psi }\right)}{\partial {\mathrm{\psi }}_j}=\sum _{i=1}^N\frac{\partial \ln P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi }\right)}{\partial {\mathrm{\psi }}_j}=\sum _{i=1}^N\frac{1}{P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi }\right)}\frac{\partial P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi }\right)}{\partial {\mathrm{\psi }}_j}\\ =\sum _{i=1}^N\frac{1}{P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi }\right)}{\int }_{{\mathrm{\theta }}_i}\frac{\partial P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)}{\partial {\mathrm{\psi }}_j}g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i\\ =\sum _{i=1}^N\frac{1}{P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi }\right)}{\int }_{{\mathrm{\theta }}_i}\left[\frac{\partial \ln P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)}{\partial {\mathrm{\psi }}_j}\right]P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i\\ =\sum _{i=1}^N{\int }_{{\mathrm{\theta }}_i}\left[\frac{\partial \ln P\left({\mathbf{y}}_i,{\mathbf{w}}_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)}{\partial {\mathrm{\psi }}_j}\right]f\left({\mathrm{\theta }}_i|{\mathbf{y}}_i,{\mathbf{w}}_i,\mathrm{\tau },\mathbf{\Psi }\right)d{\mathrm{\theta }}_i\\ =\sum _{i=1}^N{\int }_{{\mathrm{\theta }}_i}\frac{\partial }{\partial {\mathrm{\psi }}_j}\sum _{j=1}^J\left[\begin{array}{c}{Y}_{ij}{W}_{ij}\ln \left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right)+Y_{ij}\left(1-W_{ij}\right)\ln {\mathrm{\gamma }}_j\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right)\\ +\left(1-Y_{ij}\right)W_{ij}\ln {\mathrm{\varsigma }}_j{P}_j^{*}\left({\mathrm{\theta }}_i\right)+\left(1-Y_{ij}\right)\left(1-W_{ij}\right)\ln \left(1-{\mathrm{\gamma }}_j\right)\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right)\end{array}\right]f\left({\mathrm{\theta }}_i|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)d{\mathrm{\theta }}_i\\ =\sum _{i=1}^N{\int }_{{\mathrm{\theta }}_i}\left[\begin{array}{l}\frac{W_{ij}-P_j^{*}\left({\mathrm{\theta }}_i\right)}{P_j^{*}\left({\mathrm{\theta }}_i\right)\left[1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right]}\frac{\partial P_j^{*}\left({\mathrm{\theta }}_i\right)}{\partial {\mathrm{\psi }}_j}\\ +\frac{\left(Y_{ij}-{\mathrm{\gamma }}_j\right)\left(1-W_{ij}\right)}{{\mathrm{\gamma }}_j\left(1-{\mathrm{\gamma }}_j\right)}\frac{\partial {\mathrm{\gamma }}_j}{\partial {\mathrm{\psi }}_j}+\frac{W_{ij}\left(1-{\mathrm{\varsigma }}_j-Y_{ij}\right)}{{\mathrm{\varsigma }}_j\left(1-{\mathrm{\varsigma }}_j\right)}\frac{\partial {\mathrm{\varsigma }}_j}{\partial {\mathrm{\psi }}_j}\end{array}\right]\times f\left({\mathrm{\theta }}_i|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)d{\mathrm{\theta }}_i\end{array}$$

(12)

with

$$f\left({\mathrm{\theta }}_i|{\mathbf{y}}_i,{\mathbf{w}}_i,\mathrm{\tau },\mathbf{\Psi }\right)=\frac{P\left({\mathbf{y}}_{\mathbf{i}}|\mathbf{\Psi },{\mathrm{\theta }}_i\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)}{{\int }_{{\mathrm{\theta }}_i}P\left(y_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)g\left({\mathrm{\theta }}_i|\mathrm{\tau }\right)d{\mathrm{\theta }}_i},$$

(13)

$$P\left({\mathbf{y}}_i|\mathbf{\Psi },{\mathrm{\theta }}_i\right)={\displaystyle \prod }_{j=1}^J{P}_j{\left({\mathrm{\theta }}_i\right)}^{Y_{ij}}{\left(1-P_j\left({\mathrm{\theta }}_i\right)\right)}^{1-Y_{ij}}.$$

(14)

Notice that the above expression includes the latent variable $W_{ij}$ so when we take the expectation of the log-likelihood, it will be necessary to calculate $E\left(W_{ij}|Y_{ij},{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)$ . To find this, we first find the conditional distribution of $W_{ij}$ using Bayes’ rule:

$$\begin{array}{l}P\left(W_{ij}=1|Y_{ij}=1,{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)\\ =\frac{\left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right)}{\left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right)+{\mathrm{\gamma }}_j\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right)}=\frac{\left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right)}{P_j\left({\mathrm{\theta }}_i\right)},\end{array}$$

(15)

$$\begin{array}{l}P\left(W_{ij}=1|Y_{ij}=0,{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)\\ =\frac{{\mathrm{\varsigma }}_j{P}_j^{*}\left({\mathrm{\theta }}_i\right)}{{\mathrm{\varsigma }}_j{P}_j^{*}\left({\mathrm{\theta }}_i\right)+\left(1-{\mathrm{\gamma }}_j\right)\left(1-P_j^{*}\left({\mathrm{\theta }}_i\right)\right)}=\frac{{\mathrm{\varsigma }}_j{P}_j^{*}\left({\mathrm{\theta }}_i\right)}{1-P_j\left({\mathrm{\theta }}_i\right)}.\end{array}$$

(16)

Thus, we find the conditional expectation of $W_{ij}$ to be

$$\begin{array}{l}E\left(W_{ij}|Y_{ij},{\mathrm{\theta }}_i,{\mathrm{\psi }}_j\right)\\ =\frac{\left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left({\mathrm{\theta }}_i\right)}{P_j\left({\mathrm{\theta }}_i\right)}{Y}_{ij}+\frac{{\mathrm{\varsigma }}_j{P}_j^{*}\left({\mathrm{\theta }}_i\right)}{1-P_j\left({\mathrm{\theta }}_i\right)}\left(1-Y_{ij}\right).\end{array}$$

(17)

Finally, let $X_k$ , $k=1,2,\dots ,q$ , be nodes on the ability scale with associated weights $A\left(X_k\right)$ . Using the Hermite-Gauss quadrature method to approximate the integral, the first derivatives of the expected log-likelihood function for each item is:

$$E\left[\frac{\partial \ln L\left(\mathbf{Y},\mathbf{W}|\mathbf{\Psi }\right)}{\partial {\mathrm{\psi }}_j}\right]\approx \left[\begin{array}{l}\sum _{i=1}^N\sum _{k=1}^q\frac{E\left(W_{ij}|Y_{ij},X_k,\mathbf{\Psi }\right)-P_j^{*}\left(X_k\right)}{P_j^{*}\left(X_k\right)\left[1-P_j^{*}\left(X_k\right)\right]}f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)\frac{\partial P_j^{*}\left(X_k\right)}{\partial {\mathrm{\alpha }}_j}\\ \sum _{i=1}^N\sum _{k=1}^q\frac{E\left(W_{ij}|Y_{ij},X_k,\mathbf{\Psi }\right)-P_j^{*}\left(X_k\right)}{P_j^{*}\left(X_k\right)\left[1-P_j^{*}\left(X_k\right)\right]}f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)\frac{\partial P_j^{*}\left(X_k\right)}{\partial {\mathrm{\beta }}_j}\\ \sum _{i=1}^N\sum _{k=1}^q\frac{\left(Y_{ij}-{\mathrm{\gamma }}_j\right)\left[1-E\left(W_{ij}|Y_{ij},X_k,\mathbf{\Psi }\right)\right]}{{\mathrm{\gamma }}_j\left[1-{\mathrm{\gamma }}_j\right]}f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)\frac{\partial {\mathrm{\gamma }}_j}{\partial {\mathrm{\gamma }}_j}\\ \sum _{i=1}^N\sum _{k=1}^q\frac{E\left(W_{ij}|Y_{ij},X_k,\mathbf{\Psi }\right)\left(1-{\mathrm{\varsigma }}_j-{\mathrm{\gamma }}_j\right)}{{\mathrm{\varsigma }}_j\left[1-{\mathrm{\varsigma }}_j\right]}f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)\frac{\partial {\mathrm{\varsigma }}_j}{\partial {\mathrm{\varsigma }}_j}\end{array}\right]$$

(18)

with

$$\begin{array}{l}f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)=f\left(X_k|y_i,\mathrm{\tau },\mathbf{\Psi }\right)\\ =\frac{P\left(y_i|X_k,\mathbf{\Psi }\right)A\left(X_k\right)}{{{\displaystyle \sum }}_{k=1}^{q}P\left(y_i|X_k,\mathbf{\Psi }\right)A\left(X_k\right)},\end{array}$$

(19)

$$P\left(y_i|X_k,\mathbf{\Psi }\right)={\displaystyle \prod }_{j=1}^J{P}_j{\left(X_k\right)}^{Y_{ij}}{\left[1-P_j\left(X_k\right)\right]}^{1-Y_{ij}},$$

(20)

where $f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)$ is the posterior density of ${\mathrm{\theta }}_i$ at quadrature point $X_k$ . Following Bock and Aitkin (1981), we can define several “artificial data” points for our algorithm:

$$\begin{array}{l}{\overline{f}}_k={\displaystyle \sum }_{i=1}^{N}f\left(X_k|y_i,\mathrm{\tau },\mathbf{\Psi }\right),{\overline{r}}_{jk}\\ ={\displaystyle \sum }_{i=1}^N{Y}_{ij}f\left(X_k|y_i,\mathrm{\tau },\mathbf{\Psi }\right),{\overline{f}}_{jk}^{\left(W\right)}\\ ={\displaystyle \sum }_{i=1}^{N}E\left(W_{ij}|Y_{ij},X_k,\mathbf{\Psi }\right)f\left(X_k|y_i,\mathrm{\tau },\mathbf{\Psi }\right),\end{array}$$

(21)

and

$${\overline{r}}_{jk}^{\left(W\right)}\underset{i=1}{\overset{N}{{\displaystyle =\sum }}}{Y}_{ij}E\left(W_{ij}|Y_{ij},X_k\mathbf{\Psi }\right)f\left(X_k|y_i,\mathrm{\tau },\mathbf{\Psi }\right).$$

(22)

By taking their expectations, we find that these expressions have very intuitive interpretations:

1. $E\left({\overline{f}}_k\right)={\overline{f}}_k$ is the expected number of persons (out of N) with ability $X_k$ ;

2. $E\left({\overline{r}}_{jk}\right)=P_j\left(X_k\right){\overline{f}}_k$ is the expected number of persons with ability $X_k$ to answer item j correctly;

3. $E\left({\overline{f}}_{jk}^{\left(W\right)}\right)=P_j^{*}\left(X_k\right){\overline{f}}_k$ is the expected number of persons with ability $X_k$ to know the answer to item j; and

4. $E\left({\overline{r}}_{jk}^{\left(W\right)}\right)=\left(1-{\mathrm{\varsigma }}_j\right)P_j^{*}\left(X_k\right){\overline{f}}_k$ is the expected number of persons with ability $X_k$ that know the answer to item j and also answer it correctly.

Maximization Step 1: ${\mathrm{\gamma }}_j$ -Parameters

Let ${\mathrm{\lambda }}_{{\mathrm{\gamma }}_j}$ be the summation of the expected first derivatives of the log-likelihood with respect to ${\mathrm{\gamma }}_j$ and the first derivatives of the prior distribution of ${\mathrm{\gamma }}_j$ . Following usual practice (Baker & Kim, 2004), we use a beta prior with hyperparameters ${\mathrm{\upsilon }}_j^{\left(\mathrm{\gamma }\right)}$ and ${\mathrm{\rho }}_j^{\left(\mathrm{\gamma }\right)}$ . Then,

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\gamma }}_j}=\frac{\partial \ln E\left(L\left(\mathbf{Y},\mathbf{W}|\mathbf{\Psi }\right)\right)}{\partial {\mathrm{\gamma }}_j}\\ +\frac{\partial \ln g\left({\mathrm{\gamma }}_j|{\mathrm{\upsilon }}_j^{\left(\mathrm{\gamma }\right)},{\mathrm{\rho }}_j^{\left(\mathrm{\gamma }\right)}\right)}{\partial {\mathrm{\gamma }}_j}\\ \approx \frac{{\mathrm{\upsilon }}_j^{\left(\mathrm{\gamma }\right)}-1}{{\mathrm{\gamma }}_j}-\frac{{\mathrm{\rho }}_j^{\left(\mathrm{\gamma }\right)}-1}{1-{\mathrm{\gamma }}_j}\\ +{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{k=1}^q\left[\frac{\left(Y_{ij}-{\mathrm{\gamma }}_j\right)\left(1-E\left(W_{ij}|Y_{ij},X_k,\mathbf{\Psi }\right)\right)}{{\mathrm{\gamma }}_j\left(1-{\mathrm{\gamma }}_j\right)}\right]\\ f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)\\ =\frac{\left({\mathrm{\upsilon }}_j^{\left(\mathrm{\gamma }\right)}-1\right)\left(1-{\mathrm{\gamma }}_j\right)-\left({\mathrm{\rho }}_j^{\left(\mathrm{\gamma }\right)}-1\right){\mathrm{\gamma }}_j}{{\mathrm{\gamma }}_j\left(1-{\mathrm{\gamma }}_j\right)}\\ +{\displaystyle \sum }_{k=1}^q\left[\frac{{\overline{r}}_{jk}-{\mathrm{\gamma }}_j{\overline{f}}_k-{\overline{r}}_{jk}^{\left(W\right)}+{\mathrm{\gamma }}_j{\overline{f}}_{jk}^{\left(W\right)}}{{\mathrm{\gamma }}_j\left(1-{\mathrm{\gamma }}_j\right)}\right].\end{array}$$

(23)

By setting ${\mathrm{\lambda }}_{{\mathrm{\gamma }}_j}$ equal to 0 and solving for ${\mathrm{\gamma }}_j$ , we get a closed-form solution for the guessing parameter:

$${\mathrm{\gamma }}_j=\frac{{\mathrm{\upsilon }}_j^{\left(\mathrm{\gamma }\right)}-1+{{\displaystyle \sum }}_{k=1}^q\left({\overline{r}}_{jk}-{\overline{r}}_{jk}^{\left(W\right)}\right)}{{\mathrm{\upsilon }}_j^{\left(\mathrm{\gamma }\right)}+{\mathrm{\rho }}_j^{\left(\mathrm{\gamma }\right)}-2+{{\displaystyle \sum }}_{k=1}^q\left({\overline{f}}_k-{\overline{f}}_{jk}^{\left(W\right)}\right)}.$$

(24)

Standard errors will require the expected second derivative as well, which is

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\gamma }}_j{\mathrm{\gamma }}_j}=\frac{\partial }{\partial {\mathrm{\gamma }}_j}{\mathrm{\lambda }}_{{\mathrm{\gamma }}_j}\\ =-\frac{{\mathrm{\upsilon }}_j^{\left(\mathrm{\gamma }\right)}-1+{{\displaystyle \sum }}_{k=1}^q\left({\overline{r}}_{jk}-{\overline{r}}_{jk}^{\left(W\right)}\right)}{{\mathrm{\gamma }}_j^2}\\ -\frac{{\mathrm{\rho }}_j^{\left(\mathrm{\gamma }\right)}-1+{{\displaystyle \sum }}_{k=1}^q\left({\overline{f}}_k-{\overline{f}}_{jk}^{\left(W\right)}-{\overline{r}}_{jk}+{\overline{r}}_{jk}^{\left(W\right)}\right)}{{\left(1-{\mathrm{\gamma }}_j\right)}^2}.\end{array}$$

(25)

Maximization Step 2: ${\mathrm{\varsigma }}_j$ -Parameters

Let ${\mathrm{\lambda }}_{{\mathrm{\varsigma }}_j}$ be the summation of the expected first derivatives of the log-likelihood with respect to ${\mathrm{\varsigma }}_j$ and the first derivatives of the prior distribution of ${\mathrm{\varsigma }}_j$ . As with the ${\mathrm{\gamma }}_j$ -parameters, we use a beta prior with hyperparameters ${\mathrm{\upsilon }}_j^{\left(\mathrm{\varsigma }\right)}$ and ${\mathrm{\rho }}_j^{\left(\mathrm{\varsigma }\right)}$ . Then,

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\varsigma }}_j}=\frac{\partial \ln E\left(L\left(\mathbf{Y},\mathbf{W}|\mathbf{\Psi }\right)\right)}{\partial {\mathrm{\varsigma }}_j}+\frac{\partial \ln g\left({\mathrm{\varsigma }}_j|{\mathrm{\upsilon }}_j^{\left(\mathrm{\varsigma }\right)},{\mathrm{\rho }}_j^{\left(\mathrm{\varsigma }\right)}\right)}{\partial {\mathrm{\varsigma }}_j}\\ \approx \frac{{\mathrm{\upsilon }}_j^{\left(\mathrm{\varsigma }\right)}-1}{{\mathrm{\varsigma }}_j}-\frac{{\mathrm{\rho }}_j^{\left(\mathrm{\varsigma }\right)}-1}{1-{\mathrm{\varsigma }}_j}\\ +{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{k=1}^q\left[\frac{E\left(W_{ij}|Y_{ij},X_k,\mathbf{\Psi }\right)\left(1-{\mathrm{\varsigma }}_j-Y_{ij}\right)}{{\mathrm{\varsigma }}_j\left(1-{\mathrm{\varsigma }}_j\right)}\right]\\ f\left(X_k|y_i,w_i,\mathrm{\tau },\mathbf{\Psi }\right)\\ =\frac{\left({\mathrm{\upsilon }}_j^{\left(\mathrm{\varsigma }\right)}-1\right)\left(1-{\mathrm{\varsigma }}_j\right)-\left({\mathrm{\rho }}_j^{\left(\mathrm{\varsigma }\right)}-1\right){\mathrm{\varsigma }}_j}{{\mathrm{\varsigma }}_j\left(1-{\mathrm{\varsigma }}_j\right)}\\ +{\displaystyle \sum }_{k=1}^q\left[\frac{\left(1-{\mathrm{\varsigma }}_j\right){\overline{f}}_{jk}^{\left(W\right)}-{\overline{r}}_{jk}^{\left(W\right)}}{{\mathrm{\varsigma }}_j\left(1-{\mathrm{\varsigma }}_j\right)}\right].\end{array}$$

(26)

By setting ${\mathrm{\lambda }}_{{\mathrm{\varsigma }}_j}$ equals to 0 and solving for ${\mathrm{\varsigma }}_j$ , we get a closed-form solution for the slipping parameter:

$${\mathrm{\varsigma }}_j=\frac{{\mathrm{\upsilon }}_j^{\left(\mathrm{\varsigma }\right)}-1+{{\displaystyle \sum }}_{k=1}^q\left({\overline{f}}_{jk}^{\left(W\right)}-{\overline{r}}_{jk}^{\left(W\right)}\right)}{{\mathrm{\upsilon }}_j^{\left(\mathrm{\varsigma }\right)}+{\mathrm{\rho }}_j^{\left(\mathrm{\varsigma }\right)}-2+{{\displaystyle \sum }}_{k=1}^q{\overline{f}}_{jk}^{\left(W\right)}}.$$

(27)

Standard errors will require the expected second derivatives as well, which is

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\varsigma }}_j{\mathrm{\varsigma }}_j}=\frac{\partial }{\partial {\mathrm{\varsigma }}_j}{\mathrm{\lambda }}_{{\mathrm{\varsigma }}_j}=-\frac{{\mathrm{\upsilon }}_j^{\left(\mathrm{\varsigma }\right)}-1+{{\displaystyle \sum }}_{k=1}^q\left({\overline{f}}_{jk}^{\left(W\right)}-{\overline{r}}_{jk}^{\left(W\right)}\right)}{{\mathrm{\varsigma }}_j^2}\\ -\frac{{\mathrm{\rho }}_j^{\left(\mathrm{\varsigma }\right)}-1+{{\displaystyle \sum }}_{k=1}^q{\overline{r}}_{jk}^{\left(W\right)}}{{\left(1-{\mathrm{\varsigma }}_j\right)}^2}.\end{array}$$

(28)

Maximization Step 3: ${\mathrm{\alpha }}_j$ - and ${\mathrm{\beta }}_j$ -Parameters

Let ${\mathrm{\alpha }}_j^{*}=\ln \left({\mathrm{\alpha }}_j\right)$ , and ${\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}}$ and ${\mathrm{\lambda }}_{{\mathrm{\beta }}_j}$ be the summation of the expected first derivatives of the log-likelihood and the first derivatives of the prior distribution with respect to ${\mathrm{\alpha }}_j^{*}$ and ${\mathrm{\beta }}_j$ , respectively. In this maximization step, we first solve the estimates of ${\mathrm{\alpha }}_j^{*}$ , and then obtain the ${\mathrm{\alpha }}_j$ parameter via ${\mathrm{\alpha }}_j=\exp \left({\mathrm{\alpha }}_j^{*}\right)$ . For both the ${\mathrm{\alpha }}_j^{*}$ - and ${\mathrm{\beta }}_j$ -parameters, we use a normal prior with means ${\mathrm{\mu }}_{{\mathrm{\alpha }}_j^{*}}$ and ${\mathrm{\mu }}_{{\mathrm{\beta }}_j}$ and variances ${\mathrm{\sigma }}_{{\mathrm{\alpha }}_j^{*}}^2$ and ${\mathrm{\sigma }}_{{\mathrm{\beta }}_j}^2$ , respectively. Then,

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}}=\frac{\partial \ln E\left(L\left(\mathbf{Y},\mathbf{W}|\mathbf{\Psi }\right)\right)}{\partial {\mathrm{\alpha }}_j^{*}}+\frac{\partial \ln g\left({\mathrm{\alpha }}_j^{*}|{\mathrm{\mu }}_{{\mathrm{\alpha }}_j^{*}},{\mathrm{\sigma }}_{{\mathrm{\alpha }}_j^{*}}^2\right)}{\partial {\mathrm{\alpha }}_j^{*}}\\ =D\exp \left({\mathrm{\alpha }}_j^{*}\right){\displaystyle \sum }_{k=1}^q\left[X_k\left({\overline{f}}_{jk}^{\left(W\right)}-E\left({\overline{f}}_{jk}^{\left(W\right)}\right)\right)\right]-\frac{{\mathrm{\alpha }}_j^{*}-{\mathrm{\mu }}_{{\mathrm{\alpha }}_j^{*}}}{{\mathrm{\sigma }}_{{\mathrm{\alpha }}_j^{*}}^2}\\ =D\exp \left({\mathrm{\alpha }}_j^{*}\right){\displaystyle \sum }_{k=1}^q\left[X_k\left({\overline{f}}_{jk}^{\left(W\right)}-P_j^{*}\left(X_k\right){\overline{f}}_k\right)\right]-\frac{{\mathrm{\alpha }}_j^{*}-{\mathrm{\mu }}_{{\mathrm{\alpha }}_j^{*}}}{{\mathrm{\sigma }}_{{\mathrm{\alpha }}_j^{*}}^2}\end{array}$$

(29)

and

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\beta }}_j}=\frac{\partial \ln E\left(L\left(\mathbf{Y},\mathbf{W}|\mathbf{\Psi }\right)\right)}{\partial {\mathrm{\beta }}_j}+\frac{\partial \ln g\left({\mathrm{\beta }}_j|{\mathrm{\mu }}_{{\mathrm{\beta }}_j},{\mathrm{\sigma }}_{{\mathrm{\beta }}_j}^2\right)}{\partial {\mathrm{\beta }}_j}\\ =-D{\displaystyle \sum }_{k=1}^q\left[{\overline{f}}_{jk}^{\left(W\right)}-E\left({\overline{f}}_{jk}^{\left(W\right)}\right)\right]-\frac{{\mathrm{\beta }}_j-{\mathrm{\mu }}_{{\mathrm{\beta }}_j}}{{\mathrm{\sigma }}_{{\mathrm{\beta }}_j}^2}\\ =-D{\displaystyle \sum }_{k=1}^q\left[{\overline{f}}_{jk}^{\left(W\right)}-P_j^{*}\left(X_k\right){\overline{f}}_k\right]-\frac{{\mathrm{\beta }}_j-{\mathrm{\mu }}_{{\mathrm{\beta }}_j}}{{\mathrm{\sigma }}_{{\mathrm{\beta }}_j}^2}.\end{array}$$

(30)

Unlike the guessing and slipping parameters, ${\mathrm{\alpha }}_j^{*}$ and ${\mathrm{\beta }}_j$ do not have closed-form solutions for their estimates, so it is necessary to find them numerically. The method used here is the Fisher scoring algorithm, which uses the expected values of the second derivatives, along with the first derivatives. These will also be used later for standard error calculations. We find the second derivatives as follows:

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\alpha }}_j^{*}}=-D^2\exp \left(2{\mathrm{\alpha }}_j^{*}\right)\\ {\displaystyle \sum }_{k=1}^q\left\{X_k^2{P}_j^{*}\left(X_k\right)\left[1-P_j^{*}\left(X_k\right)\right]{\overline{f}}_k\right\}-\frac{1}{{\mathrm{\sigma }}_{{\mathrm{\alpha }}_j^{*}}^2},\end{array}$$

(31)

$${\mathrm{\lambda }}_{{\mathrm{\beta }}_j{\mathrm{\beta }}_j}=-D^2{\displaystyle \sum }_{k=1}^q\left\{P_j^{*}\left(X_k\right)\left[1-P_j^{*}\left(X_k\right)\right]{\overline{f}}_k\right\}-\frac{1}{{\mathrm{\sigma }}_{{\mathrm{\beta }}_j}^2},$$

(32)

and

$$\begin{array}{l}{\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\beta }}_j}=D^2\exp \left({\mathrm{\alpha }}_j^{*}\right)\\ {\displaystyle \sum }_{k=1}^q\left\{X_k{P}_j^{*}\left(X_k\right)\left[1-P_j^{*}\left(X_k\right)\right]{\overline{f}}_k\right\}.\end{array}$$

(33)

Standard Errors of Parameter Estimation in BE3M

An important index of estimation quality is the standard error (SE). However, one major criticism of the EM algorithm is that parameter estimate SEs are not a natural product of the algorithm, and so other methods have to be devised (McLachlan & Krishnan, 2007). BE3M falls prey to this criticism, just as all members of the EM algorithm family do.

It can be shown that the inverse of the negative expected value of the matrix of second derivatives of the log-likelihood (i.e., the inverse of the information matrix) is the asymptotic covariance matrix of the estimates (Stuart & Kendall, 1968). As such, the square roots of the diagonal elements of the resulting matrix are the asymptotic standard errors of the parameters. Generically, the expected second derivative matrix can be written as

$${\mathrm{\Lambda }}_j=\left(\begin{array}{cccc}{\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\alpha }}_j^{*}}& {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\beta }}_j}& {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\gamma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\varsigma }}_j}\\ {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\beta }}_j}& {\mathrm{\lambda }}_{{\mathrm{\beta }}_j{\mathrm{\beta }}_j}& {\mathrm{\lambda }}_{{\mathrm{\beta }}_j{\mathrm{\gamma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\beta }}_j{\mathrm{\varsigma }}_j}\\ {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\gamma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\beta }}_j{\mathrm{\gamma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\gamma }}_j{\mathrm{\gamma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\gamma }}_j{\mathrm{\varsigma }}_j}\\ {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\varsigma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\beta }}_j{\mathrm{\varsigma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\gamma }}_j{\mathrm{\varsigma }}_j}& {\mathrm{\lambda }}_{{\mathrm{\varsigma }}_j{\mathrm{\varsigma }}_j}\end{array}\right)$$

(34)

in which each entry is the expected partial second derivative of the log-likelihood with respect to the subscript parameters for a given item j. Convergence issues resulting from the ill-conditioning of the matrix plague the MMLE/EM for 4PLM just as it does for the 3PLM (Baker & Kim, 2004).

In contrast, due to the divide-and-conquer strategy implemented in BE3M, we contend that the resulting expected second derivative matrix does not have these issues. In BE3M, the item parameter estimation problem is separated into three smaller distinct estimation problems: (1) estimation of the guessing parameter, (2) estimation of the slipping parameter, and (3) joint estimation of the difficulty and discrimination parameters. This separation has the benefit that any instabilities present in one of the three steps will not negatively affect either of the other two steps. Statistically, this separation implies that the covariances between ${\widehat{\mathrm{\gamma }}}_j$ and ${\widehat{\mathrm{\alpha }}}_j^{*}$ , ${\widehat{\mathrm{\beta }}}_j$ , and ${\widehat{\mathrm{\varsigma }}}_j$ are zero; similarly, the covariances between ${\widehat{\mathrm{\varsigma }}}_j$ and ${\widehat{\mathrm{\alpha }}}_j^{*}$ , ${\widehat{\mathrm{\beta }}}_j$ , and ${\widehat{\mathrm{\gamma }}}_j$ are zero. Collected together, we find the item j expected second derivative matrix to be

$${\mathrm{\Lambda }}_j=\left(\begin{array}{cccc}{\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\alpha }}_j^{*}}& {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\beta }}_j}& 0& 0\\ {\mathrm{\lambda }}_{{\mathrm{\alpha }}_j^{*}{\mathrm{\beta }}_j}& {\mathrm{\lambda }}_{{\mathrm{\beta }}_j{\mathrm{\beta }}_j}& 0& 0\\ 0& 0& {\mathrm{\lambda }}_{{\mathrm{\gamma }}_j{\mathrm{\gamma }}_j}& 0\\ 0& 0& 0& {\mathrm{\lambda }}_{{\mathrm{\varsigma }}_j{\mathrm{\varsigma }}_j}\end{array}\right).$$

(35)

The standard errors of the estimates for item j, then, are the square roots of the diagonal elements of the matrix ${\mathrm{\Sigma }}_j=-{\mathrm{\Lambda }}_j^{-1}$ . Furthermore, the standard error of the ${\mathrm{\alpha }}_j$ parameter can be transformed from ${\mathrm{\Lambda }}_j$ via the delta method (Oehlert, 1992).

Simulation Studies

This section presents Monte Carlo results regarding the accuracy of the developed BE3M procedure for estimating 4PLM item parameters, compared to the BME in mirt (version 1.30) package (Chalmers, 2012) and the MCMC in R package fourPNO (Culpepper, 2016) and “bayesmh” module of Stata/SE 16. Results from two simulation studies are reported. The goal of the simulation studies was to assess the influence of various design features, including the location of item threshold, sample size, and priors in educational (Simulation 1) and psychological (Simulation 2) scenarios.

Item and Person Parameter Generation

For Simulation 1, the item parameters were generated in the same way as in Culpepper (2016). Specifically, item parameters for 20 items were generated from ${\mathrm{\alpha }}_j~N\left(2,0.5\right)I(\mathrm{\alpha }>0)$ , ${\mathrm{\beta }}_j~N\left(0,0.5\right)$ , ${\mathrm{\gamma }}_j~Beta\left(2,8\right)$ , and ${\mathrm{\varsigma }}_j|{\mathrm{\gamma }}_j~Beta\left(0,1\right)I\left({\mathrm{\varsigma }}_j<1-{\mathrm{\gamma }}_j\right)$ . For Simulation 2, we used the item parameters for the 23 items calibrated in Waller and Reise (2010) application to an adolescent self-esteem scale; these were also used in Culpepper (2016). For both testing scenarios, three sample sizes (2,500, 5,000, and 10,000) of person parameters were sampled from the standard normal distribution.

The Setting of MCMC in R Package FourPNO

Since Culpepper (2016) introduced an extra latent variable and imposed various truncation restrictions on the item parameters, the priors in the fourPNO package have obvious differences from the traditional MCMC procedures (see Culpepper (2016) for details). Thus, the current studies used the same setting in the fourPNO package as the simulation studies of Culpepper (2016) in each condition: (1) the ${\mathrm{\alpha }}_j$ parameters were sampled from the truncated normal distribution $N\left(0,2\right)I\left(\mathrm{\alpha }>0\right)$ , which implied $E\left({\mathrm{\alpha }}_j\right)=\frac{2}{\sqrt{\mathrm{\pi }}}$ in the priors; (2) the ${\mathrm{\beta }}_j$ parameters were sampled from $N\left(0,2\right)$ ; (3) the uniform truncated $B\mathrm{eta}\left(\mathrm{\upsilon },\mathrm{\rho }\right)$ priors were employed for both $\mathrm{\gamma }$ and $\mathrm{\varsigma }$ parameters $\left({\mathrm{\upsilon }}_j^{\mathrm{\gamma }}={\mathrm{\rho }}_j^{\mathrm{\gamma }}={\mathrm{\upsilon }}_j^{\mathrm{\varsigma }}={\mathrm{\rho }}_j^{\mathrm{\varsigma }}=1\right)$ ; and (4) the ${\mathrm{\theta }}_i$ parameters were sampled from $N\left(0,1\right)$ . Furthermore, since Culpepper (2016) concluded that “no more than 50,000 iterations were needed to achieve $\widehat{R}$  < 1.1 for all item parameters,” only one chain (100,000 MCMC iterations and discarded the first 50,000 as burn-in) were run in the current studies.

The Setting of BE3M, BME in R Package mirt and MCMC in Stata

Both the BE3M and mirt used 61 quadrature points to approximate the Gaussian distribution from −6 to +6 in all simulation studies and real-data examples. The priors for α and β parameters were set as follows: $\ln \mathrm{\alpha }~N\left(0,0.25\right)$ and $\mathrm{\beta }~N\left(0,2\right)$ . Furthermore, two different priors for γ and ζ were used in the BE3M, mirt, and Stata:

1. $\mathrm{\gamma },\mathrm{\varsigma }~Beta\left(2,8\right)$ , which is the default setting for the guessing parameter in the advanced IRT software flexmirt (Houts & Cai, 2015), and

2. $\mathrm{\gamma },\mathrm{\varsigma }~Beta\left(4,16\right)$ , which is the default setting for the guessing parameter in the classical IRT software BILOG-MG3 (Zimowski et al., 2003)

In sum, for both the educational and psychological testing scenarios, 21 conditions were simulated: Seven estimation method conditions (fourPNO, BE3M with prior conditions 1 and 2, mirt with prior conditions 1 and 2, and Stata with prior conditions 1 and 2) each crossed with three examinee sample size conditions (2,500, 5,000, and 10,000). For each condition, 100 replications were run to reduce sampling error. Each condition was assessed by its accuracy in recovering individual item parameter as defined by bias and the root mean squared error (RMSE) across the 100 replications:

$$bia{s}_j=\frac{{\sum }_{s=1}^{S=100}\left({\widehat{\mathrm{\psi }}}_j^{\left(s\right)}-{\mathrm{\psi }}_j\right)}{S},$$

(36)

$$RMS{E}_j=\sqrt{\frac{{\sum }_{s=1}^{S=100}{\left({\widehat{\mathrm{\psi }}}_j^{\left(s\right)}-{\mathrm{\psi }}_j\right)}^2}{S}},$$

(37)

Furthermore, since Waller and Feuerstahler (2017) pointed out that the item parameter recovery for the individual parameter cannot reflect the entire situation of the accuracy of fitting the 4PLM, the current studies also introduced the root integrated mean squared error (RIMSE; Ramsay, 1991) to assess the entire accuracy via the item response function of the 4PLM:

$$RIMS{E}_j=\frac{{\sum }_{s=1}^{S=100}\sqrt{\int {\left[{\widehat{P}}_j^{\left(s\right)}\left({\mathrm{\theta }}_i\right)-P_j\left({\mathrm{\theta }}_i\right)\right]}^{2}g\left({\mathrm{\theta }}_i\right)d{\mathrm{\theta }}_i}}{S}.$$

(38)

Results

The results of two simulation studies will be respectively reported at two levels: the individual level for each type of parameters and the whole level for the item response function of the 4PLM.

As for the individual level, the detailed results for two simulation studies are presented in tables in Supplemental Appendices B and C, respectively. Only RMSEs (Tables 1 and 2) for the condition of 5,000 examinees are summarized and presented here since those for other conditions are very similar. The results of two simulation studies can be briefly summarized as follows: (1) In general, the results of MCMC in the fourPNO package replicate the previous study in both scenarios; (2) The item parameter recovery of BE3M are very similar with MCMC although the RMSEs of the α parameters yielded by the fourPNO package and Stata are slightly larger than those from the BE3M in the first simulation study; (3) The mirt package and the BE3M have similar bias and RMSEs for α and β parameters, but the BE3M has better performance in estimating the γ and ζ parameters; and (4) Compared with the mirt package, the BE3M, and Stata can provide more stable estimates when changing priors for γ and ζ parameters.

Table 1.

The RMSEs of Item Parameter Recovery for 5,000 Examinees in the Simulation Study 1.

Item	True values		RMSEs for α							RMSEs for β
	α	β	Four PNO	BE3M (2, 8)	BE3M (4, 16)	mirt (2, 8)	mirt (4, 16)	Stata (2, 8)	Stata (4, 16)	Four PNO	BE3M (2, 8)	BE3M (4, 16)	mirt (2, 8)	mirt (4, 16)	Stata (2, 8)	Stata (4, 16)
1	1.24	.62	.36	.29	.26	.35	.31	.32	.25	.21	.16	.16	.17	.17	.20	.18
2	.84	.75	.33	.14	.12	.13	.10	.17	.18	.20	.14	.13	.14	.13	.11	.10
3	1.69	.19	.39	.32	.29	.35	.31	.27	.19	.15	.12	.12	.12	.12	.10	.09
4	.59	.04	.54	.18	.11	.14	.10	.14	.15	.28	.24	.24	.24	.23	.21	.19
5	.88	.29	.25	.17	.19	.15	.16	.19	.23	.16	.15	.15	.15	.15	.14	.13
6	.80	−.29	.40	.20	.17	.20	.15	.14	.17	.17	.15	.14	.14	.13	.16	.15
7	.51	.10	.42	.24	.19	.15	.10	.25	.22	.22	.23	.23	.23	.23	.21	.24
8	2.82	1.37	.39	.45	.37	.45	.37	.48	.41	.20	.26	.23	.25	.23	.31	.30
9	.50	.03	.59	.20	.08	.15	.13	.21	.13	.35	.31	.26	.26	.18	.30	.27
10	.94	1.27	.41	.15	.14	.17	.15	.25	.24	.29	.21	.21	.23	.23	.27	.26
11	1.81	.90	.27	.24	.25	.25	.26	.28	.32	.15	.11	.11	.11	.11	.14	.16
12	1.18	.97	.37	.20	.26	.19	.25	.24	.28	.24	.16	.16	.16	.16	.21	.20
13	.53	.72	.42	.19	.15	.13	.09	.30	.22	.20	.18	.17	.15	.14	.16	.13
14	1.76	−.77	.33	.23	.30	.26	.32	.44	.52	.11	.11	.11	.11	.11	.10	.10
15	1.13	.24	.28	.23	.23	.24	.22	.32	.35	.13	.10	.10	.10	.10	.09	.12
16	.59	.53	.34	.23	.26	.17	.22	.32	.33	.17	.12	.11	.12	.11	.09	.09
17	.67	−.17	.31	.17	.17	.12	.12	.24	.25	.15	.15	.14	.15	.14	.10	.11
18	.75	1.23	.32	.11	.13	.10	.10	.16	.17	.20	.14	.13	.13	.12	.15	.14
19	.76	.79	.33	.15	.15	.13	.12	.20	.19	.25	.17	.16	.17	.16	.13	.14
20	1.30	1.10	.35	.25	.21	.28	.24	.24	.21	.25	.16	.15	.17	.16	.17	.15
1	.25	.23	.03	.03	.03	.04	.04	.04	.03	.06	.07	.06	.08	.06	.06	.05
2	.13	.17	.03	.02	.02	.03	.02	.02	.01	.06	.06	.04	.06	.04	.05	.05
3	.25	.14	.04	.03	.03	.04	.03	.02	.02	.03	.03	.02	.03	.02	.02	.02
4	.29	.16	.07	.07	.06	.12	.10	.05	.05	.06	.05	.04	.05	.03	.04	.03
5	.21	.07	.05	.04	.04	.05	.04	.04	.03	.04	.04	.05	.04	.04	.04	.05
6	.19	.20	.06	.06	.04	.07	.04	.05	.03	.04	.03	.03	.04	.03	.03	.03
7	.26	.09	.06	.07	.06	.09	.07	.05	.04	.08	.08	.08	.05	.06	.08	.09
8	.23	.04	.01	.02	.01	.02	.01	.02	.02	.02	.01	.01	.01	.01	.01	.02
9	.25	.30	.06	.06	.04	.10	.08	.04	.04	.06	.06	.06	.12	.10	.06	.05
10	.19	.14	.02	.02	.02	.03	.02	.02	.02	.09	.05	.04	.04	.02	.07	.08
11	.07	.13	.01	.01	.01	.01	.01	.01	.01	.04	.03	.03	.03	.03	.02	.02
12	.22	.03	.02	.02	.02	.02	.02	.02	.02	.06	.05	.07	.04	.06	.05	.07
13	.07	.18	.05	.03	.03	.02	.03	.03	.03	.15	.08	.06	.08	.04	.13	.09
14	.07	.09	.04	.03	.04	.03	.04	.04	.05	.01	.01	.01	.01	.01	.01	.01
15	.08	.23	.03	.02	.02	.03	.02	.03	.03	.05	.04	.04	.05	.04	.04	.03
16	.08	.04	.06	.05	.05	.04	.05	.05	.06	.14	.13	.14	.12	.14	.16	.16
17	.20	.08	.06	.07	.05	.07	.04	.06	.06	.04	.04	.05	.03	.04	.03	.04
18	.06	.12	.02	.01	.01	.01	.01	.01	.02	.15	.08	.08	.06	.06	.12	.11
19	.18	.15	.03	.03	.02	.03	.03	.02	.02	.08	.06	.05	.05	.03	.06	.05
20	.16	.23	.02	.02	.02	.02	.02	.02	.02	.07	.07	.05	.08	.06	.06	.05

Note. Bold indicates the relatively larger RMSEs in the same conditions.

Table 2.

The RMSEs of Item Parameter Recovery for 5,000 Examinees in the Simulation Study 2.

Item	True values		RMSEs for α							RMSEs for β
	α	Β	Four PNO	BE3M (2, 8)	BE3M (4, 16)	mirt (2, 8)	mirt (4, 16)	Stata (2, 8)	Stata (4, 16)	Four PNO	BE3M (2, 8)	BE3M (4, 16)	mirt (2, 8)	mirt (4, 16)	Stata (2, 8)	Stata (4, 16)
1	1.91	−.53	.34	.29	.26	.30	.28	.30	.30	.15	.14	.16	.15	.16	.12	.11
2	1.95	−.31	.34	.28	.29	.30	.33	.30	.33	.13	.13	.15	.13	.15	.12	.13
3	1.50	.07	.22	.23	.27	.23	.27	.25	.32	.10	.10	.11	.10	.11	.10	.11
4	1.12	.07	.16	.20	.24	.20	.23	.14	.18	.11	.09	.11	.09	.11	.11	.14
5	.89	.40	.14	.16	.20	.15	.18	.22	.26	.10	.09	.09	.09	.09	.09	.09
6	1.08	−.54	.23	.20	.23	.20	.22	.34	.36	.11	.11	.11	.11	.11	.11	.12
7	1.16	−.55	.24	.19	.20	.19	.20	.15	.20	.12	.11	.12	.12	.12	.12	.12
8	1.10	.01	.18	.21	.24	.20	.23	.27	.26	.11	.09	.10	.09	.10	.11	.11
9	.78	.35	.24	.16	.12	.17	.14	.14	.12	.20	.16	.20	.18	.24	.17	.20
10	1.23	.23	.14	.18	.26	.18	.26	.21	.29	.06	.05	.06	.05	.06	.06	.07
11	1.34	.55	.12	.17	.21	.16	.21	.24	.29	.07	.07	.07	.07	.07	.06	.07
12	1.54	−.74	.30	.24	.24	.26	.26	.36	.36	.17	.17	.18	.18	.19	.20	.20
13	1.16	.21	.31	.26	.25	.28	.29	.21	.18	.18	.18	.24	.19	.27	.17	.25
14	.84	.60	.16	.15	.16	.14	.14	.15	.19	.12	.09	.10	.09	.10	.08	.10
15	1.13	.17	.24	.18	.19	.18	.19	.23	.27	.11	.11	.13	.11	.14	.13	.14
16	.79	.94	.26	.13	.11	.12	.08	.11	.14	.12	.12	.13	.12	.13	.14	.14
17	1.27	.61	.16	.16	.23	.16	.22	.24	.31	.07	.06	.07	.06	.07	.07	.09
18	.94	1.29	.26	.14	.17	.13	.15	.17	.21	.11	.13	.13	.13	.13	.16	.16
19	.84	1.21	.23	.10	.13	.08	.10	.14	.16	.11	.11	.13	.11	.12	.14	.16
20	1.14	1.73	.25	.12	.17	.12	.15	.13	.19	.09	.11	.15	.10	.14	.13	.16
21	1.10	.28	.12	.17	.25	.16	.24	.20	.26	.06	.07	.07	.06	.07	.07	.08
22	.72	.38	.22	.17	.19	.15	.16	.14	.19	.15	.14	.14	.14	.15	.17	.16
23	.88	1.37	.30	.13	.15	.10	.12	.12	.17	.14	.13	.13	.13	.12	.11	.11
1	.04	.48	.01	.01	.01	.01	.01	.01	.01	.02	.02	.02	.02	.02	.01	.01
2	.02	.52	.01	.01	.01	.01	.01	.01	.01	.02	.02	.02	.02	.02	.02	.02
3	.02	.40	.01	.01	.01	.01	.01	.01	.02	.03	.03	.03	.03	.03	.02	.02
4	.02	.37	.01	.02	.03	.02	.03	.02	.03	.04	.03	.03	.04	.03	.03	.03
5	.04	.18	.02	.02	.03	.02	.02	.02	.02	.05	.06	.06	.05	.05	.06	.06
6	.06	.17	.05	.04	.05	.04	.05	.05	.05	.02	.02	.02	.02	.02	.03	.03
7	.07	.29	.03	.03	.03	.03	.03	.01	.03	.02	.02	.02	.02	.02	.03	.02
8	.04	.27	.02	.02	.03	.02	.03	.02	.03	.04	.04	.03	.04	.03	.04	.03
9	.05	.43	.02	.02	.02	.02	.01	.02	.02	.09	.08	.08	.12	.13	.08	.08
10	.01	.10	.01	.02	.03	.02	.03	.02	.03	.03	.03	.04	.03	.04	.03	.04
11	.02	.15	.01	.01	.01	.01	.01	.01	.02	.03	.04	.04	.04	.04	.03	.04
12	.06	.41	.02	.02	.02	.02	.02	.02	.03	.02	.02	.02	.03	.02	.02	.02
13	.02	.60	.01	.01	.01	.01	.01	.01	.02	.05	.06	.07	.07	.10	.04	.07
14	.04	.25	.02	.02	.02	.01	.02	.02	.02	.07	.06	.05	.07	.05	.04	.04
15	.03	.39	.02	.02	.02	.02	.02	.02	.02	.04	.03	.03	.04	.03	.03	.03
16	.04	.27	.01	.01	.01	.01	.01	.01	.01	.09	.09	.06	.10	.07	.07	.05
17	.01	.16	.01	.01	.02	.01	.02	.02	.02	.04	.04	.04	.04	.04	.04	.04
18	.09	.06	.01	.01	.01	.01	.01	.01	.01	.13	.08	.11	.07	.10	.10	.12
19	.02	.18	.01	.01	.01	.01	.01	.01	.01	.11	.08	.07	.07	.04	.06	.05
20	.00	.18	.00	.00	.01	.00	.01	.01	.01	.12	.08	.06	.08	.06	.05	.06
21	.02	.07	.02	.02	.03	.02	.03	.03	.04	.03	.04	.04	.03	.04	.04	.05
22	.24	.05	.04	.04	.04	.04	.04	.04	.04	.06	.05	.07	.04	.06	.06	.08
23	.06	.09	.01	.01	.01	.01	.01	.01	.01	.15	.10	.11	.10	.11	.10	.12

Note. Bold indicates the relatively larger RMSEs in the same conditions.

As for the whole level, the RIMSEs for two simulation studies are presented in Figure 1. In general, with the increasing of the sample size, the RIMSEs of BE3M, fourPNO, Stata, and mirt package gradually decrease, and all of the RIMSEs across different conditions are relatively small. Furthermore, the results of RIMSEs show that there is no obvious difference between BE3M and BME in R package mirt although the BE3M has better performance in estimating the asymptotes, which is consistent with Waller and Feuerstahler (2017).

Figure 1.

The RMSEs of item parameter recovery for all simulation conditions.

Real Data Example

The 4PLM is applied to bullying items collected as part of the 2005 to 2006 Health Behavior in School-Aged Children study (Iannotti, 2005) funded by the National Institute of Child Health and Human Development. This real data example was chosen to replicate previous results (Culpepper, 2016), and to demonstrate comparability between BE3M and the MCMC approaches (fourPNO package and “bayesmh” module in Stata) in terms of estimation accuracy and BE3M’s advantage in terms of estimation time. The prior settings for the “bayesmh” module in Stata are the same as those for BE3M and BME in the simulation studies. To further compare BE3M with BME, items were also calibrated with BME in the mirt package. In addition, as we have discussed in the introduction, there is no available SEs for all item parameters in the mirt package because the Beta priors (“expbeta” in mirt) were imposed on the guessing and slipping parameters in this example.

The data includes student responses on ten items related to bullying behavior for 7,491 adolescents. The ten bullying items asked students “How often have you bullied another student(s) at school in the past couple of months in the ways listed below?” The items were dichotomized as “1→0” and “>1→1” to indicate the student has not bullied or has bullied, respectively, in the given way over the past couple of months. For this application, the latent variable $W_{ij}$ can be interpreted as

$$W_{ij}=\{\begin{array}{c}\begin{array}{l}1,\mathrm{Student}i\mathrm{bullied}\mathrm{another}\\ \mathrm{Student}\mathrm{as}\mathrm{asked}\mathrm{in}\mathrm{item}j,\end{array}\\ \begin{array}{l}0,\mathrm{Student}i\mathrm{did}\mathrm{not}\mathrm{bully}\\ \mathrm{another}\mathrm{student}\mathrm{as}\mathrm{asked}\mathrm{in}\mathrm{item}j,\end{array}\end{array}$$

so that ${\mathrm{\varsigma }}_j$ is the probability that a student who bullied over the past couple of months neglected to report bullying. On the other hand, it can be expected for ${\mathrm{\gamma }}_j$ to be zero in this application, because it is unlikely for students who did not bully to report bullying behaviors (Culpepper, 2016).

Results

Item calibration results from the three methods are summarized in Table 3. The estimates and SEs yielded by fourPNO are the same as in Culpepper’s (2016) study due to we used identical package (fourPNO) and estimation conditions. The BE3M produces point estimates comparable to both MCMC approaches, and estimates of the BE3M and Stata remain stable when changing priors. Unlike the performance of the mirt package in the simulation studies, it is interesting to note that the changing of priors has no obvious impact on the guessing and slipping parameters, but rather on the sloping parameters for the Item 1 and 2, which perhaps caused by the unstable Hessian matrix of the BME method. The implication is that careful consideration of priors should be taken if one uses the mirt package (which implements the BME method) to pursue precise estimates of the 4PLM.

Table 3.

The Estimated Item Parameters and the SEs for HBSC Data.

Item	α (SE)							β (SE)
	Four PNO	Stata		BE3M		mirt		Four PNO	Stata		BE3M		Mirt
		β (2, 8)	β (4, 16)	β (2, 8)	β (4, 16)	β (2, 8)	β (4, 16)		β (2, 8)	β (4, 16)	β (2, 8)	β (4, 16)	β (2, 8)	β (4, 16)
1	4.44 (.44)	6.16 (.72)	6.31 (.79)	5.20 (.61)	5.40 (.51)	5.75	6.32	.67 (.10)	.92 (.16)	.90 (.18)	.82 (.03)	.92 (.03)	.84	1.01
2	3.54 (.33)	3.59 (.52)	3.67 (.63)	3.54 (.10)	3.70 (.11)	3.62	4.14	.74 (.09)	.83 (.12)	.97 (.15)	.81 (.03)	.94 (.03)	.79	.94
3	1.21 (.07)	1.31 (.09)	1.39 (.09)	1.25 (.03)	1.36 (.04)	1.25	1.35	1.10 (.06)	1.15 (.08)	1.18 (.07)	1.14 (.02)	1.20 (.02)	1.12	1.18
4	1.41 (.06)	1.49 (.06)	1.55 (.09)	1.50 (.04)	1.56 (.05)	1.49	1.55	1.83 (.06)	1.91 (.05)	1.94 (.08)	1.94 (.03)	1.98 (.03)	1.91	1.96
5	1.63 (.09)	1.66 (.11)	1.74 (.10)	1.71 (.06)	1.79 (.06)	1.68	1.75	2.32 (.11)	2.35 (.13)	2.42 (.11)	2.43 (.03)	2.51 (.03)	2.38	2.46
6	1.93 (.10)	2.02 (.12)	2.06 (.11)	2.09 (.07)	2.16 (.08)	2.04	2.12	2.93 (.13)	3.04 (.15)	3.08 (.14)	3.16 (.04)	3.24 (.04)	3.09	3.17
7	2.44 (.15)	2.50 (.15)	2.58 (.14)	2.78 (.12)	2.91 (.13)	2.66	2.76	4.05 (.22)	4.15 (.21)	4.25 (.20)	4.59 (.05)	4.74 (.05)	4.39	4.52
8	1.58 (.08)	1.63 (.10)	1.68 (.09)	1.65 (.05)	1.72 (.05)	1.63	1.70	2.13 (.09)	2.16 (.11)	2.19 (.09)	2.21 (.03)	2.27 (.03)	2.18	2.24
9	2.23 (.15)	2.26 (.11)	2.36 (.16)	2.49 (.10)	2.59 (.11)	2.36	2.45	3.75 (.22)	3.79 (.17)	3.92 (.23)	4.16 (.05)	4.27 (.05)	3.95	4.07
10	2.29 (.16)	2.41 (.15)	2.44 (.19)	2.58 (.11)	2.70 (.12)	2.44	2.54	3.89 (.24)	4.07 (.22)	4.11 (.27)	4.36 (.05)	4.51 (.05)	4.13	4.27
1	.00 (.00)	.01 (.00)	.01 (.01)	.01 (.00)	.01 (.00)	.00	.01	.17 (.01)	.18 (.01)	.18 (.01)	.18 (.01)	.18 (.01)	.18	.18
2	.00 (.00)	.01 (.01)	.02 (.01)	.01 (.00)	.02 (.00)	.01	.02	.18 (.01)	.17 (.01)	.18 (.01)	.17 (.01)	.17 (.01)	.17	.17
3	.01 (.01)	.01 (.01)	.02 (.01)	.01 (.00)	.02 (.00)	.01	.02	.03 (.02)	.05 (.02)	.07 (.02)	.03 (.00)	.06 (.01)	.03	.06
4	.00 (.00)	.00 (.00)	.00 (.00)	.00 (.00)	.00 (.00)	.00	.00	.01 (.01)	.02 (.02)	.04 (.02)	.01 (.00)	.03 (.01)	.01	.03
5	.01 (.00)	.01 (.00)	.01 (.00)	.01 (.00)	.01 (.00)	.01	.01	.01 (.01)	.02 (.01)	.04 (.02)	.01 (.00)	.03 (.01)	.01	.03
6	.00 (.00)	.00 (.00)	.00 (.00)	.00 (.00)	.01 (.00)	.00	.00	.01 (.01)	.02 (.01)	.04 (.02)	.01 (.00)	.03 (.01)	.01	.03
7	.00 (.00)	.00 (.00)	.00 (.00)	.00 (.00)	.00 (.00)	.00	.00	.01 (.01)	.02 (.01)	.04 (.02)	.01 (.01)	.03 (.01)	.01	.03
8	.01 (.00)	.01 (.00)	.01 (.00)	.01 (.00)	.01 (.00)	.01	.01	.01 (.01)	.02 (.02)	.04 (.02)	.01 (.00)	.03 (.01)	.01	.03
9	.01 (.00)	.00 (.00)	.01 (.00)	.01 (.00)	.01 (.00)	.01	.01	.01 (.01)	.02 (.01)	.04 (.02)	.01 (.01)	.03 (.01)	.01	.03
10	.00 (.00)	.00 (.00)	.00 (.00)	.00 (.00)	.01 (.00)	.00	.00	.01 (.01)	.02 (.01)	.04 (.02)	.01 (.01)	.03 (.01)	.01	.03

Note. Bold indicates an obvious difference in estimates when changing priors.

In terms of SEs, the two priors do not seem to influence the quality of estimation with similar values for the SEs in general. We also find that the SEs for the MCMC approaches are similar to those obtained with BE3M, with the exceptions of the sloping parameters.

The main finding is that BE3M is much more computationally efficient than the MCMC. The Gibbs sampler in the fourPNO package required approximately 47 minutes to complete 100,000 iterations and the MH sampler in Stata required approximately 3 hours to complete 10,000 iterations with N = 7,491 using a 2.4 GHz processor and 6 GB of RAM while BE3M only needed about 1 minute. We hope that this finding in particular will help to promote further applications of the 4PLM in practice.

Discussion

In recent years, interest in the four-parameter IRT model has been on the rise in measurement research. The current study presented a new Bayesian formulation of the 4PLM model. Through two simulation studies and a real data analysis, we demonstrated that the BE3M combines the strengths of these two approaches; the BE3M can produce estimates as accurately as the Gibbs sampling method as quickly as the Bayes model estimation approach using mirt. To facilitate calibration, the authors have made available an R package IRTBEMM (will be published on CRAN soon). The authors believe this package offers a better alternative for 4PLM calibration, helping to promote the use of the 4PLM in educational and psychological research.

Several opportunities exist for future research related to the 4PLM. First, future research can help to identify other substantive research areas and large survey datasets that could benefit from modeling the slipping effect (Beghetto, 2019). Second, for successful, high-efficient integration into a CAT testing system, being able to implement online calibration of item parameters is important (Zheng, 2016). Since the method proposed at present is a member of the EM family to which most online calibration algorithms belong, there is little barrier to developing online calibration methods based on the BE3M.

Supplemental Material

sj-docx-1-sgo-10.1177_21582440211052556 – Supplemental material for Fast Bayesian Estimation for the Four-Parameter Logistic Model (4PLM)

Supplemental material, sj-docx-1-sgo-10.1177_21582440211052556 for Fast Bayesian Estimation for the Four-Parameter Logistic Model (4PLM) by Chanjin Zheng, Shaoyang Guo and Justin L Kern in SAGE Open