Bayesian Methods for Scalable Multivariate Value-Added Assessment

Specification of Variable Persistence Model for a Single Subject

McCaffrey et al. (2004) provided a general longitudinal model for T years of annual test score data for a single subject, where throughout this article subject refers to an academic content area such as mathematics or reading rather than an individual student. A special case of this model that includes teacher effects but not school effects is the following:

The test scores are denoted Y _it for student i’s score in year t, t = 1, . . . , T. The model includes an overall mean μ _t for each year, and covariate vectors x _it that can include both time invariant and time varying background variables. The model also includes teacher effects θ _t for each year. For consistency with the literature and for simplicity in presentation we use the term teacher effects when describing the random components included at the classroom levels. The effects of interest are not necessarily causal effects or intrinsic characteristics of teachers. Rather, they account for unexplained heterogeneity at the classroom level. Ideally, they provide information about teacher performance, but there might be many sources of this heterogeneity, including omitted student characteristics (McCaffrey et al., 2004; McCaffrey, Lockwood, Koretz, & Hamilton, 2003). Furthermore, the estimated teacher effects are limited in scope to the subject matter aligned to the examination that generated the observed scores and to the appropriateness and validity of the vertical linking procedures employed by the test developer (Reckase, 2004).

The teacher effects are linked to students by the vectors φ _it . Elements of these vectors measure the share of the students’ instruction provided by each teacher and take on values in the interval from 0 to 1 (inclusive) with Σ _j φ _itj = 1. Zeros indicate a student had no instruction from the corresponding teacher, exactly one element equal to 1 indicates the corresponding teacher provided the student with all of his or her instruction, and fractional values indicate a student received instruction from multiple teachers either because of transfer during the year or team teaching. The validity of the chosen fractions would need to be verified for each application, but the model can easily accommodate multiple instructors. For ease of presentation and because our data did not require fractional linkages, for the remainder of the article, we assume that exactly one element of each φ _it equals 1 and the rest are 0.

The value of α _{tt ^*} for t ^* = t is defined to equal 1 and is used only to make the notation in Equation 1 more compact. The values of α _{tt ^*} for t ^* < t moderate the covariances of scores in year t among students who shared a teacher at a prior time t ^* by modeling the persistence of current teacher effects in subsequent years of testing. A special case of particular interest to the current article is when α _{tt ^*} ≡ 1 for all t ^* ≤ t, which matches the assumptions made by the complete persistence model. The residual error terms (ɛ _{i 1}, . . . , ɛ _iT ) are assumed to be multivariate normal with mean vector 0 and an unstructured variance-covariance matrix, allowing for different variances at each time point and possibly nonzero and nonconstant correlation of scores from different years (grades), thus making efficient use of the repeated measures at the student level.

Multivariate Specification

The model in Equation 1 is for longitudinal outcomes for a single academic subject such as reading or mathematics. More commonly, students are tested in multiple subjects per year, and it may be of interest to examine teacher effects for these different subjects. Because outcomes for the same student are likely to be positively correlated across subjects, a model for the joint distribution of all outcomes is substantively more appropriate than independent models for different subjects, and exploiting the covariances among all outcomes may lead to reduced bias and improved efficiency in estimates of model parameters (Thum, 1997, 2003).

We generalize the model in Equation 1 to multiple tested subjects per test administration as follows. Y _i denotes the vector of test scores for student i. Y _i is of length ST, the number of subjects (S) times the number of years (T), and is organized by subject, and then by test administration (year) within subject. X _i denotes the (ST × p) design matrix of both time-invariant and time-varying student background variables for the p-dimensional vector of regression coefficients

where the subscripts on the vectors denote subject and year. That is, the multi-subject generalization of the model is constructed as parallel copies of the single-subject model, with all parameters from the single-subject model being expanded to include separate components for each subject. The teacher effects are expanded along with the parameters of the mean structure so that each teacher has different effects for each subject. We organize the teacher effects by subject and year as θ′ = (θ′₁₁, . . . , θ′_{1 T} , θ′₂₁, . . . , θ′_{2 T} , . . . , θ′ _{S 1}, . . . , θ′ _ST ) of length n _θ, where θ _st provides the teacher effects for subject s in year t. The matrices Φ _i generalize to multiple subjects the assignment vectors φ _it by specifying the linkages of students to teachers by subject. Specifically, Φ _i is (ST × n _θ) with only 0 or 1 entries and row sums equal to 1, with the nonzero element in each row corresponding to student i’s teacher for a given year and subject. The contribution of teacher effects to the outcomes for student i is then given by A Φ _i θ, where A is a (ST × ST) block diagonal matrix consisting of S distinct (T × T) lower triangular blocks corresponding to subjects. The (t,t ^*) element of the block for subject s is α _{s,tt ^*} for t ^* ≤ t and 0 otherwise, where α _{s , tt ′} denotes the teacher effect persistence parameters for subject s.

