Introduction to Multilevel Item Response Theory Analysis

Why use IRT?

IRT is a probabilistic framework for the development and analysis of a multi-item instrument that has the potential to jointly shed light on items and person characteristics by conjointly linking item parameters and latent trait values on the same measurement scale. IRT models can be used for several aims, specifically to (a) choose items which effectively discriminate across individuals with a different magnitude of the latent trait; (b) select informative items to build up a briefer instrument and remove redundant items; (c) quantify the informative power of each item and an instrument to assess which segments of the latent trait are poorly measured; (d) measure individuals and item/item-categories (or item-categories for polytomous item response) location along the latent trait; (e) place individuals and item (or item-categories for polytomous data) characteristics on the same scale; (f) assess the minimum level of individual latent trait needed to likely provide a positive answer to an item or endorse a specific category of response to a certain item; (g) provide a measure of reliability that varies along the latent trait. Although a full description of IRT models is beyond the scope of this article and a listing of all the papers that have discussed this topic is impossible, interested learners are referred to De Ayala (2009), Edelen and Reeve (2007), Lord (1980), and Toland (2014) for further information.

Basic concepts and common dichotomous IRT models

Even though there are a variety of IRT models, for simplicity this study focuses on two IRT models for dichotomously scored item responses. The two most well-known IRT models for dichotomous item responses are the one-parameter logistic (1PL; Rasch, 1960) model and the two-parameter logistic (2PL; Birnbaum, 1968) model. The models are so called because of the number of item parameters each model contains. The 2PL model is the more general model and can be described by the following mathematical expression (Lord, 1980):

P ( y j i = 1 | θ j , b i , a i ) = e a i ( θ j − b i ) 1 + e a i ( θ j − b i ) ,

where P is the probability that a randomly chosen person j with latent trait value θ will answer item i correctly or endorse it. The values of b_i and a_i are parameters characterizing item i, and e is the mathematical constant 2.71828. The relationship between a correct item response and the latent trait can be modeled using a logistic (S-shaped) function known as an item response function (IRF). This function specifies that as the level of the latent trait (math) increases, the probability of a correct response on an item will increase nonlinearly.

The two item parameters, a_i and b_i, are called item “discrimination” (a) parameter and item “location” or “difficulty” (b) parameter. The discrimination parameter (a_i) provides information on the steepness of the logistic function signaling the power of an item to differentiate across individuals with different levels of the latent trait. The item location parameter b_i signals the position of the item alongside the latent trait; specifically the minimum amount of latent trait required to have a chance of 0.50 to answer it correctly. So, greater values of this parameter indicate that higher values of the latent trait are required to provide a positive answer to the item.

The values of the person parameter (θ), which identifies the individual position along the latent trait, and the item parameters are measured on the same metric. This links both parameters on the same continuum. If items have the same estimated discrimination power, then the parameter a_i becomes constant parameter (a) and the 2PL model reduces to the 1PL model. If the a parameter is set or fixed at 1 (not estimated), then the 1PL model is known as the Rasch model (Rasch, 1960).

Figure 1 plots four hypothetical IRFs. Items 1 and 2 depict items with equal discrimination as expected by a 1PL model, whereas Items 2 and 3 depict items with varying discrimination in accordance with a 2PL model. Notice that as one moves from low to high along the θ continuum, the probability to respond positively to an item increases, but is dependent on person location (θ)

OPEN IN VIEWER

Figure 1.

Example one-parameter (Rasch) and two-parameter logistic models item response functions for four hypothetical items dichotomously scored.

Statistical frameworks

The basic IRT model can be viewed from two different perspectives: namely as a multivariate regression model (e.g., multiple items can be viewed as a set of outcome variables) or as a multilevel regression model (e.g., the set of items can be viewed as repeated measurements on the same student). In the multivariate regression framework, the basic IRT models can be set as a single-factor measurement model with categorical responses (StataCorp, 2013). The relationships between the observed responses to each item and the latent factor are described by I logistic regression equations, each of which is a function of an item intercept and a loading. Thus, IRT models can be specified with software which fits confirmatory factor analysis with categorical items such as Mplus or the package for generalized structural equation models such as the Stata program SEM. Within this tutorial, we follow the multilevel regression approach, which is also the approach adopted in the Stata program GLLAMM. In both frameworks, IRT item parameters (a_i and b_i) are usually specified as fixed effects, whereas the person parameter θ, which is based on responses provided by the same person, can be specified as a fixed or random term, which is usually assumed to follow a normal distribution, N ( 0 , σ 2 ) . Readers interested in replicating the analysis within a multivariate regression model or SEM framework are invited to read the supplementary materials containing data and command files that use the Stata program SEM, Mplus, and flexMIRT.

Fixed- or random effects approach

The nesting of units in clusters (e.g., the same person responds to several items, or students belonging to the same class) is usually modeled by introducing fixed or random effects among the predictors of a regression model. In a fixed-effects approach, cluster affiliation is considered as a categorical covariate and a number of dummy variables, equal to the number of clusters minus 1 (i.e., a value of 1 is used for those units belonging to the same cluster on a given dummy variable and 0 elsewhere), are introduced as predictors. The effect of belonging to a cluster on the expected response is summarized by the related parameter estimate. As the number of clusters increases, the number of parameters that have to be estimated increases and the fixed-effects approach becomes inefficient. The random effects approach (more commonly known as a Multilevel approach) overcomes this issue by treating cluster effects as random terms (usually it is assumed that they are generated from a normal distribution) which assume the same values for units in the same cluster. The degree of similarity in responses of individuals in the same cluster (which are not explained by individual covariates) is considered by adding to the other predictors of the regression a random intercept which is shared by units belonging to the same cluster (e.g., responses of the same student or students in the same class). In a random effects framework, likely sources of dependence within higher level units can be easily taken into account by adding further random terms at higher levels of the data structure (e.g., classes in the same school). In this way, units belonging to the same cluster are treated as more alike than those which do not and therefore students in the same class and in the same schools share two common random terms. Moreover, a random effects formulation of the latent trait person parameters is recommended when a brief instrument is under analysis, when the number of individuals increases, or when the analysis focuses on person-level covariates. Ultimately, we adopt a random effects approach to specify the person parameter and to account for the nesting of items within students and students within classes.

