Multilevel Higher-Order Item Response Theory Models

Abstract

In the social sciences, latent traits often have a hierarchical structure, and data can be sampled from multiple levels. Both hierarchical latent traits and multilevel data can occur simultaneously. In this study, we developed a general class of item response theory models to accommodate both hierarchical latent traits and multilevel data. The freeware WinBUGS was used for parameter estimation. A series of simulations were conducted to evaluate the parameter recovery and the consequence of ignoring the multilevel structure. The results indicated that the parameters were recovered fairly well; ignoring multilevel structures led to poor parameter estimation, overestimation of test reliability for the second-order latent trait, and underestimation of test reliability for the first-order latent traits. The Bayesian deviance information criterion and posterior predictive model checking were helpful for model comparison and model-data fit assessment. Two empirical examples that involve an ability test and a teaching effectiveness assessment are provided.

Keywords

item response theory higher-order latent trait multilevel model Markov chain Monte Carlo (MCMC) estimation

In some large-scale educational and psychological testing programs, the test takers have a multilevel structure and the measured latent traits have a hierarchical structure. For example, in the Programme for International Student Assessment (or PISA), approximately 150 schools are first randomly selected from a country, and approximately 30 to 50 students are then randomly sampled from each sampled school. Such a two-stage sampling creates a so-called multilevel data structure. In addition, the latent traits measured in the Programme for International Student Assessment can be regarded as having a hierarchical structure in which the domains of a subject are the first-order latent traits (specific domain abilities), and the subject is the second-order latent trait (the overall ability). Consider the subject of mathematics as an example. We can treat the four domains of Quality, Space and Shape, Change and Relationship, and Uncertainty as the four first-order latent traits, which are governed by the second-order latent trait of “mathematical problem solving.” Multilevel data structure and hierarchical latent traits also occur in other large-scale testing programs, such as the Trends in International Mathematics and Science Study (TIMSS), the Progress in International Reading Literacy Study, the International Civic and Citizenship Education Study, and the National Assessment of Educational Progress.

It is very likely that subjects (e.g., students) who are sampled from the same cluster (e.g., school) are more homogenous than subjects who are sampled from different clusters because those from the same cluster share the same environment and experience (Hill & Rowe, 1996; Sellstrom & Bremberg, 2006). Failing to consider multilevel structures results in biased parameter estimates and inadequate standard errors (Raudenbush & Bryk, 2002). To resolve this problem, researchers have developed multilevel models to account for multilevel structures in categorical and continuous data (Fox, 2010; Raudenbush & Bryk, 2002; Wang & Qiu, 2013), but these models fail to consider hierarchical structures in latent traits. On the other hand, researchers have developed higher-order item response theory (IRT) models to account for hierarchical structures in latent traits (de la Torre & Hong, 2010; de la Torre & Song, 2009; H.-Y. Huang & Wang, 2013; H.-Y. Huang, Wang, Chen, & Su, 2013; Sheng & Wikle, 2008), but these models fail to consider multilevel structures. Given that data have a multilevel structure and latent traits have a hierarchical structure and that multilevel IRT models for single-order latent traits and higher-order IRT models for single-level data are available, it is of great value to integrate both types of models to create multilevel higher-order IRT (MHIRT) models in such a way that the data can be analyzed appropriately, which is the major purpose of this study.

We introduce the new class of MHIRT models, Markov chain Monte Carlo (MCMC) methods for parameter estimation, and Bayesian model-data fit checking techniques. We then outline the results of a series of simulations that were conducted to assess the parameter recovery and efficiency of model-data fit indices by using the freeware WinBUGS (Spiegelhalter, Thomas, & Best, 2003). We give two empirical examples of achievement testing and psychological inventory to illustrate the new models. Finally, we draw conclusions about the development of the new models and provide suggestions for future studies.

Multilevel Higher-Order Item Response Theory Models

There are three components in MHIRT models: an item response function, a higher-order structure of latent traits, and a multilevel structure of parameters. First, consider the item response functions. Let θ_n denote a vector of latent traits for person n; let f(X_ni| θ_n) denote the item response function, which is conditional on θ_n; and let X_ni be the item response to item i for person n. In a general approach, the item response function f(X_ni| θ_n) in the MHIRT models should not be limited to a specific type; instead, they can be standard IRT models, such as the one-, two-, and three-parameter logistic models for dichotomous items (1PLM, 2PLM, 3PLM; Birnbaum, 1968; Rasch 1960), and the rating scale models (RSMs; Andrich, 1978), the partial credit model (PCM; Masters, 1982), the generalized partial credit model (GPCM; Muraki, 1992), and the graded response model (GRM; Samejima, 1969) for polytomous items; or other less common IRT models, such as the linear logistic test model (Fischer, 1983), the facets model (Linacre, 1989), testlet (bifactor) response models (Wainer, Bradlow, & Wang, 2007), speeded IRT models (Goegebeur, De Boeck, Wollack, & Cohen, 2008), and unfolding IRT models (Andrich, 1996). More importantly, users should be allowed to create their own customized IRT models.

Second, consider higher-order structures of latent traits. Let $θ_{nv}^{(p)}$ be the pth-order (p = 1, . . ., P) vth latent trait (v = 1, . . ., V) for person n (n = 1, . . ., N). The relationship between a lower-order and its higher-order latent traits is assumed to be

θ_{n v}^{(p)} = β_{0 v}^{(p)} + β' (p) θ_{n}^{(p + 1)} + ε_{n v}^{(p)}