The distribution for a single student’s score vector, Y _i , conditional on the model parameters, teacher effects, and all covariate and linkage information is

where N _ST denotes the ST dimensional multivariate normal distribution and Σ is a (ST × ST) unstructured positive definite covariance matrix. Outcomes for different individuals are assumed to be conditionally independent given all of these parameters. The teacher effects for each subject s and timepoint t are treated as random effects. We assume that given the variance τ θ , s t 2 , the components of θ _st are independent and identically distributed N ( 0 , τ θ , s t 2 ) . The teacher effects are assumed to be independent across subjects and timepoints. More general covariance structures for the random effects (e.g., allowing correlation within teachers across subjects) are possible but are not explored here.

Computational Challenges

Fitting either the single subject or multisubject model to test score data presents considerable computational challenges. Standard mixed-model routines (e.g., those available in R, S-plus, SAS, HLM, and MLWin) are in general not equipped to estimate the persistence parameters. Moreover, even if the persistence parameters are fixed, as in the complete persistence model, the multiple membership structure makes likelihood estimation difficult for realistically sized data sets. Unlike hierarchical models, cross-classified or multiple membership structures prevent reduction of the marginal covariance matrix of the scores into the simple forms that allow for fast computation (Rasbash & Browne, in press). Without these simple forms, likelihood estimation requires inversion of a sparse and generally large matrix. Although specialized likelihood methods for crossed random effects have been proposed (Clayton & Rashbash, 1999; Rasbash & Goldstein, 1994) and additional methods are in development (Bates & DebRoy, 2004; DebRoy & Bates, 2003; Rasbash & Browne, 2001), to date, only the work of Sanders and colleagues (Ballou, Sanders, & Wright, 2004; Sanders et al., 1997) has been able to scale likelihood estimation of the complete persistence model to large data sets. A scalable implementation of the variable persistence model does not exist.

Recently, Bayesian methods have been proposed as a fruitful alternative approach to the computational challenges presented by crossed random effects (Browne & Draper, 2006; Browne, Draper, Goldstein, & Rasbash, 2002; Browne et al., 2001; Rasbash & Browne, in press; Simonite & Browne, 2003). The utility of Bayesian methods for these problems derives from using Markov Chain Monte Carlo (MCMC) sampling algorithms, particularly successive substitution sampling or “Gibbs” sampling (Carlin & Louis, 2000; Gelman, Carlin, Stern, & Rubin, 1995; Gilks, Richardson, & Spiegelhalter, 1996). The key feature of MCMC algorithms that makes them well suited to handling models with complex relational structures is that sampling the joint posterior distribution of all unknown parameters requires only repeatedly sampling from the conditional posterior distributions of certain parameters given the data and values of all other parameters. Conditioning on random effects reduces the complex covariance matrices to simple, computationally tractable block diagonal forms. Similarly, conditioning facilitates the estimation of the persistence parameters, which can be viewed as unknown regression coefficients of known predictors conditional on the random effects.

Bayesian Formulation of the Variable Persistence Model

To capitalize on these advantages, we developed a Bayesian version of the multisubject (or single-subject) persistence model in Equation 3. Given the likelihood, the Bayesian model requires a joint prior distribution for all unknown parameters. The unknown parameters for the model in Equation 3 are μ, θ, α, Σ as well as the variance components τ 2 = { τ θ , s t 2 } . We used the following distribution with fixed hyperparameters μ _μ, V _μ, {δ_{θ, st} }, {ν_{θ, st} }, μ _α, V _α, q _Σ, Q _Σ:

Here N _d again denotes a multivariate normal distribution of dimension d, W denotes the Wishart distribution, and U(· |ν, δ) denotes the uniform distribution on the interval (ν, δ). Also the counts n _{θ, st} denote the lengths of the vectors θ _st . The joint distribution of the data and parameters is given by the product of Equation 3 across all students and the prior in Equation 4. The posterior distribution of the parameters is proportional to this product. Note that this specification uses uniform distributions on the standard deviation components rather than gamma distributions on the inverse variance components, the latter leading to closed form full conditional distributions, which are gamma distributions (Gelman et al., 1995). We used the uniform priors because they allow for a somewhat more natural introduction of substantive prior information into the model and also because “noninformative” gamma priors have infinite density at zero, which can lead to undesirable properties of the posterior distribution and associated sampling algorithm.