Considering items as fixed parameters and adopting a random effects formulation of the person parameter, the most common IRT models for dichotomous and polytomous items can be classified as generalized linear or nonlinear mixed effects models (De Boeck & Wilson 2004; Fisher & Molenaar, 1995). This approach is fruitful as it can be used to extend IRT models in order to deal with a hierarchical or nested data structure, a latent variable that has a multidimensional structure (multidimensional IRT), and with covariates at person level (De Boeck & Wilson, 2004; Fisher & Molenaar, 1995; Pastor, 2003; Rijmen & Briggs, 2004).

Link between IRT and multilevel logistic model

At the turn of the 21st century, a number of pioneering papers (see, for example, Kamata, 2001; Singer, 1998) show that the class of Rasch models (i.e., where a is set at 1 and not estimated) for dichotomous and polytomous item response data can be parameterized as two-level logistic regression models. In this framework, a random intercept is specified at the person level (Level 2) and item indicator (dummy) variables are specified as predictors at the item response level (Level 1).

In order to adopt such an approach to estimate IRT parameters, the data have to be converted from wide (multivariate, person-level) format (see Appendix Table A1) to long (stacked, person-period) format (see Appendix Table A2) wherein the responses to the 25 items (m1-m25) of the 1,141 individuals are stacked in a unique variable m and each individual has a number of records/lines/rows equal to the number of items (see INVALSIMATHLONG.csv). The variable m (which has length or rows = 28,525 = 1,141 by 25 and assumes values 1 or 0) is used as the dependent variable in a multilevel logistic regression model containing as predictors 25 binary variables (x_i, for i in 1,…, 25; see Appendix Table A3). Each binary variable has a value 1 when the response in the m variable is associated to item M _i , otherwise 0. When the model is estimated, the related coefficient parameters of the x _i variables are item intercepts.

In the formulation of the IRT model as a two-level multilevel logistic regression model (Adams et al., 1997; Kamata, 2001), the person parameter θ _j is treated as a random term, which is made up of responses to items provided by the same person. This term captures the unexplained (residual) variability in item responses attributable to differences in person latent trait after controlling for item characteristics (i.e., item location parameters).

The Rasch model can be specified as a two-level multilevel logistic regression model with item responses (Level 1 units) nested within students (Level 2 units). In this model, item indicators are covariates, which are the x_i binary variables that identify responses to the same item and expressed as follows (Jeon & Rabe-Hesketh, 2012; Kamata, 2001):

P ( m i j = 1 | θ j , α i ) = e α i x i + θ j 1 + e α i x i + θ j .

In Equation 2, the probability that student j provides a correct response (m_ij = 1) to item i depends on an item intercept α _i , and a common latent trait (random term), which takes the same value for all the responses related to the same person θ j . Equation 2 can be obtained from Equation 1 by setting − α i = b i and a = 1 . The item intercept α _i can be interpreted as an easiness parameter, the larger its value the greater the probability to get a correct response. The item loading λ can be interpreted as a discrimination parameter, the larger the value the greater the relationship between the item and latent trait. Note, item loading is not estimated in this formulation as the value is fixed at 1 as is done in Rasch modeling. The model in Equation 2 can be estimated with any software which allows users to specify a generalized linear mixed effects model, such as Stata, Mplus, MLwiN, flexMIRT, SPSS, and SAS procedures PROC NLMIXED or PROC GLIMMIX.

The 2PL model in Equation 1 can be expressed in terms of an item intercept α i and loading λ i unique for each item as (Skrondal & Rabe-Hesketh, 2004)

P ( m i j = 1 | θ j , α i , λ i ) = e α i x i + λ i θ j 1 + e α i x i + λ i θ j ,

where α i = − a i b i and λ i = a i (McDonald, 1999). The item intercept α i can be expressed in terms of “location/difficulty parameter” as b i = − α i / λ i . ¹ As the discrimination-location form does not generalize well to the truly multilevel IRT model, we adopted the intercept-loading parameterization for this model and all remaining IRT models. Equation 3 is an IRT model with a unidimensional latent trait expressed in intercept-slope parametrization. The 2PL model is also a generalized nonlinear model since in Equation 3 the loading parameter is multiplied by the person parameter (Jeon & Rijmen, 2015).

To identify Equation 3, it is required to fix the variance of the random term θ j equal to 1 or a loading equal to 1. If we suppose that all items have the same discrimination power of 1, then Equation 3 simplifies to Equation 2.

The 2PL model can be expressed in the generalized nonlinear mixed model framework and estimated using routines for mixed models that allow varying loadings, such as Stata programs GLLAMM and SEM, the PROC NLMIXED procedure in SAS as illustrated by Jeon and Rabe-Hesketh (2012), Rijmen, Tuerlinckx, De Boeck, and Kuppens (2003), and StataCorp (2013). Another option is to use software which relies on Bayesian methods such as WinBUGS (Cho & Suh, 2012).

Why use MLIRT?

Within the MLIRT framework, we gain all the potential of IRT along with four additional features of being able to (a) assess how much an item discriminates between classes and within individuals; (b) assess how much of the variability in the responses is due to individuals and to classes; (c) estimate a value of the latent trait at class level; (d) use the estimates of the latent trait at class level as indicators (value-added measures) to measure effectiveness removing the effects of the heterogeneity in students’ characteristics (Goldstein & Spiegelhalter, 1996); (e) assess the amount of information (reliability) of the latent trait at both class and student levels.

Dichotomous MLIRT models

By extending the 2PL model in Equation 3, we can consider the nesting of students in the K classes to obtain a multilevel logistic regression model and call it a multilevel 2PL (MLM-2PL) model. Now, the MLM-2PL model can be used to express the probability to respond correctly to item i as

P ( m i j k = 1 | θ , α , λ ) = e α i x i + λ i ( s ) θ j k + λ i ( c ) θ . k 1 + e α i x i + λ i ( s ) θ j k + λ i ( c ) θ . k ,