where $β_{0 v}^{(p)}$ is the intercept for the pth-order vth latent trait, $β^{' (p)}$ is a vector of regression parameters, $θ_{n}^{(p + 1)}$ is a vector of (p+ 1)th-order latent traits that govern $θ_{nv}^{(p)}$ , and $ε_{nv}^{(p)}$ is the residuals and is assumed to be normally distributed. As in multiple regression, the equation can be polynomial. In addition, a lower-order latent trait can be governed by more than one higher-order latent trait, and higher-order latent traits can have an interaction effect on lower-order latent traits (H.-Y. Huang et al., 2013). For illustrative simplicity, let there be two orders of latent trait, including $θ_{nv}^{(1)}$ and $θ_{n}^{(2)}$ , and let their relationship be

θ_{n v}^{(1)} = β_{0 v} + β_{1 v} θ_{n}^{(2)} + ε_{n v}^{(1)},

where $β_{0 v}$ , $β_{1 v}$ , and $ε_{nv}^{(1)}$ are the intercept, slope, and residual for the vth first-order latent trait, respectively; $ε_{nv}^{(1)}$ is assumed to be normally distributed with a mean of zero and a variance of $σ_{ε}^{2}$ and independent of $θ_{n}^{(2)}$ .

Third, consider multilevel structures of parameters. For simplicity, let there be two levels (e.g., schools are sampled from a population and students are sampled from schools). Let $θ_{ngv}^{(1)}$ and $θ_{ng}^{(2)}$ denote the vth first-order latent trait and the second-order latent trait, for student n sampled from school g, respectively. At Level 1 (student level), the relationship between these two orders of latent traits is assumed to be

θ_{n g v}^{(1)} = β_{0 g v} + β_{1 v g} θ_{n g}^{(2)} + ε_{n g v}^{(1)},

where $β_{0 gv}$ is the intercept and $β_{1 vg}$ is the slope (factor loading) for school g; $ε_{ngv}^{(1)}$ is the within-school residual and is typically assumed to be homogeneous across schools and normally distributed with a mean of zero and independent of the other ε and θ variables. When appropriate, Level-1 predictors (e.g., a student’s gender and ethnicity) can be added to Equation (3):

θ_{n g v}^{(1)} = β_{0 g v} + β_{1 v g} θ_{n g}^{(2)} + \sum_{k = 2}^{K} β_{k g v} W_{k n g} + ε_{n g v}^{(1)},

where $W_{kng}$ is the kth Level-1 predictor for person n in school g and $β_{kgv}$ is its corresponding regression weight for the vth latent trait and school g.

At Level 2 (the school level), the parameters in Equation (4) can be modeled as

β_{0 g v} = γ_{00 v} + u_{0 g v},

β_{1 g v} = γ_{10 v} + u_{1 g v},

β_{k g v} = γ_{k 0 v} + u_{k g v}, k = 2, \dots, K,

where the residuals $u_{g}^{T} \equiv (u_{0 gv}, u_{1 gv}, u_{2 gv}, \dots, u_{Kgv})$ are assumed to follow a multivariate normal distribution with a mean vector of zero and a covariance matrix of $∑_{u}$ . In Equations (3) or (4), different schools have different slopes, which suggests that the metric is not invariant across schools. The variance of $u_{1 gv}$ across schools depicts the magnitude of the metric variance on the vth first-order latent trait: The larger the variance, the larger is the metric variance. When the variance is sufficiently small, one can constrain $β_{1 gv} = γ_{10 v}$ , to make the metric invariant. It is also possible to set constraints on Equation (5). For example, if the variance of $u_{0 gv}$ across schools is sufficiently small, then one can constrain: $β_{0 gv} = γ_{00 v}$ , which makes the scalar invariant. In the following simulation studies, we assume that the metric is invariant, whereas the scalar is not. Note that such settings do not result in any loss in the generality of the MHIRT models.

When appropriate, Level-2 predictors (e.g., school type, school size) can be added to Equation (5):

β_{0 g v} = γ_{00 v} + \sum_{k = 1}^{K} γ_{k v} Z_{k g} + u_{0 g v},

where $Z_{kg}$ is the kth Level-2 predictor and $γ_{kv}$ is its corresponding regression weight for the vth latent trait.

Note that $β_{0 gv}$ (g = 1, . . ., G; v = 1, . . ., V) is a random intercept, and there are a total of V random variables. In Equations (5) or (8), it is assumed that $u_{0 gv}$ follows a multivariate normal distribution with a mean vector of zeros that have length V and a full covariance matrix with the size of V by V (denoted as $∑_{v}$ ). No structural relationship is specified on the covariance matrix, which makes the model more general than when the structural relationship is specified. Within the framework of multilevel but single-order factor analysis, Muthén and Muthén (2012, pp. 274-276) propose a so-called cluster-level latent trait to account for the covariance matrix:

β_{0 g v} = μ_{v} + λ_{v} ξ_{g},

where $μ_{v}$ is the mean of the intercepts across clusters for the vth latent trait, $λ_{v}$ is the factor loading for the vth latent trait, $ξ_{g}$ is the factor score for cluster g, and the residuals are fixed at zero. Apparently, Equation (9) is more parsimonious than the full covariance matrix with any structural relationship.

Several procedures are required to identify MHIRT models. For each test that measures a first-order latent trait, either the mean of the item difficulties in the test should be set to zero, in such a way that all the random intercepts $β_{0 gv}$ can be freely estimated, or $γ_{00 v}$ (the grand mean of $β_{0 gv}$ across schools) should be set to zero, in such a way that all the item difficulties in the test can be freely estimated. If the IRT model has a slope parameter (e.g., the 2PLM, 3PL, GPCM, and GRM), in addition to setting $θ^{(2)}$ to have a unit variance, the slope of an item (e.g., Item 1) in each test should be set to unity, or the error variance of $ε_{ngv}^{(1)}$ for each test should be set to some value (one minus $β_{1 v}^{2}$ was set in this study).

Markov Chain Monte Carlo Estimation and Bayesian Model-Data Fit Assessment

MHIRT models involve both high dimensionality (many random effects) and a complex multilevel structure, in such a way that marginal maximum likelihood estimation procedures become inefficient. We thus recommend Bayesian estimation with MCMC methods as an alternative and used the freeware WinBUGS to calibrate the model parameters. WinBUGS not only readily furnishes an easy estimation procedure for users who are not familiar with the complicated algorithms in Bayesian estimation but also enables users to fit a wide range of IRT models and even customized models. The large amount of flexibility is attainable by changing a few lines of previously written code for an item response function (e.g., the 1PLM) when a new item response function is implemented (e.g., the 3PLM).

In Bayesian estimation, the specifications of a statistical model, prior distributions of model parameters, and observed data are required to produce a joint posterior distribution for the model parameters. MCMC methods provide alternative and simple ways to simulate the joint posterior distribution of the unknown quantities and obtain simulation-based estimates of posterior parameters of interest. The next section details the prior distributions for the model parameters according to previous Bayesian IRT analyses. Although only WinBUGS was used to assess the parameter recovery and produce a fit to the empirical data, the general conclusions would not change when other computer programs are used. For example, the JAGS program (Just Another Gibbs Sampler; Plummer, 2003) is found to be computationally equivalent to WinBUGS (Curtis, 2010). In addition, the parameter estimates that are obtained from WinBUGS are comparable to those obtained by using marginal maximum likelihood estimation procedures (Jiao, Wang, & He, 2013), although Bayesian methods are computationally more efficient in highly complicated and multidimensional models. The computer program Mplus (Muthén & Muthén, 2012), although very flexible for either multilevel models or higher-order models, is not feasible for multilevel higher-order models. Furthermore, the item response functions that are applicable in Mplus are cumulative logits (Agresti, 2010), which means that only the 2PLM for dichotomous items and the GRM for polytomous items are applicable; other types of IRT models (e.g., the 3PLM, PCM, GPCM, and customized models) are not.

Bayesian model-data fitting can be conducted by posterior predictive model checking (PPMC), which assesses the plausibility of posterior predictive replicated data against observed data and has the advantage of a strong theoretical basis and an intuitively appealing simplicity that can be applied as numerical evidence (Gelman, Meng, & Stern, 1996). A test statistic is chosen to detect the systematic discrepancy between the observed and replicated data given the model parameters. An extreme p value (close to 0 or 1) indicates a model-data misfit.

Three statistics are chosen in this study to assess the model-data fit based on PPMC to satisfy the nature of the complicated MHIRT models. The first statistic is the Bayesian chi-square test (Sinharay, Johnson, & Stern, 2006), which assesses the overall model-data fit. The second statistic is chosen to assess the factor structure and overall model-data fit by comparing the reproduced correlation matrix with the original correlation matrix. The reproduced correlation between the first-order latent traits can be computed through β_vβ_v′ ( $v \neq v'$ ). If the difference between the two correlation matrixes is large, then the model-data fit is poor. The third statistic is obtained by comparing the standard deviation (SD) of the observed group’s mean raw scores with that of the replicated group’s mean raw scores for each test. If the multilevel structure is ignored and analyzed with a single-level or lower level model, then a systematic discrepancy between the two standard deviations should suggest poor model-data fit.

To conduct a model comparison, the Bayesian deviance information criterion (DIC) can be used, which simultaneously accounts for the model fit and model complexity (Spiegelhalter, Best, Carlin, & van der Linde, 2002): $DIC = \bar{D} + PD$ , where $\bar{D}$ is the posterior expectation of the deviance and PD is defined as the difference between the posterior mean of the deviance and the deviance at the posterior mean of the parameters and is referred to as the effective number of parameters. A smaller DIC means a better model-data fit.

Simulation Studies

Simulation Design

Two simulation studies were performed. Simulation Study 1 focused on the parameter recovery and model-data fit indices for MHIRT models without covariates in three dichotomous IRT models and four polytomous IRT models and on the consequences of ignoring the multilevel structure during parameter estimation. Simulation Study 2 focused on the parameter recovery for MHIRT models with covariates when the 3PLM and the GPCM were used as the IRT model. In Simulation Study 1, there were three first-order latent traits and one second-order latent trait. Each test had either 20 dichotomous items or 20 four-point polytomous items. In each test, the 1PLM, 2PLM, and 3PLM were used to generate responses to 20 dichotomous items, and the GPCM, PCM, GRM, and RSM were used to generate responses to 20 polytomous items. The second-order latent traits were generated from a standard normal distribution. The factor loadings that represented the correlations between the second-order latent trait and the three first-order latent traits were set at 0.9, 0.8, and 0.7, and their corresponding variances of the within-school residuals were set at $1 - 0 . 9^{2},$ $1 - 0 . 8^{2}$ , and $1 - 0 . 7^{2}$ , respectively, for the three first-order latent traits. For the 60 dichotomous items, their difficulty parameters were generated from a uniform distribution U(−2, 2), their discrimination parameters were generated from N(1, 0.25) and truncated above zero, and their pseudo-guessing parameters were generated from a uniform distribution U(0, 0.3). For the 60 four-point polytomous items, the location parameters were generated from a uniform distribution U(−2.5, 2.5) for the GPCM and PCM, and the overall difficulty parameters were generated from a uniform distribution U(−1.5, 1.5) with negative 1, 0, and 1 as the category threshold parameters, respectively, for both the GRM and RSM. The procedure to generate the discrimination parameters for both the GPCM and GRM was the same as in the dichotomous-item condition. The specifications of item parameters were in line with those commonly found in practice.

A sample of 2,000 persons for the 1PLM, 2PLM, and all polytomous response models and 5,000 examinees for the 3PLM were generated; all the examinees were divided into 100 schools, with each school having 20 and 50 persons when the sample size was 2,000 and 5,000, respectively. A large sample size of 5,000 was used in the 3PLM because the 3PLM often requires a large sample size to yield stable parameter estimates. The variance–covariance matrix ( $∑_{v}$ ) at the school level between latent traits was set to have variances of one and covariances of 0.6. The data sets were generated according to the MHIRT models. The generated data sets were analyzed with two models: one model was the data-generating model (i.e., the MHIRT models), which assessed the parameter recovery, and the other model was the single-level HIRT, which assessed the consequences of ignoring the multilevel structure.

In Simulation Study 2, item responses were generated from the 3PLM and GPCM. The simulation design was identical to that in Simulation Study 1, except that there was a school-level predictor (i.e., MHIRT models with covariates). For the 50 schools that had higher latent traits, the value of the predictor was 0, whereas for the 50 schools that had lower latent traits, the value was 1. The school-level regression weight was set to negative 1; in other words, the overall mean intercept for the higher proficiency schools was 0, and for the lower proficiency schools, it was negative 1.

A Matlab computer program was written to generate responses for each condition. A total of 20 replications were conducted for each condition, mainly because each replication can take more than 3 days of computer time in Simulation Study 1 and more than 1 week in Simulation Study 2. Fortunately, in our experience, 20 replications appeared to be sufficient to gain reliable inferences in Bayesian IRT models because the empirical sampling variations across the replications were small. Many simulation studies with Bayesian estimation use only 10 replications (e.g., Bolt & Lall, 2003; Cao & Stokes, 2008).

Analysis

WinBUGS 1.4 was used to calibrate the model parameters. A normal prior with a mean of 0 and a variance of 4 was set for all the location, threshold, and regression weight parameters; a lognormal prior with a mean of 0 and a variance of 4 was set for the discrimination parameters; a beta prior with both hyperparameters equal to 1 was set for the pseudo-guessing parameters; and a normal prior with a mean of 0.5 and a variance of 10 was set for the factor loadings. A Wishart distribution with a diagonal matrix equal to 0.1 and degrees of freedom of 3 (i.e., the number of tests) was set as the prior for the inverse of the variance–covariance matrix.

A total of 15,000 iterations were conducted, with the first 5,000 iterations treated as burn-in after monitoring the convergence diagnostic according to the multivariate potential scale reduction factor (Brooks & Gelman, 1998), with three parallel chains for the first simulated data set across all analysis models. In addition, 200 MCMC samples obtained from the remaining 10,000 iterations with a thinning factor of 50 were used to compute the test reliability (the precision of the person measures).

For each estimator, we computed the bias and the root mean square error (RMSE). It was expected that (a) the parameter recovery of the MHIRT models would be satisfactory, (b) the multilevel structure by fitting single-level HIRT models, when ignored, would yield poor parameter estimates, (c) the Bayesian model-data fit checking techniques (i.e., PPMC and Bayesian DIC) would be able to identify the generating MHIRT models, and (d) the parameter recovery of the MHIRT models with covariates would be satisfactory.

Results

Simulation Study 1: Parameter Recovery and Model-Data Fit Checking for MHIRT Models

Because of space constraints, only the mean, standard deviation, and range of the bias and RMSE across parameters are reported. Table 1 summarizes the parameter recovery for dichotomous items generated from the 1P-, 2P-, and 3P-MHIRT, when the data-generating model and analysis model were identical. With respect to the bias, it was found that the bias values appeared to be close to 0 for all the models. In the 1P-MHIRT, the mean RMSE was 0.106 for the difficulty parameters, 0.021 for the common slope parameters, and 0.149 for the variance–covariance matrix; and the RMSE was 0.030, 0.021, and 0.027 for the three factor loadings. The parameter recovery appeared to be good. In the 2P-MHIRT, the parameter recovery was satisfactory, although slightly worse than that in the 1P-MHIRT. The mean RMSE was 0.116 for the difficulty parameters, 0.064 for the slope parameters, and 0.163 for the variance–covariance matrix; the RMSE was 0.033, 0.018, and 0.014 for the three factor loadings. The parameter recovery in the 3P-MHIRT was also satisfactory, although slightly worse than that in the 2P-MHIRT. The mean RMSE was 0.207 for the difficulty parameters, 0.083 for the slope parameters, 0.049 for the pseudo-guessing parameters, and 0.123 for the variance–covariance matrix; the RMSE was 0.019, 0.018, and 0.028 for the three factor loadings.

Table 1.

Parameter Recovery of Dichotomous Items in Simulation Study 1.

Model	1PLM				2PLM				3PLM
Multilevel	MHIRT		HIRT		MHIRT		HIRT		MHIRT		HIRT
Parameter	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
Variance–covariance
Mean	0.041	0.149	—	—	0.033	0.163	—	—	0.022	0.123	—	—
SD	0.039	0.029	—	—	0.031	0.046	—	—	0.025	0.016	—	—
Max	0.079	0.209	—	—	0.091	0.254	—	—	0.068	0.153	—	—
Min	−0.035	0.116	—	—	0.005	0.105	—	—	−0.017	0.110	—	—
Difficulty
Mean	0.007	0.106	0.002	0.110	0.004	0.116	0.040	0.342	0.056	0.207	0.076	0.412
SD	0.023	0.020	0.021	0.024	0.030	0.037	0.366	0.162	0.128	0.130	0.424	0.190
Max	0.079	0.155	0.075	0.160	0.088	0.210	0.573	0.584	0.747	0.855	0.851	0.885
Min	0.040	0.063	−0.039	0.059	−0.056	0.050	−0.574	0.081	−0.121	0.081	−0.595	0.091
Slope
Mean	0.000	0.021	0.414	0.418	−0.003	0.064	0.395	0.409	0.014	0.083	0.417	0.436
SD	0.010	0.001	0.015	0.015	0.019	0.018	0.125	0.127	0.030	0.031	0.133	0.137
Max	0.012	0.023	0.435	0.439	0.043	0.109	0.656	0.674	0.084	0.208	0.686	0.720
Min	−0.013	0.019	0.401	0.405	−0.048	0.040	0.136	0.147	−0.066	0.038	0.170	0.200
Pseudo-guessing
Mean	—	—	—	—	—	—	—	—	0.013	0.049	0.010	0.052
SD	—	—	—	—	—	—	—	—	0.031	0.029	0.035	0.029
Max	—	—	—	—	—	—	—	—	0.107	0.130	0.115	0.133
Min	—	—	—	—	—	—	—	—	−0.043	0.009	−0.048	0.010
Factor loading
β₁	0.009	0.030	−0.066	0.073	−0.004	0.033	−0.085	0.096	0.005	0.019	−0.057	0.070
β₂	0.005	0.021	0.018	0.032	0.005	0.018	−0.016	0.034	−0.007	0.018	−0.029	0.039
β₃	−0.002	0.027	0.020	0.036	0.001	0.014	0.050	0.054	0.005	0.028	0.033	0.045

Note. HIRT = higher-order item response theory; MHIRT = multilevel higher-order item response theory; RMSE = root mean square error; 1PLM = one-parameter logistic model; 2PLM = two-parameter logistic model; 3PLM = three-parameter logistic model.

When the multilevel structure was ignored by fitting single-level HIRT models, the resulting parameter recovery was much worse, as shown in Table 1. Take the 3P-HIRT as an example. The mean RMSE was 0.412 for the difficulty parameters, 0.436 for the slope parameters, and 0.052 for the pseudo-guessing parameters; the RMSE was 0.070, 0.039, and 0.045 for the three factor loadings. Regarding the bias, it appeared that the slope parameters were considerably overestimated, as is the case for the 2P- and 1P-HIRTs.

Table 2 summarizes the parameter recovery for polytomous items generated from the PC-, GPC-, RS-, and GR-MHIRT when the MHIRT models were fit to the MHIRT data. In the PC-MHIRT, the mean RMSE was 0.138 for the location parameters, 0.021 for the common slope parameters, and 0.099 for the variance–covariance matrix; the RMSE was 0.013, 0.017, and 0.011 for the three factor loadings. In the GPC-MHIRT, the mean RMSE was 0.140 for the location parameters, 0.046 for the slope parameters, and 0.169 for the variance–covariance matrix; the RMSE was 0.027, 0.019, and 0.017 for the three factor loadings. The parameter recovery for the GPC-HHIRT was slightly worse than that for the PCM-MHIRT. In the RS-MHIRT, the mean RMSE was 0.114 for the location parameters, 0.024 for the common slope parameters, and 0.120 for the variance–covariance matrix; the RMSE was 0.012, 0.020, and 0.018 for the three factor loadings. In the GR-MHIRT, the mean RMSE was 0.107 for the location parameters, 0.043 for the slope parameters, and 0.221 for the variance–covariance matrix; the RMSE was 0.018, 0.021, and 0.014 for the three factor loadings. Compared with the RS-MHIRT, slightly worse parameter estimates were obtained with the GR-MHIRT. In summary, all the models appeared to have a good parameter recovery, although the GPC-MHIRT was slightly inferior to the other three models.

Table 2.

Parameter Recovery of Polytomous Items in Simulation Study 1.

Model	PCM				GPCM				RSM				GRM
Multilevel	MHIRT		HIRT		MHIRT		HIRT		MHIRT		HIRT		MHIRT		HIRT
Parameter	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
Variance–covariance
Mean	−0.033	0.099	—	—	0.054	0.169	—	—	0.017	0.120	—	—	0.131	0.221	—	—
SD	0.016	0.024	—	—	0.017	0.035	—	—	0.029	0.026	—	—	0.022	0.016	—	—
Max	−0.010	0.132	—	—	0.078	0.229	—	—	0.047	0.166	—	—	0.177	0.251	—	—
Min	−0.056	0.070	—	—	0.025	0.111	—	—	−0.030	0.080	—	—	0.112	0.206	—	—
Location
Mean	−0.074	0.138	−0.075	0.147	−0.035	0.140	−0.037	0.365	−0.059	0.114	−0.065	0.116	0.019	0.107	0.042	0.309
SD	0.029	0.025	0.033	0.027	0.034	0.034	0.393	0.188	0.024	0.020	0.023	0.016	0.028	0.014	0.353	0.194
Max	0.016	0.224	0.001	0.245	0.094	0.293	0.710	0.718	−0.011	0.161	−0.017	0.155	0.078	0.135	0.771	0.775
Min	−0.147	0.086	−0.163	0.089	−0.118	0.071	−0.701	0.072	−0.116	0.077	−0.120	0.082	−0.055	0.072	−0.742	0.059
Slope
Mean	0.010	0.021	0.418	0.420	−0.012	0.046	0.377	0.386	0.000	0.024	0.407	0.409	−0.008	0.043	0.418	0.424
SD	0.002	0.003	0.016	0.015	0.013	0.016	0.116	0.120	0.002	0.005	0.007	0.006	0.014	0.014	0.130	0.131
Max	0.013	0.024	0.437	0.439	0.016	0.103	0.686	0.719	0.001	0.028	0.416	0.417	0.034	0.080	0.712	0.720
Min	0.008	0.017	0.399	0.402	−0.043	0.020	0.120	0.125	−0.004	0.017	0.399	0.401	−0.042	0.016	0.144	0.149
Loading
β₁	−0.003	0.013	−0.077	0.082	0.001	0.027	−0.058	0.066	0.003	0.012	−0.042	0.055	0.001	0.018	−0.059	0.067
β₂	0.000	0.017	−0.009	0.035	−0.010	0.019	−0.028	0.044	−0.008	0.020	−0.015	0.038	−0.011	0.021	−0.012	0.042
β₃	0.003	0.011	0.029	0.050	−0.011	0.017	0.023	0.041	−0.006	0.018	0.010	0.021	0.006	0.014	0.049	0.056

Note. HIRT = higher-order item response theory; MHIRT = multilevel higher-order item response theory; RSM = rating scale model; RMSE = root mean square error; PCM = partial credit model; GPCM = generalized partial credit model; GRM = graded response model.

When the multilevel structure was ignored by fitting single-level HIRT models, as shown in Table 2, the resulting parameter recovery was much worse. Take the GPC-HIRT as an example. The mean RMSE was 0.365 for the location parameters and 0.386 for the slope parameters; and the RMSE was 0.066, 0.044, and 0.041 for the three factor loadings. Regarding the bias, the slope parameters were considerably overestimated, as is the case for the PC-HIRT, RS-HIRT, and GR-HIRT.

Table 3 lists the mean RMSE of the person measures and the mean and standard deviation of the test reliability across 20 replications for dichotomous and polytomous items. With respect to the RMSE of the person measures in the first- and second-order latent traits, fitting the true MHIRT models always yielded a much smaller RMSE than that yielded from the corresponding single-level HIRT models; and simpler models (e.g., the family of Rasch models) yielded a smaller RMSE than that yielded from complicated models (e.g., the 3PLM and GPCM). With respect to the test reliability for the dichotomous items, the 1PL-MHIRT had the highest test reliability, followed by the 2PL-MHIRT and then the 3PL-MHIRT. Treating the test reliability obtained from the MHIRT models as the gold standard (because MHIRT models were the data-generating models), one finds that the test reliability for the second-order latent traits was considerably overestimated by single-level HIRT models; in contrast, the test reliability for the first-order latent traits was slightly underestimated by single-level HIRT models. This biased estimation of the test reliability in single-level HIRT models occurred because the multilevel structure was ignored; thus, items were not locally independent, and the school-level random effect was mistakenly treated as a person-level variance.

Table 3.

Mean RMSE and Mean and SD of Test Reliability for the Latent Trait Estimate in Simulation Study 1.

		Test 1			Test 2			Test 3			Second-Order
Model	Multilevel	RMSE	Reliability	SD	RMSE	Reliability	SD	RMSE	Reliability	SD	RMSE	Reliability	SD
1PLM	MHIRT	0.466	0.895	0.007	0.478	0.894	0.008	0.492	0.886	0.006	0.533	0.740	0.007
	HIRT	0.640	0.871	0.008	0.666	0.873	0.006	0.651	0.868	0.006	0.809	0.804	0.013
2PLM	MHIRT	0.489	0.886	0.008	0.462	0.900	0.010	0.508	0.881	0.006	0.548	0.734	0.011
	HIRT	0.668	0.854	0.009	0.634	0.881	0.010	0.664	0.859	0.005	0.829	0.795	0.018
3PLM	MHIRT	0.529	0.869	0.007	0.510	0.873	0.009	0.562	0.848	0.009	0.555	0.710	0.005
	HIRT	0.718	0.826	0.009	0.679	0.842	0.009	0.726	0.809	0.008	0.817	0.772	0.017
PCM	MHIRT	0.323	0.954	0.001	0.324	0.953	0.002	0.335	0.953	0.002	0.460	0.810	0.007
	HIRT	0.531	0.948	0.001	0.516	0.947	0.002	0.552	0.949	0.002	0.792	0.824	0.011
GPCM	MHIRT	0.337	0.950	0.003	0.311	0.958	0.003	0.337	0.947	0.003	0.464	0.808	0.013
	HIRT	0.527	0.943	0.003	0.515	0.953	0.003	0.517	0.943	0.002	0.794	0.830	0.016
RSM	MHIRT	0.317	0.955	0.002	0.318	0.954	0.002	0.333	0.952	0.002	0.458	0.811	0.010
	HIRT	0.524	0.949	0.002	0.525	0.949	0.001	0.522	0.948	0.002	0.794	0.837	0.009
GRM	MHIRT	0.336	0.951	0.003	0.306	0.960	0.001	0.336	0.950	0.003	0.465	0.807	0.008
	HIRT	0.548	0.944	0.003	0.539	0.955	0.001	0.551	0.945	0.003	0.802	0.837	0.017

Note. HIRT = higher-order item response theory; MHIRT = multilevel higher-order item response theory; RMSE = root mean square error; RSM = rating scale model; PCM = partial credit model; GPCM = generalized partial credit model; GRM = graded response model; 1PLM = one-parameter logistic model; 2PLM = two-parameter logistic model; 3-PLM = three-parameter logistic model.

With regard to the polytomous items, as shown at the bottom of Table 3, the same findings for dichotomous items applied. Under most conditions, the RS-MHIRT had the highest test reliability, followed by the PC-MHIRT and then both the GPC-MHIRT and GR-MHIRT. Ignoring the multilevel structure led to underestimation of the test reliability for the first-order latent traits and overestimation of the test reliability for the second-order latent trait.

Table 4 presents the PPMC with the Bayesian chi-square, reproduced correlation, and SD of the group mean scores for the data generated from the MHIRT models and analyzed with the MHIRT models and single-level HIRT models. When the analysis models were the data-generating MHIRT models, the Bayesian chi-square and SD of the group mean scores always indicated a good model-data fit (the Type I error rate might be too conservative), whereas the reproduced correlation declared a poor model-data fit from 14 to 20 out of 20 replications (the Type I error rate was very high). When the analysis models were single-level HIRT models, only the SD of the group mean scores was able to reject the false model with a power of 100%, whereas the other statistics had zero power. Thus, the SD of the group mean scores can be recommended because it had a conservative Type I error rate and a high power for the multilevel structure. Note that the Type I error rate and power should be interpreted with caution because only 20 replications were performed.

Table 4.

Number of Misfits in Posterior Predictive Model Checking Across 20 Replications.

Model	MHIRT	HIRT
1PLM
Bayesian chi-square	0	0
Reproduced correlation	16	0
SD of the group mean scores	0	20
2PLM
Bayesian chi-square	0	0
Reproduced correlation	16	0
SD of the group mean scores	0	20
3PLM
Bayesian chi-square	0	0
Reproduced correlation	19	0
SD of the group mean scores	0	20
PCM
Bayesian chi-square	0	0
Reproduced correlation	20	0
SD of the group mean scores	0	20
GPCM
Bayesian chi-square	0	0
Reproduced correlation	16	0
SD of the group mean scores	0	20
RSM
Bayesian chi-square	0	0
Reproduced correlation	14	0
SD of the group mean scores	0	20
GRM
Bayesian chi-square	0	0
Reproduced correlation	16	0
SD of the group mean scores	0	20

Note. HIRT = higher-order item response theory; MHIRT = multilevel higher-order item response theory; RSM = rating scale model; PCM = partial credit model; GPCM = generalized partial credit model; GRM = graded response model; 1PLM = one-parameter logistic model; 2PLM = two-parameter logistic model; 3-PLM = three-parameter logistic model.

In addition to PPMC, the Bayesian DIC was computed to compare the MHIRT models and single-level HIRT models. The results (which are not shown because of space constraints but are available on request) suggested that the MHIRT models always had a smaller Bayesian DIC value than the corresponding single-level HIRT models. Thus, it was recommended that the Bayesian DIC be used for model selection and the PPMC be used for model-data fit checking.

Simulation Study 2: Parameter Recovery for MHIRT Models With Covariates at the School Level

Table 5 summarizes the parameter recovery for dichotomous and polytomous items generated from the 3P-MHIRT and GPC-MHIRT, with covariates at the school level (i.e., the 3PL- and GPC-MHIRT-C). In the 3P-MHIRT-C, the mean RMSE was 0.210 for the difficulty parameters, 0.086 for the slope parameters, 0.037 for the pseudo-guessing parameters, and 0.165 for the variance–covariance matrix; the RMSE was 0.013, 0.014, and 0.017 for the three factor loadings and 0.184, 0.147, and 0.180 for the three school-level regression weights. The parameter recovery in the MHIRT-C, although acceptable, appeared to be slightly worse than that in the MHIRT without covariates.

Table 5.

Parameter Recovery for the MHIRT With Covariates in Simulation Study 2.

Model	3P-MHIRT-C		GPC-MHIRT-C
Parameter	Bias	RMSE	Bias	RMSE
Variance–covariance
Mean	−0.040	0.165	0.010	0.144
SD	0.020	0.012	0.055	0.064
Max	−0.005	0.190	0.129	0.267
Min	−0.071	0.154	−0.031	0.059
Difficulty (location)
Mean	0.038	0.210	0.036	0.140
SD	0.123	0.126	0.039	0.040
Max	0.832	0.871	0.146	0.298
Min	−0.131	0.102	−0.047	0.068
Slope
Mean	0.015	0.086	−0.010	0.046
SD	0.031	0.033	0.015	0.018
Max	0.100	0.210	0.025	0.107
Min	−0.049	0.039	−0.054	0.017
Pseudo-guessing
Mean	0.007	0.037	—	—
SD	0.023	0.023	—	—
Max	0.120	0.125	—	—
Min	−0.034	0.006	—	—
Factor loading
β₁	0.006	0.013	−0.006	0.012
β₂	−0.008	0.014	−0.003	0.013
β₃	−0.007	0.017	−0.006	0.016
Regression weight
Test 1	0.092	0.184	0.058	0.193
Test 2	0.009	0.147	0.043	0.123
Test 3	0.124	0.180	0.007	0.206

Note. 3P-MHIRT-C = three-parameter multilevel higher-order item response theory model with covariates; GPC-MHIRT-C = generalized partial credit multilevel higher-order item response theory model with covariates; RMSE = root mean square error.

The same findings applied to the GPC-MHIRT-C. As shown on the right-hand side of Table 5, the mean RMSE was 0.140 for the location parameters, 0.046 for the slope parameters, and 0.144 for the variance–covariance matrix; the RMSE was 0.012, 0.013, and 0.016 for the three factor loadings and 0.193, 0.123, and 0.206 for the three school-level regression weights. In summary, the MHIRT-C for both dichotomous and polytomous items had an acceptable parameter recovery.

Two Empirical Examples

Ability Test With Dichotomous Items

The mathematics assessment of TIMSS 2007 administered to fourth-grade students in Taiwan was demonstrated using the multilevel HIRT as the analysis model. In TIMMS 2007, there were three content domains to measure the mathematics ability: Number, Geometric Shapes and Measurement, and Data Display, each of which was treated as measuring a first-order latent trait. Thus, the three first-order latent traits were treated as being governed by a second-order latent trait, which was referred to as “mathematics ability.” All the multiple-choice items were included and analyzed with different types of dichotomous IRT models. As a result, there were 48 items in Number, 32 items in Geometric Shapes and Measurement, and 14 items in Data Display. A total of 4,131 examinees were recruited to respond to the mathematics assessment, and these examinees were sampled from 150 schools in Taiwan. Whether schools offer enrichment in mathematics to students was used as the school-level covariate; the school was coded as 1 if enrichment was offered, whereas the school was coded as 0 if enrichment was not offered.

We were especially interested in the following questions:

Did the items share a common slope parameter for each test? Was it necessary to incorporate the pseudo-guessing parameters for all the items?

Was there a substantial multilevel structure?

Was the school-level predictor useful?

To answer these three questions, five models were fit to the data. The upper panel of Table 6 lists the posterior expectation of the deviance ( $\bar{D}$ ), the difference between the posterior mean of the deviance and the deviance at the posterior mean of the parameters (PD), and the DIC for these models. To answer Question 1, we compared the 1P-, 2P-, and 3P-HIRT and found that the 2P-HIRT had the smallest DIC of 47510.300; this result indicates that these items had different slope parameters but that the pseudo-guessing parameters were not necessary. To answer Question 2, we compared the 2P-HIRT with the 2P-MHIRT and found that the 2P-MHIRT had a smaller DIC of 47444.200, which suggests that the multilevel structure should be considered. To answer Question 3, we compared the 2P-MHIRT with the 2P-MHIRT-C and found that the 2P-MHIRT-C had a smaller DIC of 47427.000 but a larger PD of 3136.800, which suggests that the 2P-MHIRT-C was not parsimonious. In addition, the school-level regression weight estimates were so small that these values were not significantly different from zero. Thus, the 2P-MHIRT was the final model of choice. According to the PPMC computed by the SD of the group mean scores, the 2P-MHIRT had a good fit because the p value was .41 for tests measuring Number, .37 for tests measuring Geometric Shapes and Measurement, and .48 for tests measuring Data Display, all of which were far different from both 0 and 1.

Table 6.

Comparison of Model-Data Fit in the Two Empirical Examples.

Model	$\bar{D}$	PD	DIC
Example 1 (dichotomous items)
1P-HIRT	44948.300	3108.110	48056.400
2P-HIRT	44386.700	3123.620	47510.300
3P-HIRT	44391.400	3393.070	47784.500
2P-MHIRT	44367.900	3070.160	47438.000
2P-MHIRT-C	44288.600	3136.800	47425.400
Example 2 (polytomous items)
PC-HIRT	257,642	19812.50	277,454
GPC-HIRT	250,444	20879.10	271,323
RS-HIRT	259,560	19801.70	279,361
GR-HIRT	255,472	20381.60	275,854
GPC-MHIRT	250,399	19709.00	270,108
GPC-MHIRT-C	250,404	19708.40	270,112

Note. $\bar{D}$ = the posterior expectation of the deviance; PD = the difference between the posterior mean of the deviance and the deviance at the posterior mean of the parameters; DIC = the Bayesian deviance information criterion; HIRT = higher-order item response theory model; MHIRT = multilevel higher-order item response theory model; MHIRT-C = multilevel higher-order item response theory model with covariates.

Under the 2P-MHIRT, the estimates were 0.29 to 5.21 (M = 1.23) for the slope parameters; −3.78 to 1.75 (M = −1.42) for the difficulty parameters; and 0.98, 0.98, and 0.96 for the factor loadings of Number, Geometric Shapes and Measurement, and Data Display, respectively. The school-level variance was estimated as 0.12, 0.11, and 0.11 for the three respective tests, and the covariance was 0.10 between Number and Geometric Shapes and Measurement, 0.09 between Number and Data Display, and 0.09 between Geometric Shapes and Measurement and Data Display. In summary, a two-parameter HIRT was needed, the multilevel structure that resulted from the school effects was substantial, and the offer of mathematical enrichment by schools had little impact on the performance of examinees in mathematics.

Teaching Assessment With Polytomous Items

Student evaluation of teaching effectiveness is a common practice in higher education. C.-J. Huang (2001) used the Chinese version of Students’ Evaluation of Educational Quality, which was originally developed by Marsh (1987), to investigate and assess teaching quality at the Changhua University of Education in Taiwan. The revised Students’ Evaluation of Educational Quality had 29 nine-point rating scale items, which were clustered into eight scales: Value, Instructor Enthusiasm, Organization, Group Interaction, Individual Rapport, Breadth of Coverage, Examinations, and Assignments. These eight scales can be assumed to measure the “overall teaching effectiveness.” A total of 3,566 students were sampled from 120 classes. Students were treated as Level 1 and classes as Level 2. The teacher’s gender was treated as the Level-2 covariate (male coded as 1 and female as 0).

We were interested in the following three questions:

Did the items have a common slope parameter for each test? Was a common set of threshold parameters sufficient for all the items?

Was there a substantial multilevel structure?

Was the gender predictor useful?

Six models were fit to the data. The results are shown in the lower panel of Table 6. To answer Question 1, we compared the PC-, GPC, RS-, and GR-HIRT and found that the GPC-HIRT had the smallest DIC, which indicates that different items had different slope parameters and their own set of threshold parameters. To answer Question 2, we compared the GPC-HIRT and GPC-MHIRT and found that the GPC-MHIRT had a smaller DIC, which suggests that the multilevel structure was needed. This result was expected because teachers in different classes perform their teaching differently. To answer Question 3, we compared the GPC-MHIRT and GPC-MHIRT-C and found that the GPC-MHIRT had a smaller DIC, which suggests that the predictor of the teacher’s gender was not useful at the class level. The implication was that the teaching performance was not related to the teacher’s gender. Thus, the GPC-MHIRT was the final model of choice. According to the PPMC computed by the SD of the group mean scores, the GPC-MHIRT provided a good fit because the p value was substantially beyond 0 and 1 for the eight scales.

Under the GPC-MHIRT, the estimates were 0.57 to 3.05 (M = 1.54) for the slope parameters and −4.49 to 2.73 (M = −0.92) for the location parameters. The factor loadings were 0.82, 0.87, 0.92, 0.77, 0.79, 0.85, 0.83, and 0.78 for the eight scales. The class-level variance estimates were 0.38 to 0.79, and the covariance estimates were 0.25 to 0.62. In summary, a multilevel HIRT was needed to fit the data to account for both the higher-order latent traits and the multilevel data, and the effect of the teachers’ gender was not sufficient to account for the interclass variation in teaching performance.

Conclusions

Conventional IRT models have been developed to either accommodate hierarchical latent traits or fit multilevel data. In practice, both hierarchical latent traits and multilevel data can occur simultaneously. In this study, we developed a general class of MHIRT models to accommodate both hierarchical latent traits and multilevel data. The underlying IRT models are not limited to a specific function; instead, they can be any IRT model. We used the freeware WinBUGS for the parameter estimation; as a result, no effort was required to derive the parameter estimation procedures and to develop the computer programs. A series of simulations were conducted to evaluate the parameter recovery of the MHIRT and the consequences of ignoring the multilevel structure. The results indicated that the parameters were recovered fairly well by using WinBUGS and that ignoring the multilevel structure led to poor parameter estimation, overestimation of the test reliability for the second-order latent trait, and underestimation of the test reliability for the first-order latent traits.

Three test statistics were investigated for PPMC. Although the findings were not decisive due to a small number of replications, it was clear that using the SD of the group mean scores was helpful in the assessment of the model-data fit. The implications and applications of the MHIRT were demonstrated with two empirical examples, one for ability tests and the other for a psychological inventory. Both the empirical examples showed that the MHIRT models can furnish a better model-data fit than single-level HIRT models and that multilevel structures should be considered.

The simulations were conducted using personal computers with 2.66-GHz Intel Core i5 processors. The simulations took approximately 3 to 10 days per replication. The computation time was feasible for most real data analyses but might not be feasible for comprehensive simulations, which was the main reason why only 20 replications were conducted in this study. To facilitate the parameter estimation of the MHIRT models, one can adopt standard IRT computer programs, such as BILOG-MG, MULTILOG, or PARSCALE, to fit each test at the first order and then use the parameter estimates as starting values for MHIRT models while using WinBUGS. Future studies are encouraged to develop more efficient and user-friendly computer programs.

For simplicity of explanation and simulation, this study focused on two orders of latent traits and two levels of data. In fact, MHIRT models can be easily generalized to more than two orders and two levels. Additionally, a lower latent trait can be governed by more than one higher-order latent trait, and the relationship between a higher- and lower-order latent trait can be nonlinear; these topics are of interest and are left for future study. In this study, the metrics are assumed to be invariant across schools, and the items are assumed to be invariant across different groups of test takers (i.e., no differential item function). Measurement invariance at the item and latent trait levels is an important issue in test development. How MHIRT models can be extended to address this concern requires further investigation. Recently, many cognitive diagnosis models have been proposed to provide profile information about test takers’ binary attributes (Rupp, Templin, & Henson, 2010). Higher-order latent continuous traits and latent binary attributes can be integrated to form the so-called higher-order cognitive diagnosis models (de la Torre & Douglas, 2004). Future studies can be conducted to develop a class of multilevel higher-order cognitive diagnosis models.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The first author was sponsored by the National Science Committee, Taiwan (No. 101-2410-H-133-001) and the second author was sponsored by the Public Policy Research Funding Scheme, Hong Kong (No. 8012-PPR-10).

References

Agresti

(2010). Analysis of ordinal categorical data (2nd ed.). New York, NY: Wiley.

Andrich

(1978). A rating formulation for ordered response categories. Psychometrika, 43, 561-573.

Andrich

(1996). A hyperbolic cosine latent trait model for unfolding polytomous responses: Reconciling Thurstone and Likert methodologies. British Journal of Mathematical and Statistical Psychology, 49, 347-365.

Birnbaum

(1968). Some latent trait models and their use in inferring an examinees’ ability. In Lord

F. M.

Novick

M. R.

(Eds.), Statistical theories of mental test scores (pp. 395-479). Reading, MA: Addison-Wesley.

Bolt

D. M.

Lall

V. F.

(2003). Estimation of compensatory and noncompensatory multidimensional item response models using Markov chain Monte Carlo. Applied Psychological Measurement, 27, 395-414.

Brooks

S. P.

Gelman

(1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7, 434-455.

Cao

Stokes

S. L.

(2008). Bayesian IRT guessing models for partial guessing behaviors. Psychometrika, 73, 209-230.

Curtis

S. M.

(2010). BUGS code for item response theory. Journal of Statistical Software, 36, 1-34.

de la Torre

Douglas

J. A.

(2004). Higher-order latent trait models for cognitive diagnosis. Psychometrika, 69, 333-353.

10.

de la Torre

Hong

(2010). Parameter estimation with small sample size a higher-order IRT model approach. Applied Psychological Measurement, 34, 267-285.

11.

de la Torre

Song

(2009). Simultaneously estimation of overall and domain abilities: A higher-order IRT model approach. Applied Psychological Measurement, 33, 620-639.

12.

Fischer

G. H.

(1983). Logistic latent trait models with linear constraints. Psychometrika, 48, 3-26.

13.

Fox

J.-P.

(2010). Bayesian item response modeling: Theory and applications. New York, NY: Springer.

14.

Gelman

Meng

X.-L.

Stern

H. S.

(1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6, 733-807.

15.

Goegebeur

De Boeck

Wollack

J. A.

Cohen

A. S.

(2008). A speeded item response model with gradual process change. Psychometrika, 73, 65-87.

16.

Hill

P. W.

Rowe

K. J.

(1996). Multilevel modelling in school effectiveness research. School Effectiveness and School Improvement, 7, 1-34.

17.

Huang

C.-J.

(2001). Is student evaluation of teaching performance effective? (NSC Research Report NO. NSC 90-2413-H-018-017). Taipei, Taiwan: National Science Council.

18.

Huang

H.-Y.

Wang

W.-C.

(2013). Higher-order testlet response models for hierarchical latent traits and testlet-based items. Educational and Psychological Measurement, 73, 491-511.

19.

Huang

H.-Y.

Wang

W.-C.

Chen

P.-H.

C.-M.

(2013). Higher-order item response models for hierarchical latent traits. Applied Psychological Measurement, 37, 619–637.

20.

Jiao

Wang

(2013). Estimation methods for one-parameter testlet models. Journal of Educational Measurement, 50, 186-203.

21.

Linacre

J. M.

(1989). Many-faceted Rasch measurement. Chicago, IL: MESA.

22.

Marsh

H. W.

(1987). Students’ evaluations of university teaching: Research findings, methodological issues and directions for future research. International Journal of Educational Research, 11, 253-388.

23.

Masters

G. N.

(1982). A Rasch model for partial credit scoring. Psychometrika, 47, 149-174.

24.

Muraki

(1992). A generalized partial credit model: Application of an EM algorithm. Applied Psychological Measurement, 16, 159-176.

25.

Muthén

L. K.

Muthén

B. O.

(2012). Mplus user’s guide (7th ed.). Los Angeles, CA: Muthén & Muthén.

26.

Plummer

(2003, March). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. Paper presented at the Third International Workshop on Distributed Statistical Computing, Vienna, Austria.

27.

Rasch

(1960). Probabilistic models for some intelligence and attainment tests (Expanded ed., 1980). Chicago, IL: University of Chicago Press.

28.

Raudenbush

S. W.

Bryk

A. S.

(2002). Hierarchical linear models: Applications and data analysis methods (2nd ed.). Thousand Oaks, CA: Sage.

29.

Rupp

A. A.

Templin

J. L.

Henson

R. A.

(2010). Diagnostic assessment: Theory, methods, and application. New York, NY: Guilford Press.

30.

Samejima

(1969). Estimation of latent ability using a response pattern of graded scores. Psychometrika Monograph Supplement (No. 17). Richmond, VA: Psychometric Society.

31.

Sellstrom

Bremberg

(2006). Is there a “school effect” on pupil outcomes? A review of multilevel studies. Journal of Epidemiology and Community Health, 60, 149-155.

32.

Sheng

Wikle

C. K.

(2008). Bayesian multidimensional IRT models with a hierarchical structure. Educational and Psychological Measurement, 68, 413-430.

33.

Sinharay

Johnson

M. S.

Stern

H. S.

(2006). Posterior predictive assessment of item response theory models. Applied Psychological Measurement, 30, 298-321.

34.

Spiegelhalter

D. J.

Best

N. G.

Carlin

B. P.

van der Linde

(2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, Series B, Methodological, 64, 583-616.

35.

Spiegelhalter

D. J.

Thomas

Best

(2003). WinBUGS version 1.4 [Computer Program]. Cambridge, England: MRC Biostatistics Unit, Institute of Public Health.

36.

Wainer

Bradlow

E. T.

Wang

(2007). Testlet response theory and its applications. New York, NY: Cambridge.

37.

Wang

W.-C.

Qiu

X.-L.

(2013). A multidimensional and multilevel extension of a random-effect approach to subjective judgment in rating scales. Multivariate Behavioral Research, 48, 398-427.