We implemented this model in WinBUGS (Spiegelhalter, Thomas, & Best, 1999), free software for Bayesian model fitting available at http://www.mrc-bsu.cam.ac.uk/bugs/welcome.shtml. The WinBUGS implementation was sufficiently fast to estimate the variable persistence model with small data sets and the complete persistence model with small or large data sets. However, fitting the variable persistence model with large data sets or multiple subjects was prohibitively slow. Therefore, we developed an MCMC algorithm (specified in detail in the appendix) for sampling the posterior distribution and programmed it in C, resulting in an efficient and scalable implementation of the multisubject variable persistence model. We use this algorithm for all results reported in this article.

We conducted an extensive simulation study to verify the correctness of the algorithm and its scalability to many teachers and students. We simulated hundreds of data sets from the multisubject persistence model using various configurations of the model parameters, including the persistence parameters, various amounts of missing test scores, and varying numbers of student-level covariates. In all cases the algorithm was able to recover reasonable estimates of the model parameters. Computing times on a standard desktop PC required for the MCMC algorithm to converge sufficiently and to obtain samples for inference were extremely reasonable even for large problems. Depending on the complexity of the simulated data set (i.e., the numbers of students, teachers, years, subjects and regressors, the extent of missing data, and whether or not the α parameters are estimated), the computing times ranged from on the order of several minutes for problems involving a few hundred teachers and several thousand students to a few hours for a huge problem of fitting the complete persistence model for five subjects, 5 years, considerable missing test score data, and following 50,000 students across 50,000 teacher effects (2,000 teachers per subject per year). Fitting the variable persistence model in such a large data set, with five subjects, increases the computing time by about a factor of 6, and we are working on ways to refine the updating of the persistence parameters to reduce this burden. In all cases, RAM requirements are quite low, even for data sets with tens of thousands of students and thousands of teachers.

Identifiability of Persistence Parameters

A technical consideration of the variable persistence model is whether the multiplication of the unknown persistence parameters α and teacher effects θ makes them unidentified. The issue is further complicated in the Bayesian framework where parameters may be identified only by the prior distribution. Fortunately, assuming that classroom groupings change over time, the data appear to provide sufficient information to uniquely identify α and θ. From a theoretical standpoint, the fact that α _tt is defined to equal 1 for every t (that is, the teacher effects in year t show up in the model for the year t scores without an unknown coefficient) implies that the signs and scales of the elements of θ can be estimated without confounding from α and vice versa. Our empirical investigations support this notion. In particular, in the simulations described earlier, the model recovered values of α near 1 when the true values were 1 and values distinctly less than 1 when the true values were less than 1. Moreover, the persistence parameter estimates were not sensitive to prior specification or arbitrary rescaling of the data; this was true in both the simulated data and the actual data used in this study. Finally, we revisited the estimates of persistence obtained via likelihood methods by McCaffrey et al. (2004) and were able to recover those estimates using the Bayesian methods developed here.

Data Description

The data contain vertically linked mathematics and reading scale scores on a norm-referenced standardized test administered during the spring of the years 1998 to 2002. For this analysis, we focused on the cohort of students who under normal grade promotion would have been in Grade 1 during the 1997-1998 school year and Grade 5 during the 2001-2002 school year. We thus estimated effects for teachers of Grade 1 during the 1997-1998 school year, Grade 2 during the 1998-1999 school year, Grade 3 during the 1999-2000 school year, Grade 4 during the 2000-2001 school year, and Grade 5 during the 2001-2002 school year. Teacher links to students did not vary by subject; the same teachers were linked to student mathematics and reading scores.

A total of 10,332 students in our data linked to these teachers. However, some of these students had no valid test scores or other problems such as unusual patterns of grades across years that suggested incorrect linking of student records. We deleted records for these students. The final data set includes 9,295 students with 168 unique observation patterns (patterns of missing and observed test scores for both subjects over time), and only about 20% of the students had fully observed scores.

For our analyses we standardized the test scores by subtracting 400 and dividing by 40. We did this to make the variances approximately one and to keep the scores positive with a mean that was consistent with the scale of the variance. Keeping the variance near one made the choice of priors less complex. Keeping the scores positive allowed us to easily code missing values as the only negative values.

Prior Distributions

