A straightforward (single level) regression model, showing the connection between the reading scores of Year 5 students and their corresponding scores when they were in Year 3. Might be formulated like so:

$$y_{i,t} = \beta_0 + \beta_1*y_{i,t-1} +e_i \tag{1}$$

Where $y_{i,t}$ represents the reading test score of the $i^{th}$ student at time $t$ (e.g., in Year 5, 2021), $yi,t-1$ is the reading test score of the $i^{th}$ student at a previous time t - 1 (e.g., in Year 3, 2019).

I am tasked with developing a multi-level model to estimate the 'value add' of individual schools that the students are nested within. The value add model is a common approach for quality improvement in education. In theory, allow a fair comparison of schools by controlling for student-level factors and school-level contextual factors, and estimate the random-effects component for each school to provide a ranking of positive/negative impact on the 'typical' students' performance.

This is given by

$$ y_{i,j} = \beta_{0} + \beta_{1}*X_{1ij} + ... + \beta_{n}*X_{nij} + \beta_{n+1}*W_{1j} + ... \beta_{n+k}*W_{kj} + u_{0j} + e_{ij} \tag{2} $$

where $X_{1ij}$ to $X_{nij}$ are $n$ student level controlling factors, which include the student's prior score. $W_{1j}$ to $W_{kj}$ are $k$ school level controlling factors and $\beta_{0}$ is the average performance of all students, conditioned on the student and school factors.

In the context of school effect analysis, $u_{0j}$ is of primary interest as it signifies the contribution each school makes to its students' learning beyond what can be predicted from student (e.g., background, prior academic achievement) and school characteristics (e.g., demographic and academic composition). $u_{0j}$ is assumed to follow a normal distribution with a mean of zero and variance of $\sigma_{\mu0}^2$ across the population of schools.

Typically, a value added model will estimate random effects for each school (usually random intercept) for the focal year. For example, with it being 2024 we would fit the model on this year's results to derive random effect estimates for this year to then rank order schools and differentiate between schools that benefit student learning growth from those that do not. But what if we want to know if the random effect for School X has improved over time? Or if the schools, within the school system, are improving in general, given by the VA values, stable or declining in general? Typically VA values are mean-centred or scaled to the focal year and therefore our 'average' reference school belongs to the focal year.

I am now considering a 'benchmarking' approach by fitting the linear model in equation $(1)$ to historical student performance data (2017), and using the estimated fixed-effect coefficients (eg for prior scores) to 'set' the fixed-effects component in equation $(2)$ for the focal year (2024). I am able to do this using this approach: How to manually set coefficients for variables in linear model? [duplicate] but just updating the function like:

