mirt : A Multidimensional Item Response Theory Package for the R Environment

Item response theory (IRT) is widely used in assessment and evaluation research to explain how participants respond to item level stimuli. Several R packages can be used to estimate the parameters in various IRT models, the most ﬂexible being the ltm (Ri-zopoulos 2006), eRm (Mair and Hatzinger 2007), and MCMCpack (Martin, Quinn, and Park 2011) packages. However these packages have limitations in that ltm and eRm can only analyze unidimensional IRT models eﬀectively and the exploratory multidimensional extensions available in MCMCpack requires prior understanding of Bayesian estimation convergence diagnostics and are computationally intensive. Most importantly, multidimensional conﬁrmatory item factor analysis methods have not been implemented in any R package. The mirt package was created for estimating multidimensional item response theory parameters for exploratory and conﬁrmatory models by using maximum-likelihood meth-ods. The Gauss-Hermite quadrature method used in traditional EM estimation (e.g., Bock and Aitkin 1981) is presented for exploratory item response models as well as for conﬁrmatory bifactor models (Gibbons and Hedeker 1992). Exploratory and conﬁrma-tory models are estimated by a stochastic algorithm described by Cai (2010a,b). Various program comparisons are presented and future directions for the package are discussed.


Introduction
Item response theory (IRT) is widely used in educational and psychological research to model how participants respond to test items in isolation and in bundles (Thissen and Wainer 2001). It is a general framework for specifying the functional relationship between a respondent's underlying latent trait level (i.e., commonly known as 'ability' in educational testing, or 'factor score' in the factor analysis tradition 1 ) and an item level stimulus. IRT methodology attempts to model individual response patterns by specifying how the underlying latent traits interact with the item's characteristics -such as an item's easiness or discrimination ability -to form an expected probability of the response pattern. As such, a major goal of IRT is to separate the item parameters and population sampled characteristics from manifest data so that both may be understood and studied separately. This parameter separation often requires advanced numerical analysis techniques for effective estimation and can become computationally burdensome as the model complexity increases.
The simplest and most popular IRT models are those that specify a single (i.e., unidimensional) latent trait. Unidimensional IRT models have been predominant across social science and educational research mainly because of historical traditions, but also because multidimensional parameter estimation procedures were not fully developed or studied (Baker and Kim 2004;Reckase 2009). While unidimensional models are often simpler and can have various interesting and important measurement properties (e.g., Rasch models), many psychological constructs are unavoidably multidimensional in nature. For instance, unobservable constructs might be understood as a combination of sub-scale components nested within -or along-side -a more general construct, or as compensatory or noncompensatory factors that combine to influence the item response probabilities. A major impediment when deciding to utilize these models has been that the estimation of the item parameters in higher dimensional space (due to increasing the number of factors) is computationally difficult for standard numerical integration techniques. However, with recent advances in estimation theory, coupled with the advances in computational power of personal computers, multidimensional IRT research is finally beginning to blossom as a feasible statistical analysis methodology (Edwards 2010;Reckase 2009;Wirth and Edwards 2007).
Several R (R Development Core Team 2012) packages can be used to fit IRT models, such as: the ltm package (Rizopoulos 2006), which can handle the Rasch, general latent trait, threeparameter logistic, and graded response models; the eRm package (Mair and Hatzinger 2007), which can fit the rating scale and partial credit models; and the MCMCpack package (Martin et al. 2011), which can estimate k-dimensional unconstrained two-parameter item response models (normal, heteroscedastic, and robust estimation) using a Markov chain Monte Carlo (MCMC) approach. While useful in their own right, these packages have limitations in that ltm and eRm can only analyze unidimensional IRT models effectively while the multidimensional extensions available in MCMCpack require prior understanding of Bayesian estimation diagnostics, are computationally demanding, can require a large amount of memory storage, and are only available for dichotomous item response sets.

