Introduction

Entire volumes (e.g., Baker & Kim, 2004) have been dedicated to the discussion of statistical parameter estimation techniques for item response theory (IRT) models. There has also been much recent development in the technical literature on improved methods for estimating complex IRT models (e.g., Cai, 2010a, 2010b; Edwards, 2010; Rijmen & Jeon, 2013). We offer here a map to help researchers and graduate students understand the fundamental challenges of IRT parameter estimation, and appropriately contextualize the underlying logic of some of the proposed solutions. We assume that the reader is familiar with elementary probability concepts such as prior, posterior, and likelihood, as well as the equations for describing statistical models for categorical observed data, for example, logistic regression. For methodologically inclined readers interested in studying IRT parameter estimation and in trying out some of the approaches discussed here, the combination of conceptual sections and more technical sections should be sufficient as a basis of software implementation.

We do not discuss limited-information estimation methods derived from categorical factor analysis or structural equation modeling (see, e.g., Bolt, 2005), but not because estimators based on polychoric correlation matrices and weighted least squares are not useful. They certainly may be, when conditions are appropriate (see, e.g., Wirth & Edwards, 2007), and for a long time, limited-information methods provided the only practically viable means for conducting formal model appraisal; although on that latter point, the situation has changed dramatically in the past few years (see, e.g., Maydeu-Olivares, 2013). In choosing to focus exclusively on full-information approaches that are based on either likelihood or Bayesian derivations, we believe that we provide readers with insight that tends to be obscured by the technicality that tends to accompany limited-information approaches. That is, the latent variables in IRT models are missing data, and had the latent variable scores been available, estimation for IRT models would have been a rather straightforward task. For the sake of simplicity, we contextualize our discussion with unidimensional logistic IRT models for dichotomously scored outcomes, but the missing data formulation applies generally across a far wider range of statistical modeling frameworks, IRT modeling included.

Univariate Logistic Regression

Item Response Theory Model as Multivariate Logistic Regression

Some Notation

Suppose a hypothetical assessment is made up of i = 1, …, I dichotomously scored items. An item score of one indicates a correct or endorsement response, and zero otherwise. Furthermore, suppose that the assumption of unidimensionality holds for this set of items. Let us use the standard notation of θ _j to denote the latent variable score for individual j. The two-parameter logistic (2PL) item response model specifies the conditional response probability curve (also known as the traceline) of a correct response or endorsement as a function of the latent variable and the item parameters: where α _i and β _i are the item intercept and slope parameters. The parentheses in T _i (θ; α _i, β _i)highlight the fact that the response probabilities are conditional on θ _j, and that they also depend on the item parameters. Let Y _ij be a Bernoulli (0–1) random variable representing individual i’s response to item j, and let y _ij be a realization of Y _ij. This suggests a formulation of the conditional probability of the event Y _ij = y _ij similar to Equation (3.4),

Under the assumption of unidimensionality, the latent variable θ alone explains all the observed covariations among the items. In other words, conditionally on θ, the item response probabilities are independent for an individual, that is, the probability of response pattern y _j = (y _1j , …, y _Ij) factors into a product over individual item response probabilities: where on the left-hand side we collect all item intercept and slope parameters into γ = (α ₁, …, α _I, β ₁, …, β _I), a 2I-dimensional vector. The joint probability of the observed and latent variables is equal to the product of the conditional probability of the observed variables given the latent variables, times the prior probability of the latent variables: where h(θ) is prior (population) distribution of the latent variable θ. In IRT applications, it is customary to resolve the location and scale indeterminacy of the latent variable by assuming that the θ’s are standard normal, so h(θ) does not contain free parameters.

From Equation (3.14), a natural derived quantity is the marginal probability of the response pattern, after integrating the joint probability over θ: Unfortunately Equation (3.15) is already the simplest form that we can obtain, given the combination of the IRT model and normally distributed latent variable. Note that the marginal probability does not depend on the unobserved latent variable scores; it is a function solely of the observed item response pattern and the item parameters.

