Multilevel Structural Equation Models for the Analysis of Comparative Data on Educational Performance

One-Dimensional Latent Variable Models for Student Performance

The standard psychometric procedure for the modeling of test item responses is commonly known as item response theory (Lord, 1980). A simple, basic, latent-trait model of this type relates the responses on a set of test items to one or more underlying latent abilities. A basic version can be expressed as follows.

For a student (i) who responds to item (r), the probability (π _ri ) of a correct response is given by

where g is a link function, typically the logit, and the response, y _ri , is 1 if the item is correctly answered and 0 if not and the y _ri are mutually independent. This is just a binary factor model with a single factor (θ) and a set of loadings (λ _r ). We refer to the term β_{0 r} , often referred to as the “facility” for Item r, as belonging to the fixed part of the model, which will later be augmented with further predictors. In the PISA data, we have some graded or partial-credit items where the correct responses are either partially correct (coded 1) or fully correct (coded 2). In this case, for such an item, the first line of the model can be written, for a response coded s (s =0, 1), as

where π r i f is the probability of a response in Category f . We here model the cumulative probabilities where cumulation is from Category 0 (incorrect) rather than in the binary case where, by convention, the first category is the correct response. Thus, in the fixed part of the model, the probability of an incorrect response is simply g ⁻¹(β_{0 r} ), and the probability of an incorrect or partially correct response is g ⁻¹(β_{0 r} + α₁).

An important assumption in Models 1 and 2 is that the responses y _ri are conditionally independent. Because some of the items involve responses to the same text or figure, it is possible that this assumption will be violated, as the (conditional) probability of a correct response to one item may depend on the outcome with respect to an earlier item. Thus, for example, Items R104Q01, R104Q02, R104Q05, and R104Q06 all refer to a passage about telephone use and feature the same person (Pedro) in each question.1 To avoid this problem, an alternative is to consider the complete set of item responses, treating the number of correct responses to this set of four binary items as a single-item-graded response. Thissen, Steinberg, and Mooney (1989) discussed this, and Steinberg and Thissen’s (1996) study patterned item combinations as elementary response units (testlets). In the above example, we could, for instance, compute a total score ranging from 0 to 4. Scott and Ip (2002) considered an alternative approach to what they term the “item clustering” effect. For each specified item cluster, they added an individual-level random effect designed to identify an individual’s additional response to each item as a member of that cluster (see below). We have not pursued these possibilities here, but they are an interesting area for further work.

A special case of Model 1 is the so-called Rasch model, where the loadings λ _r are constrained to be equal. This is the model, with a logit link, used in the PISA analyses.

In this article, we use a probit link rather than a logit (see, e.g., Lord & Novick, 1968, chap. 16). The two link functions are, in fact, very similar so that we can expect resulting estimated probabilities to be very close. The probit, however, has certain advantages computationally and also has a useful interpretation in terms of an underlying normal “propensity” distribution for the responses. Thus, for a binary response, we can suppose that there is an underlying continuous response for an item with a threshold value (X) such that responses above that value are correct and those below are incorrect. Formally, we write the probability of a correct response as

where φ is the standard normal density. This model is discussed, among others, by Fox and Glas (2001) and by Goldstein and Browne (2004) and is here extended to the ordered category case given by Model 2 (see the appendix for a complete specification).

The above model can be extended in several ways. First, we can add further fixed-part explanatory variables such as gender, country, age, and so on. Second, we can make the model multilevel by recognizing between-school variation and explicitly incorporating school-level latent variables or factors. Third, we can add further factor dimensions along which student responses can vary. Fourth, we can allow the loadings to be functions of explanatory variables. Finally, we can allow the factor values or scores also to be functions of explanatory variables. It is possible to extend the model to consider more general structural equation models, but we shall not pursue this here. Steele and Goldstein (2006) gave a further discussion and an application to women’s status data.