Item response models
Item response models typically follow a monotonically increasing probability form with respect to the underlying latent traits. Two well known and commonly used logistic response models for dichotomous and polytomous item responses are the Birnbaum (1968) three-parameter model (3PL; which can be reduced to a 1PL or 2PL model) and the Samejima (1969) ordinal response model, respectively. Although first introduced as unidimensional modeling functions, both models can readily generalize to incorporate more than one factor. Let i = 1, . . . , N represent the distinct participants, j = 1, . . . , n the test items, and suppose that there are m latent factors θ i = (θ i1 , . . . , θ im ) with associated item slopes α j = (α 1 , . . . , α m ). For the multidimensional 3PL model, the probability of answering a dichotomous item correctly is where d j is the item intercept, γ j is the so-called 'guessing' parameter, and D is a scaling adjustment (usually 1.702) used to make the logistic metric more closely correspond to the traditional normal ogive metric (Reckase 2009).
For Samejima's (1969) multidimensional ordinal response model, suppose that there are C j unique categories for item j, with intercepts d j = (d 1 , . . . , d (C j −1) ). Here we define the boundary of response probabilities as These boundaries lead to the conditional probability for the response x ij = k to be Note that the 3PL is in fact a special case of the ordinal response model (with the inclusion of a lower asymptote parameter) and can be defined with boundaries in the same fashion, where (2) would consist of only two possible values: Recognizing this, and letting Ψ be the collection of all item parameters, allows us to declare the likelihood equations more concisely.
Expressing the data in indicator form, where Assuming a distributional form g(θ) (most often a multivariate normal) the marginal distribution becomes where there are m-fold integrals to be evaluated. Finally, this brings us to the observed data likelihood function. Letting X represent the complete N × n data matrix, the observed likelihood equation is

Exploratory and confirmatory item analysis
IRT can be applied in a way that is analogous to exploratory and confirmatory factor analysis for continuous variables (McDonald 1999). Historically, IRT models began in a confirmatory spirit by modeling the item response probabilities as a function of a single underlying factor, with varying degrees of how the item slopes (e.g., Rasch versus 2PL) and intercepts (e.g., ordinal, nominal, or partial credit) were related. But IRT can also be applied in an exploratory manner, whereby the number of dimensions are not assumed known beforehand, and are instead estimated empirically by comparing nested models (Bock and Aitkin 1981) or by rotating the factor loadings matrix to find a more parsimonious structure (Bock, Gibbons, and Muraki 1988). The TESTFACT program (Wood et al. 2003) was specifically designed for this approach, but other software exist that use different methods of estimation, such as NO-HARM (Fraser 1998) and Mplus's various WLS estimators (Muthén and Muthén 2008), which use limited-information algorithms, and BMIRT (Yao 2008) which uses Bayesian MCMC estimation techniques.
Confirmatory item analysis is useful when more than one factor is thought to be present in the data but various constraints (such as zero slopes) should be imposed. One of the first approaches in this spirit was the bifactor method (Holzinger and Swineford 1937) explicated by Gibbons and Hedeker (1992) for dichotomous data. The inspiration for bifactor models is that a single factor is believed to be present in all items, but with additional clusters of local dependencies formed by other independent specific factors. This approach was later generalized to polytomous data (Gibbons et al. 2007) and further expanded to accommodate more than one local dependency caused by specific factors (Cai 2010c).
A more general approach that accommodates linear constraints and missing data can be found in stochastic estimation techniques, such as Bayesian MCMC estimation with Metropolis-Hastings sampling (Metropolis, Rosenbluth, Teller, and Teller 1953;Hastings 1970), Gibbs sampling (Casella and George 1992), or by employing the Metropolis-Hastings Robbins-Monro (MH-RM) algorithm (Cai 2010b). The MH-RM algorithm is explained in more detail below since it is implemented in the mirt package.

