BUGS Code for Item Response Theory

I present BUGS code to fit common models from item response theory (IRT), such as the two parameter logistic model, three parameter logisitic model, graded response model, generalized partial credit model, testlet model, and generalized testlet models. I demonstrate how the code in this article can easily be extended to fit more complicated IRT models, when the data at hand require a more sophisticated approach. Specifically, I describe modifications to the BUGS code that accommodate longitudinal item response data.


Introduction
In this paper, I present BUGS (Gilks, Thomas, and Spiegelhalter 1994) code to fit several models from item response theory (IRT). Several different software packages are available for fitting IRT models. These programs include packages from Scientific Software International (du Toit 2003), such as PARSCALE , BILOG-MG (Zimowski, Muraki, Mislevy, and Bock 2005), MULTILOG (Thissen, Chen, and Bock 2003), and TESTFACT (Wood, Wilson, Gibbons, Schilling, Muraki, and Bock 2003). The Comprehensive R Archive Network (CRAN) task view "Psychometric Models and Methods" (Mair and Hatzinger 2010) contains a description of many different R packages that can be used to fit IRT models in the R computing environment (R Development Core Team 2010). Among these R packages are ltm (Rizopoulos 2006) and gpcm (Johnson 2007), which contain several functions to fit IRT models using marginal maximum likelihood methods, and eRm (Mair and Hatzinger 2007), which contains functions to fit several variations of the Rasch model (Fischer and Molenaar 1995). Volume 20 of the Journal of Statistical Software is devoted to "Psychometrics in R" (de Leeuw and Mair 2007) and contains articles on how to fit a multilevel Rasch model with the lme4 package (Doran, Bates, Bliese, and Dowling 2007;Bates and Maechler 2010), how to fit multilevel polytomous item response models using Markov chain Monte Carlo (MCMC) methods (Fox 2007), and how to fit item response models that account for response times (Fox, Entink, and van der Linden 2007). The SCORIGHT software (Wang, Bradlow, and Wainer 2005) uses MCMC methods to estimate parameters in testlet models.
The use of BUGS software to estimate IRT models, however, allows the user to alter existing code to fit new variations of current models that cannot be fit in existing software packages. For example, longitudinal or multilevel data can easily be accommodated by small changes to existing BUGS code. The BUGS software takes care of the "grunt work" involved in estimating model parameters by constructing an MCMC algorithm to sample from the posterior distribution. As one anonymous reviewer stated it, in BUGS "the user does not have to bother thinking about how extensions of a model can be estimated". Thus, using BUGS frees the user to experiment with different models that may be more appropriate for specialized data than the models that can currently be fit in other software packages.
Of course, more complicated models involve more parameters than simpler models, and the analyst must specify prior distributions for these new parameters. This is a small price to pay, however, for the flexibility that the Bayesian framework and BUGS software provide.
Throughout this article, I assume that the reader has some basic understanding of IRT models and working knowledge of a software implementation of the BUGS language. However, if this is not the case, I give some references in Section 2 to sources that discuss IRT models and references to sources that contain tutorials and descriptions of BUGS.

Item response theory
I assume that the reader has working knowledge of basic IRT models; however, to establish notation, I briefly discuss each IRT model for which I provide code. The reader is encouraged to consult other sources for more detailed descriptions of the models discussed here. Excellent sources for learning IRT are Baker and Kim (2004), who provide a mathematically detailed introduction to IRT; Hambleton, Swaminathan, and Rogers (1991), who give an intuitive introduction to the topic; and Wainer, Bradlow, and Wang (2007), who provide an introduction to testlet models.