Generalizing the notation of Goldstein and McDonald (1988) and McDonald and Goldstein (1989), a basic multidimensional two-level model for a continuous normal response is given by Model 3. In the case of binary or ordered responses, this models the underlying propensity as defined above, and our exposition thus will be expressed in terms of the following basic model.

where h indexes the fixed-part explanatory variables, R is the number of responses, F the number of Level 2 factors, and G the number of Level 1 factors. The v f j ( 2 ) , v g i j ( 1 ) are respectively sets of common factors at Levels 2 and 1 with corresponding uniquenesses u _rj , e _rij . Note that correlations between the factors are allowed, although we do not fit correlated factors in the present article. Where the response is binary or ordered, then, for the underlying propensity distribution, we have e _rij ∼ N(0, 1). The subscript r refers to the item, i to the student, and j to the school. In this article, we shall assume that the diagonal terms of Ω₁, Ω₂ are 1. The alternative is to fix one or more loadings and allow Ω₁, Ω₂ to be general covariance matrices to be estimated. Goldstein and Browne (2004) gave a further discussion, and the steps needed for both formulations are described in the appendix.

Several other authors have studied Model 3 and extensions to it. Zhu and Lee (1999) and Fox and Glas (2001) used Markov Chain Monte Carlo (MCMC) estimation, the former based on Gibbs sampling for a single-level factor model, whereas the latter authors consider the binary response two-level model with a single factor at Level 1 and use Gibbs sampling with a probit link function.

A number of authors have extended Model 3, using maximum likelihood procedures, to include categorical responses and more general structural equation formulations. Thus, Muthen (1997) considered applications to latent growth curve modeling, and a more general discussion of these models is also given by Muthen (1989, 2002). More recently, a very general framework for multilevel structural equation modeling is provided by Rabe-Hesketh, Skrondal, and Pickles (2004) that includes most of the models to be discussed below. These authors obtain maximum likelihood estimates, typically based on quadrature, although other authors (e.g., Raudenbush, 1995) use an expectation maximization (EM) algorithm. Song and Lee (2004) fit this model for mixtures of normal, binary, and ordered responses using a mixture of Gibbs and Metropolis-Hastings sampling.

An advantage of carrying out the estimation using MCMC methods, as in the present article, is the ability to incorporate prior information in a Bayesian sense and to provide exact interval estimates for parameters or functions of parameters. Also, because of the modularity of the algorithm, it is possible to add additional complexity relatively straightforwardly, including the possibility of incorporating distributional assumptions other than normality. The algorithm described in the appendix assumes diffuse priors but is readily extended to incorporate informative prior distributions. It extends previous work in particular by allowing for missing data and parameter constraints among fixed coefficients. It also proposes the use of the Deviance Information Criterion (DIC; Spiegelhalter, Best, Carlin, & Van der Linde, 2002) that provides a measure of model complexity and can be used for comparing nonnested models. In the single-level case, the DIC is analogous to the Akaike Information Criterion (AIC) and can be considered a generalization of this. Unlike other model fit procedures such as Bayes factors, it does not require improper (diffuse) priors such as are used for some of the parameters.

The procedures have been implemented using MATLAB (Mathworks, 2004) and MLwiN software (Rasbash, Browne, & Steele, 2004). Although MLwiN has some basic facilities for multilevel factor modeling, the algorithm as written in MATLAB is more flexible, although computationally slow.

In the next section, we describe the analysis models used in PISA and conduct some simple comparisons using the above formulation, before describing our more detailed analyses.

The OECD PISA Models