We used the following specifications of the hyperparameters of the prior distributions. μ _μ was set to the means of the national norming sample for the test, standardized by subtracting 400 and dividing by 40, for consistency with our rescaling of the data. However, the variance V _μ was set to a diagonal matrix of 1,000s to allow the means to be driven by the data. We used independent uniform distributions on the range of [ν_{θ, st} , δ_{θ, st} ] = [0, 2] for the priors for the teacher standard deviation components. Given our rescaling of the data, this range conservatively avoided restricting the teacher variance components. We set μ _α = 1 because previous researchers (Ballou et al., 2004; Rowan et al., 2002; Sanders et al., 1997) have assumed this value for their analyses. However, V _α was set to a diagonal matrix of 1,000s so these values are determined essentially entirely from the data. Q _Σ was determined by the inverse of a covariance matrix that assumed a correlation of .7 among student residuals from the same subject, a correlation of .7 between mathematics and reading scores from the same year, and .49 correlation between mathematics and reading scores from different years. These correlations were stylized from analyses we have conducted with numerous student test score databases in our previous research. The standard deviations were set to the values from the national norming sample divided by 40. The degrees of freedom parameter q _Σ was set to the number of years times the number of subjects plus one (11 in this case). This specification makes the prior minimally informative while maintaining a bounded density over all positive definite covariance matrices.

In addition to these prior specifications, we generated the posterior distribution for the parameters from one model specification using various more informative prior distributions for all parameters as well as overdispersion of the initial teacher effects. The posterior distributions were unaffected as measured by virtually perfect correspondence of the posterior means and standard deviations of all parameters. Finally, as noted previously, gamma priors on the inverse teacher variance components would lead to gamma full conditional distributions for these parameters. Some of our earlier implementations of the model code used this specification rather than the uniform priors on the teacher standard deviations reported here. The posterior distributions of the model parameters were unaffected by this choice.

Estimation

The algorithms previously described and presented in the appendix provide a sample from the posterior distribution for the teacher effects θ and other model parameters. Any number of features of these distributions such as means, standard deviations, quantiles, and probabilities assigned to different sets might be used to summarize the distributions and compare the parameter estimates across different models. For our comparisons of teacher effects, we focus on two primary summaries of the posterior distributions. We examine the posterior means of the teacher effects E(θ| y ), which are the analogs to the best linear unbiased predictors (BLUPs) produced by classical mixed-model estimation procedures. We examine the correlations of the posterior means of the teacher effects across different models, with high correlations indicating that the models are ordering the teacher effects nearly the same. To encompass both central tendencies of the effects and their uncertainties, we also examine Pr(θ > 0| y ), which can be used to identify teacher effects that are “extreme” (either positive or negative) relative to the mean effect of zero. For the comparisons presented in this article, we characterize an effect as extreme and positive if Pr(θ > 0| y ) ≥ 0.95 and extreme and negative if Pr(θ > 0| y ) ≤0.05 (note that given our model for the teacher effects, Pr(θ = 0| y ) = 0). These cutoffs may or may not be appropriate in policy applications, depending on the desired inferences, but are sufficient to quantify the sensitivities that we are exploring. The sensitivity of a teacher effect to a model choice is defined for our purposes as the teacher effect being flagged as extreme under one model and either being flagged as extreme in the opposite direction or not flagged as extreme under another model. These discrepancies should be interpreted with caution as the classifications can be sensitive to small changes in the posterior probability distributions under different models.

For each model that we examined, we used our MCMC algorithm to generate a sequence or chain of parameters sampled from the posterior distribution from that model. We “burned in” each chain for 5,000 iterations and based our inferences on 10,000 post-burn-in iterations. We diagnosed convergence of the chains using the Gelman-Rubin diagnostic (Gelman & Rubin, 1992) implemented in the coda package (Best, Cowles, & Vines, 1995) for the R statistics environment ( R Development Core Team, 2005). Using five parallel chains, 5,000 burn-in iterations were clearly sufficient for the ostensible convergence of all the model parameters, including the teacher effects. The posterior means and standard deviations that we report are based on all 10,000 post-burn-in iterations, whereas other inferences regarding teacher effects are based on 1,000 evenly spaced parameter vectors from the sample of 10,000. Replication of these calculations on different posterior samples obtained from independent chains with dispersed initial values indicated that these posterior sample sizes were sufficient for stability of the inferences.

Implementation Issues

The complexities of the data forced us to address several challenging issues in implementing the model. We discuss four of these: incomplete test score data for a majority of students, missing student-teacher links, nonzero means for estimated teacher effects, and the availability of both reading and mathematics data.

Incomplete Test Score Data

As noted, only 20% of the students had complete testing data for both reading and mathematics over the 5 years for which they were followed. Fortunately, the Bayesian framework provides a particularly simple method for implementing a missing at random (MAR; Little & Rubin, 2002) missing data model known as data augmentation (Schafer, 1997; Tanner & Wong, 1987; van Dyk & Meng, 2001). In essence, missing data are treated as any other unknown parameter, being sampled from their conditional distributions given the observed scores and the values of model parameters (details are provided in the appendix). The inferences about the parameters of interest thus automatically account for the additional uncertainty arising from the missing test score data, provided that the MAR assumption is reasonable.

