1 Current methodology for measuring distance to SDG targets

A key argument of this article is that the progress countries are making to achieve their SDG targets generally, and in education particularly, can be conceptualized as a rate of change. However, the major policy reports on country-level progress to the SDG targets presents the problem as one of measuring the distance to the target of interest. A review of the extant reports on measuring the distance to the SDG targets have revealed differences in approaches depending on whether the report emanates from the OECD or the UN. In this section, we review the approaches taken by both organizations along three lines (a) selection of indicators, (b) specification of target levels, (c) calculating distance measures, and (d) estimating trends. Our review of the methodologies used by the OECD is taken from OECD (2019). Our review of the methodologies used by the UN are taken from Lafortune et al. (Reference Lafortune, Fuller, Moreno, Schmidt-Traub and Kroll2018)).

2 Critique of the OECD and UN methodologies

The importance of the OECD and UN reports notwithstanding, there are several problems with their methodology for computing the distance measures that need to be raised.

3 Bayesian growth curve modeling: Specification, identification, and estimation

In this section, we outline the specific form of the growth model that we will employ in our empirical example. In anticipation of our example, we will contextualize the problem as one of predicting the pace of progress in country-level mathematics proficiency as measured by PISA. To begin, we write the intra-country model of the pace of progress in mathematics outcomes as

(1) $$ \begin{align} y_{ci} = \pi_{0i} + \pi_{1i}C_{i} + r_{ci}, \end{align} $$

where $y_{ci}$ is mathematics outcome for country i ( $i = 1, \ldots n)$ at PISA assessment cycle c ( $c=1 \ldots C$ ), $\pi _{0i}$ is the intercept capturing country i’s proficiency level at a fixed point in time (usually the first cycle—for this article, 2009), $\pi _{1i}$ is the slope (linear pace of progress over time) in mathematics proficiency for country i at cycle c, and $r_{ci}$ is the residual term.Footnote ¹ Together $\pi _{0}$ and $\pi _{1}$ are random effects, typically referred to as growth parameters in the growth curve modeling literature (e.g., Bollen & Curran, Reference Bollen and Curran2006; Grimm et al., Reference Grimm, Ram and Estabrook2017). For this article, $\pi _{0i}$ is the estimated 2009 country-level achievement for country i, which we will refer to simply as the starting point, and $\pi _{1i}$ is an estimate of the pace at which country i is increasing or decreasing in mathematics proficiency. Our goal is to obtain optimal estimates of these parameters to be used to predict future mathematics proficiency scores. The linear growth curve model in Equation (1) can be extended to include time-varying predictors as

(2) $$ \begin{align} y_{ci} = \pi_{0i} + \pi_{1i}C_{ci} + \beta_{ci}z_{ci} + r_{ci}, \end{align} $$

where $z_{ci}$ is a predictor that also changes over time with the outcome and is hypothesized to predict the outcome at cycle c, but is not an outcome of interest in itself. The parameter $\beta _{ci}$ describes the size and sign of the prediction. An example of a time-varying predictor might be the number of math teachers in country i at cycle c.

It is well known that higher order terms can be specified, such as quadratic pace of progress, which can be written as

(3) $$ \begin{align} y_{ci} = \pi_{0i} + \pi_{1i}C_{ci} + \pi_{2i}C^2_{ci}+ r_{ci}, \end{align} $$

and indeed, the latent variable modeling approach to growth curve modeling allows for considerable flexibility in non-linear curve fitting (see Bollen & Curran, Reference Bollen and Curran2006; Grimm et al., Reference Grimm, Ram and Estabrook2017), and our approach to stacking growth models can be easily extended to allow the analyst to focus on any shape component (linear or non-linear) of substantive interest.

An important flexibility in growth curve modeling allows for the estimation of non-linear trajectories using latent basis methods. This specification requires that some of the time points be fixed to constant values while allowing the remaining time points to be estimated from the data directly. Latent basis modeling yields data-based estimates of the time points and often provides better fit of the model to the empirical trajectories of change over time compared to forcing, say, a linear trend on the data. In the context of latent basis methods, the pace parameter $\pi _{1i}$ is best conceived of as a shape parameter, but we will continue to refer to this as a pace of progress parameter. These, and other extensions, are discussed in Bollen & Curran (Reference Bollen and Curran2006) and Grimm et al. (Reference Grimm, Ram and Estabrook2017). For this article, we will examine differences in predictive performance across three latent curve models: (a) a simple linear growth curve model, (b) a latent basis model with that last time point (2022) freed to reflect a decrease in scores due to the COVID-19 pandemic, and (c) a latent basis model in which the basis functions from 2015 to 2022 are freed, reflecting the idea that the trend downward had started prior to 2022.