The unidimensionality assumption was used by PISA to determine the structure of the tests as well as the subsequent modeling. The analysis essentially involves two stages (Adams & Wu, 2002). The first stage consists of fitting Models 1 and 2 to the complete data set using equal loadings. This provides estimates of the intercept parameters β_{0 r} . Treating these estimates as known parameter values, a further unidimensional Rasch model is fitted to the responses, but this time including as fixed explanatory variables (conditioning variables) a set of scales formed from a principal components analysis of the data in the student questionnaire. These data include items related to both the school and social background of the students. This is done separately for each country. To take account of the uncertainty in the factor scores, an estimate of the posterior distribution of the factor score for each student is obtained, and five values are sampled randomly from this distribution for each student. For example, if MCMC were used for the analysis, these could be five (approximately) independent values from the chain of factor scores for each student, obtained by selecting values suitably far apart in the chain. Alternatively, values can be selected from parallel chains. These plausible values are then used in subsequent analyses to compare countries and so on. Essentially, five analyses are performed, each one using just one plausible value for each student, and the analysis estimates are then combined to obtain inferences. Mislevy (1991) described the procedure.

Certain problems arise with this approach. The first is that although the plausible values may be expected to perform reasonably well for models that use a subset of the conditioning data, this will not generally be true for variables not in the conditioning data set. This includes school-level variables and in particular applies to multilevel analyses where school is the higher level unit (see Mislevy, 1991). The second problem is that the uncertainty attached to the intercept parameter estimates is not taken into account, although it is not clear whether this is a serious problem. A related issue is that the intercept parameter estimates (item difficulties) are based on a model that does not include the conditioning variables. The PISA analyses use sampling weights that reflect the achieved sample characteristics. In the present article, we shall ignore these because they appear to make only small differences.

In subsequent sections, we show how a fully efficient analysis, avoiding the use of plausible values, can be performed. Our first simple analysis, however, looks at the assumption of equal loadings. We fitted Model 1 across the whole data set with and without the equal loadings constraint, and the results are given in Tables 1 and 2, with the loading and intercept estimates together with their standard errors. For these and subsequent analyses, we used a “burn in” of 1,000 and a subsequent 5,000 iterations for the chain. For the main parameters of interest, the fixed-part coefficients and the loadings, the mixing of the chains is reasonable, as it is for the residual variance estimates. The mixing for the threshold parameters is poor using the Gibbs sampling approach but is satisfactory when Metropolis-Hastings sampling is used (see appendix, Step 1, Ordered Responses of the algorithm). Factor scores are estimated by computing the mean of the chain values for each student. The variances of these means are also estimated from the chain, and the inverses of these variance estimates are used as weights in Table 3 to provide a valid analysis comparing the French and English mean scores; standard errors of parameters are computed using sandwich estimators. This analysis and that in Table 4 are thus just weighted least squares regression models where we have followed the PISA procedure of computing a factor score for each student and then using these in subsequent analyses.

It is evident from inspecting the estimates and their standard errors that there is strong evidence for differential loadings. This is confirmed by comparing the values of the DIC. In the present case, there is a reduction of 664 for the model with unconstrained loadings. The two assumptions for the loadings also provide somewhat different fixed-part parameter estimates for the difference between France and England, although in both cases the 95% interval includes zero. The point estimates for the equal loadings model have a variance about two thirds that for the unequal loadings model, although the results are similar in terms of effect sizes to those obtained using PISA procedures (ONS, 2002). Table 4 extends the comparison by including a gender effect and an interaction between gender and country. This shows a slightly smaller difference in favor of girls in England than in France (although not reaching the formal 5% significance level) with similar results in terms of effect sizes for both models.

We now compare this procedure with a single-model approach where the explanatory variables are included in the factor model. In this and subsequent analyses, we use the model with unconstrained loadings. Thus, for a model with the single additional explanatory variable country (denoted by the dummy variable x ₁), we have the single-level model

Note that if we were to adopt the item cluster effect model of Scott and Ip (2002), Model 4 would become, for cluster k,

Instead of a constant β₁, we could allow a different country coefficient for each response (β_{1 r} ), but interest here lies in an overall comparison between countries so that these are constrained to be equal. We discuss how to interpret this parameter below.