setCoeffs <- function(frml, benchmark_weights, len){
  el <- paste0("offset(", benchmark_weights[-1], "*", 
           unlist(strsplit(as.character(frml)[-(1:2)], " +\\+ +"))[-(2)], ")")
el <- c(paste0("offset(rep(", benchmark_weights[1], ",", len, "))"), el)                       
as.formula(paste(as.character(frml)[2], "~", 
                   paste(el, collapse = " + "), " + (1|school_number) + -1"))

First fitting the linear model ($1$) to historical student performance data:

mod <- lm(score ~ score_previous_year
                       , data = data
                       , subset=completed_year==2017

Reserve the fixed-effect coefficients from $(1)$:

benchmark_weights <- mod$coef

Then fitting the mixed-effects model $(2)$ to the focal year:

frml = <- score ~ score_previous_year + (1|school)
mod_mixed <- lmer(frml 
                       , data = data
                       , subset=completed_year==2024
                       , REML = FALSE)

And then update the fixed-effects component with the 'benchmark' fixed effects for the 'typical' student in 2017:

benchmark_mod_mixed   <- update(mod_mixed , setCoeffs(frml, benchmark_weights, nrow(data)))

And this gives two mixed-effects models that yield different value-add values per school:

va_values <- extract_random_effects(mod_mixed)
va_values_fixed <- extract_random_effects(benchmark_mod_mixed)

Taking the very first school as an example, the fitted 'value add' estimate for School A using 'mod_mixed' where the fixed effects are estimated by the focal year (2024), is -13. Whereas with the 'benchmark_mod_mixed' models, which borrow the fixed effect estimates for the 'average or typical student' from 2017 would be -27.4.

Given the first value, and ignoring for a moment the CIs (both negative, under zero line), the takeaway would be that the school, with reference to all other schools in the focal year (2024) has had a negative impact on student growth. But given the second value, we would say that this school's value add for a typical student's growth, benchmarked to 2017 is worse. Furthermore, if we had had the value add for the same school in 2017, we would see whether or not the school has improved or declined in performance over time. But we wouldn't be able to know this with the original approach that uses the focal year's student performance data to estimate the 'typical student' performance.

My main question about this is whether this is a reasonable approach to evaluate the quality of an individual school's value add over time. I would also like to know if there are some considerations regarding update.merMod and if the function is appropriate given that it fixes the intercept and fixed effects coefficients from the LM to the MLM. Also, suggestions for alternative approaches would be welcome. I am currently reviewing modernVA package which considers time series in the value add modelling approach as well as other interesting avenues to explore further. All cautions, critiques and considerations warmly welcomed.

Use case

The following results may help further illustrate the intended use case of a value add model that 'benchmarks' to 'typical' student performance set by a historical cohort period (in this case 2015-2017, year 3 to year 5 results). Here I set the MLM fixed effect coefficients to the Benchmark model or allowed the MLM to re-estimate the typical student performance each year:

enter image description here

As expected the centre of the distribution of VA values for the Relative model are zero, whereas the centre of VA values for the Benchmark model varies above and below the zero line. I'm left wondering if the results validly indicate that there has been systemic improvement or decline in student performance over time (eg based on state-wide investment or effects of the pandemic) and if the Benchmark model results for VA values (2019,2021) are appropriately estimated.

Does this approach address the critique provided in this paper (page 13):

"while value-added measures are likely to be useful for individual schools aiming to improve their performance, all value is added relative to the average school. This means that it is generally not possible for all schools to achieve positive values, or to use the measures to track changes in systemic performance over time. Alternative formulations of the value-added measure (such as measuring value-added relative to a school in a baseline year, or measuring the change in schools’ value-added scores) could be used to create a benchmark of performance, but further work is required to validate these models."

Or is there a better way?

  • $\begingroup$ Thanks for providing all this detail. I am not sure about using a lm model to benchmark the 2017 data. In theory, you could run a value added model on the 2017 data and the 2024 data and compare the estimates, e.g., by running correlations across time, but you need to "reflate" the EB predictions (or run a multivariate model). See this article by Leckie (2018) in JEBS (doi.org/10.3102/1076998618755351). There is a Stata program for doing this (reffadjust; see stata-journal.com/article.html?article=st0327), but it should be doable with a little coding in R. $\endgroup$
    – Erik Ruzek
    Feb 9 at 14:40
  • $\begingroup$ Thanks Erik, probably what I'm looking for is more closely linked to the criticism of the value add modelling approach on how it might not track systemic improvement over time as VA values always wrap around zero (see page 13) education.nsw.gov.au/content/dam/main-education/about-us/… $\endgroup$
    – mkrasmus
    Feb 9 at 22:32
  • $\begingroup$ I see. From your question, I believe you are aware that the 0 the model is wrapping around is the weighted population mean ($\gamma_{00}$) based on your multilevel model. So you can add the model-estimated population mean to the EB prediction. Then you have a non-zero mean. At any rate, it's not clear to me why the Leckie (2018) multivariate approach doesn't suit your needs. $\endgroup$
    – Erik Ruzek
    Feb 9 at 22:43
  • $\begingroup$ Hi Eric, just extended with a use case for what I'm trying to achieve in an attempt to follow on from the suggestion in the paper noted earlier. Hope this clarifies more on what I am hoping to gain from the approach. Maybe a comparison of correlations between Relative and Benchmark VA values could help to determine if the approach is appropriate (i.e., would at least expect a strong positive correlation). I will look at Leckie for support with that. $\endgroup$
    – mkrasmus
    Feb 10 at 6:57

1 Answer 1


Note that to estimate the random effect, you need the fixed effects and the variance components. In particular, for a random intercepts model, the formula is $$\hat{b}_i = \frac{n_i \sigma_b^2}{\sigma^2 + n_i \sigma_b^2} \frac{1}{n_i} \sum\limits_{j = 1}^{n_i}(y_{ij} - x_{ij}^\top \beta),$$ where $n_i$ is the number of measurements for group $i$, $\sigma_b^2$ the variance of the random intercepts, $\sigma^2$ the measurement error variance, $y_{ij}$ the data at group $i$, and $x_{ij}^\top \beta$ the fixed-effects part.

The following code illustrates how this can be done using lme4:

fm <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)

# data for a subject
data_i <- sleepstudy[sleepstudy$Subject == 308, ]
# outcome for the subject
y_i <- data_i$Reaction
# number of measurements
n_i <- length(y_i)
# fixed effects part for subject
mu_i <- predict(fm, newdata = data_i, re.form = NA)
# variance components
vc <- VarCorr(fm)
# random intercepts variance
sigma2_b <- vc$Subject
# measurement error variance
sigma2 <- attr(vc, "sc")^2

# random effect estimate
mean(y_i - mu_i) * n_i * sigma2_b / (n_i * sigma2_b + sigma2)

# check that is works using the ranef() that only works 
# for subjects in the dataset used to fit the model
ranef(fm)$Subject[1L, ]

Hence, you can fit a model for a specific focal year (e.g., 2017) and, based on this model, estimate the random effects using the 2024 data.

  • $\begingroup$ Thank you @Dimitris especially for outlining those useful steps in the code, a really great learning experience. Extending a bit more to get CIs I use L95_i - b_i - 1.96*(sqrt(sigma2_b * sigma2 / (sigma2_bn_i + sigma2))) and U95_i <- b_i + 1.96*(sqrt(sigma2_b * sigma2 / (sigma2_bn_i + sigma2))) $\endgroup$
    – mkrasmus
    Feb 14 at 5:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.