Missing Student-Teacher Links

Although the Bayesian approach naturally handles missing test score data via data augmentation, the model as specified requires teacher links for all students for all subjects and all years. However, records missing test scores were typically missing teacher information, most likely because the students were not in the school district at the time of testing. How we treated missing teacher links depended on the missingness pattern of the student. Because the accumulation of the teacher effects assumed by the model implies that future teacher effects are not related to past scores, the potential future teacher effects on students who drop out are not required for estimation. We thus assume that all future teacher effects for students who drop out are zero. More care has to be taken with teacher effects that are missing prior to the last observed score for a student. The model for observed test scores depends on these missing links through the persistence of prior teachers, and thus the link must be specified to estimate the model parameters.

We implemented three procedures for dealing with the missing link information. The first and most naive method, which we call M1 throughout the remainder of the article, assumes that like unobserved future effects for students who drop out, unobserved effects prior to dropout are all zero. The other two methods assign a teacher effect to each unobserved predropout link that is distinct from the actual teacher effects in the data and occurs at the level of subject within year within student. We refer to these as pseudo-teacher effects. The pseudo-teacher effects are treated in the same manner as real teacher effects (real referring to actual teachers rather than actual values of the unknown effects), accumulating over time to determine the mean structure for the outcomes. Like the real teacher effects, the pseudo-teacher effects are treated as independent random effects with variance components that vary as a function of year and subject. The two methods that use pseudo-teacher effects differ in their assumptions about these variance components. The simpler method (M2) assumes that the pseudo-teachers for a particular year and subject share the variance component of the observed teachers for that year and subject. This is notionally similar to assuming that the unobserved teacher link is drawn from the population of appropriate teacher effects in the school district, but the pseudo-teacher effect is independent of the real effects. The other method (M3) allows the pseudo-teacher effects to be drawn from their own distributions by year and subject, and the variance components for these distributions are estimated separately from those of the real teachers. This relaxes the assumption that the unobserved effect comes from the same distribution as the real teachers. It reflects the fact that we know nothing about the type of education the students are exposed to when we do not observe them and acknowledges the possibility that these students were in another school district with teachers that differed markedly from the current district. It is important to note that the pseudo-teacher effects are not intrinsically interesting but rather are used only as a tool to obtain more valid estimates of the effects for the real teachers.

Nonzero Means for Estimated Teacher Effects

An additional complication worth noting is that missing teacher links and missing test score data can combine to cause the estimated teacher effects (i.e., the posterior means) to be centered at a small positive value for each grade rather than at zero. This was observed in the actual data, and we reproduced the effect using simulated data as well. The basic intuition behind this phenomenon is that missing test score records, which as noted also tend to be missing teacher links, come from students who tend to score lower. When including the missing teacher effects as part of the total population of teacher effects, which by definition has mean zero by grade, the mean effect for the real teachers is ostensibly positive. This mean shift has no impact on the ordering of real teacher effects but can affect inferences about whether or not effects are positive. We judged that practitioners most often would be interested in knowing whether particular teacher effects were larger than the average teacher effect in the population of teachers being modeled, rather than larger than zero, as the latter question can never really be answered (because the mean structure is always based on the population of students and teachers in the data). Thus, we based our inferences on teacher effects using the posterior distributions of the centered effects θ^* = (θ − θ̄ _t ), where θ̄ _t is the average effect of real teachers in year t. For example, E(θ^*| y ) is the posterior mean teacher effect relative to the average effect for that grade, and Pr(θ^* > 0| y ) is the posterior probability that a teacher effect is better than the average effect for that grade.

The Availability of Both Reading and Mathematics Data

A final implementation issue arises from the multivariate structure of the data. We allow each teacher in the data to have a separate effect for mathematics and reading, but these effects can be estimated either by separate univariate longitudinal models (e.g., using Equation 1) fit to each subject or by a joint model (using Equation 3) that makes use of the cross-covariances among the subjects and estimates the effects simultaneously. An interesting and important empirical question is to what extent inferences about the effects may be sensitive to this modeling choice.

Taken in summary, although the primary comparison of interest regards how assumptions about persistence might affect inferences, the empirical issues about how missing teacher links are handled and whether or not the full multivariate structure of the data is exploited are also of interest. We explore these issues by estimating both mathematics and reading effects for all teachers under 12 different settings obtained by a full crossing of the following three factors:

Sensitivity to Methods for Missing Teacher Links