The OECD analyses essentially fit the equivalent of the following model:

albeit with equal loadings, where the term (X ^{( c )}β^{( c )}) _ij represents the effect of the conditioning variables. However, the two distributional assumptions in Model 5 will not in general both be satisfied; one case where they are satisfied is where {ν, x, ν^*} have a joint multivariate normal distribution. Instead, the following model avoids this problem:

which on substitution gives

This is a structural equation model, and Model 6 is not equivalent to Model 4. It is an example of a Multiple Indicator Multiple Indicator Cause (MIMIC) model (Muthen, 1989). More generally, the structural part of the model will include further variables of interest at individual and school levels. A choice between Models 6 and 4 in any particular case can be determined by the model that has the better fit to the data. We note, however, that the interpretation of Model 4 is not completely straightforward because the responses have different variances; only the Level 1 residual variances are equal. Thus, although the coefficients are equal across responses in Model 4, the response distributions themselves are not standardized. We therefore use models with structural predictors based on Model 6 in the following analyses.

In the present case, fitting Model 4 to the data gives an estimate for β₁ of 0.009 (0.019); with equal loadings, this is −0.002 (0.011). Fitting Model 6 gives an estimate for α₁ of 0.014 (0.026) with a DIC value of 89,145.5, compared with a rather similar value of 89,142.8 for Model 4.

Note that in the (Rasch) case of equal loadings, Models 4 and 6 do become equivalent. We also note that for the case of more than one factor at a level and also for factors at more than one level, we may wish to model each factor score as a function of explanatory variables, in which case Model 4 is inadequate and we must use models that are an extension of Model 6.

We now look at the case of the two-level model, where Model 4 becomes

and the PISA analyses fit essentially the (variance components) model,

where again we have the problem that the distributional assumptions cannot in general simultaneously be satisfied. An alternative two-level model is

which becomes, on substitution,

Comparing Model 10 with Model 4, we see that effectively the Level 1 and Level 2 loading vectors are constrained to be equal and there is no additional Level 2 residual term in Model 10. Thus, Model 10 becomes a highly constrained model. We show below, however, that an extension of Model 10 does provide a useful model. Thus, Table 5 shows a simple two-level model fit from which it is clear that the Level 1 and Level 2 loading vectors are very different. Furthermore, even with equal loadings, Models 8 and 10 are not equivalent.

In the next section, we develop Model 7, introducing further explanatory variables, and in a later section, we look at the dimensionality structure of the data.

One-Dimensional Models With Several Explanatory Variables

In view of the problems of differential lengths of schooling discussed in the introduction, we fit in our initial models the interactions between age and country. Age is categorized as a dummy variable January to August versus September to December births. The additional use of grade is problematic. For example, if we compare the September to December births in the two countries, all the English are in Grade 10, whereas the repeating French are in Grade 9. If we only compare Grade 10 pupils, then the French will tend to do relatively better because of the strong negative association between performance and repetition. For the children born in January to August, all the English are in Grade 11, whereas the French nonrepeaters are in Grade 10, and the repeaters are in Grade 9. For these reasons, we have not used grade in our comparisons, but a special analysis of the effect of grade in the 2003 PISA survey in France will be reported elsewhere.

Table 6 extends Model 6 by including age-group, gender, and the first-order interactions of age-group, gender, and country; interactions between age trends and country are negligible and not displayed.

We see that there is an advantage to the older pupils and a negative trend with month of birth for the period January to August, with the older children scoring higher and a gender effect in favor of girls. There is little evidence for a trend for the period September to December, or for a country difference, or for any interactions, except possibly for country by gender.

Two-Level Models

We now fit the full two-level factor model with structural predictors and a single factor at each level together with just gender and age terms. No interactions are significant, and fitting a coefficient for country also gave a high standard error for the England–France difference, as well as a rather badly mixing chain with very high serially correlated values. Running the chain for longer provided no evidence that the coefficient was significant. The results are presented in Table 7.