where θ represents the vector of the latent trait values at the student level ( θ j k ~ N[0, σ θ ( s ) 2 ] ) and at the class level, ( θ 0 k ~ N [ 0 , σ θ ( c ) 2 ] ) . Also, λ is the vector containing the loadings at the student level ( λ i ( s ) ) and class level ( λ i ( c ) ) and as before α i is the item intercept. The student-level loading measures item i’s capability to differentiate across students (between-students within-classes) with different achievement in math conditional upon the class level of the latent trait, whereas the loading at the class level provides “value-added” information on item i’s capacity to discriminate between classes on the basis of students’ achievement in math. It is worth highlighting that in MLM, the sample size of units at each level of analysis (e.g., item responses [Level 1], students [Level 2], and classes [Level 3]) is relevant in order to assess the power of statistical tests related to fixed and random parameters. Moreover, the appropriate sample size at higher levels (Level 2 and Level 3) depends on the parameter(s) being tested (e.g., fixed vs. random) and the main focus of the research question(s) (Snijders, 2005). If the aim of an analysis is to test fixed effects, then the sample size of higher level units is not relevant. A guide and statistical software program for determining statistical power and sample size in MLM is provided by Browne, Golalizadeh Lahi, and Parker (2009). However, to our knowledge, there is not a program readily available for directly testing statistical power in the MLIRT context.

To identify the MLM-2PL model, a loading at each level of analysis (students and classes) has to be constrained to 1 and then variances of the related random terms are freely estimated. Or, the model can be identified by constraining the variances of the random terms to 1. In our application of the MLIRT models for comparing the relative contribution of each level of analysis to the overall variability and when comparing competing models we allowed the variances of the random terms to be freely estimated. When conducting MLIRT solely for measurement description purposes, the variances of the random terms were fixed to 1 to standardize the metric and improve interpretability.

If loadings λ i ( c ) at class level are supposed to have the same capacity to discriminate across classes, then all loadings are set equal to 1 at the class level and the MLM-2PL model reduces to the multilevel 2PL constrained (MLM-2PLC) model which can be expressed as

P ( m i j k = 1 | θ , α , λ i ( s ) ) = e α i x i + λ i ( s ) θ j k + θ . k 1 + e α i x i + λ i ( s ) θ j k + θ . k .

To test if the MLM-2PLC model works, a comparison needs to be made with the model described in Equation 6 using a LRT. If the specification of loadings at the class level do not improve the model fit, then it means that the loadings at the class level do not differ in a significant way from 1 and the loadings are the same for the i items. The model expressed in Equation 7 may possibly be realistic when classes perform differently at the instrument level, but not at the item level. Finally, if items have the same power to discriminate across students, then loadings at the student level are constant across items and equal to 1 (i.e., λ i ( s ) = λ ( s ) = 1 ). Thus, Equation 7 reduces to the multilevel Rasch (MLM-Rasch) model expressed as

P ( m i j k = 1 | θ , α ) = e α i x i + θ j k + θ . k 1 + e α i x i + θ j k + θ . k .

Equation 8 represents the probability of answering correctly an item as a function of two random terms, namely θ _jk and θ _k , which respectively capture the difference between the k class’s level of the latent trait and the overall mean, and the difference in math achievement of student j in class k with respect to the overall class k achievement. Classes with higher values of θ _k tend to have students with higher math achievement, whereas lower values of θ _k signal that students have on average lower levels of math achievement. The item intercept ( − α i = b i ) signals the difference in the probability of a correct answer due to a difference in item difficulty.

Notice that in the MLM-2PLC model and MLM-Rasch model that in order to assess the relative contribution of both levels to the overall variability, the variance of the random terms, which capture students and class variability, are estimated and the values of the loading parameter are fixed equal to 1. Note, Equation 5 differs from Equation 8 in that it has an overall intercept across all items, whereas Equation 8 allows each item to have its own unique intercept.

The MLM-2PL model in Equation 6 is estimated in the Stata program GLLAMM as a three-level model with loadings at student (Level 2) and class level (Level 3) which represent the within-class between-students loadings and the between-classes loadings (see Rabe-Hesketh et al., 2004b; Skrondal & Rabe-Hesketh, 2004).

Similar to traditional dichotomous IRT models, dichotomous MLIRT models make similar assumptions. One assumption is that the IRT model is modeling the correct number of latent traits that exist in the data at both student- and class level. A second assumption is conditional independence. This implies that conditional upon the value of the latent trait at student- and class-level responses to the item are due only to the latent trait(s) being modeled and nothing else. A third assumption is that the model appropriately represents the functional form of the data. Multilevel dimensionality can be assessed using exploratory and confirmatory factor analytic methods via Stata and Mplus. Unfortunately, direct tests of the second and third assumptions are not standard features available in Stata and Mplus (see Jiao, Kamata, Wang, & Jin, 2012; Ravand, 2015). Nonetheless, testing these assumptions is beyond the scope of this article since it would add unnecessary length. It is suggested that an extension of Yen’s (1984) Q3 index be used to examine the conditional independence assumption. Q3 is the correlation between the residuals for a pair of items where the residual is the difference between observed response (0/1) and the predicted response probability from the MLIRT model. So, the residual for student j in class k on item i is as follows: d_jik = y_jik – p_ijk(θ _ik θ _k ), whereas for the same person the residual for item m is d_jmk = y_jmk – p_jmk(θ _ik θ _k ). The closer Q3 is to zero the stronger the evidence that the assumption is tenable. Yen suggested using a cut value of 0.2. The standard has remained in place for more than 20 years, and although not perfect, continues to be shown as the optimal criterion for assessing this assumption with Yen’s Q3 index.

MLIRT results

Table 3 provides a summary of the characteristics of the five IRT models we focus on in this section. We investigated the advantages of fitting an MLIRT model in the analysis of students’ achievement in math by answering the following question: Are differences in mathematics achievement between classes relevant? In order to answer this question, several MLIRT models (MLM-Rasch, MLM-2PLC, and MLM-2PL models) were fit to each item on the mathematics achievement test and results compared with models that have similar proprieties, but ignore the multilevel data structure (Rasch, 2PL). Table 4 summarizes the item intercept and factor loading parameters for each of the five models and the model statistics. The models were identified by setting the first item loading to 1.

OPEN IN VIEWER

Table 3.

Characteristics of the Five Descriptive IRT Models.

Model characteristics	Rasch	2PL	MLM-Rasch	MLM-2PLC	MLM-2PL
Item parameter
Intercept	Yes	Yes	Yes	Yes	Yes
Factor loading
@ student level		Yes		Yes	Yes
@ class level					Yes
Random effect
@ student level	Yes	Yes	Yes	Yes	Yes
@ class level			Yes	Yes	Yes