The posterior means E(θ^*| y ) for the teacher effects were extremely robust to the assumptions about missing teacher links. For the variable persistence model, the correlations for the three sets of estimates exceeded .99 for every combination of grade, subject, and modeling approach (marginal or joint). Moreover, graphically the estimates align nearly perfectly along the 45 degree line (not shown), indicating the estimates also share a common variance. Estimates for the other model parameters were also quite robust to the assumptions about missing links, and the models provided very similar estimates of Pr(θ^* > 0| y ). For the complete persistence model, the minimal correlation of estimated teacher effects across missing link methods also was .99. However, there was an additional sensitivity for this model in that missing link method M2 yielded substantially smaller estimates of the teacher variance components than the alternatives for Grade 2 teachers for both mathematics and reading. Fortunately this had a minimal impact on the quantities of interest such as Pr(θ^* > 0| y ).

Given these findings of robustness to the methods for missing teacher links, the remaining results are all based on missing link model M2, which offers a reasonable compromise between the restrictiveness of model M1 and the tenuous information about the variance components of unobserved teacher effects required by model M3. However, the restriction to model M2 did not affect the substantive conclusions of the subsequent results.

Sensitivity to Joint Versus Marginal Modeling of Subjects

Jointly modeling the test scores for both subjects versus modeling them marginally with independent replicates of the single-subject model has some effect on the parameter estimates, but the resulting teacher effects are quite similar. The joint model reduces the estimates of the persistence parameters α and the estimates of the teacher variance components τ 2 = { τ θ , s t 2 } for both subjects. The estimated persistence parameters from the marginal models range from about 10% to 80% larger and average 40% larger than those from the joint model. The differences in persistence parameters tend to be largest for the first grade parameters (e.g., α _{m ,21}, α _{m ,31}, etc.); otherwise no clear patterns exist in the differences. Similarly, the teacher variance components τ ² are uniformly larger for the marginal models than for the joint model, with the variance components for the marginal variable persistence model between 10% and 30% larger than those for the joint variable persistence model. The discrepancies in the complete persistence model are somewhat larger, with the ratios of τ ² from the marginal models to those of the joint model decreasing from about 2.7 for Grade 1 down to 1.1 for Grade 5. We suspect that the differences between these parameters when estimated by marginal rather than joint models results from the joint model exploiting the additional information on the students to more fully separate student effects from teacher effects.

Fortunately, the estimated teacher effects show very little sensitivity to the choice of joint or marginal modeling. The correlation between E(θ^*| y ) under the joint and marginal models exceeds .99 for the variable persistence model and .97 for the complete persistence model for both subjects and all grades. As we would expect, the posterior standard deviations V^1/2(θ^*| y ) of the teacher effects are nearly always smaller for the joint model than for the marginal model because pooling information across subjects for students provides more precision in estimating teacher effects. However, the increase in precision for joint modeling is generally modest, particularly for teachers in later grades. For the variable persistence model, V^1/2(θ^*| y ) from the marginal model average 7%, 6%, 5%, and 3% larger than those from the joint model for math in Grades 2 through 5, respectively; the analogous percentages for reading are 8%, 4%, 4%, and 2%. The differences in the complete persistence model are more pronounced, with average percentage increases of 20%, 14%, 7%, and 3% for math in Grades 2 through 5 and 16%, 11%, 6%, and 2% for reading in Grades 2 through 5.

The increased precision of the teacher effects under the joint model along with its impacts on the estimates of the other model parameters combine to make the classification of teachers as extreme weakly sensitive to the choice of joint or marginal modeling. Joint models flag somewhat fewer teachers as extreme than do marginal models, and the complete persistence model is slightly more sensitive than the variable persistence model. However, for Grades 2 to 5 these differences are rather small.

Sensitivity to Assumptions About Teacher Effect Persistence

The primary goal of our application was to gather additional information about the persistence parameters and sensitivity of results to these parameters. We report our comparison of the complete persistence model (α ≡ 1) to the variable persistence model (α estimated from the data) using the joint reading and mathematics model with missing teacher link model M2 (pseudo-teacher effects sharing variance components with real teacher effects). However, as noted, the substantive conclusions hold for all models.

Figure 1 presents posterior means and 95% posterior intervals for the components of α estimated from the variable persistence model. The most important finding is that the α parameters imply long-term persistence of past effects with significant decay in the strength of that persistence over time. All estimates are positive and substantially less than one, with posterior probabilities of being between zero and one approximately equal to one. Thus, the empirical results do not support the assumptions of the complete persistence model. This finding of small, positive estimates of α is consistent with the empirical results presented in McCaffrey et al. (2004) obtained with a different data set.