In addition to adding time-varying predictors to the model as in Equation (2), it is common practice to model the starting point and pace parameters as a function of time-invariant predictors. For this article, the inter-country model for the starting point $\pi _{0i}$ and pace $\pi _{1i}$ , respectively, can be written as

(4) $$ \begin{align} \pi_{0i} = \beta_{00} + \sum_{q = 1}^{Q}\beta_{0q}x_{qi} + \epsilon_{\pi_{0i}}, \end{align} $$

and

(5) $$ \begin{align} \pi_{1i} = \beta_{10} + \sum_{q = 1}^{Q}\beta_{1q}x_{qi} + \epsilon_{\pi_{1i}}, \end{align} $$

where $x_{qi}$ are values on Q predictors for country i, $\beta _{00}$ , $\beta _{01}$ , $\beta _{10}$ , and $\beta _{1q}$ are the intercept and slopes associated with the time-invariant predictors of the starting point and pace of progress, and $\epsilon _{\pi _{0i}}$ and $\epsilon _{\pi _{1i}}$ are errors. An example of a time-invariant predictor might be the type of political system of a country.Footnote ² For this article, we will develop a set of distinct models for predicting $\pi _{0}$ and $\pi _{1}$ .Footnote ³

The Bayesian framework for the growth curve model in Equations (1) through (5) requires specifying a probability model for the outcome and placing priors on all model parameters (Kaplan, Reference Kaplan2023). Following Kaplan & Huang (Reference Kaplan and Huang2021), the priors for our growth curve models will be non-informative or weakly-informative (see, e.g., Gelman et al., Reference Gelman, Simpson and Betancourt2017). This choice allows the data to “speak” while stabilizing the analysis without impacting inferences, particularly in the presence of small sample sizes (Kaplan, Reference Kaplan2023).

4 Bayesian probabilistic prediction

In the previous section, we discussed the specification, identification, and estimation of the growth curve modeling framework that we will use in our example. Again, our goal is to optimize the prediction of the pace at which countries are trending toward or away from the SDGs for the purposes of accurately predicting future outcomes of country-level mathematics competencies, and thus, the focus of this article is on prediction in the longitudinal context. We argue that a central characteristic of statistics is to develop accurate predictive models, and, all other things being equal, a given model is to be preferred over other competing models if it provides better predictions of what actually occurred (Dawid, Reference Dawid1984). Indeed, it is hard to feel confident about inferences drawn from a model that does a poor job of predicting the extant data. For our problem, the question is how to develop accurate predictive models of country-level pace, and, importantly, how to evaluate the accuracy of the predictions. Only then may we feel comfortable using optimal predictions of pace to predict future observable outcomes. We argue that the evaluation of Bayesian predictive models is best situated in the context of Bayesian decision theory.

Bayesian decision theory (see e.g. Berger, Reference Berger2013; Good, Reference Good1952; Lindley, Reference Lindley1991) provides a natural and intuitive approach to evaluating Bayesian predictive models. Specifically, as will be expanded on below, Bayesian decision theory casts the problem of predictive evaluation in the context of minimizing expected loss—that is, the penalty that is accrued from using a particular model to predict future observations. The less the expected loss, the better the model is at predictive performance in comparison to other models.

4.2 Loss functions for evaluating predictions

The goal of predictive modeling is to minimize the loss associated with taking an action a among a set of actions in the action space $\mathcal {A}$ . A number of loss functions exist, but common loss functions rest on the negative quadratic loss function

(6) $$ \begin{align} L(a,\tilde{y} ) = (\tilde{y} - a)^2. \end{align} $$

The optimal action $a^*$ is the one that minimizes the posterior expected loss, written as

(7) $$ \begin{align} a^* = \operatorname*{\mbox{arg min}}_{a \in \mathcal{A}} \int_{\Omega}L(\omega, a)p(\omega|D)d\omega. \end{align} $$

The idea here is to take an action a that minimizes the loss L when the future observation is $\tilde {y}$ . Clyde & Iversen (Reference Clyde and Iversen2013) show that the optimal decision obtains when $a^* = E(\tilde {y}|D)$ , which is the posterior predictive mean of $\tilde {y}$ given the data D. Under the assumption that the true model exists and is among the set of models under consideration, this can be expressed as

(8) $$ \begin{align} E(\tilde{y}|D) = \sum_{k=1}^{K}E(\tilde{y}|M_k,D)p(M_k|D)=\sum_{k=1}^{K}p(M_k|D)\hat{\tilde{y}}_{M_k}, \end{align} $$

where $\hat {\tilde {y}}_{M_k}$ is the posterior predictive mean of $\tilde {y}$ under $M_k$ and $p(M_k|D)$ is the posterior model probability (PMP) associated with model k. The PMP can be expressed as