Note. IRT = item response theory; 2PL = two-parameter logistic; MLM = multilevel modeling; MLM-2PLC = multilevel 2PL constrained; MLM-2PL = multilevel 2PL.

OPEN IN VIEWER

Table 4.

Summary of Item and Model Characteristics for Five IRT Model Solutions.

Item	Rasch		2PL				MLM-Rasch		MLM-2PLC				MLM-2PL
	α i	SE	α i	SE	λ i ( s ) λi(s)	SE	α i	SE	α i	SE	λ i ( s )	SE	α i	SE	λ i ( s )	SE	λ i ( c )	SE
1	1.95	0.09	2.10	0.13	1.00		1.93	0.13	2.06	0.15	1.00		2.13	0.181	1.00		1.00
2	0.69	0.08	0.60	0.07	0.43	0.07	0.68	0.12	0.64	0.12	0.40	0.09	0.58	0.080	0.42	0.09	0.369	0.10
3	−0.10	0.07	−0.09	0.06	0.30	0.06	−0.12	0.12	−0.10	0.12	−0.03	0.07	−0.10	0.099	0.00	0.07	0.617	0.13
4	1.85	0.09	1.87	0.11	0.85	0.13	1.82	0.13	1.84	0.14	0.78	0.14	1.84	0.143	0.80	0.15	0.758	0.18
5	0.31	0.07	0.32	0.07	0.77	0.11	0.29	0.12	0.29	0.12	0.59	0.11	0.30	0.128	0.59	0.11	0.853	0.16
6	1.27	0.08	1.32	0.09	0.89	0.12	1.25	0.13	1.27	0.13	0.80	0.14	1.29	0.141	0.77	0.13	0.865	0.18
7	−0.50	0.07	−0.50	0.08	0.99	0.13	−0.52	0.12	−0.51	0.12	0.67	0.12	−0.52	0.152	0.74	0.13	1.053	0.19
8	−0.08	0.07	−0.05	0.08	1.08	0.14	−0.10	0.12	−0.11	0.12	1.03	0.16	−0.10	0.136	1.05	0.16	0.893	0.17
9	−1.00	0.08	−1.04	0.09	1.04	0.13	−1.02	0.13	−1.09	0.13	1.01	0.16	−1.08	0.132	0.99	0.16	0.790	0.15
10	1.97	0.09	1.86	0.10	0.65	0.11	1.95	0.13	1.86	0.13	0.50	0.12	1.87	0.145	0.50	0.12	0.812	0.19
11	0.33	0.07	0.33	0.07	0.76	0.10	0.31	0.12	0.31	0.12	0.68	0.12	0.30	0.112	0.66	0.12	0.701	0.14
12	−0.15	0.07	−0.13	0.08	0.91	0.12	−0.17	0.12	−0.18	0.12	0.95	0.15	−0.17	0.119	0.93	0.15	0.741	0.14
13	0.05	0.07	0.08	0.08	0.95	0.12	0.03	0.12	0.04	0.12	0.88	0.14	0.04	0.129	0.88	0.14	0.847	0.17
14	−0.32	0.07	−0.28	0.06	0.42	0.07	−0.34	0.12	−0.30	0.12	0.09	0.07	−0.30	0.108	0.12	0.07	0.698	0.14
15	−0.14	0.07	−0.12	0.08	1.19	0.15	−0.16	0.12	−0.17	0.12	1.05	0.16	−0.16	0.148	1.11	0.17	0.995	0.18
16	−0.24	0.07	−0.21	0.07	0.75	0.10	−0.25	0.12	−0.25	0.12	0.78	0.13	−0.25	0.102	0.73	0.12	0.596	0.12
17	−0.70	0.08	−0.71	0.08	1.00	0.13	−0.72	0.13	−0.72	0.12	0.74	0.13	−0.73	0.155	0.76	0.13	1.066	0.19
18	−0.21	0.07	−0.19	0.09	1.30	0.16	−0.22	0.12	−0.22	0.12	0.74	0.12	−0.06	0.286	0.91	0.15	2.145	0.40
19	0.89	0.08	0.93	0.08	0.90	0.12	0.87	0.13	0.95	0.13	1.09	0.17	0.90	0.106	1.03	0.16	0.523	0.12
20	0.21	0.07	0.24	0.08	0.91	0.12	0.20	0.12	0.21	0.12	0.87	0.14	0.21	0.123	0.85	0.14	0.785	0.15
21	0.87	0.08	0.91	0.08	0.90	0.12	0.85	0.13	0.85	0.12	0.73	0.13	0.91	0.154	0.73	0.13	1.036	0.21
22	0.18	0.07	0.20	0.07	0.81	0.11	0.17	0.12	0.18	0.12	0.98	0.16	0.16	0.095	0.97	0.16	0.483	0.11
23	1.47	0.08	1.76	0.12	1.24	0.17	1.45	0.13	1.91	0.17	1.73	0.27	1.86	0.164	1.67	0.26	0.675	0.16
24	2.32	0.10	2.69	0.17	1.21	0.18	2.29	0.14	2.85	0.21	1.57	0.26	2.79	0.205	1.54	0.25	0.641	0.18
25	1.45	0.08	1.36	0.09	0.64	0.10	1.43	0.13	1.37	0.13	0.53	0.11	1.35	0.123	0.51	0.11	0.715	0.16
σ θ ( s ) 2	1.16	0.06	1.51	0.31			0.57	0.04	1.08	0.27			1.09	0.27
σ θ ( c ) 2							0.62	0.12	0.58	0.12			0.95	0.36
#pars	26		50				27		51				75
LL	−16,149.84		−16,008.65				−15,909.32		−15,750.19				−15,648.45					j
VPC(s)	0.26		0.31				0.13		0.22				0.20
VPC(c)							0.14		0.12				0.18

Note. IRT = item response theory; 2PL = two-parameter logistic; MLM = multilevel modeling; MLM-2PLC = multilevel 2PL constrained; MLM-2PL = multilevel 2PL; LL = log likelihood; VPC = variance partitioning coefficient.

Table 4 shows that across the models there are slight differences in the estimates of the item intercept. Also, the five models provide similar rankings of the items in terms of item intercept as the Spearman correlation coefficient estimates (r_s) between pairs of rankings of the estimates of intercept parameters obtained with different models ranged between .99 (2PL vs. MLM-2PL) to 1.00 (Rasch vs. ML Rasch). Recall, the relationship between item intercepts and location parameters can be used to interpret the item intercept (α) results in Table 4 in terms of “item difficulty parameters” rather than “item intercept” by multiplying each item intercept by −1.