The disparity between the degree of persistence assumed by the complete persistence model and that estimated by the variable persistence model affects the estimated teacher effects, variance components, and inferences about teachers. Figure 2 plots E(θ^*| y ) for reading as estimated by the variable persistence model against the corresponding estimates from the complete persistence model, separately by grade. The corresponding comparisons for the mathematics effects are similar. As shown in the figure, the estimated effects from the alternative models of persistence are moderately to strongly positively correlated, but with notable differences especially in earlier years. Moreover, the data do not lie on the diagonal, except for Grade 5, as the estimates from the variable persistence model are substantially more dispersed than the corresponding estimates from the complete persistence model. The difference is most pronounced in Grades 2 and 3, and it is also evident in Figure 3, which presents “caterpillar” plots of third-grade reading teacher effects estimates from the two alternative persistence models. The figures provide a quick, high-level view of the strength of the signal at the teacher level and the rough proportion of teachers whose effects are likely to be different from zero. As clearly shown by the plots, the estimates from the variable persistence model are more variable, and more effects are distinct from the mean.

The increased dispersion in the estimated teacher effects from the variable persistence model results from dramatic differences between the estimated teacher variance components τ ² from the two models. The posterior means for τ ² for the variable persistence model are .35, .26, .22, .13, and .17 for mathematics and .41, .16, .18, .13, and .13 for reading, which range from a minimum of 17% to a maximum of 47% of the corresponding posterior means of the residual variance (i.e., the diagonal elements of Σ) for each year and subject. The analogous estimates for the complete persistence model are .06, .05, .06, .08, and .12 for mathematics and .06, .05, .05, .07, and .10 for reading, which range from 5% to 20% of the corresponding residual variance for each year and subject. The posterior standard deviations of τ ² are all .04 or less, so the differences are significant. The complete persistence model requires smaller τ ² in general to maintain consistency with the marginal variance of the data (which is approximately constant over time). The small values of α allow τ ² in the variable persistence model to be large in general without conflicting with the data. The elements of Σ also differ between the two models. In particular, the posterior means of the residual error variances for each subject are substantially larger for the complete persistence model in Grade 1 (0.94 vs. 0.75 for mathematics, 1.15 vs. 0.91 for reading). By Grades 4 and 5, the estimates from the variable persistence model tend to exceed those from the complete model by about 0.04; however, these differences tend to be within the posterior standard deviations, which are about 0.02.

Table 1 summarizes how these differences affect teachers’ effects identified as extreme. If the two models perfectly corresponded in their identifications, all counts would fall on the diagonals of the (3 × 3) tables. This is not the case. Generally, the variable persistence model flags substantially more teachers as extreme (both positive and negative), which is consistent with the larger estimated teacher variance components. Particularly important is the fact that a small number of teacher effects (given by the upper right-most cells in each table) are flagged as extremely positive by the variable persistence model and extremely negative by the complete persistence model.

The differences in the estimated teacher effects arise because the two models use the data in markedly different ways. As is evident from the model specification, α _{tt ^*} ≡ 1 implies that gain scores (i.e., Y _it – Y _{i , t – 1}) in a given year are a function of educational inputs only through the current year teacher, and thus gain scores provide most of the information about the teacher effects for the complete persistence model. Alternatively, if α _{tt ^*} ≡ 0, then past teacher effects do not persist into the future, and current year level scores provide most of the information about current year teacher effects. When α _{tt ^*} takes on other values, both level scores and gain scores depend on both current and past teachers, and the model makes more complicated use of the data to estimate effects; however, if α _{tt ^*} is very small, then current year levels score are likely to be the primary source of information about teacher effects.

Figure 4 demonstrates how these heuristic notions of the differences between the complete and variable persistence models manifest in our application. The figure is based on reading teacher effects, but the results for math teacher effects are similar. For each teacher, we calculated the average annual reading score for the students linked to that teacher so that each teacher in each grade is associated with five average reading scores summarizing the performance of his or her students, not only in the grade g that the teacher instructs but in each of Grades 1 to 5. We then regressed the estimated teacher effects on these annual average scores. The figure plots the standardized coefficients by grade and model. The R ² of each regression is presented on the horizontal axis. Although the regression on aggregate scores is substantially less complex than the actual model because it uses only the observed test scores and makes no adjustment for the effects of teachers in other years, the fact that all of the R ² values exceed .8 indicates that the estimated effects are largely identified by simple contrasts. For the complete persistence model, the coefficients for the grade g regression are dominated by a large positive value in grade g and large negative value in grade g − 1, consistent with the intuition that the complete persistence model identifies teacher effects via average gain scores. On the other hand, the coefficients for the variable persistence model are dominated by the current grade average score. This is consistent with the small positive values of α _{tt ^*} estimated from the model.

Empirical Findings