Parameter estimation
IRT parameter estimation has been a progressive science over that past 60 years, moving from heuristic estimation techniques to more advanced Bayesian MCMC methods (Baker and Kim 2004). The early focus was on estimating the item specific parameters for unidimensional models, and until Bock and Aitkin (1981) introduced an EM based estimation solution IRT, applications were largely limited to small testing situations (Baker and Kim 2004). The EM algorithm, which was introduced by using fixed Gauss-Hermite quadrature, appeared to be a reasonable solution for lower dimensional models without compromising numerical accuracy. Unfortunately this technique quickly becomes inefficient as the number of dimensions increases, since the number of quadrature points required for estimating the 'E-step' increases exponentially and must be accommodated for by decreasing the number of quadrature. A partial solution for a moderate number of dimensions was described by Schilling and Bock (2005), where the authors demonstrated that adaptive quadrature could be used for better accuracy when a smaller number of quadratures per dimension is used, but the problem of high-dimensional solutions still remained.
More recently, a solution to the high-dimensionality problem has been to employ stochastic estimation methods for exploratory and confirmatory item analysis. Bayesian MCMC methods have been explored by Edwards (2010) and Sheng (2010), and both authors have released software to estimate the parameters for polytomous and dichotomous response models, respectively. These methods are not implemented in the mirt package, so they will not be discussed further, but see Bolt (2005) and Wirth and Edwards (2007) for more thorough reviews of using these full-information estimation methods and for item response model estimation in general.
The two methods that are implemented in the mirt package are the fixed quadrature EM method for exploratory (Bock et al. 1988;Muraki and Carlson 1995) and bifactor (Gibbons and Hedeker 1992;Gibbons et al. 2007) models, and the Metropolis-Hastings Robbins-Monro method for exploratory (Cai 2010a) and confirmatory (Cai 2010b) polytomous models.
3.1. Estimation using the expectation-maximization algorithm Bock and Aitkin (1981) were the first to propose a feasible method for estimating the item parameters of large scale tests using a method similar to the Expectation-Maximization (EM) algorithm (Dempster, Laird, and Rubin 1977). As explained in Bock et al. (1988) and Muraki and Carlson (1995) this method is appropriate for low to moderate factor solutions, so long as the number of quadratures per dimension decreases as the number of factors increases. For the following EM estimation methods we will explore only the special case when all the data are dichotomous. To begin, approximate (4) for each unique response vector by using m-fold Gauss-Hermite quadraturẽ From this result the observed likelihood based on the u unique response patterns with r u individuals with identical patterns becomes Differentiating with respect to an arbitrary item parameter within item j and integrating out the m-dimensions of θ gives Approximating (8) by using quadratures gives The 'E-step' of the EM algorithm consists of finding (9) and (10) by treating Ψ as provisionally known when computing L (K). The 'M-step' then consists of finding the 0 root of (11) independently for each item. The EM process is repeated until the change between iterations falls below some pre-specified tolerance.

A special case for EM estimation: The bifactor model
The full-information bifactor model (Gibbons and Hedeker 1992;Gibbons et al. 2007) combines a unidimensional model with the so-called 'simple structure' item loadings model (Thurstone 1947). The purpose is to estimate a common latent trait alongside independent components for each item so that local dependencies are accounted for properly. The bifactor model can specify only one additional item specific factor (although see Cai 2010c), but is not limited in the number of factors estimated since the quadratures remain fixed regardless of the number of specific factors extracted.
Define P as The gradient with respect to an arbitrary parameter from item j, expressed in quadrature form, is then where andN As before, (15) and (16) are computed by treating Ψ as provisionally known, and (13) is then solved for each item independently. This model has the benefit of having a fixed number of quadratures regardless of the number of specific factors estimated, and is closely related to 'Testlet' response models (Wainer, Bradlow, and Wang 2007).