There remains a large gender difference in favor of girls and an advantage to the older pupils. There is a negative trend with month of birth for the period January to August and also evidence for a positive trend for the period September to December. The latter seems difficult to explain.

Exploring Dimensionality

We now explore the dimensionality structure of the data. We have performed a series of analyses at a single level that establishes the existence of at least two dimensions. In the PISA analyses, items were a priori selected for membership of the three separate proficiencies, with each item identified with just one proficiency. The items used in our analyses of the Retrieving Information proficiency subscale are therefore unique to that scale. In the present analyses, we fit orthogonal factors so that we can detect dimensions along which countries may differ meaningfully (see Steele & Goldstein, 2006, for an example with correlated factors).

In common with all factor models, we have choices to make to ensure identifiability. In the models of this article where Ω₁, Ω₂ are identity matrices, a simple procedure at any given level is to set, for the jth factor (j > 1), λ _k = 0, k =1, …j − 1 (Goldstein & Browne, 2004). In Table 8, we do this to fit two Level 1 factors, setting the first item of the second factor to zero. We start with a model including just the intercepts and a fixed effect for country varying across responses. This model (DIC =87,028.9) provides a better fit and somewhat different loading estimates compared with a model fitting intercepts only (DIC =87,687.3) and a much better fit than the basic model for one factor (DIC =88,471.3).

We see clear evidence in Table 8 for at least two factors. If we fix all the loadings below 0.2 to 0 and reestimate, we obtain the results in Table 9, with a somewhat higher value of DIC (87,312.6).

The interpretation of factors estimated in this way is problematic because a different choice of zero loading will, in general, lead to different loading patterns. In fact, using different starting values, we find that the loadings are not stable, moving from one factor to the other. To explore the various possibilities is time-consuming, and we have not done this because our principal aim is to see whether more than one dimension exists; further factors can be fitted in similar ways, however. Another approach would be to fit simple-structure models where each item loads on only one factor at each level, but the factors are allowed to be correlated. This involves choosing appropriate subsets of items, and we have not pursued this. An exploration of the factor space will need to make choices about the loadings to be fixed based on substantive considerations of item formats, positioning, and content. In addition, when carrying out such an exploration, we should fit a two-level model and fit explanatory variables such as age and country and also allow for the possibility that factor structures may vary across countries. It may also be useful to carry out exploratory analyses separately at each level based on a separate modeling of estimated Level 1 and Level 2 residual covariance matrices (see, e.g., Rowe, 2003).

A Constrained Multilevel Structural Model

Rather than fitting the following Model 11 and estimating the loadings for each change in model parameters, we can consider fitting a standardizing model such as Model 6 and subsequently treating the posterior mean estimates of the loadings as fixed in further analyses. The advantage of this approach is that we are dealing with essentially the same factors, as defined by the loadings, in each analysis. Reestimating the loadings for each fitted model will complicate interpretation.

Clearly, various choices for the standardizing model are possible—for example, fitting a 2-level structure with the loadings for each response constrained to be equal across levels. In practical applications, sensitivity analyses can be performed to see whether inferences are strongly affected by different choices.

We present here only the results from fitting a single factor, but the model can be extended by fitting structural parameters in the case of more than one factor.

where the structural predictors z _k are distinct from the fixed-part predictors x _h and the Level 2 random effect is incorporated into the Level 1 structural model. We may allow different variances for different groups at both Level 1 and Level 2, and in the present case, we fit different country variances at both Level 1 and Level 2. The second line of Model 11 becomes

When Models 11 and 12 are combined, because the Level 1 loadings are assumed known, we have the random coefficient factor model

We fit Model 13 with the Level 1 loadings of Table 5 and the intercept and structural predictors of Table 7 without the gender and country interactions (which are not significant at the 5% level), and we obtain the results in Table 10.