The empirical findings presented in this article using our Bayesian models have several important implications. Our estimates of persistence of teacher effects suggest that scores remain correlated for students who shared a teacher in the past. This finding is consistent with the hypothesis of persistent teacher effects. However, similar to the results obtained by McCaffrey et al. (2004), the estimates are significantly less than 1 and thus do not support the assumption of complete persistence made by the Sanders et al. (1997) layered model and the Raudenbush and Bryk (2002) cross-classified model. Although the teacher effects from the different models are relatively strongly correlated, the assumptions about persistence have substantive implications for teacher effects being identified as extreme. Further study is necessary to examine which of the models provides “better” estimates of teacher effects as they use the data in notably different ways.

Fortunately, the inferences from the models were largely insensitive to the other implementation choices examined here. Somewhat surprisingly given the degree of missingness in the data, the inferences from the models were almost invariant to the different models for missing teacher links. This robustness may help to alleviate some concerns about the viability of VAA but should be examined in other empirical settings before general conclusions can be made.

Joint versus marginal modeling of the subjects similarly does not appear to exert strong influence on the estimated effects, though the joint models do result in a slight increase in precision for the estimated effects. The fact that the joint model reduces the estimated teacher variance components and persistence is interesting and warrants further consideration; a possible interpretation is that part of what the marginal models are identifying as teacher variance is actually unmeasured student heterogeneity that is more completely removed by the joint models.

Our prior specification for teacher effects assumed that a teacher’s effects for reading and mathematics were independent, and alternative priors might provide more precision in the estimated effects. This choice to treat effects as independent prevents potential bias that would result from shrinking the teacher’s reading and mathematics effects toward one another when they are disparate, for example, when a teacher effect is significantly larger for one subject than the other. However, the correlation of the posterior means of the effects were .56, .53, .55, and .64 in Grades 2 to 5, respectively, for the complete persistence model and systematically larger at .71, .73, .70, and .76 for the variable persistence model. Clearly the data suggest that teacher effects are correlated and a prior that captures that correlation could allow for sharing information across reading and mathematics, increasing precision of the estimates for both subjects.

There are a few important limitations to the empirical results presented here. First, we do not include any student background variables in our models. As described in McCaffrey et al. (2003), omitted covariates can result in bias, and hence our results might be biased by such variables. We do not adjust for covariates because the current usage of the layered model excludes these variables, and we felt it important to contrast models in a manner consistent with current usage. Also, many data sets that might be used for estimating teacher or school effects might have very limited student measures; for example, the only background variable included in our data is a race-ethnicity indicator. Future work should include studies that compare the models when covariates are included.

Another limitation of our study is the MAR assumption. There is some evidence (Dunn, Kadane, & Garrow, 2003) that students who miss testing will tend to score lower than others even after controlling for prior test scores. We might find that different models for missing links or the persistence of effects might be more or less sensitive to this assumption. Future work should explore alternative assumptions about missing data.

Advantages of the Bayesian Approach

Although the empirical questions surrounding the persistence of teacher effects ultimately motivated our exploration of Bayesian methods, as shown earlier and as previously noted by Thum (2003), the Bayesian paradigm provides several important benefits for VAA in general. First, as noted, the sequential sampling methods used to make inferences from Bayesian models allow complex cross-classified and multiple membership models for teacher effects to be taken to scale. This advance in capability is practically useful as the demand to perform VAA estimation with large, complicated data sets continues to grow. Second, the Bayesian framework deals naturally with rich latent parametric structures. In a manner similar to the way in which Bayesian estimation methods simplify the treatment of complicated covariance structures, the methods can leverage the sequential sampling approaches to deal successfully with potentially complex functions of multiple latent parameters (e.g., the persistence parameters α).

The benefits of the Bayesian framework go beyond its computational advantages. The inferential structure of posterior probability distributions for all parameters, including teacher effects, provides estimates of the uncertainty of the effects of interest that account for uncertainty about all unknown parameters, including variance components. Furthermore, the Bayesian framework facilitates communication of results to policy makers and other consumers of VAA in intuitive terms and straightforwardly allows inferences to be tailored to specific decision analyses that might be required in a given application of VAA. For example, Pr(θ^* > 0| y ) provides an accessible and simple summary of teacher effects, but inferences about more complex facets of teacher performance, such as Thum’s (2003) teacher productivity indices, are similarly easily obtained.

Our future work will include exploring more efficient algorithms for estimating the persistence parameters. In addition, the ease with which the Bayesian framework handles latent parameters provides promise for extending the class of VAA models to include models with richer latent parametric structure. For example, models could include not only latent teacher effects but also latent student effects (akin to the cross-classified model of Raudenbush & Bryk, 2002). The additional latent structure could be used to explore more realistic missing not at random models for student test scores and to explore models that parameterize teacher effects as functions of latent student effects.