(9) $$ \begin{align} p(M_k|D) = \frac{p(D|M_k)p(M_k)}{\sum_{l=1}^Kp(D|M_l)p(M_l)}, \qquad{l \ne k}. \end{align} $$

where $p(M_k)$ the prior probability for model k. The typical default uniform prior mass over the model space and non-informative priors for the parameters of each model, but software programs such as BMS (Zeugner & Feldkircher, Reference Zeugner and Feldkircher2015) allows for other prior choices over both the model space and the parameters. Equations (8) and (9) define BMA.

It is important to note that when considering the selection of a single model, one might be tempted to choose the model with the highest PMP. In the case of only two models, the model with the largest PMP will be the closest to the BMA solution. However, for more than two models, Clyde & Iversen (Reference Clyde and Iversen2013) point out that the model closest to the BMA solution might not be the one with the largest PMP.

4.3 Scoring rules for probabilistic prediction

A critical part of building ensemble prediction models is to have a method for assessing the quality of an ensemble’s predictive performance, sometimes referred to as a model’s predictive skill. Popular methods used in economic forecasting and weather forecasting, among other areas, for assessing predictive skill are referred to as scoring rules (see, e.g., Bernardo & Smith, Reference Bernardo and Smith2000; Gneiting & Raftery, Reference Gneiting and Raftery2007; Jose et al., Reference Jose, Nau and Winkler2008; Merkle & Steyvers, Reference Merkle and Steyvers2013; Winkler, Reference Winkler1996). Scoring rules provide a measure of the accuracy of probabilistic predictions, and a prediction can be said to be “well-calibrated” if the assigned probability of the outcome matches the actual proportion of times that the outcome occurred (Dawid, Reference Dawid1982). For this article, we focus on one strictly proper scoring rule that is commonly used to evaluate predictive accuracy—namely, the Kullback–Leibler Divergence score (Kullback & Leibler, Reference Kullback and Leibler1951; Kullback, Reference Kullback1959, Reference Kullback1987).

4.4 Kullback–Leibler divergence score

Consider two distributions, $p(y)$ and $g(y|\theta )$ , where $p(y)$ could be the distribution of observed mathematics proficiency scores, and $g(y|\theta )$ could be the prediction of these mathematics scores based on a model. The KLD between these two distributions can be written as

(10) $$ \begin{align} \mbox{KLD}(f,g) = \int p(y)\mbox{log} \left(\frac{p(y)}{g(y|\theta)}\right)dy, \end{align} $$

where $\mbox {KLD}(f,g)$ is the information lost when $g(y|\theta )$ is used to approximate $p(y)$ . For example, the actual mathematics outcome scores might be compared to the predicted outcome using BMA along with different choices of model and parameter priors. The model with the lowest KLD measure is deemed best in the sense that the information lost when approximating the actual mathematics outcome distribution with the distribution predicted on the basis of the model is lowest.

5 The tripartite $\mathcal {M}$ -framework

We noted earlier that BMA rests on a very restrictive assumption, namely that there is a true (or correct) model for predicting the pace of progress, denoted as $M_T$ , and that this true model for the pace of progress $\pi _1$ is in the set of models that is being averaged. If this assumption does not hold, then conventional BMA does not makes sense because the priors on the model space are elicited to reflect the analyst’s belief about the existence of the true model within the full set of models under consideration. Nevertheless, if we assume that $M_T$ is in the space of models under consideration, this is referred to as the $\mathcal {M}\text{-}closed$ framework, introduced as one of three modeling frameworks ( $\mathcal {M}$ -frameworks) by Bernardo & Smith (Reference Bernardo and Smith2000) and further discussed in Clyde & Iversen (Reference Clyde and Iversen2013).

The $\mathcal {M}\text{-}{closed}$ framework for BMA may be especially difficult to justify in the social and behavioral sciences. However, as pointed out by Bernardo & Smith (Reference Bernardo and Smith2000), there may be cases in which it is reasonable to act as though there is a true model. For example, we may wish to act as though $\mathcal {M}\text{-}{closed}$ holds when a model has demonstrated good predictive capabilities under a wide variety of situations, but that under a new situation, new uncertainties arise. Such justification might be reasonable in cross-sectional studies, but may be particularly difficult in longitudinal studies where information from previous longitudinal studies may be hard to come by. Still, as long as the analyst is comfortable assigning model priors, then the $\mathcal {M}\text{-}{closed}$ framework can be adopted. Nevertheless, the truth or falsity of the $\mathcal {M}\text{-}{closed}$ framework notwithstanding, it is important to reiterate that conventional BMA takes place under the $\mathcal {M}\text{-}{closed}$ framework and, indeed, readily available BMA software typically employ a non-informative prior to the space of models as a default, with the idea that the true model lies in the model space.

6 Ensemble prediction using Bayesian stacking