Two parameter logistic model
The 2PLM is used for data collected on n individuals who have each given responses on p different items. The items have binary outcomes, i.e., the items are scored as 1 if correct and 0 if not. The i-th individual in the sample is assumed to have a latent ability θ i , and the i-th individual's response on the j − th item is a random variable Y ij with a Bernoulli distribution. The probability that the i − th individual correctly answers the j − th item (i.e., the probability that Y ij = 1) is assumed to have the following form where α j is called the discrimination parameter and δ j is called the difficulty parameter for where logit(x) = log x 1−x . Each latent ability θ i is assumed to come from a standard normal distribution. Additionally, in a Bayesian analysis, item parameters are given prior distributions where m α , s α , m δ , and s δ are constants specified before the data analysis, N m, s 2 denotes a normal distribution with mean m and variance s 2 , and N + m, s 2 denotes the normal distribution truncated to the positive real line. Table 1 contains BUGS code to fit the 2PLM. In general, the code speaks for itself; however, I list a few comments below that may clarify some aspects of the code.
Line 5 uses the logit function to specify the logistic ogive as the link function for the model. A normal ogive could be used by changing line 5 to The BUGS language parametrizes the normal distribution in terms of the precision-the inverse of the variance. Thus, standard deviations for prior distributions on the item parameters need to be converted to precisions in lines 14 and 15.
A truncated normal distribution is specified for the discrimination parameters in line 12 by using the I(0, ) operator. Strictly speaking, in WinBUGS, there are no truncated distributions; the I(.,.) operator is used only to denote censored observations (Spiegelhalter et al. 2003). However, when all parameters in a prior distribution are observed, the I(.,.) operator will mimic the behavior of a truncated distribution (Lunn et al. 2009, p. 3061).
JAGS and OpenBUGS remove this ambiguity between truncation and censoring by introducing the truncation operator T(.,.). The truncation operator, however, is not available for all distributions in OpenBUGS, and it is unclear from the JAGS manual whether the truncation operator can be used on all distributions implemented in JAGS.
OpenBUGS still accepts the I(.,.) operator, but JAGS does not. Therefore, if the code in Table 1 is to run in JAGS, line 12 should be changed to alpha[j]~dnorm(m.alpha, pr.alpha) T(0, ) Some authors have used a log-normal distribution as a prior for the discrimination parameters α j (e.g., Patz and Junker 1999). The log normal distribution can be specified for the discrimination parameters in BUGS with the code alpha[j]~dlnorm(m.alpha, pr.alpha) Be aware, however, that the log-normal distribution has two parameters and that the mean and variance of the log-normal distribution are functions of both parameters. Therefore, changing only one parameter (e.g., m.alpha in the previous code) in the log-normal distribution will change both the mean and the variance of the log-normal distribution, which makes prior specification tricky.
The Rasch model is a special case of the 2PLM where each discrimination parameter alpha[j] is set equal to one (see, for example, Baker and Kim 2004, Chapter 5

Three parameter logistic model
The 3PLM is often used to analyze data from multiple choice tests where subjects try to choose the correct answer from a list of possible answers and may end up choosing the correct answer just by chance. The 3PLM is similar to the 2PLM except that the probability that the i-th individual will respond positively to the j-th item is dependent on a parameter η j that is constrained to lie in the unit interval: The parameter η j is sometimes called a "guessing" parameter because it represents the probability that an individual of extremely low ability could guess the correct answer on the j-th item just by chance.  Equation (3) can also be written as Not all items are required to have a guessing parameter. For items with no guessing parameter, η j ≡ 0.0. For items with a guessing parameter, the η j parameters are assigned beta prior distributions η j ∼ Beta(a η , b η ) .
All remaining model parameters are assigned priors as in the 2PLM of Section 3.  Table 2. However, specifying each item probability by hand leaves the code more cluttered than in Table 2 and does not allow the code to be easily reused on other data sets.

Graded response model
As with the 2PLM and 3PLM, the GRM is used for data collected on n subjects who have responded to each of p items. However, each item can have more than 2 ordered response categories. Thus, the response Y ij of the i-th individual to the j-th item can take values in the set {1, . . . , K j }, where K j is the largest category of the j-th item. The probability that the i-th subject will select the k-th category on the j-th item is constructed by first considering the cumulative probabilities where κ jk is a threshold, and F L (·) is the CDF of the logistic distribution. Each item has K j − 1 thresholds κ j1 , . . . , κ j,K j −1 that must satisfy the order constraint κ j1 < · · · < κ j,K j −1 .
The probability p ijk that the i-th subject will select the k-th category on item j can now be written as Priors on item parameters α j and δ j are the same as the priors for the 2PLM. The priors for the threshold parameters must account for the order constraint κ j1 < · · · < κ j,K j −1 . A prior on the threshold parameters can be induced by defining unconstrained auxiliary parameters κ * j1 , . . . , κ * j,K j −1 such that κ * jk ∼ N m κ , s 2 κ for k = 1, . . . , K j − 1. Prior distributions on the thresholds for the j-th item are obtained by setting κ jk equal to the k-th order statistic of the auxiliary variables κ * j,1 , . . . , κ * j,K j −1 for the . . .
where κ j, [k] denotes the k-th order statistic of κ * j1 , . . . , κ * j,K j −1 . This approach to modeling thresholds is recommended by Plummer (2010a, p. 36). Table 3 contains BUGS code to fit the GRM. Below are some comments on the code.
Because the responses are no longer have binary, the dcat function must be used to specify the distribution of the data. The dcat function defines a categorical distribution with more than two categories. The dcat function requires that the data for item j are coded as 1, 2, ..., K[j], i.e., the data cannot contain any zeros.
The ranked function on line 31 returns the k-th smallest value in the vector kappa.star[j, Some tests may not have the same numbers of categories for all items. Therefore, not all values in the matrix kappa[,] will be defined. WinBUGS and OpenBUGS run without any problems with these undefined parameters in the GRM. However, JAGS crashes when it tries to set monitors for the defined kappa parameters in the presence of undefined kappa parameters. The GRM can be fit in JAGS by setting the undefined item-step parameters to some arbitrary fixed value (e.g., zero) by passing those values to JAGS as data. For example, if the numbers of categories for a 5 item test are 2, 3, 4, 4, and 4, then a 5 × 3 matrix of the following format should be passed to JAGS as data for the matrix kappa.
The JAGS implementation of BUGS does not have the ranked function used on line 31; instead, JAGS has a sort function. The sort function in JAGS takes a vector as its input and returns the vector sorted in ascending order. Thus, lines 28-33 in Table 3 should be replaced with for the code to work in JAGS.

Generalized partial credit model
The GPCM (Muraki 1992) is an alternative to the GRM for ordinal item responses. In the GPCM, the probability that the i-th individual selects category k on item j is where β j1 ≡ 0 for all j. The parameters β jk are often called "item-step" parameters. Unlike the threshold parameters in the GRM, the item-step parameters do not need to satisfy any order restrictions and can be given normal prior distributions with no additional constraints . . , p and = 2, . . . , K j .
Even though, mathematically, the item-step parameters do not need to satisfy order constraints, if parameter estimates of item-step parameters for a certain item do not satisfy the ordering β j1 < · · · < β jK j , then it might be more appropriate to model the offending item with less categories (see, for example, Reckase 2009, page 34).
The discrimination parameters α j are assigned a truncated normal prior distribution as in the previous models. As usual, the latent abilities θ i are assumed to follow a standard normal distribution. Table 4 contains BUGS code for the GPCM. Below are some comments on this code.
As with the threshold parameters in the GRM, not all of the item-step parameters are defined when the items have differing numbers of response categories. WinBUGS and OpenBUGS run without problems in spite of the undefined parameters. JAGS, however, crashes with an error when running this model. The GPCM can be fit in JAGS by setting the undefined item-step parameters to some arbitrary fixed value (e.g., zero) by passing those values to JAGS as data. For example, if the data contain 5 items with numbers of categories 2, 2, 3, 3, and 4, then a 5 × 4 matrix with the following structure should be passed to JAGS as data for the matrix beta [,].
The partial credit model (Masters 1982, PCM) is a special case of the GPCM where each discrimination parameter alpha[j] is set equal to one. The code in Table 4 can be adapted to fit the PCM by changing line 21 to read or by removing all references to alpha[] from the code.

Testlet model
The testlet model is used for tests that are structured into groups of items that share some common feature. These groups of items are called "testlets". For example, many tests require a subject to read a certain passage and then answer two or more questions about the passage. Responses to items from the same testlet will tend to be more highly correlated than items from different testlets after accounting for the latent ability.
In a testlet model for binary responses, the probability that the i-th subject answers the j-th item correctly is assumed to have the following form where i = 1, . . . , n and j = 1, . . . , p.
Each subject has n T testlet effects γ id(j) , which can be thought of as testlet-specific abilities. The testlet effects explicitly model the correlation among items in a testlet after accounting for the latent ability. The function d(·) maps values in {1, . . . , p} to {0, 1, . . . , n T }. In other words, the function d(·) denotes which items belong to which testlets, so if d(5) = 3, then the 5-th item belongs to the 3 rd testlet. When d(j) = 0, the j-th item does not belong to any testlet and, therefore, has a testlet effect equal to zero for all subjects, i.e., γ i0 ≡ 0 for i = 1, . . . , n.
Testlet-effect parameters are assumed to come from normal distributions, and each testlet is given its own testlet-specific variance σ 2 The testlet variances are assigned inverse-gamma prior distributions Priors for all other model parameters are assigned as in the 2PLM. Table 5 contains BUGS code for the testlet model. Below are some comments on this code.
As with the BUGS code for other models, the terms n and p denote the number of subjects and items, respectively, and must be passed to BUGS as data. Additionally, the user must also specify the number of testlets n.t and the testlet-identifier vector The specification of d[j] in the BUGS code differs somewhat from the mathematical formulation outlined in the beginning of this section. BUGS does not allow subscripts equal to 0, thus, we cannot use gamma[i, 0] <-0.0 to set γ i0 ≡ 0. Instead, we let gamma[i, n.t + 1] <-0.0. Users must account for this specification in their data by assigning a value of n.t + 1 to d[j] if the j-th item does not belong to a testlet. For example, in a ten item test, if items 1-4 are part of the first testlet, item 5 is not a part of any testlet, items 6-9 are part of the second testlet, and item 10 is not a part of any testlet, then n.t=2 and the vector d[] should be d = c(1, 1, 1, 1, 3, 2, 2, 2, 2, 3) The BUGS language has no direct method of specifying an inverse-gamma prior distribution. The typical method of specifying an inverse-gamma prior for a parameter is to specify a gamma prior on the inverse of the parameter. This "trick" is used in lines 24 and 25 to specify an inverse-gamma prior on the sigsq.gamma[k] parameters.
Inverse-gamma priors have typically been used on variance parameters for reasons of convenience. The inverse-gamma prior is the conjugate prior for a variance parameter from a normal likelihood, so the update in an MCMC algorithm is a simple random draw from an inverse-gamma distribution. However, since BUGS frees the user from having to directly code the MCMC algorithm, other priors for variance parameters may be preferable. Gelman and Hill (2007) (2006) for further discussion on priors for variance parameters in hierarchical models.
In this section, I have defined the testlet model only for binary responses. However, the code in Table 5 can be adapted for other models (e.g., the GRM) by combining the code in Table 5 with code for models in the earlier sections. As mentioned previously, this ability to quickly adapt BUGS for new models is one of the major advantages to using BUGS.

Generalized testlet model
In the testlet model, the testlet effects γ id(j) are forced to have the same discrimination parameter α j as the latent ability. The generalized testlet model relaxes this restriction by introducing p new discrimination parameters α 21 , . . . , α 2p for the testlet effects. In the generalized testlet model, the probability that the i-th subject correctly answers the j-th item is modeled as .
Both sets of discrimination parameters are given truncated normal prior distributions and the difficulty parameters ζ j are given normal prior distributions  The ability parameters are assumed to come from a standard normal distribution. As in the testlet model, the testlet effects γ d(j) are assumed to come from normal distributions. However, for identifiability, the testlet effects are no longer permitted to have testlet specific variances. The variances of the testlet effects are restricted to be one γ id(j) ∼ N(0, 1) . Table 6 contains BUGS code for the generalized testlet model. The code uses the same "trick" to define a zero testlet effect as in the testlet model, so the comments of Section 7 apply to the code here as well.
9. Extending the BUGS code In many studies, subjects are administered the same set of test items on multiple occasions. Thus, instead of having a single vector of responses, the i-th individual has a matrix of responses composed of T response vectors one response vector Y it at each time point t, and an ability vector where parameter θ it is the ability of the i-th subject at time t.
Because the ability parameters θ i1 , . . . , θ iT are measured on the same subject over a period of time, latent abilities on the same subject will be more correlated than latent abilities among different subjects. A statistical model for this scenario should account for this correlation structure.

Population-averaged covariance models
In longitudinal data analysis, there are two main modeling strategies for modeling the covariance of repeated observations (Davidian 2005). The first strategy is called "population averaged", which, in the context of item response theory, involves assuming a covariance structure directly on the population distribution of all ability vectors θ i . Ability vectors θ i are assumed to come from a normal distribution with some mean structure µ θ and covariance Σ θ . Popular structures for the covariance are the "uniform" or "compound symmetric" covariance and the Markov structure for unequally spaced time points where σ 2 > 0, −1 < ρ < 1, d ij = |t i − t j |, and t k is the time at the k-th observation.
Each of the above models implies a constant variance of the latent ability across time. Less restrictive versions of each of the above covariance models can be obtained by allowing the variance to be different at each time point. These "heterogeneous" covariance structures can be modeled by pre and post multiplying the correlation matrices R in the above models by a diagonal matrix with unique standard deviations along the diagonal.
For identifiability, the mean and variance at one time point must be specified to establish the location and scale of the latent ability distribution (Tavares and Andrade 2006, pp. 105, 106). Setting the first element of µ θ to zero and the first variance σ 2 1 in Σ θ to one is one possibility. Analogous to the univariate normal distribution, the multivariate normal distribution in BUGS is parametrized in terms of the precision matrix, which is the inverse of the covariance matrix. Therefore, the covariance matrix Sigma.theta[,] must be inverted in line 30. The resulting precision matrix Pr.theta [,] is then used in the dmnorm function.

BUGS code for AR(1) models
The syntax for matrix operations is more flexible in JAGS than in WinBUGS or Open-BUGS. For example, line 30 in Table 7 can be replaced with the following cleaner code in JAGS

Pr.theta <-inverse(Sigma.theta)
Because the variance of the latent ability is restricted to be constant across time and the variance must be specified exactly for at least one time period to establish identifiability, sigsq.theta is set equal to 1.0 in line 21.
WinBUGS and OpenBUGS do not experience problems when fitting the AR(1) model. However, the default MCMC algorithm constructed by JAGS does not reach convergence. Table 8 contains BUGS code for the AR(1) covariance structure with heterogeneous variances. Below are some comments on this code. The only major substantive difference between the code in Table 7 and in Table 8 is that sigsq.theta[] is now a vector where the first element of the vector is set to one and the remaining elements are given gamma prior distributions. There is also a slight difference in how the heterogeneous AR(1) structure is constructed (see line 27 of Table 8) relative to the AR(1) structure with constant variance (see line 26 in Table 7). Line 27 of Table 8 implements the pre and post multiplication of the correlation matrix R by the standard deviation matrix D θ .

BUGS code for unstructured covariance
Researchers may be reluctant to specify a particular structure for the covariance matrix of the latent abilities. Wishart distributions are typically used to specify priors for unstructured covariance matrices. However, a Wishart prior cannot be used for the covariance of latent abilities, because at least one variance must be set to a constant for identifiability.
One approach to solve this problem is to parametrize the covariance matrix Σ θ in terms of its Cholesky decomposition where L θ is a lower triangular matrix with positive entries on the diagonal and unrestricted entries below the diagonal. Setting the first element of L θ equal to one ensures that the first   element of Σ θ is also equal to one. Putting priors on the elements of L θ is straightforward, because the only elements that must satisfy restrictions are the diagonal elements, which need only be positive. Table 9 contains BUGS code for fitting a model with an unstructured covariance for the population distribution of the ability vectors θ i . Below are some comments on this code.
Gamma priors are used to restrict the diagonal elements of L.theta to be positive, but other similar priors may be used (e.g., dlnorm).
WinBUGS and OpenBUGS do not have a built-in matrix multiplication function. Therefore, the matrix multiplication in Equation (5) is performed explicitly using loops and the inprod function. The inprod function computes the inner product of two vectors.
JAGS offers more support for vector and matrix calculations than WinBUGS and Open-BUGS. In JAGS, lines 29-33 can be replaced by the single line Sigma.theta <-L.theta %*% t(L.theta) In general, the mixing of the Markov chains for this model is poor in all of the BUGS packages, and the chains must be run for a long time (several hundred thousand iterations) to achieve a good sample from the posterior.

Subject-specific covariance models
The other modeling strategy for longitudinal data is sometimes called "subject specific", which, in the context of item response theory, involves assuming that the ability θ it at time t is some function (usually a linear combination) of random coefficients and time. For example, one subject-specific model is a linear function over time In this model, each individual in the data is assumed to have two "random coefficients" that determine the individual's linear trajectory of ability over time.
Like the population averaged models, the subject-specific models must incorporate restrictions on the location and scale of the latent abilities to establish identifiability. In populationaveraged models, the mean and variance of the ability population at some time point t can be set to constants-usually zero and one, respectively. In subject-specific models, the mean and variance of the population of one of the random coefficients can be set to constants. For example, in the linear, subject-specific model specified above, setting µ γ 0 = 0 and σ 2 γ 0 = 1 establishes the location and scale of the latent abilities by setting the distribution of the abilities at t = 0 to be equal to the standard normal distribution.

Example
In this section I demonstrate how to use the BUGS code in the previous sections to analyze a real data set. I use the Environment data set from the R package ltm. The data are a subset of responses from the 1990 British Social Attitudes Survey and contain responses of 291 individuals to six questions involving environmental issues. The possible answers to each question are "not very concerned", "slightly concerned", and "very concerned". The ordinal nature of the responses implies that the GRM or the GPCM should be used to analyze these data. In what follows, I use the GRM.
It is often convenient to call WinBUGS (or OpenBUGS, or JAGS) from some other dataprocessing software as this allows the user to compute posterior summaries and create plots of posterior draws that are not available directly in BUGS software. For this example, I use the openbugs function in the R2WinBUGS package to call OpenBUGS from R.
In the data analysis of this section, I follow the approach of Rizopoulos (2006) by fitting two versions of the GRM-a constrained model, where the discrimination parameters are forced to be equal (α 1 = · · · = α 6 = α) for each of the six survey questions, and an unconstrained model, where all six discrimination parameters are freely estimated. The constrained model can be fit by removing the subscript [j] from the parameter alpha in line 10 of and removing the prior on line 24 of Table 3 from the loop and the subscript [j] from the parameter alpha as in alpha~dnorm(m.alpha, pr.alpha) I(0, ) The relevant R packages and the data can be loaded in the R session with the following code: R> library("R2WinBUGS") R> library("ltm") R> library("mcmcplots") R> data("Environment") The data set Environment is a data.frame with six variables each of class factor. Because the data are inherently ordinal, the variables in the Environment data set must first be converted to class ordered by calling the ordered function on each variable and setting the levels attribute. This can be done by using the lapply function to apply the ordered function with argument levels = c("not very concerned", "slightly concerned", "very concerned") to each of the variables in Environment.
Additionally, in order to be processed by OpenBUGS, the data must be coerced to a matrix of numeric values that correspond to each level of the ordinal response: 1 → "not very concerned" 2 → "slightly concerned" 3 → "very concerned" The above operations can be executed in one line of R code: R> Y <-data.matrix(data.frame(lapply(Environment, ordered, levels = + c("not very concerned", "slightly concerned", "very concerned")))) Several constants-sample size (n), number of survey questions (p), number of response categories per question (K), and parameters for prior distributions (m.alpha, s.alpha, m.kappa, and s.kappa)-need to be passed to OpenBUGS as data. These constants are first stored as variables in the R session, then the names of these variables and the data matrix (Y) are stored in the character vector data, which will be passed to the openbugs function.
R> monitor <-c("alpha", "theta", "kappa") R> n.burn <-4000 R> n.thin <-10 R> n.sim <-500 * n.thin + n.burn I use a burn-in period of 4,000 iterations to ensure that OpenBUGS can complete the adaptive phase of the MCMC algorithm. The number of parallel chains is kept at the default value of 3.
The openbugs function can be used to call OpenBUGS and generate draws from the posterior distribution of the parameters for the constrained model. The BUGS model file is contained in the grmeq.bug file of the bugs subdirectory of the supplemental material. The model.file option in the code below will need to be modified if the working directory has not been set to the base directory of the supplemental material.
R> eq.alpha.out <-openbugs(data = data, inits = NULL, + parameters.to.save = monitor, + model.file = file.path(getwd(), "bugs/grmeq.bug"), + n.iter = n.sim, n.thin = n.thin, n.burnin = n.burn)) The openbugs function returns an object of class bugs. One way to access posterior draws in a bugs object is via the sims.matrix component of the bugs object. The sims.matrix component is a matrix object with number of rows equal to the number of saved iterations of the MCMC simulation and number of columns equal to the number of monitored parameters. Posterior draws of a single parameter can be accessed by passing the name of the parameter to the subset function for matrices, as in eq.alpha.out$sims.matrix [, "theta[162]"] or eq.alpha.out$sims.matrix[, "alpha"] The grep function provides a convenient way to access groups of parameters. For example, the code R> parnames <-colnames(eq.alpha.out$sims.matrix) R> eq.alpha.out$sims.matrix [, grep("theta", parnames)] can be used to return a matrix of the posterior draws for all the ability parameters θ 1 , . . . , θ 291 .
A call to the generic function plot produces some summary plots for the MCMC simulation.
R> plot(eq.alpha.out) Among the several summary plots is a plot of the Gelman-Rubin (GR) convergence diagnostics (Gelman and Rubin 1992) for a subset of the parameters. The GR diagnostics are all close to 1, which gives some assurance that the chains have converged.
Additionally, standard, MCMC diagnostic plots (such as trace and autocorrelation plots) can be created in an HTML file and viewed in a browser with a call to the mcmcplot function in the mcmcplots package (Curtis 2010).
R> mcmcplot(eq.alpha.out, random = 20) The fit of this model can be appraised using posterior predictive model checks (see, for example, Gelman, Carlin, Stern, and Rubin 2003, Chapter 6). This method of model checking involves simulating several new data sets using the posterior distribution of the model parameters. Summary statistics from the simulated data sets are compared with the same summary statistics from the observed data. If major discrepancies exist between the distribution of the summary statistics from simulated data and the summary statistics from the observed data, then the model is deemed to have poor fit.
In many latent variable procedures, the goal of model fitting is to find values of model parameters that give a model implied correlation or covariance matrix that is close to the raw correlation matrix computed from the sample data (see, for example, Bollen 1989). This idea can be used in the posterior predictive model checking procedure in the current example by comparing correlation matrices computed on simulated data sets with the correlation matrix from the observed data Y. Because the data in this example are ordinal, a nonparametric measure of correlation (e.g., Kendall's tau) must be used. (b) Compute Kendall's correlation matrix for the new data set generated in the previous step, and save the values for later analysis.
2. Compute Kendall's correlation matrix for the raw data.
3. Compare the distribution of simulated values of Kendall's correlation to the observed values of Kendall's correlation.
The details of the procedure are implemented in the function ppktau contained in the functions.R file of the supplemental material to this article. The function ppktau requires the plyr package (Wickham 2009). In this example, the ppktau function returns a matrix with M rows (where M = 1, 500 is the number of posterior simulations-500 iterations from 3 parallel chains) and 15 columns for the p(p − 1)/2 correlations among the 6 variables. Appropriate column names are given to the posterior predictive output via the makeNames utility function (also included in the functions.R file). Kendall's tau values for the real data are computed with the cor function using the method="kendall" option.
R> pp.eq <-ppktau(eq.alpha.out$sims.matrix) R> colnames(pp.eq) <-makeNames("ktau", 1:p, 1:p, + symmetric.matrix = TRUE, diag = FALSE) R> ktau <-cor(Y, method = "kendall")[lower.tri(diag(p))] The simulated values of Kendall's tau can be compared to the observed values of Kendall's tau by producing density plots of the simulations and marking on the density plots the corresponding value of the observed Kendall's tau. The function plotpppval (in functions.R) automates the process of producing density plots for each column of the simulated output. Additionally, the plotpppval function shades the area under the density from the observed Kendall's tau to the nearest tail. This tail area is known as the posterior predictive p value or Bayesian p value (Rubin 1984;Meng 1994;Gelman, Meng, and Stern 1996). The Bayesian p values are included on the left-hand side of the density plots produced by plotpppval.
R> plotpppval(pp.eq, ktau) The caterplot function in the mcmcplots package can be used to create "caterpillar" plots of 95% predictive intervals for the simulated Kendall's correlations in pp.eq. The caterpoints function can be used to overlay the observed values of the Kendall's tau on the same plot. Observed values of Kendall's tau that fall outside the prediction intervals indicate lack of fit.
The results of the posterior predictive model check provide some assurance that the unre-

Summary
In this article, I have presented several snippets of BUGS code to fit many of the common IRT models found in the literature. I have also shown how to extend these models to other modeling situations. I hope that researchers who use IRT models will find it useful to have BUGS code to fit these basic models in one convenient source and will find that the code can quickly be adapted for use in novel settings.