The specification or not of the multilevel structure does not seem to have a relevant impact on the estimates of the item intercept parameters. However, the drawback of omitting the multilevel structure in the IRT specification is greater with respect to the estimates of the item loading or discrimination parameters as the Spearman correlation estimates ranged between .75 and .85 between the 2PL and MLM-2PL models. The specification of between-classes factor loadings in the MLM-2PL (class level or Level 3 in GLLAMM) does not have a relevant impact on the within individuals loadings ranking (student level or Level 2 in GLAMM) for models that account for the multilevel data structure. For instance, the Spearman correlation coefficient estimate for item loadings estimates in MLM-2PLC and MLM-2PL was .97.

Now, focusing on the variances of the random terms at class- and student level, the main findings can be summarized as follows:

To find the best fitting IRT model, a series of LRTs were conducted and the results showed that the MLM-2PL model had the best fit to the data. For instance, the LRT for the comparison of the MLM-2PLC model to the MLM-2PL model was −2(−15,750.19 − −15,648.45) = 203.48 with df = 75 – 51 = 24 which has an associated p value < .001. Therefore, the LRTs show that including the class-level loadings improves the explanation of the variability in item responses. In conclusion, an MLM-2PL model is recommended to appropriately represent the data.

MLIRT IRF

Within the MLIRT framework, the probability of responding correctly to an item is a function of both the class parameter θ_k and the person parameter θ_jk. The class parameter quantifies the effect of belonging to class k on student math achievement, whereas the latter represents the effect of the student component on the observed math achievement score. To jointly represent the item response probabilities as a function of θ_k and θ_jk, a tridimensional plot can be used. However, such a plot can be difficult to interpret and is not easily created, especially when the number of nested levels increases and dimensionality increases. Furthermore, there currently is not a standard approach for graphically representing multilevel IRT functions (Vector Psychometric Group, 2015). Therefore, we believe it is helpful to grasp differences in function behavior in two dimensions. To do so we create IRFs at student- and class level conditioned at critical points at the class level or student level of the latent trait (math achievement) continuum; we call these conditional item response functions (CIRFs). To create the CIRFs for any item, it is convenient to use Equation 6 and the MLIRT estimates where the model is identified by setting the variances of the latent trait to 1. The estimates of the model with variance of the random term fix equal to 1 can be approximated by multiplying each loading estimate in the MLM-2PL model in Table 5 by the standard deviation of the respective random term (Rabe-Hesketh, Skrondal, & Pickles, 2004a, 2004b). A CIRF that is built up conditioning on the student level (fixing θ_jk at a critical point and allowing θ_k free to vary) of the latent trait represents the probability to provide a correct response would vary across classes if they were composed of students with equal achievement in math for a given critical point. On the contrary, a CIRF that is built up conditioned on the class level (fixing θ_k at a critical point and allowing θ_jk free to vary) of the latent trait describes the probability of responding correctly to an item if students would attend classes which provide them the same contribution to their overall math achievement (i.e., classes with equal value of θ_k). The CIRFs for item 18 at student level with θ_k set equal to −1, 0, and 1 are plotted in Figure 2, whereas the CIRFs at class level with θ_jk set equal to −2, 0, and 2 are plotted in Figure 3. An inspection of Figures 2 and 3 clearly shows how Item 18 has greater discrimination power at the class level ( λ 18 ( c ) = 2.09; see Figure 3) versus the student level ( λ 18 ( s ) = 0.95; see Figure 2) given the steeper CIRFs. In order to show the implications of allowing the loadings to vary at student level and class level and to highlight practical differences across MLIRT models with different loading constraints, we plot in Figure 4 the CIRFs at student level (fixing θ ( c ) = 0 ; plots a, c, and e in Figure 4) and class level (fixing θ ( s ) = 0 ; plots b, d, and f in Figure 4) of the MLM-Rasch, MLM-2PLC, and MLM-2PL model for items m2, m10, m11, m18. Comparing the MLM-2PL (plots e and f) versus the MLM-Rasch (plots a and b) shows that the latter model inflates the steepness of item m2 loading and deflates the steepness of item m18 loading at the class level. Thus, the MLM-Rasch model provides misleading information about the contribution of the items to the class-level latent trait values. A similar pattern is observed with MLM-2PLC model (plots c and d) versus the MLM-2PL model (plots e and f).

OPEN IN VIEWER

Table 5.

Summary of MLM-2PL With Covariates.

Item	α i	SE	λ i ( s )	SE	λ i ( c )	SE
1	2.20	0.19	1		1.00
2	0.61	0.08	0.43	0.10	0.37	0.09
3	−0.10	0.10	0.00	0.07	0.59	0.12
4	1.89	0.15	0.83	0.15	0.72	0.17
5	0.35	0.13	0.63	0.12	0.82	0.16
6	1.35	0.14	0.86	0.15	0.78	0.17
7	−0.46	0.15	0.79	0.14	1.01	0.19
8	−0.02	0.14	1.10	0.18	0.85	0.16
9	−1.01	0.13	1.04	0.17	0.73	0.14
10	1.91	0.15	0.54	0.12	0.76	0.18
11	0.35	0.11	0.69	0.12	0.66	0.13
12	−0.11	0.12	0.96	0.15	0.72	0.14
13	0.11	0.13	1.02	0.17	0.74	0.15
14	−0.29	0.11	0.13	0.08	0.68	0.14
15	−0.08	0.15	1.15	0.18	0.94	0.18
16	−0.20	0.10	0.76	0.13	0.56	0.12
17	−0.67	0.16	0.77	0.13	1.02	0.19
18	0.02	0.28	0.97	0.17	2.03	0.38
19	0.98	0.11	1.12	0.18	0.46	0.11
20	0.27	0.13	0.88	0.14	0.76	0.15
21	0.96	0.15	0.78	0.14	0.96	0.20
22	0.24	0.10	1.07	0.17	0.42	0.10
23	1.97	0.18	1.75	0.28	0.61	0.15
24	2.86	0.21	1.58	0.26	0.58	0.17
25	1.39	0.12	0.57	0.12	0.67	0.15
Gender	−0.13*	0.07
Avg_pedc	0.10	0.02
σ θ ( s ) 2	0.88	0.23
σ θ ( c ) 2	1.04	0.39
#pars	77
LL	−15,605.54
VPC(s)	0.17
VPC(c)	0.20