The method of stacking was originally developed in the machine learning literature by Wolpert (Reference Wolpert1992) and Breiman (Reference Breiman1996) and brought into the Bayesian paradigm by Clyde & Iversen (Reference Clyde and Iversen2013). A review of Bayesian stacking applied to large-scale educational assessments can be found in Kaplan (Reference Kaplan2021) and extensions of Bayesian stacking applied to multilevel models can be found in Huang & Kaplan (Reference Huang and Kaplan2024). The basic idea behind stacking is to enumerate a set of $K (k=1,2,\ldots K)$ models and then create a weighted combination of their predictions.

In what follows, we describe the process of ensemble modeling via Bayesian stacking for the pace of progress parameter, $\pi _{1}$ noting that the same approach was used for the starting point, $\pi _{0}$ . Returning to our example, we can specify a set of ensemble member models of the pace of progress in mathematics proficiency as

(11) $$ \begin{align} \pi_{1} = f_{k}(x) + \epsilon, \end{align} $$

were $f_{k}(x)$ are different models for the pace of progress $\pi _{1}$ conditional on a vector of predictors of the pace—e.g., some models may include only demographic predictors, others may include various combinations of attitudes and behaviors related to mathematics, and still others may be highly complex functional forms for the prediction of rates of change. To begin, define a set of weights as a simplex,

(12) $$ \begin{align} \mathcal{W}_1^K = \Biggl\{w \in [0,1]^K : \sum_{k=1}^K w_k =1\Biggr\}. \end{align} $$

The stacking problem can be written in terms of either minimizing the divergence d as

(13) $$ \begin{align} \min_{w\in \mathcal{W}_1^K} d \Biggl(\sum_{k=1}^K w_kp(\tilde\pi_1|\pi_1,M_k),p_t(\tilde\pi_1|\pi_1)\Biggr), \end{align} $$

or maximizing the log score

(14) $$ \begin{align} \max_{w\in \mathcal{W}_1^K} S \Biggl(\sum_{k=1}^K w_kp(\tilde\pi_1|\pi_1,M_k),p_t(\tilde\pi_1|\pi_1)\biggr). \end{align} $$

To approximate the full predictive distribution $p(\tilde \pi _1|\pi _1,M_k)$ Yao et al. (Reference Yao, Vehtari, Simpson and Gelman2018) use the leave-one-out (LOO) predictive distribution

(15) $$ \begin{align}{} \hat{p}_{k,-1}(\pi_{1i}) = \int p(\pi_{1i}|\theta_k, M_k)p(\theta_k|\pi_{1,-i},M_k)d\theta_k. \end{align} $$

The stacking weights using the log score are the solution to

(16) $$ \begin{align} \max_{w\in\mathcal{W}_1^K}\frac{1}{n}\sum_{i=1}^n \text{log}\sum_{k=1}^K w_k\hat{p}(\pi_{1i}|\pi_{1,-i},M_k). \end{align} $$

Optimization of Equation 16 uses algorithms in Stan (Stan Development Team, 2021). Inspection of Equations (15) and (16) reveal that the the LOO predictive density is being used twice and so could leave to overally optimistic conclusions. To remedy this, a full Bayesian approach referred to as Bayesian hierarchical stacking (BHS) could be used, but was not implemented in this article (see Huang & Kaplan, Reference Huang and Kaplan2024, for an application of BHS to large-scale assessments).

6.1 LOO cross-validation

We see from Equation (16) that a method is needed to estimate $\pi _{1i}$ based on $n-1$ observations leaving the ith country out, and the most common approach is referred to as leave-one-out cross validation. LOO-cross-validation (LOO-CV) is a special case of q-fold cross-validation (q-fold CV) when $q=n$ . In q-fold CV, a sample is split into q groups (folds) and each fold is taken to be the validation set with the remaining $q-1$ folds serving as the training set. For LOO-CV, each observation serves as the validation set with the remaining $n-1$ observations serving as the training set. Leave-one-out cross-validation is available in the R software program loo (Vehtari et al., Reference Vehtari, Gabry, Yao and Gelman2019).Footnote ⁵

Following Vehtari et al. (Reference Vehtari, Gelman and Gabry2017), let $\pi _{1i} (i=1,\ldots ,n)$ be an n-dimensional vector of pace parameters following a distribution conditional on parameters $\theta $ - viz. $p(\pi _{1}|\theta ) = \prod _{i=1}^{n}p(\pi _{1i}|\theta )$ . Given a prior distribution on the parameters, $p(\theta )$ , we can obtain the posterior distribution, $p(\theta |\pi _{1})$ as well as a posterior predictive distribution of predicted values $\tilde {\pi }_{1}$ written as $p(\tilde {\pi }_{1}|{\pi }_{1}) = \int p(\tilde {\pi }_{1}|\theta )p(\theta |{\pi }_{1})d\theta $ . The Bayesian LOO-CV rests on the derivation of the expected log point-wise predictive ity (ELPD) for new data defined as