Estimation via the Metropolis-Hastings Robbins-Monro algorithm
The Metropolis-Hastings Robbins-Monro (MH-RM) algorithm estimates the item parameters by using a stochastically imputed complete-data likelihood with an assumed population distributional form (typically, multivariate normal) For exploratory item factor analysis the population mean vector µ is usually assumed to be a fixed m × 1 vector of zeros, and Σ is assumed to be a fixed m × m 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 (17) the complete-data log-likelihood is composed of two additive components: the log-likelihood for a multivariate ordinal regression and the log-likelihood relationship between the factor scores.
The MH-RM algorithm deals with the integration problem in a different way than the traditional EM approach. For the EM algorithm θ is treated as set of 'nuance' parameters with a known distribution and are then integrated out of the likelihood equation (using numerical quadrature) so that the marginal frequencies (r j andN j ) can be approximated. The marginal frequencies are then used to update the item level parameters, and with the newly updated parameters the process is repeated. The MH-RM and related methods (such as the stochastic approximated EM) use stochastic methods to temporarily fill in the missing θ instead, and once filled in the values are treated as if they were 'known'. Given the newly augmented data the item parameters can then be updated by using conventional root-finding methods that use the complete-data log-likelihood function directly. Imputation methods are not exact or determinate but often allow for easier and more convenient evaluation of higher dimensional integrals than their numerical quadrature counterparts. The MH-RM is a more recent attempt to control for the inaccuracies borne out of using stochastic imputations to approximate the θ parameters.
Cai (2010a) demonstrated that when using a stochastically imputed complete-data model, and properly accounting for the error in the stochastic imputations, maximum-likelihood estimation and observed parameter standard errors can be calculated. The estimation begins by computing θ (d) , given initial start values in Ψ, by using a Metropolis-Hastings sampler (Metropolis et al. 1953;Hastings 1970). This allows the complete-data gradient vector and hessian matrix to be calculated and utilized in updating the initial parameters 2 . The initial parameters are then updated for the next iteration by using a single Newton-Raphson correction and these new parameters are then used to sample θ (d+1) . This process is repeated for several cycles and constitutes the so-called 'burn-in' period. After sufficient burn-in iterations, a set of parameter estimates are reserved to compute the approximate starting values for the MH-RM algorithm. Finally, the last set of parameter updates are conducted using the same method as before, but are now controlled by using the Robbins-Monro (Robbins and Monro 1951) root finding algorithm, which slowly converges as the decreasing gain constant approaches zero. In this way, the inaccuracies borne from the Metropolis-Hastings sampler are properly accounted for when attempting to maximize the parameter estimates. The MH-RM algorithm is useful for estimating both exploratory and confirmatory item response models. For specifying confirmatory models, several linear restrictions can be imposed on single parameters (e.g., α 1j = 0) or between parameters (e.g., α 1j ≡ α 2j ). This is accomplished by defining a matrix L, which selects the parameters that are to be estimated, and a vector c, which contains the fixed values for the parameters. For example, is a vech stacked vector containing all of the possible group and item parameters. This example shows that the slope and intercept are estimated for item 1, the guessing parameter is fixed at 0.1, and the latent mean and variance are fixed at 0 and 1, respectively. The usefulness of this formulation lies in how to manipulate the gradient and hessian given the user defined restrictions, since the restricted gradient and hessian can be easily expressed. This result implies that one can first estimate the complete-data gradient and hessian given no restrictions, and simply apply (19) and (20) to update only a subset of the parameters that are to be estimated.
The MH-RM algorithm offers a flexible way to specify both confirmatory and exploratory item models but has three main shortcomings: the estimation times will be much larger for lower dimensional problems when compared to an EM solution, the observed data log-likelihood is not calculated automatically and must be estimated use Monte Carlo integration, and parameter estimates will not be identical between different estimations. It is recommended to use the MH-RM algorithm only when the dimensionality increases to the point where the EM solution becomes difficult due to the large number of quadratures required (approximately around four or more factors), or when a multidimensional confirmatory model is required. However, the user should keep in mind that even the MH-RM method will become slower as the dimensions increase, so testing parameters in, for example, a bifactor model with eight specific factors may take an extended amount of time.

Implementation in the mirt package
Four separate datasets are used to demonstrate how each major estimation function typically behaves. These are: the well known 5-item LSAT7 data (Bock and Lieberman 1970); the SAT12 data set, which is available as an example from the TESTFACT (Wood et al. 2003) manual; a simulated data set constructed from the item parameters with orthogonal factors found in Reckase (2009, p. 153); and a simulated set derived from the parameters shown in Appendix B. All simulated data were constructed by using simdata() from the mirt package.
There are four main functions for estimating MIRT models: mirt(), bfactor(), polymirt(), and confmirt(), where the latter two employ the MH-RM algorithm. Individual item plots can be generated using itemplot(), factor scores (with MAP or EAP estimation) can be estimated using the fscores() function, and MIRT data with known parameters can be simulated using simdata(). The subroutines that do not directly relate to model estimation are demonstrated in the appendices.
The data that are passed to all the estimation functions must be either a matrix or data.frame consisting of numeric values for the item responses. For example, coding an 'ability' test as 0 for incorrect and 1 for correct, or coding a Likert type format with 1 representing the lowest category and 5 as the highest category are conceptually the preferred layout, although technically the choice of direction is arbitrary. Responses that are omitted must be coded as NA.
Only complete data-sets can be passed to these functions, so if the data are in a tabulated format (see below) the use of expand.

An example with the mirt() function
The LSAT7 data found in Bock and Lieberman (1970) initially were presented in tabulated form, with the number of individuals with identical response pattern in the rightmost column. This type of input can be modified easily with the expand.table(), and here we see the default use of mirt() for a one-factor model R> library("mirt") R> data <-expand. The values is the first column (a_1) reflect the item slopes for factor one, while the values is the second column (d_1) correspond to the item intercept. The names of these columns also reflects their relationship to equation (1), although the γ parameter has been renamed to guess instead to avoid confusion.