As in the case of logistic regression, we assume the individuals are independent, with latent variable scores sampled independently from the population distribution. Let Y be a matrix of all observed item responses. If we treat the item responses as fixed once observed, the marginal likelihood function for all the item parameters in γ , based on observed item response data, can be expressed as: Because the marginal likelihood L ( γ|Y ) does not depend on the unobserved θ values, it may be referred to as the observed data likelihood.

Under some circumstances, this likelihood function can be optimized directly, again using Newton Raphson or Fisher Scoring-type algorithms (see, e.g., Bock & Lieberman, 1970), but those circumstances are rather limited. In particular, Bock and Lieberman (1970) noted that this direct approach does not generalize well to the case of many items and many parameters because of computing demands. We would add that even as computers have become faster and storage cheaper, what the direct approach glosses over is a missing data formulation of latent variable models that is central to our understanding of IRT and of other modern statistical techniques such as random effects regression modeling, or modeling of survey nonresponse. This missing data formulation was made transparent by Dempster, Laird, and Rubin’s (1977) classical paper that coined the term Expectation-Maximization (EM) algorithm.

Bock-Aitkin EM Algorithm

Bock and Aitkin (1981) began with the insight that the marginal probability can be approximated to arbitrary precision by replacing the integration with a summation over a set of Q quadrature points over θ: where X _q is a quadrature point, and W _q is the corresponding weight. In the simplest case, one may take the quadrature points as a set of equally spaced real numbers over an interval that captures sufficiently the probability mass of the population distribution, for example, from −6 to +6 in increments of 0.1, and the corresponding weights as a set of normalized ordinates of the quadrature points from the population distribution .

Another important insight of Bock and Aitkin (1981) is that the height of the posterior distribution at quadrature point X _q can be approximated to arbitrary precision as well:

Ignoring constants involving the prior distribution h(θ _j) from Equation (3.17), the complete data log-likelihood for the item parameters can be written as: Following the logic inspired by the Fisher Identity, the conditional expected complete data likelihood given provisional item parameter values approximated by quadrature, case by case, as follows:

The third and arguably most important insight from Bock and Aitkin (1981) is that by interchanging the order of summation, they realized that the posterior probabilities can be accumulated over individuals first: where is understood as the conditional expected proportion of individuals that respond positively/correctly to item i, and is the conditional expected proportion of individuals that respond negatively/incorrectly to item i, at quadrature point X _q. Taken together, let be the conditional expected proportion of individuals at quadrature point X _q, then we have:

Equation (3.25) highlights the fact that the conditional expected complete data log-likelihood is a set of I independent logistic regression log-likelihoods, with the quadrature points X _q serving as the predictor values, and weights given by n _iq, and r _iq serving as the positive outcome “frequency” at X _q. The inner summation over the quadrature points bears striking similarity to the log-likelihood given in Equation (3.5). The only difference is that in standard logistic regression, the weights n _j and number of successes y _j are integers, whereas in the case of Bock-Aitkin EM, n _iq and r _iq may be fractional and will change from cycle to cycle, given different item parameter values. With the Fisher Scoring algorithm developed in Section 2, optimization of Q(γ|Y;γ *) is straightforward, which leads to updated parameter estimates that may be used in the next cycle.

In general, Bock-Aitkin EM (or any EM algorithm) alternates between the following two steps from a set of initial parameter estimates, say γ ⁽⁰⁾, and it generates a sequence of parameter estimates γ ⁽⁰⁾,…,γ ^(k),…, where , that converges under some very general conditions to the MLE of γ as the number of cycles k tends to infinity (Wu, 1983):