(17) $$ \begin{align} \mbox{ELPD} = \sum_{i=1}^{n}\int p_t(\tilde{\pi}_{1i})\mbox{log } p(\tilde{\pi}_{1i}|{\pi}_{1})d\tilde{\pi}_{1i}, \end{align} $$

where $p_t(\tilde {\pi }_{1i})$ represents the distribution of the true but unknown data-generating process for each country’s pace of progress $\tilde {\pi }_{1i}$ and where Equation (17) is approximated by cross-validation procedures. The ELPD provides a measure of predictive accuracy for the n data points taken one at a time (Vehtari et al., Reference Vehtari, Gelman and Gabry2017). From here, the Bayesian LOO estimate can be written as

(18) $$ \begin{align} \mbox{ELPD}_{loo} = \sum_{i=1}^{n}\mbox{log}\,p({\pi}_{1i}|{\pi}_{1-i}), \end{align} $$

where

(19) $$ \begin{align} p({\pi}_{1i}|{\pi}_{1-i}) = \int\,p({\pi}_{1}|\theta)p(\theta|{\pi}_{1-i})d\theta, \end{align} $$

which is the LOO predictive distribution using the log predictive score to assess predictive accuracy. It is useful to note that an information criterion based on LOO (LOO-IC) can be easily derived as

(20) $$ \begin{align} \mbox{LOO-IC} = -2\,{\widehat{\mbox{ELPD}}_{loo}}. \end{align} $$

which places the LOO-IC on the “deviance scale” (see Vehtari et al., Reference Vehtari, Gelman and Gabry2017 for more details on the implementation of the LOO-IC in loo). Among a set of competing models, the one with the smallest LOO-IC is considered best from an out-of-sample point-wise predictive point of view.

As pointed out by Vehtari et al. (Reference Vehtari, Gelman and Gabry2017), it can be time-consuming to calculate exact LOO-CV and this may be a reason why LOO-CV is not widely adopted. To remedy this, Vehtari et al. (Reference Vehtari, Gelman and Gabry2017) developed a fast and stable approach to obtaining LOO-CV referred to as Pareto-smoothed importance sampling (PSIS-LOO) (see Vehtari et al., Reference Vehtari, Gelman and Gabry2017, for more details). The PSIS approach is implemented in loo (Vehtari et al., Reference Vehtari, Gabry, Yao and Gelman2019).

6.2 Other types of stacking weights

In addition to stacking weights based on the ELPD, we will also examine the performance of two alternative stacking weights: pseudo-BMA (PBMA) and pseudo-BMA+ weights.

6.3 Pseudo-BMA weights

Pseudo-BMA (PBMA) weights were proposed by (Geisser & Eddy, Reference Geisser and Eddy1979, see also; Gelfand, Reference Gelfand, Gilks, Richardson and Spiegelhalter1996; Yao et al., Reference Yao, Vehtari, Simpson and Gelman2018). The basic idea behind PBMA is as follows. First, as discussed in Yao et al. (Reference Yao, Pirš, Vehtari and Gelman2021)), LOO-CV has connections to other types of weights that can be used for stacking. For example, in the case of maximum likelihood estimation, LOO-CV weights are asymptotically equivalent to Akaike information criterion (AIC) weights (Akaike, Reference Akaike, Petrov and Csaki1973) that are used in frequentist model averaging applications (see also; Burnham & Anderson, Reference Burnham and Anderson2002; Fletcher, Reference Fletcher2018; Yao et al., Reference Yao, Vehtari, Simpson and Gelman2018). As a method of model selection, earlier work by Geisser & Eddy (Reference Geisser and Eddy1979), see also; Gelfand, Reference Gelfand, Gilks, Richardson and Spiegelhalter1996) criticised the underpinnings of Bayes factors and suggested substituting the marginal likelihood of the kth model, $p(y|M_k)$ , used in the calculation of Bayes factors with Bayesian leave-one-out cross-validation predictive densities, defined as $\prod _{i=1}^n p(\pi _{1i}|\pi _{1-i},M_k)$ . Yao et al. (Reference Yao, Vehtari, Simpson and Gelman2018) refer to AIC weighting using LOO-CV predictive densities as pseudo-BMA weighting.

6.4 Pseudo-BMA+ weights