Note. Gender = student gender (boys = 0, girls = 1); Avg_pedc = average parents’ education centered. MLM = multilevel modeling; 2PL = two-parameter logistic; LL = log likelihood; VPC = variance partitioning coefficient.

*

p < .05.

OPEN IN VIEWER

Figure 2.

Student-level item response function conditioned at given class levels for Item 18 on the 25-item math achievement test fit by the MLM-2PL model.

OPEN IN VIEWER

Figure 3.

Class-level item response function conditioned at given student levels for Item 18 on the 25-item math achievement test fit by the MLM-2PL model.

OPEN IN VIEWER

Figure 4.

IRFs at student level (setting theta(c) = 0) and at class level (setting theta(s) = 0) for MLM-Rasch (a, b), MLM-2PLC (c, d), and MLM-2PL (e, f).

MLIRT information and reliability

Within the framework of classical test theory (CTT), reliability is the consistency of a group of persons’ scores within one administration of a given multi-item instrument. A constant value is given to represent the degree of reliability associated with any raw score on a given instrument. Within this framework, a constant standard error of measurement value can be calculated to represent the precision or degree of error associated with any raw score on a given instrument.

Similarly, a single reliability index can be found within an IRT framework. A single IRT reliability index reflects, in general, the average or marginal reliability across the latent trait continuum (De Ayala, 2009, p. 177). Moreover, this IRT reliability index does not necessary reflect the amount of reliability in scores at the group average. For instance, an instrument may provide maximum reliability at θ = −1 and less elsewhere, but the group average latent trait score could be above the location of maximal reliability. Thus, the IRT reliability index could over- or under-estimate the reliability reflected in the group average trait score. For this reason, it is not recommended to report a single IRT reliability index unless the total information function (TIF) is uniform (De Ayala, 2009). Information in IRT terms is similar to the concept of reliability, but information is allowed to vary based on the location along the latent trait continuum. As such, it is important and useful for researchers conducting MLIRT analyses to know the degree of conditional information at each level of analysis (e.g., class level, student level) and when appropriate the amount of information associated with various discrete and continuous covariates used within an MLIRT framework to explain or predict differences on the latent trait. We can use this information to identify the amount of reliable information at various points along the latent trait continuum for each discrete and continuous covariate.

TIF is a function of the cumulative information available in each item and as such each item has its own item information function (IIF). Drawing from classical IRT wherein item information at a given location is

IIF i ( θ ) = ( P i ′ ) 2 P i ( 1 − P i ) = λ i 2 P i ( 1 − P i ) ,

where P i (as defined in Equation 3) is the probability of a correct response for item i and P i ′ is its first derivative with respect to the student parameter θ .

We can generate the IIF at the student level by deriving P_i (as defined in Equation 6) with respect to the person parameter and is expressed as

IIF i ( s ) ( θ ( s ) , θ ( c ) ) = ( P i ( θ ( s ) ) ′ ) 2 P i ( 1 − P i ) = λ i ( s ) 2 P i ( 1 − P i ) .

Similarly, an IIF can be derived at the class level by deriving P_i with respect to the class parameter and is expressed as

IIF i ( c ) ( θ ( s ) , θ ( c ) ) = ( P i ( θ ( c ) ) ′ ) 2 P i ( 1 − P i ) = λ i ( c ) 2 P i ( 1 − P i ) .

The IIFs for Item 18 at student level with θ _k set equal to −1, 0, and 1 are plotted in Figure 5, whereas the IIFs at class level with θ _jk set equal to −2, 0, and 2 are plotted in Figure 6. An inspection of Figures 5 and 6 clearly shows how Item 18 has more information at the class level (see Figure 6) versus the student level (see Figure 5) given the higher peaks of the IIFs. An inspection of IIFs in Figure 5 shows that when conditioned at a class level of math achievement, maximal student-level information occurs at different points along the student-level latent continuum. For instance, given θ^(c) = 1, the maximal student-level information occurs around θ^(s) −2.5, whereas when θ^(c) = −1, maximal student-level information occurs around 2.5. The IIFs at student level and class level provide information about the contribution of each item for assessing students and class achievement. The inspection of the IIFs at both levels helps to identify the amount of information each item provides to detect differences in mathematics achievement across students and classes. In practice, IIFs at both levels can be used to identify items that are helping to build up a test information function and which are not contributing much and need to be tossed or refined. Plus, the IIFs can be used to identify exemplary items that behave well at both levels for writing future items.

OPEN IN VIEWER

Figure 5.

Student-level item information function conditioned at given class levels for Item 18 on the 25-item math achievement test fit by the MLM-2PL model.

OPEN IN VIEWER

Figure 6.

Class-level item information function conditioned at given student levels for Item 18 on the 25-item math achievement test fit by the MLM-2PL model.

Broadly speaking, when class-level discrimination is zero, the IIFs conditioned at class level will be located at the same point, but when the ratio between class-level to student-level loading increases, the location of the IIFs will depart more and look similar to IIFs plotted in Figure 4. Interestingly, when student-level discrimination is zero, the IIFs conditioned at student level will be located at the same point, but when the ratio between student level to class level increases, the location of the IIFs will depart more.

The related TIFs are

TIF .. ( s ) ( θ ( s ) , θ ( c ) ) = ∑ i IIF i ( s ) ( θ ( s ) θ ( c ) )

and

TIF .. ( c ) ( θ ( s ) , θ ( c ) ) = ∑ i IIF i ( c ) ( θ ( s ) θ ( c ) ) .

Figure 7 shows the student-level TIFs conditioned at three points at the class level (−1, 0, and 1). Similar to the trends observed for the student-level IIFs, the shape of the TIFs show that when class-level achievement is high (1), the test is more informative for lower level students whereas when the class-level achievement is low (−1), the test is more informative for higher level students. Figure 8 shows the class-level TIFs conditioned at three points at the student level (−2, 0, and 2). Similar to the trends observed for the class-level IIFs, the shape of the TIFs show that when student-level achievement is high (2), the test tends to provide more information at lower points at the class level, whereas when student-level achievement is low (−2), the test tends to provide more information at higher points at the class level. The TIFs at student level and class level provide information about how well the test precisely assesses students and class achievement. The inspection of the TIFs at both levels helps to identify the amount of information a test provides to detect differences in mathematics achievement across students and classes. In practice, TIFs at both levels can be used to identify how well a test is measuring points along the latent trait continuum and if more items are needed to better measure points on this continuum at any given level of the data structure.