We have also fitted Model 13 where the predictors are in the fixed part of the model rather than the structural part (Table 11).

We note that for the structural model, the coefficients tend to be smaller compared to their standard errors than for the fixed-part predictor coefficient model, and the latter is also a better fit with a DIC of 88,559.2, compared with 88,573.1 in the structural model.

We see that the ratio of Level 2 to Level 1 plus Level 2 factor variances, the variance partition coefficient (VPC; Goldstein, Browne, & Rasbash, 2002), is 21% for England and 49% for France. These values are similar to those presented in PISA (Adams & Wu, 2002). Goldstein (2004) suggested that the explanation for the high value for France is that the data contain a mixture of pupils from Grades 9 and 10. As pointed out above, repetition implies a greater variation among schools. We have therefore conducted an analysis for France only using Model 13 and fitting only intercept terms, and we find that for Grade 9 (collège) pupils, the Level 2 variance estimate is 0.19, and for Level 1, the variance estimate is 0.79, giving a VPC of 19%. For lycée pupils, the variances are 0.14 and 0.54 with a VPC of 21%; thus, the VPC estimate for each school type is close to the English estimate. This explanation for the apparently high between-school variation also accounts for results from the Trends in International Mathematics and Science Study (TIMSS; Mullis et al., 2001), which show similar values for the two countries and where the sampling for France was carried out only in collège.

Discussion

Our analyses and discussion have shown that comparisons between two educational systems with different pupil progression structures are complex. The combination of different ages of starting school and different allocation to year groups on the basis of birth date and repetition of grades makes any meaningful comparison extremely difficult. Although we have here compared only England and France, our view is that the same problems occur when nonrepetition systems are compared with those that have important percentages of repetition, such as those of Spain, Portugal, and Belgium.

We have demonstrated that, even within a single proficiency domain, the data structure appears to contain at least two dimensions, although we have not conducted a full multilevel analysis of the dimensionality structure, nor have we attempted to identify and label factors as such. Nevertheless, even in the one-dimensional case, the (Rasch) assumption of equal item loadings is not supported by the data.

We have shown how a valid multilevel factor model can be fitted and, in particular, how to structure the factor variances at each level in order to properly study between-school variability. Model 13 is an example of a random coefficient factor model that can readily be extended to include further explanatory variables such as gender or age and also, for example, to cross-classifications. Thus, the full range of multilevel modeling procedures can be incorporated into these analyses, and such analyses will often lead to inferences that differ from those based on single-level models. The procedure is also much simpler than the plausible value procedure proposed by the OECD because it requires only a single fitting of a multilevel model. An issue with this approach is the choice of loadings to use. In our case, we have chosen a set of loadings from a two-level model with a single factor at each level. Other choices are possible, such as including fixed predictors in the initial model. In general, it would be useful to perform sensitivity analyses to determine whether such choices substantially affect inferences. Once the loadings are chosen, they effectively define the latent factors, and it is meaningful to make comparisons across subgroups only if we then use the same set of loadings in all analyses. We have not taken into account the uncertainty in the estimates of the loadings. Rather, we take the view that the first stage that determines the values of these loadings provides a practically useful metric for further analysis. Nevertheless, it is important to have a suitably large sample to ensure that sampling variability is small. If necessary, we can incorporate prior information, for example, from previous studies, into the estimation of these parameters.

Finally, although the main thrust of this article is to present a methodology for handling complex multilevel data in comparative studies, we should not ignore the serious drawback of a lack of longitudinal data in surveys such as PISA and other similar surveys such as TIMSS. Without such measures of prior performance on the same sample of students, it is not possible to overcome the comparability problems that arise from the different ways in which educational systems are organized, as we have described. Likewise, without such prior measures, it is not possible to attribute any observed differences between systems or subgroups to the education systems per se rather than, for example, social, cultural, or economic factors. Goldstein (2004) discussed this issue in more detail in the context of the stated aims of the PISA study.