The difficulty with PBMA weights is that they do not take into account uncertainty in the LOO estimation of the weights. To address this Yao et al. (Reference Yao, Vehtari, Simpson and Gelman2018) proposed an approach that combines the Bayesian bootstrap (see Rubin, Reference Rubin1981) with the ELPD defined earlier. They refer to this approach as pseudo-BMA+ (PBMA+). Following Yao et al. (Reference Yao, Vehtari, Simpson and Gelman2018), the essential idea behind PBMA+ is that the posterior distribution of the realizations of a random variable Z, that is $z_{i}, (i = 1,\ldots ,n)$ , follows a Dirichlet(1,…,1) distribution—i.e., a uniform distribution. Taking samples from this distribution yields Bayesian bootstrap samples from which parameters from this distribution can be calculated. Specifically, let

(21) $$ \begin{align} \{\delta\}_{i=1}^n \sim \overbrace{\text{Dirichlet}(1,\ldots,1)}^n \end{align} $$

be a set of posterior probabilities for all $z_i$ representing one Bayesian bootstrap replication. From here, a parameter of interest represented as a function of Z, $\phi (Z|\delta )$ can be obtained as

(22) $$ \begin{align} \hat{\phi}(Z|\delta) = \sum_{i=1}^n\delta\phi(z_i). \end{align} $$

Repeated sampling from $\pi $ then results in an estimate of $\phi (Z)$ .

With regard to stacking, Yao et al. (Reference Yao, Vehtari, Simpson and Gelman2018) note that the ELPD based on LOO can be highly skewed and argue that the Bayesian bootstrap might be an improvement over the usual Gaussian approximation. The PBMA+ weighting follows essentially the same line of argument as the conventional Bayesian bootstrap. That is, define for each model k,

(23) $$ \begin{align} \{z\}_{i=1}^k = \bigl\{{\widehat{{\text{ELPD}}}_{loo}}\bigr\}_{k=1}^K. \end{align} $$

Taking B bootstrap samples $(\delta _{1b},\ldots ,\delta _{nb}$ ), $b=1,\ldots ,B$ from $\overbrace {\text {Dirichlet}(1,\ldots ,1)}^n$ allows us to calculate the weighted means as

(24) $$ \begin{align} \bar{z}_{b}^k = \sum_{i=1}^n\delta_{ib}z_{i}^k. \end{align} $$

From here, a Bayesian bootstrap sample of the stacking weight for model k based on bootstrap samples of size B can be obtained as

(25) $$ \begin{align} w_{kb} = \frac{\text{exp}(n\bar{z}_{b}^k)}{\sum_{k=1}^K \text{exp}(n\bar{z}_{b}^k)}, \qquad b=1,\ldots,B \end{align} $$

leading to the final PBMA+ weight for model k

(26) $$ \begin{align} w_k = \frac{1}{B}\sum_{b=1}^B w_{k,b}. \end{align} $$

Of importance to this article, Yao et al. (Reference Yao, Vehtari, Simpson and Gelman2018) showed that PBMA+ performs better than BMA and PBMA in $\mathcal {M}$ -open settings, but not as well as stacking using the log score. This article adds to the existing literature by comparing stacking based on the ELPD $_{loo}$ to PBMA and PBMA+ weights in the context of growth curve models applied to large-scale assessments.

7 Example: Stacking growth curve models of PISA mathematics proficiency

This article will apply Bayesian stacking to data from 53 countries that have participated in PISA from 2009 to 2022. Launched in 2000 by the OECD, PISA is a triennial international survey that aims to evaluate education systems worldwide by testing the skills and knowledge of 15-year-old students and is, arguably, the most important policy-relevant international survey currently operating. In 2022, 690,000 15-year-old students attending educational institutions in lower secondary education grades or higher from 81 countries and economies took an internationally agreed-upon two-hour test (OECD, 2023b). Students were assessed in reading, mathematics, science, collaborative problem solving and financial literacy. Available country-level results already account for the complex sampling design and plausible value methodology for obtaining mathematics literacy scores. A detailed account of the PISA design can be found in the overview by Kaplan & Kuger (Reference Kaplan, Kuger, Kuger, Klieme, Jude and Kaplan2016).

We will focus on country-level longitudinal outcomes in mathematics proficiency using data from PISA 2009 to PISA 2022. Although longer time points are available, it was decided to only include countries with complete data on the mathematics outcome so that predictive analyses were not dependent on imputation of missing data. The list of these countries is shown in Table 1.

Table 1

Fifty-four countries and economies with complete mathematics assessment data from 2009 to 2022

Note: Luxembourg did not have data for PISA 2022 and the Russian Federation was excluded from the PISA 2022 assessment due to the war in Ukraine.

8 Proposed workflow