Moreover, if it was important to report a single value to capture the precision of the entire range of the latent trait, where the TIF is uniform, a marginal reliability coefficient can be estimated by extending Green, Bock, Humphreys, Linn, and Rechase (1984) concept of marginal reliability to MLIRT. Marginal reliability at the student level ≅ 1 – 1/Information(θ^[s]|θ^[c]) and at the class level ≅ 1 – 1/Information(θ^[c]|θ^[s]). For instance, if we fix θ^(c) = 0, then marginal reliability at the student level is 1 – 1 / 4.3 = .77 at θ^(s) = −0.8. Similarly, if we fix θ^(s) = 2, then marginal reliability at the class level is 1 – 1 / 3.4 = .71 at θ^(c) = −1.1.

OPEN IN VIEWER

Figure 7.

Student-level total information function conditioned at given class levels on the 25-item math achievement test fit by the MLM-2PL model.

OPEN IN VIEWER

Figure 8.

Class-level total information function conditioned at given student levels on the 25-item math achievement test fit by the MLM-2PL model.

Explanatory MLIRT

Once the MLIRT model has been built-up for descriptive purposes, researchers can attempt to explain the variability in the value of the latent trait (math achievement) at both student level and class level by including relevant covariates. The latent traits can be related to student and class characteristics expressed as

P ( m i j k = 1 | θ , α , λ , β s , β c ) = e α i x i + λ i ( s ) ( ∑ s = 1 S β s z j k s + θ ~ j k ) + λ i ( c ) ( ∑ c = 1 C β c w . k c + θ ~ k ) 1 + e α i x i + λ i ( s ) ( ∑ s = 1 S β s z j k s + θ ~ j k ) + λ i ( c ) ( ∑ c = 1 C β c w . k c + θ ~ k ) ,

where θ ~ j k is the adjusted student-level value of the latent trait after having removed the effect of student-level covariates z (S = number of student-level covariates), whereas θ ~ k is the adjusted class-level value of the latent trait after having removed the effect of the class-level covariates w (C = number of class-level covariates). The model expressed in Equation 14 is an extension of a latent regression model (Zwinderman, 1991). The model in Equation 14 includes two random terms. The variability of the random term θ ~ j k ( N [ 0 , σ ( θ ~ ( s ) ) 2 ] ) captures the between-students within-classes (student level) variance after having accounted for any student-level covariates, whereas the random term θ ~ k ( N [ 0 , σ ( θ ~ ( c ) ) 2 ] ) captures the between-class variance after adjusting for any class-level covariates.

Explanatory MLIRT results

Using Equation 14, we can add covariates of interest and test the following research question: What is the value added by fitting an explanatory MLIRT model with covariates at the student level to the analysis of students’ achievement in math? Specifically, in statistical terms, what is the nature and significance of the slope between student-level covariates (gender and parents education) as they predict scores on the latent variable (math)?

We answered this question by specifying an MLM-2PL model with covariates (MLM-2PL WC) which models the variability in the latent trait at student level as a function of gender (gender) and average level (years) of education of the students’ parents centered (avg_pedc) at the grand mean (M = 10.74, SD = 3.1). Observed values for avg_ped ranged from 5 to 17. A comparison of the MLM-2PL WC model to the model without covariates (MLM-2PL) shows it had better fit to the sample data. Specifically, the LRT was −2(−15,648.45 – −15,605.54) = 85.82 with df = 77 – 75 = 2, which has an associated p value < .001. Therefore, the LRT shows that including the student-level covariates improves the explanation of the variability in item responses.

Table 5 summarizes the results for the MLM-2PL WC model. An inspection of the student-level covariates shows gender and avg_pedc have significant slopes with the latent trait at the student level. Specifically, female students have on average lower achievement in math than male students, whereas more educated parents have on average a positive slope with the latent trait. The coefficient parameters highlight that the 577 female students have an average achievement in math that is 0.13 lower than the 564 males conditional upon the parents’ level of education. The 30 male students with avg_ped = 17 have an overall achievement that is 0.40 higher than the 101 male students who have avg_ped = 13 is about 0.63 higher than male students whose parents have a level of education equal to the average (avg_ped: M = 10.73). If we compare two students with the highest and lowest profile of both explanatory predictors in terms of change to score correctly (i.e., specifically the 30 “Male students with avg_ped = 17” and the 11 “Female students with avg_ped = 5”), the former have an achievement in math test that is on average about 1.33 higher than the latter. Translating this difference in terms of standard deviations from the average math achievement on the latent variable (that is equal to zero), this means that the two students’ profiles record on average an expected difference in the achievement equal to 1.55 z scores. This is akin to an effect size of 1.55 standard deviation differences between these groups. This can happen because the latent trait metric is standardized by default to mean of 0 and unit variance. It is worth highlighting that the adjustment for the two students’ covariates shift up the lower bound of the range of variation of the latent trait (math) and this clearly arises by comparing the range of variation of prediction of the latent trait at the student level [−2.60 to 2.49] with the range of variation of the adjusted student level [−2.39 to 2.48]. These results suggest that the socio-cultural characteristics of the family, summarized by the variable avg_ped, partially explain differences in students’ achievement in math achievement. This is highlighted also by comparing the variability of the random terms of models MLM-2PL versus MLM-2PL WC, from which it can be stated that the adjustment for covariates reduces the unexplained variability of the random terms at the student level and there is an increase of the variability at the class level. Thus, the omission of relevant students’ covariates leads to an overestimation of individual students’ differences. This means that the adjustment for covariates allows us to better understand how much of the divergences in math achievement are due to the two student-level characteristics, as it clearly arises by a comparison between the values of the variances of the two random terms, namely σ θ ~ w 2 = 0.88 and σ θ ~ B 2 = 1.03. The two student-level covariates considered in the analysis account for approximately 3% of the variability in math achievement (see Table 4, MLM-2PL VPC(s) = 0.20 and Table 5, MLM-2PL WC VPC(s) = 0.17). The MLM-2PLWC model highlights how accounting for students’ covariates increases the shared amount of unexplained variability due to differences across classes.