R> summary(mod1)
Unrotated factor loadings: The object returned by mirt() has various diagnostic tools available to determine where the model may be malfunctioning. By default, residuals() computes the local dependence (LD) pairwise statistic between each pair of items, which is very similar to a signed χ 2 value (Chen and Thissen 1997). Also, a standardized version of the LD statistic (Cramer's V) is printed above the diagonal to aid in interpretation when items contain more than two response options and hence more degrees of freedom (i.e., objects return by polymirt() and confmirt() ). The residuals can also be in the form of the marginal expected frequencies for each response pattern by specifying the input the option restype = "exp". Estimating more than one factor with mirt() is performed simply by changing the second numeric input value. There are several areas that should be considered when increasing the number of dimensions to extract. To begin, by default the number of quadrature values used during estimation decreases so that estimation time is lower, and so that there are not any memory leaking problems. However, while this means that solutions using mirt() do not take as long to estimate, it does mean that the accuracy of estimating the parameters will suffer. For moderate to high-dimensional solutions it may be better to use the polymirt() and confmirt() functions (see below).

R> (mod2 <-mirt(data, 2))
Call: mirt(data = data, nfact = 2) Full-information factor analysis with 2 factors Converged in 11 iterations using 15 quadrature points. Log-likelihood = -2652.096 AIC = 5334.192 BIC = 5407.808 G^2 = 21.82, df = 17, p = 0.1915, RMSEA = 0.017 Again, the coefficients can be extracted as above, but now summary() holds a different purpose. Since the orientation of the factor loadings is arbitrary the initially extracted solution should be rotated to a simpler structure to better facilitate interpretation. The default rotation method is the varimax criterion, but many other rotations available in the GPArotation package (Bernaards and Jennrich 2005)  Nested model comparison can be performed using the anova() generic, returning a likelihoodratio χ 2 test as well as returning the difference in AIC and BIC values. As seen below, the difference between mod1 and mod2 is marginally-significant at the α = 0.05 cut-off, while the AIC and BIC decrease indicate that the extra parameters estimated likely do not contribute additional useful information.

R> anova(mod1, mod2)
Chi-squared difference: X2 = 9.884, df = 4, p = 0.0424 AIC difference = -0.116 BIC difference = -24.655 Finally, mirt() contains plotting features used for didactic illustration (see Appendix A), but in general are not as flexible as the plink package (Weeks 2010). The generic functions demonstrated in this section are applicable for the remaining estimation methods as well, and essentially perform the same behavior. For more detailed information refer to the documentation found in the mirt package.

An example with the bfactor() function
Next we examine a bifactor model for the SAT12 data. First, we must change the raw data in SAT12 into dichotomous (correct-incorrect) form by using the key2binary() function. From here we declare which specific-factor affects which item in numeric sequence vector, where each element corresponds to its respective item. Here, for example, item one is a function of the general factor and of specific factor 2, while item two is influenced by specific factor 3 and the general factor, etcetera. As an added feature not mentioned before, fixed values for the lower asymptotes may be specified a priori for all estimation functions.
If number of factors is greater than one then multivariate discrimination and intercepts are also available. The multivariate discrimination and intercepts are defined as respectively. These indices are potentially useful for determining an item's overall utility across factors (Reckase and McKinley 1991), and are printed for all mirt objects when using coef().
R> b_mod2 <-bfactor(data, specific, guess, par.prior = list(int = c(0, 4), + int.items = 32)) Intercept prior for item ( As we can see the intercept parameter appears to be pulled slightly towards 0 which helps to add numerical stability to situations where the intercepts appear to be approaching −∞ or ∞.

4.
3. An example with the polymirt() function polymirt() and confmirt() both estimate the model parameters with the MH-RM algorithm. polymirt() is applicable for exploratory item analysis for dichotomous and polytomous data, and the object returned has many commonalities with objects returned by mirt(). There are a few pros and cons to using these stochastic functions, the pros being that parameter standard errors are automatically computed as a by-product of estimation, the models stay accurate even with higher dimensionality, lower asymptotes may be estimated (with β priors automatically added), and in confmirt() various item constraints can be imposed. The cons are that the time to estimate these models will be longer than mirt() or bfactor() for low-dimensional models since the actual estimation of the parameters takes more time, computation of such useful statistics as the observed log-likelihood must be estimated by Monte Carlo methods, and the parameters will vary slightly between independent estimations. However the added estimation time is not a major concern since the overall execution times often fall well within reasonable limits (see below).
Specification of a three-dimensional exploratory factor analysis model for the first simulated data-set is The behavior of summary(), anova(), residuals(), and plot() function the same as before, with this addition of logLik() that computes a Monte Carlo estimated integral for the observed log-likelihood. By default, the log-likelihood is approximated with 3000 draws at the end of both estimation methods, but can be suppressed by specifying calcLL = FALSE in the function calls.

An example with the confmirt() function
For this example we assume that the form of the loadings, and the relationships among the factors, are known or suspected a priori. Here we will try to recover the parameters used to simulate the data (which can be found in Appendix B). To begin, we must declare where the factors load, the relationships among the loadings, the relationships among the factors, as well as any additional parameter constraints. A model is specified by indicating which latent factors affect which numerically labeled item and by utilizing a select few keywords (e.g., COV, MEAN, INT, SLOPE, etc.) for additional parameter relations. For example, the following code declares a two-factor confirmatory model where the first factor (F1) affects items 1 to 4 while the second factor (F2) affects items 4 to 8 and the COV option allows the covariance between F1 and F2 to be freely estimated. Comparison of these two models suggests that the added restrictions do not significantly make the model fit any worse than the less restricted one.

Program comparisons
As is useful with all new software, comparing results with previously established programs to check the accuracy and potential benefits of the new software is beneficial for front-end users.
Here we compare the estimation results of mirt() and polymirt() with those obtained from TESTFACT, MCMCpack using MCMCirtKd() 3 , and ltm using ltm() (however, ltm() cannot estimate more than two factors). Two-, three-, and four-factor models are extracted from the first simulated data set, with the three-factor solutions rotated with the varimax criterion that was used for subsequent comparisons. Note that all computations were performed on a desktop workstation with an AMD Phenom 9600 Quad-Core 2.31 GHz processor with 4-GB of RAM, and each subroutine was run five times to obtain the average computation time.
Deterministic methods were set to terminate when all parameters changed less than 0.001 between consecutive iterations, and the number of quadratures used during estimation were 20, 10, and 7, respectively. polymirt() was set to have 100 burn-ins, 50 cycles to find approximate starting values, and once the RM stage was implemented, the estimation was terminated when all parameters changed less than 0.001 between iterations on three consecutive cycles. For MCMCirtKd(), the burn-in iterations were set at 10000, with 25000 MCMC draws, thinning every 5 draws, and storing only the item parameters. Finally, for the stochastic algorithms, the first model estimated was selected for subsequent comparisons.
As can be seen in Table 1, mirt() and polymirt() were much more efficient compared to the ltm() and MCMCirtKd() functions, and while mirt() was consistently more efficient than TESTFACT, polymirt() did not become more efficient than TESTFACT until there were more than 2 factors. Also note that the estimation time for MCMCirtKd() was quite long, spanning between 35-41 minutes on average per model.

Discussion
Several useful applications of the mirt package are possible that were not demonstrated in this article. For instance, confmirt() can estimate Rasch-type models stochastically by simply constraining all of the slope parameters to be equal (or exactly to 1/1.702, for the traditional Rasch model). This may be a useful strategy if the number of participants is large or the number of test items is large, since the MH-RM is well equipped to handle both of these situations. Non-linear factor combination and noncompensetory item response relationships may also be included when specifying confmirt() models. Finally, factor scores and information plots (and surfaces) for individual items or the whole test are available, and simulation functions are also readily at the user's disposal (see Appendix A).
As it stands, the mirt package appears to be a useful tool for researchers interested in exploratory and confirmatory item response models, and improves upon the overall estimation time and parameter recovery efficacy compared to various previously published software. The package is actively being developed, and some of the future developments may include: adding limited-information model fit statistics providing standard errors for EM solutions performing multiple-group estimation, and utilizing nominal and rating scale intercept methods for polytomous data These are only a few of the potential development areas, and user interest will largely guide which features will be developed. Popular options that will be available in the IRTPRO software (Cai, du Toit, and Thissen 2011) also may be given precedence depending on user feedback and interests in using open-source software along with proprietary software in their item analysis work.

A. Additional features
Additional features in plot(), itemplot() and fscores() are illustrated in the following.
Unidimensional and two-factor test information plots for the LSAT7 data can be produced with the plot() method (see Figure 1).
itemplot(mod1, combine = 5, auto.key = list(space = 'right')) The following code returns either a table of EAP or MAP factor scores of each unique response pattern, or appends the factor scores to the last column of the input data matrix for each response pattern.