For this article, we use the R program blavaan, which is a lavaan-type SEM interface to rstan (Merkle et al., Reference Merkle, Fitzsimmons, Uanhoro and Goodrich2021; Rosseel, Reference Rosseel2012; Stan Development Team, 2023) for the estimation of the starting point and pace of progress. The program mice (van Buuren & Groothuis-Oudshoorn, Reference van Buuren and Groothuis-Oudshoorn2011) was used for predictive mean matching imputation of missing data. Weight calculations under EPLD $_{loo}$ , PBMA, and PBMA+ used the program loo (Vehtari et al., Reference Vehtari, Gabry, Yao and Gelman2019), and model evaluation used LaplacesDemon (Statisticat & LLC., 2021). The full analysis with eight models, including growth curve modeling, missing data imputation, weight calculation, model evaluation, and knitted with rmarkdown (Allaire et al., Reference Allaire, Xie, Dervieux, McPherson, Luraschi, Ushey and Iannone2024), took forty-five minutes using twenty cores on a Dell laptop. The code for this analysis is available at https://bise.wceruw.org/publications.html.

The workflow for our example is as follows.

• Step 1: Estimate a latent growth curve model without predictors. Different growth curve models should be specified and compared via LOO-CV.

  • Step 1a: If each model gives rise to very different growth rates, it may be useful to perform Bayesian stacking on these models first. The stacked predictive distributions could then be used below in a stacking that involves models with different predictors. We refer to this as super-ensemble modeling.

• Step 2: Create a list of ensemble members, each one being a substantively distinct model predicting the growth rate. For each member model, estimate the pace of progress from Equation (1) and compare models based on the ELPD $_{loo}$ values, and, relatedly the LOO-IC. Compare ELPD $_{loo}$ weights to PBMA and PBMA+ weights. In practice, the analyst would likely choose just one of these weights, but there is no harm in comparing outcomes using the other weights.

  • 2a: (Optional, but recommended): Here, one could check on the relative distinctiveness of the individual member models via posterior predictive checking of each model separately (see Nold et al., Reference Nold, Meinfelder and Kaplan2024; Yao et al., Reference Yao, Pirš, Vehtari and Gelman2021).Footnote ⁶

• Step 3: Combine the predictive distributions of the growth rates weighted by the stacking weights. The code we use for stacking predictive distributions is available at (blinded for review)

• Step 4: Evaluate the stacked predictive distribution of linear growth against the linear growth rate in Step 1 via the KLD scoring rule in Equation (10). In theory, the stacked predictive distribution should have a smaller KLD than that of any of the model members.

• Step 5: Based on the average pace of progress from the stacked predictive distribution, produce prediction plots.

  • Step 5a (Optional but recommended): For this step, pseudo out-of-sample predictions may be desirable. For pseudo out-of-sample prediction, the stacked pace of progress is estimated on all but the Tth-1 time point. Then, with the stacked pace of progress in hand, one simply estimates the proficiency scores for the Tth time point and compares the predictions to the actual values at time T using KLD. One can also compare the different weight methods described above, and the weighting method that yields the lowest KLD would be chosen for the prediction task at hand. Once the analyst has settled on a weighting method that showed the best performance using the KLD, then the analyst can move on to out-of-sample prediction using the intercept and stacked pace of progress.

9 Summary and discussion

The purpose of this article was to demonstrate an approach to combining predictions of the pace of progress from a set of growth curve models. In line with a long tradition of multi-model inference, we argued that combining predictions from multiple models into an ensemble prediction yields overall better predictive skill than what could be achieved from the selection of a single model.

For this article, we presented a workflow where the first step was to examine predictive skill resulting from relatively minor changes in a growth curve model without adding predictors. We argued that this was similar, though not identical to, single-model ensemble methods found in weather forecasting. In our particular example, we did not find any important differences in the predictive performance of the three changes to the growth models. Had we found important differences, it would have been advisable to use stacking methods to create an ensemble prediction. Not having to create the ensemble, we decided to obtain Bayesian estimates of country-level pace of progress in mathematics proficiency from a latent basis model that allowed for some non-linearity due to the impact of COVID-19.

For the next step of our workflow, we specified a number of different models for the pace of progress based on a collection and merging of data sources containing country-level indicators of education. Six models were specified, and again this phase of our workflow is similar, though not identical to, multi-model ensemble prediction, also found in the weather forecasting domain. As mentioned earlier, we did not judge the predictive quality of the individual models or the ensemble based on models for the starting point. We felt it was more important to stack models for the pace of progress, and it remains a future area of research on how to combine decisions based on related parameters of interest. That said, our results showed that although the boys started off in 2009 with higher mathematics proficiency scores on average across the countries, they have been declining at a noticeably faster rate than girls, which we also observed from the growth rates in Table (2). Of course, gender inequity in the decline in mathematics performance as seen in this article, is unacceptable under any circumstances, but it should also be noted that in real terms, these declines are not very large, resulting, on average, in about two score points every three years. Nevertheless, we argue that the methodologies and workflow provided in this article could be informative to national educational systems as they consider policies to reverse these trends. We hasten to add, however, that although this article provides a promising approach to policy-relevant prediction in longitudinal settings, where obtaining optimal estimates of the pace of progress is the focus, in no way should the results of this specific study be interpreted as informing specific policy decisions.