• E-step. Given γ(^k), evaluate the conditional expected complete data log-likelihood Q( γ|Yγ (^k), which is taken to be a function of γ.
• M-step. Maximize Q( γ|Yγ (^k) to yield updated parameter estimates γ ^(k+1). Go back to E-step and repeat. The cycles are terminated when the estimates from adjacent cycles stabilize.

The application of the EM algorithm to IRT epitomizes the elegance of the missing data formulation in statistical computing. Finding MLEs in logistic regression analysis is a task that statisticians already know how to do. The goal of the E-step, then, is to replace the missing data with conditional expectations that depend on values of θ, represented using a set of discrete quadrature points. Once the missing data are filled in, complete data estimation can be accomplished with tools that are already available. Leveraging the conditional independence built into the IRT model, the M-step logit analyses can even be run in parallel and the overall demand on computing resources is rather low. Although the EM algorithm is only first-order (linearly) convergent, and may be slow (by optimization researchers’ standard), the statistical intuition is simply too elegant to ignore. Thissen (1982) extended the unconstrained Bock-Aitkin EM to handle parameter restrictions and used it to estimate the Rasch IRT model.

Metropolis-Hastings Robbins-Monro Algorithm

Implementing the Metropolis-Hastings Sampler

At this point, a critical missing link is a method to draw random values of from its posterior predictive distribution. Cai (2006) proposed the use of the Metropolis-Hastings method, for several reasons. First, we see from Equation (3.18) that while the posterior predictive distribution is analytically intractable in that it does not belong to any “named” distribution family, it is proportional to the joint probability of observed item responses and latent variables: The Metropolis-Hastings method is ideally suited to the task of sampling a posterior when the normalization constant is not readily available. In addition, the right-hand side of Equation (3.32) is the complete data likelihood at γ(^k), which is evaluated in any event to compute the item gradients and information matrices required in the approximation step of MH-RM. Furthermore, the sampling of the θ values can be accomplished in parallel, as the individual P(θ _j|y _j; γ ^(k))’s are fully independent. Finally, the Monte Carlo approximation in Equation (3.28) remains unbiased even if the draws are not independent, for example from a Markov chain.

Implementing the Metropolis-Hastings method is straightforward. For each individual j, we begin with some initial value of θ _j, say, , and let us call it the current state of θ _j. We now draw a random increment from an independent normal sampler, with mean 0 and standard deviation equal to σ. Let this increment value be denoted δ _j . By adding the increment to the current state, we have produced a proposal for a new state of . We now evaluate the right-hand side of Equation (3.32) at both current and proposal states, and form the following likelihood ratio: If is larger than 1.0, meaning that the proposed move to a new state increased the likelihood relative to the current state, we accept the move and set the proposal state as the new current state. If is smaller than 1.0, meaning that the proposed move decreased the likelihood, we accept the move with probability equal to the likelihood ratio. This can be accomplished by drawing, independently, a uniform (0,1) random number, u _j and comparing it to . If u, _j is smaller than , we accept the proposed move and set the proposal state as the new current state. If u _j is larger than the likelihood ratio, we reject the proposal, and remain at the current state. Iterating this sampler will produce a Markov chain that converges to P(θ _j|y _j;γ ^(k)) in distribution.

As the chain evolves, dependent samples from this chain can be regarded as samples from the target distribution. To avoid excessive dependence on the initial state, one can drop the samples in the so-called burn-in phase of the chain. For the IRT model, experience suggests that this burn-in phase typically amounts to not more than 10 iterations of the Metropolis-Hastings sampler. Of course, this assumes that the chain is appropriately tuned by monitoring the rate of acceptance of the proposed moves and scaling the increment density standard deviation σ up (for decreased acceptance rate) or down (for increased acceptance rate). Roberts and Rosenthal (2001) discussed the statistical efficiency of Metropolis-Hastings samplers, and its relationship to optimal scaling. Asymptotically efficient chains can be obtained by tuning the acceptance rate to around 25 percent.

Application

We analyze a well-known data set (Social Life Feelings), analyzed by Bartholomew (1998), among others, to illustrate the Bock-Aitkin EM and MH-RM algorithms. The data set contains responses from J = 1,490 German respondents to five statements on perceptions of social life. The responses were dichotomous (endorsement vs. non-endorsement of the statements). Table 3.1 presents the 2⁵ = 32 response patterns and their associated observed response frequencies.

Let us first examine the most frequently encountered response pattern (0, 1, 1, 0, 0), wherein 208 respondents endorsed items 2 and 3 and none of the others. Following the logic of Bock-Aitkin EM, we must first choose a set of quadrature points for

approximating the E-step integrals. Here we use a set of 49 quadrature points equally spaced between −6 and +6. Next we must also choose a set of initial values for the item parameters. For the sake of variety, we let the initial values of the item intercepts be α ¹ = −1.5, α ²= −1, α ³ = 0, α ⁴ = 1, α ⁵ = 1.5, and let all initial slopes be equal to 1.0. We are now ready to begin our first E-step.

Figure 3.1 contains a set of three plots showing the relationship between the prior distribution, the likelihood function for response pattern (0,1,1,0,0) evaluated at the initial values of the item parameters, and the implied posterior distribution—formed by multiplying the likelihood and the prior, point by point over the quadrature points, and then normalized to sum to one. The prior and the posterior distributions are shown as discrete probability point masses over the quadrature points. The ordinates of the normalized prior distribution have been multiplied by the observed sample size (1,490), and those of the posterior distribution have been multiplied by the observed frequency associated with the response pattern (208).

For each item, depending on the response (0 or 1), the posterior probabilities are accumulated as per Equation (3.24). For instance, item 1’s response is zero, which means

Figure 3.1   Multi-panel plot showing the relationship among the prior (population) distribution, the likelihood (products of tracelines shown) for response pattern (0,1,1,0,0), and the posterior

that the current set of posterior probabilities must be added into the values over the Q quadrature points. Similarly, because item 2’s response is one, for that item, the posterior probabilities are added into the r _2q values for all q. Regardless of the response, the posterior probabilities are added into n_iq for all items and quadrature points.

For each response pattern, there is a set of corresponding three-panel plots that generate the posterior probabilities over the same set of quadrature points. These posterior probabilities are accumulated into the item-specific r _iq and n _iq values, depending on the item response. At the end of the E-step, the weights n _iq and (artificial) response frequencies n _iq are submitted to the M-step for logit analyses.

Figure 3.2 presents the current and updated tracelines for item 1 after one cycle of E- and M-step. The current tracelines (dashed curves) are at their initial values of α ¹ = −1.5 and β ² = 1. The ordinates of the solid dots are equal to r _1q /n _1q , representing the conditional expected probability of the endorsement response for item 1 at each of the quadrature points. The size of each of the solid dots is proportional to the conditional expected number of respondents at each of the corresponding quadrature points. The updated tracelines (solid curves) correspond to α ¹ = −2.12 and β ¹ = −1.11. It is obvious that the updated tracelines are much closer approximations of the “data” generated by the E-step conditional expectations. Other items can be handled similarly. Thus iterating the E- and M-steps leads to a sequence of item parameter estimates that eventually converges to the MLE. At the maximum, the following item parameters are obtained: α ₁ = −2.35, α ₂ = 0.80, α ₃ = 0.99, α ₄ = −0.67, α ₅ = 1.10, β ₁ = 1.20, β ₂ = 0.71, β ₃ = 1.53, β ₄ = 2.55, β ₅ = 0.92.

Let us now turn to the application of the MH-RM algorithm. The MH-RM algorithm also requires the characterization of the posterior distribution of θ, but it uses it differently than Bock-Aitkin EM: The Metropolis-Hastings sampler is used to generate dependent draws from this posterior, given provisional item parameter values and the

Figure 3.2   Current and updated tracelines for item 1 after one cycle of E- and M-step

Figure 3.3   Multi-panel plot showing the relationship among the prior (population) density, the likelihood (products of tracelines shown) for response pattern (0,1,1,0,0), and the posterior density approximated in two ways. The solid posterior curve is found by numerically evaluating the normalized posterior ordinates over a range of θ values. The dashed posterior curve is found by plotting the empirically estimated density of the posterior from the random draws produced by a Metropolis-Hastings sampler for θ

samples are used in complete data estimation with the Robbins-Monro method. Figure 3.3 plots the relationship among the prior density (standard normal), the likelihood function for response pattern (0,1,1,0,0) evaluated at the initial values of the item parameters (α ¹ = −1.5, α ² = 1, α ³ = 0, α ⁴ = 1, α ⁵ = 1.5, β ¹ = … β = 1.0), and the implied normalized posterior distribution. The prior and the posterior are represented as smooth solid curves.

There are 208 individuals associated with this response pattern. For each individual, we iterate the Metropolis-Hastings sampler 10 times and take the last draw as the posterior sample. We then empirically estimate a density function from the 208 posterior samples. The estimated density is shown as the dashed curve superimposed on the true implied posterior density. The two are obviously quite close, indicating the Metropolis-Hastings method can generate adequate posterior samples. For our sampler, the tuning constant (proposal dispersion σ) is equal to 2.0. The starting value of θ is equal to the standardized total score associated with response pattern (0,1,1,0,0). The total score for this response pattern is 2.0. The sample average total score over all response patterns is 2.17, and

Figure 3.4   The iteration history of the slope parameter estimate for item 1 from MH-RM. The solid line is the MLE

the sample standard deviation of the total score is 1.32, so the standardized total score is −0.13.

Together with the observed item responses, the posterior samples for all 1,490 individuals form the complete data set, with the posterior draws serving the role of predictor values. Complete data derivatives are evaluated and the item parameters are updated according to Equation (3.31) with the Robbins-Monro method. Figure 3.4 shows a sequence of parameter estimates from 100 iterations of the MH-RM algorithm for the slope parameter of item 1. The solid line is the MLE of that parameter from the Bock-Aitkin EM algorithm (β ¹ = 1.20). The MH-RM estimates contain random error initially, but as the number of cycles increases, the Robbins-Monro method filters out the error to achieve convergence near the MLE.

Discussion and Conclusion

The key emphasis of our discussion of IRT and IRT parameter estimation is on a missing data formulation: The unobserved latent variable θ amounts to missing data. Had the missing data been observed, IRT parameter estimation would be standard logit analysis. Motivated by this missing data formulation, we described estimation algorithms that augment the observed data by replacing the missing data either deterministically by their posterior expected values or stochastically by multiple imputations from the posterior predictive distribution of θ. The former approach (Bock-Aitkin EM) requires numerical integration with quadrature. The latter approach (MH-RM) requires the combination of elements of Markov chain Monte Carlo (Metropolis-Hastings sampler) with stochastic approximation (Robbins-Monro method). In both approaches, it is revealed that an insight due to Fisher (1925) provided the key equation that connects the complete data and observed data models. We illustrated the estimation algorithms with an empirical data set.

This presentation has been restricted to parameter estimation for unidimensional IRT models for dichotomous responses, to keep the focus on the essential ideas. The procedures described here straightforwardly generalize to either multidimensional IRT models, or IRT models for polytomous responses, such as those used in the PROMIS® measures (Reeve et al., 2007), or both. We have alluded to the generalization to multidimensional IRT; that simply adds multidimensional quadrature grids, or vector-valued random draws, to the procedures described in the previous sections. Parameter estimation for IRT models for polytomous responses requires that the computations described in this chapter for each of the two dichotomous responses be carried out for each of the several polytomous responses, and that the values of partial derivatives be calculated for each parameter of the model. The necessary partial derivatives for most commonly used IRT models are available from a variety of sources, and are brought together in the book-length treatment of this topic by Baker and Kim (2004).

Author’s Note: Li Cai is supported by grants from the Institute of Education Sciences (R305D140046 and R305B080016) and National Institute on Drug Abuse (R01DA026943 and R01DA030466). David Thissen has been supported by a PROMIS® cooperative agreement from the National Institutes of Health (NIH) Common Fund Initiative (U01AR052181). The views expressed here belong to the authors and do not reflect the views or policies of the funding agencies or grantees.