Table 5 provides the item parameters from the MLM-2PLWC model. If we rank the items according to the value of the item intercept, we observe that items “m₂₄, m₁, and m₁₀” are the ones for which the probability to answer correctly is highest, whereas items “m₉, m₁₇, and m₇” are the ones for which the probability to answer correctly is lowest. It is worth highlighting that the adjustment for student-level covariates did not substantively change the variation of the item intercept parameters in the MLM-2PL model (see Table 4) versus the MLM-2PL WC model (see Table 5).

The within-class between-student loadings (Level 2) and the between-classes factor loadings (Level 3) ranged from 0 to 1.75 and 0.37 to 2.03, respectively (see Table 5). The correlation between the loadings at class level and the loadings at student level was low ( ρ = 0.06), which indicates that the same item may have high discrimination power at student level and lower at class level or vice-versa. In essence, item loadings are not always the same at all levels of analysis. For instance, item m₂₃, which is the item with the greatest power to discriminate between students ( λ 23 ( s ) = 1 . 75 ) , but has a lower power to discriminate between classes ( λ 23 ( c ) = 0 . 72 ) .

Predictions of the latent traits

Once item parameters and variances of the random terms (if requested) are estimated, with or without covariates, random person (student) and cluster (class) parameters are estimated (predicted) in a second step (see Muthén & Muthén, 1998-2015; Skrondal & Rabe-Hesketh, 2004, 2009; StataCorp, 2013). In general, the information or details about any given latent student and class parameters is found by taking the resulting predictions for θ ~ j k and θ ~ k with the respective measures of uncertainty (which can be interpreted as a standard error). For instance, student-level θ ~ j k is an adjusted indicator of student math achievement estimated by subtracting the effect of factors z on the parameter θ _jk . The literature on factor models refers to these measures as “factor scores,” whereas the literature on multilevel models calls them “empirical Bayes predictions,” “posterior prediction,” or “higher level residuals” (Rabe-Hesketh et al., 2004b; Skrondal & Rabe-Hesketh, 2009). Finally, comparisons of each class to the average class-level achievement can be done using the predicted random term at the class-level θ _k with the respective measure of uncertainty S D k θ . Figure 9 plots classes according to their latent level of math achievement based on 95% CIs. If we want to make a paired-comparison between two classes using a 95% CI, then we need to construct the caterpillar plot by taking the standard 1.96 /√2 (Goldstein & Healy, 1995; Goldstein & Spiegelhalter, 1996).

OPEN IN VIEWER

Figure 9.

Expected posterior predictions at class level with 95% CI.

Estimation methods

The likelihood function of IRT models which assumes a continuous distribution of the random effects (as the likelihood of the models just described) cannot be maximized using analytical methods because it does not have a closed form solution. The integrals involved in the functions need to be first approximated by sums before to start the maximization process. The Stata procedures GLLAMM and SEM and the program Mplus have the option of using a maximum likelihood estimation method that marginalizes the likelihood function by integrating out the random terms using a numerical quadrature method. The numerical integrated marginal likelihood is then maximized using the Newthon-Raphson algorithm. In all models fitted in this article with Stata, we adopted marginal maximum likelihood (MML) with adaptive quadrature, which ensures good estimates of the parameters compared with other quadrature methods (e.g., Gauss Hermite is less accurate also in cases where standard quadrature work poorly, see for more details Rabe-Hesketh, Skrondal, & Pickles, 2004a). Furthermore, adaptive quadrature works faster and needs fewer quadrature points [nodes]. In flexMIRT, we used both Bock-Aitkin EM (BAEM) algorithm (Bock & Aitkin, 1981) and the Metropolis Hastings-Robbins Monro (MH-RM; see Cai, 2010) algorithm for estimating all descriptive IRT and MLIRT models. In flexMIRT, MLIRT models can be estimated using also the default BAEM; however, it can take noticeable time to converge, whereas the MLIRT WC models can be estimated only with the MH-EM algorithm. In our experiences, the BAEM in flexMIRT takes considerably longer with MLIRT models estimating more than one latent trait dimension, but are efficient with MLIRT without covariates as considered herein.

The estimates of the MH-EM algorithm in flexMIRT for Rasch and 2PL models almost match (Cai, 2009, 2010) with the results which rely on the other marginal maximum likelihood algorithms available in flexMIRT, Stata and Mplus, whereas the estimates of the MLIRT models with the MH-EM show some slight differences with respect to the results observed with the other estimation methods. In our experience the biggest divergences between the MLM IRT model estimates obtained with marginal maximum likelihood algorithms in Stata, Mplus and flexMIRT and the MH-EM results are observed in the estimation of the variance of the class-level random term. Specifically, with respect to the algorithms which rely on the likelihood, the MH-EM algorithm increases the estimates of the class-level random term in the MLM-Rasch and downward the class-level random term in the MLM-2PL and MLM-2PL WC.

In general, Stata, Mplus, and flexMIRT provide similar results for the variety of models considered in this article, but the main drawbacks of Stata GLLAMM are related to the computational time required to approximate the likelihood which increases as a function of the number of latent variables, the number of levels, the number of quadrature points, the number of the observations in the data set, and the square of the number of parameters involved in the estimation. It is interesting to highlight that for models that can be specified with both Stata routines, it is more convenient to use the Stata program SEM that is more efficient in terms of computational time. Specifically, Stata program SEM is faster than GLAMM, but for more complex models (e.g., multidimensional MLIRT with more than three levels) GLAMM can take considerably longer to converge (Rabe-Hesketh et al., 2004b; StataCorp, 2013). From our experiences, this efficiency is also true when using Mplus and flexMIRT. For the models considered herein, the most efficient (fastest) software, in rank order, were flexMIRT, Mplus, Stata program SEM, and Stata program GLLAMM.

We presented the results using the Stata program GLAMM because it is the only routine that estimates MLIRT using a data set in long format within a generalized multilevel or three-level multilevel model, whereas the Stata program SEM and Mplus estimates MLIRT models using a data set in wide format within a two-level multilevel multivariate model. Appendix B provides a summary comparison of some of the software capabilities and features with respect to the MLIRT models for binary responses considered in this article. In general, Stata programs GLAMM and SEM, Mplus, and flexMIRT have similar capabilities and limitations, except Stata programs GLAMM and SEM allow users to model more than three levels of data.