Although we found Bayesian stacking to provide better predictive performance than any single model in the ensemble, there are open issues with the Bayesian stacking that set the stage for future research. First, the performance of Bayesian stacking is, of course, highly dependent on the set of the member models to be ensembled. We believe that the eight categories of models that we specified are defensible in terms of our review of the extant literature on indicators of global education systems, but naturally other specifications are possible, and we can’t even be certain we have captured functional forms correctly. Nevertheless, the tools available for assessing the quality of predictive models allow for alternative models to be specified and compared in terms of their predictive performance.

Second, even with well developed ensemble members, Bayesian stacking is argued to work well when models are as distinct as possible (Breiman, Reference Breiman1996; Clarke, Reference Clarke2003). However, as pointed out by Yao et al. (Reference Yao, Pirš, Vehtari and Gelman2021), model distinctiveness is an ideal and there is presently not much guidance on how to quantify distinctiveness among models, or any existing knowledge as to how serious a problem this might be for predictive performance. That said, Figure 2 could provide useful information for assessing the capacity of the ensemble member models to capture broad features of the empirical pace of progress. Specifically, Figure 2 showed that the predictive distributions from the individual models did not seem to capture the lower tail of the baseline empirical pace of progress and this could indicate that other model specifications that capture flatter growth rates should be specified and included in the stack. Using a plot such as Figure 2 to assess the extent to which ensemble member distributions cover the range of possible values of the distribution is beyond the scope of this study but should constitute future research.

Finally, a third area for future research in the context of Bayesian stacking concerns missing data. With regard to this article, we originally considered taking the difference in the predictors from 2022 to 2009 as a measure of their stability over time. However, two important problems emerged. First, for many countries, data on important indicators were simply not reported by the countries. Second, for those that were reported, many were not measured at each time point coincident with the cycles of PISA used in this study (2009–2022). As such, we felt that there were too many missing data points on a number of important indicators and too few countries to begin with to feel comfortable using multiple imputation methods. Thus, we only examined predictors measured in 2009 under the assumption that these predictors have been relatively stable over time. An inspection of some basic descriptive statistics suggests evidence of stability in many, but not all of the predictors. For missing data in 2009, one predictive mean matching imputation (van Buuren & Groothuis-Oudshoorn, Reference van Buuren and Groothuis-Oudshoorn2011) was used.

We recognize that it is better to analyze many (e.g., ¿20) imputed data sets. One approach to analyzing multiply imputed data sets in a Bayesian analysis was proposed by Zhou & Reiter (Reference Zhou and Reiter2010) who recommended analyzing each imputed data set separately and then mixing and summarizing the posterior draws. They find this approach to yield less biased parameter estimates than averaging the parameter estimates. Combining the approach suggested by Zhou & Reiter (Reference Zhou and Reiter2010) along with Bayesian stacking was felt to be beyond the scope of this article, but it does offer an interesting area for future research. Nevertheless, we recognize that our results could change if we had made other decisions regarding missing data.

To conclude, we recognize that the domain of weather forecasting has certain advantages compared to monitoring trends in educational outcomes. In particular, the non-linear dynamics of weather, as well as the near-continuous collection of data, provide rich information on which to build and stack predictions from complex models. Indeed, as noted by the World Climate Service, modern weather ensemble forecasts have been known to have between 12 and 51 model members (Dutton, Reference Dutton2021). Moreover, weather forecasting models are flexible enough to handle exogenous shocks to weather systems, such as the 2024 Icelandic volcano eruption. In the case of monitoring trends in educational outcomes, as we have seen, often simple linear models will suffice, and theory as to why trends in proficiency outcomes have been, in our example, declining for many countries, are not well developed. Still, the methods we are proposing are also flexible enough to handle exogenous shocks to the educational system, such as the disruption to schooling caused by the 2019 global pandemic. It would be interesting to examine whether the workflow that we propose in this article would be applicable to other education targets or even targets associated with different SDGs.

An additional contribution of this article was the demonstration of how ILSAs generally, and PISA in particular, can be leveraged to provide information relevant to tracking progress to the educational sustainable development targets. We believe that this contribution is important because PISA, in particular, has been used by organizations such as the World Bank to monitor and forecast the global impact of the pandemic on educational trends (Azevedo et al., Reference Azevedo, Hasan, Goldemberg, Iqbal and Geven2020). These endeavors have contributed greatly to our understanding of the impact of the pandemic on schooling. This article adds to the growing literature on monitoring and forecasting educational trends, as well as exogenous shocks to those trends, by combining models that are designed for estimating rates of change, along with Bayesian approaches to optimizing prediction, which we maintain can provide important insights into country-level progress in education.