From Many to One: Consensus Inference in a MIP

Increasing our knowledge of the behavior of geophysical processes in the past, at the present, and in the future requires a scientific understanding of the processes, of the parameters that govern them and of the observing systems measuring them. That understanding can be quantified through a geophysical model and the mathematical and statistical relationships embedded in it. The presence of uncertainty in the model formulation is generally recognized, although where and how it enters is usually a topic of more or less disagreement (e.g., Knutti et al.,, 2010). Reconciliation and consensus estimation in these settings is a topic of central interest in many geophysical studies of global import, including the Coupled Model Intercomparison Projects (CMIP, Meehl et al.,, 2000), and the Ice Sheet Mass Balance Inter-comparison Exercises (IMBIE, Shepherd et al.,, 2018).

This article considers the Model Intercomparison Project (MIP) approach, commonly used to understand the similarities and differences in competing geophysical models. Here a number of teams produce outputs that estimate the underlying geophysical process in question, but the outputs generally differ due to modeling choices or assumptions. The collection of outputs, which we simply call the ‘MIP’s outputs,’ is sometimes also called a ‘multi-model ensemble’. When analyzing and summarizing a MIP’s outputs, it is common to use the ensemble average (‘one vote per model’; see Flato et al.,, 2014, Section 9.2.2.3), and either the range from the ensemble’s lowermost to its uppermost or the ensemble standard deviation, to put bounds on how the process did, does, or will behave. However, it is recognized that this practice of treating all outputs equally is frequently inappropriate (Tebaldi and Knutti,, 2007), as some outputs reproduce observational data more faithfully than others and should therefore receive more attention (Knutti,, 2010). At the same time, similarities between the methodologies of certain teams suggest that these team’s outputs should not be treated as completely independent and therefore each one should receive less attention (Tebaldi and Knutti,, 2007; Abramowitz et al.,, 2019). The most common paradigm (and the one used in this article) for converting these ideas into a method for combining outputs is to weight the outputs according to given criteria. Under this paradigm, the ensemble is typically summarized through a weighted mean rather than the unweighted mean (i.e., the average), and the uncertainty quantification also takes into account the weights (e.g., Tebaldi et al.,, 2005).

A key question is how to construct the ensemble weights. The methods used to do so depend on the details of the application. It is natural to use any available observational data to help with the weighting and, broadly speaking, applications fall onto a spectrum depending on how much verification data are available. At one end of the spectrum, when the underlying geophysical process is contemporaneous, such as in the prediction of meteorological fields, weighting methods can typically take advantage of large amounts of verification data (e.g., Krishnamurti et al.,, 2000). At the other end of the spectrum, when the geophysical process is extremely difficult to measure directly, little-to-no direct verification data may be available; this situation is the focus of this article.

In the middle of the spectrum, limited verification data are available. The archetypal case of this is the projection of future climate, where verification data are available for a number of years up to the beginning of the projections. For example, see Figure 11.9 on page 981 of the Intergovernmental Panel on Climate Change Fifth Assessment Report (Kirtman et al.,, 2013), which includes 42 projections from the CMIP5 ensemble members. For 1986–2005, the projections are based on historical observed greenhouse gas concentrations. Afterwards, they are based on a pre-specified path of future emissions, called a representative concentration pathway (RCP, in this case RCP 4.5). Verification data are available for the years 1986–2005, and they are used to evaluate the performance of the models. In this and similar settings, there are many reasons why historical performance may not guarantee future accuracy (Reifen and Toumi,, 2009) but, nonetheless, one should not discount this historical information entirely (Knutti,, 2010).

For MIPs like CMIP5, where limited verification data are available, Giorgi and Mearns, (2002) present a framework called the Reliability Ensemble Average (REA), in which heuristic weights are derived that reward ensemble members that are simultaneously (i) in consensus and (ii) in agreement with verification data. Tebaldi et al., (2005) formalize REA in a statistical framework for analyzing model projections of future global temperatures across 22 regions. In their approach, the regions are treated independently, and the individual outputs are modeled as random quantities that center on the true climate state. The variances differ between the outputs, which imply different weights, and the estimated weights balance the REA criteria (i) and (ii). There are several extensions of the Tebaldi et al., (2005) method: Smith et al., (2009) treat the regions jointly and describe a cross-validation technique to validate the model fit, while Tebaldi and Sansó, (2009) allow for joint modeling of both temperature and precipitation.

A key limitation of the REA framework and its associated methods is that they do not account for dependence between outputs (Abramowitz et al.,, 2019). More recent developments have attempted to address this by modifying the weights, typically by down-weighting sub-groups of outputs that are considered to be dependent (Bishop and Abramowitz,, 2013). However, there is disagreement on what it means for certain MIP outputs to be dependent, and approaches to understanding this include the distance between outputs (e.g., Sanderson et al.,, 2017), the historical ‘genealogy’ of climate models (e.g., Knutti et al.,, 2013, 2017), component-wise similarity between models (e.g., Boé,, 2018), dependence between errors (e.g., Bishop and Abramowitz,, 2013; Haughton et al.,, 2015), and probabilistic/statistical dependence (e.g., Chandler,, 2013; Sunyer et al.,, 2014).

We now consider the case that is our main focus, where little to no direct verification data are available. An example of this, described in detail in Section 5, is the estimation of surface sources and sinks (fluxes) of CO₂ in a MIP; see Figure 1, which shows the nine teams’ estimates of global land and global ocean non-fossil-fuel fluxes, as well as the unweighted ensemble average. Fluxes of CO₂ occur over large temporal and geographical scales and are difficult to measure directly. One way their distribution is inferred indirectly is through their consequences on atmospheric concentrations of CO₂, for which observations are available (Ciais et al.,, 2010). However, there are complications in using such indirect data taken anywhere in the atmosphere for verification, because issues such as overfitting and model misspecification mean that an accurate reproduction of the observed concentrations may not imply accurate fluxes (e.g., Houweling et al.,, 2010). Therefore, in this case and in cases like it, one is not applying REA’s criterion (ii), of agreement with flux data. However, according to REA criterion (i), one can infer the consensus between the MIP’s outputs (while taking into consideration model dependence, if needed). To apply the first criterion, we present a framework that enables statistically unbiased prediction and estimation (SUPE) based on an analysis of variance (ANOVA) statistical model, originally used for analyzing agricultural field experiments (Fisher,, 1935; Snedecor,, 1937). The framework promotes models that are in agreement and discounts those that are in disagreement, but it does so selectively in regions and seasons thought to have their own geophysical characteristics. For each of these regions/seasons, it provides a consensus summary and accompanying uncertainty quantification. It can also cater for model interdependencies when generating the summary, provided these have been, or can be, identified from a priori or expert knowledge.

Section 2 presents the ANOVA model in general, by recognizing that different factor combinations (e.g., spatial regions, seasons, geophysical-model configurations) will lead to agreement in differing degrees. We model the teams’ outputs for each factor combination as having different variances and, in Section 3, obtain statistically optimal weights based on these variances that allow optimal SUPE and its uncertainty quantification. In Section 4, we show how to estimate these weights using likelihood methods derived from the statistical distributions in the ANOVA model. Section 5 gives an illustration of our statistical-weighting methodology on the MIP published by Crowell et al., (2019), where nine teams produced outputs of global carbon-dioxide fluxes that were compared and summarized. A discussion of our framework, which we call SUPE-ANOVA, and conclusions, are given in Section 6.

The area of statistics known as Experimental Design was born out of a need to control variations in agricultural field experiments, in order to compare and then select crop treatments based on replicate observations. The statistical ANOVA model was devised to analyze these experiments, and much has been written about it, starting with Fisher, (1935), Snedecor, (1937), Snedecor and Cochran in its eight editions (1957, …, 2014), and many other texts. It has been adopted in a number of disciplines driven by experimentation, including industrial manufacturing (e.g., Wu and Hamada,, 2000), and computer experiments where the so-called treatments are complex algorithms and code (e.g., Santner et al.,, 2010). We propose using the ANOVA framework for the analysis of outputs from important geophysical experiments performed within a MIP. It enables the output from many to be weighted in a statistically efficient manner to obtain one consensus estimate and a quantification of its uncertainty.

Let ‘$f$’ represent a combination of factors (e.g., season, region) that the underlying geophysical process is expected to depend on, and let ‘$j$’ index the output of the various teams participating in the MIP. By stratifying on each factor combination, we can write the output for the $f$th factor combination ($f=1,\ldots,F$) and the $j$th team ($j=1,\ldots,J$) as

where $I(f)$ represents the number of replicates for factor combination $f$. Each combination $(f,i)$ corresponds to a particular region/time-period in the MIP, and note that $I(f)$ may differ for different factor combinations, reflecting the ability to capture complex MIP designs. In the analysis of the CO₂-flux MIP given in Section 5, $I(f)$ is a constant number (equal to six) for all regions and seasons.

The notation recognizes that each output is split into $F$ blocks (factor combinations) and, within the $f$th block, there are $(I(f)\cdot J)$ outputs that are observing the same part of the geophysical process, corresponding to the $I(f)$ replicates and the $J$ teams participating in the MIP. For example, suppose the MIP’s goal is to compare $20$ climate models’ projections of North American mean temperature between the years $2021$ and $2070$, inclusive, at yearly intervals. Say the factor $f$ corresponds to decade, so $f\in\{1,2,\ldots,5\}$. Further, $i\in\{1,2,\ldots,10\}$ corresponds to the $I(f)=10$ years within a decade, and $j\in\{1,2,\ldots,20\}$ corresponds to the $20$ climate models. For climate model 17 (say) and decade 4 (say, which corresponds to the years between $2051$ to $2060$, inclusive), the output at yearly intervals within the decade is $\{Y_{4,i}^{(17)}:i=1,\ldots,10\}$.

The statistical ANOVA model proposed for (1) is, for $f=1,\ldots,F$, $i=1,\ldots,I(f)$, and $j=1,\ldots,J$,

where $\{\eta_{f,i}^{(1)},\ldots,\eta_{f,i}^{(J)}\}$ each have mean zero, $\mathrm{var}(\eta_{f,i}^{(j)})=(\sigma_{f}^{(j)})^{2}$, and $\mathrm{cov}(\eta_{f,i}^{(j)},\eta_{f,i}^{(j^{\prime})})=\rho_{f}^{(j,j^{\prime})}\sigma_{f}^{(j)}\sigma_{f}^{(j^{\prime})}$. Here, $Y_{f,i}$ is the MIP’s consensus of the geophysical process that each of the MIP outputs is attempting to estimate, and $\eta_{f,i}^{(j)}$ is the deviation from the consensus of the $j$th output, assumed to have team-specific variances and covariances and to be independent of $Y_{f,i}$. Marginally, $\eta_{f,i}^{(j)}\sim\mathrm{Dist}(0,(\sigma_{f}^{(j)})^{2})$, where $\mathrm{Dist}(\mu,\sigma^{2})$ denotes a generic distribution with mean $\mu$ and variance $\sigma^{2}$. Dependence between teams’ outputs is captured with correlations $\{\rho_{f}^{(j,j^{\prime})}:j,j^{\prime}=1,\ldots,J\}$. The statistical ANOVA model has variances that may differ from team to team. A mathematical-statistical result called the Gauss–Markov theorem (e.g., Casella and Berger,, 2006) implies that the hereto-referred weights are inversely proportional to these variances, and there is a more general version of the theorem that includes the dependencies between outputs; see Appendix B.

For $f=1,\ldots,F$ and $i=1,\ldots,I(f)$, we write the consensus as,

where $\{\alpha_{f,1},\ldots,\alpha_{f,I(f)}\}$ each have mean zero. Here, $\mu_{f}$ is the unknown mean of the $I(f)$ replicates, which we refer to as the climatological mean, and $\{\alpha_{f,1},\ldots,\alpha_{f,I(f)}\}$ are random effects that capture variability in the replicates, and whose distribution we denote as $\mathrm{Dist}(0,\tau^{2}_{f})$. In the illustrative example of North American temperatures, $f$ indexes decade, $\mu_{f}$ corresponds to the consensus temperature in the $f$th decade, while $\alpha_{f,i}$ is the deviation from $\mu_{f}$ for the $i$th year within the decade. The variability of temperatures can differ between decades, reflected by $\tau^{2}_{f}$ for $f=1,\ldots,5$. A temporal trend within the $f$th decade could also be included in the model for $\{\alpha_{f,i}:i=1,\ldots,I(f)\}$. The team-specific variances $\{(\sigma_{f}^{(j)})^{2}:j=1,\ldots,J\}$ and correlations $\{\rho_{f}^{(j,j^{\prime})}:j,j^{\prime}=1,\ldots,J\}$ defined through (2) may also be different for different decades. In the North American temperature example, a team’s output with higher weight in one decade may not have higher weight in another.

The statistical ANOVA model in (2) and (3) allows statistically unbiased prediction and estimation (SUPE) on the consensus for the climatological means, $\{\mu_{f}:f=1,\ldots,F\}$, and on the consensus for the underlying geophysical process, $\{Y_{f,i}:i=1,\ldots,I(f);\enskip f=1,\ldots,F\}$. Notice the heterogeneity of variation in (2) and (3), where the variances (and correlations) depend on both $f$ and $j$. In the next section, we explain how knowledge of these variance (and correlation) parameters leads to statistically optimal weighting of the $J$ MIP outputs (1) for each $f$ when making inference on $\{\mu_{f}\}$ and $\{Y_{f,i}\}$.

In the simple example from Section 2, of combining 20 climate models, not all models are expected to show the same variability. Different assumptions are made, say, about the model’s base resolution at which energy and water exchange takes place, how nonlinear dynamics are approximated, and so forth. This possibility is taken into account through (2), where the outputs from each team receives a different variance. A consequence of this is that, for optimal estimation and prediction of the consensus, different weights should be assigned to the different teams. This section discusses how to construct the optimal weights.

To illustrate optimal weighting in a simpler setting, without factor combinations, consider each climate model’s 2070 projected North American mean temperature (Tmp), $Y_{\mathrm{Tmp}}^{(1)},\ldots,Y_{\mathrm{Tmp}}^{(20)}$. In this case, there is only one replication, so we omit the subscript $i$ on all quantities. To start with, assume a naïve statistical model wherein the 20 climate models all have the same variability:

where independently, $\eta_{\mathrm{Tmp}}^{(j)}\sim\mathrm{Dist}(0,\sigma^{2}_{\mathrm{Tmp}})$. For simplicity of exposition, the correlations in this example are assumed to be zero. From (4), $Y_{\mathrm{Tmp}}^{(j)}\sim\mathrm{Dist}(\mu_{\mathrm{Tmp}},\sigma^{2}_{\mathrm{Tmp}})$ for $j=1,\ldots,20$. Then the best linear unbiased estimator (BLUE) of $\mu_{\mathrm{Tmp}}$ is the unweighted mean, and its variance is straightforward to calculate. That is,

To justify using the word ‘best’ in the acronym ‘BLUE,’ one needs to build a statistical model. Assuming the model (4), the estimate (5) is the most precise linear estimator, in that $\mathrm{var}(\bar{\mu}_{\mathrm{Tmp}})\leq\mathrm{var}(\tilde{\mu}_{\mathrm{Tmp}})$, for any other linear unbiased estimator, $\tilde{\mu}_{\mathrm{Tmp}}$, of $\mu_{\mathrm{Tmp}}$. However, (4) does not allow for different variabilities in the different MIP outputs.

Now, for the sake of argument, suppose that the first 15 out of the 20 climate models have an underlying resolution of $0.5\degree\times 0.5\degree$, while the remaining five models have a resolution of $1\degree\times 1\degree$. Consequently, an assumption of equal variance in the statistical distributions in (4) will need modification. Assume variances $(\sigma_{\mathrm{Tmp}}^{(\text{lo})})^{2}$ and $(\sigma_{\mathrm{Tmp}}^{(\text{hi})})^{2}$ are known, and

In Section 4, we show how these variances can be estimated. The parameter $\mu_{\mathrm{Tmp}}$ is still the target of inference but, under the statistical model (6), the Gauss–Markov theorem (e.g., Casella and Berger,, 2006) can be used to establish that the estimate $\bar{\mu}_{\mathrm{Tmp}}$ in (5) is no longer the BLUE. We define the weights for estimating the projected mean temperature in the year 2070, $\mu_{\mathrm{Tmp}}$, as follows:

Under the statistical model (6), the Gauss–Markov theorem can be used to establish that the BLUE (the most precise linear unbiased estimator) and its variance are, respectively,

That is, the statistically optimal weights for estimation of $\mu_{\mathrm{Tmp}}$ are given by (7), and the best estimator of $\mu_{\mathrm{Tmp}}$ is obtained by weighting each datum inversely proportional to its variance.

A more general statistical model, analogous to what we consider in (3), involves the random effect, $\alpha_{\mathrm{Tmp}}$, as follows:

where $\alpha_{\mathrm{Tmp}}\sim\mathrm{Dist}(0,\tau_{\mathrm{Tmp}}^{2})$, and $\tau_{\mathrm{Tmp}}^{2}$ is assumed known. When $\alpha_{\mathrm{Tmp}}=0$, we recover the model (6), however the random effect $\alpha_{\mathrm{Tmp}}$ plays an important role in recognizing the stochasticity in the projected temperature. Hence, for $j=1,\ldots,20$,

Notice that there are between-model correlation terms induced by the random effect; for $j\neq j^{\prime}$,

where recall that in this simple example $\eta_{\mathrm{Tmp}}^{(j)}$ and $\eta_{\mathrm{Tmp}}^{(j^{\prime})}$ are uncorrelated. Prediction of $Y_{\mathrm{Tmp}}$, which is the appropriate inference in the presence of the random effect, is made through the predictive distribution of $Y_{\mathrm{Tmp}}$ given all the data from all the climate-model outputs; the optimal predictor of $Y_{\mathrm{Tmp}}$ is the mean of this predictive distribution. We use the best linear unbiased predictor (BLUP) to make inference on $Y_{\mathrm{Tmp}}$. In the simple example of North American temperatures, the BLUP is

where

Notice that, for both optimal prediction (BLUP) and optimal estimation (BLUE), outputs with larger variability around the consensus receive lower weights. Additionally, for the BLUP given by (11), the weight on $\mu_{\mathrm{Tmp}}$ depends on how variable the underlying process is; that is, it depends on $\tau_{\mathrm{Tmp}}^{2}$. Thus, the more variable the geophysical process is (captured through the parameter $\tau^{2}_{\mathrm{Tmp}}$), the less weight its mean $\mu_{\mathrm{Tmp}}$ receives in the optimal predictor, $Y_{\mathrm{Tmp}}^{*}$, of $Y_{\mathrm{Tmp}}$.

To justify the word ‘best’ in the acronym ‘BLUP,’ an uncertainty measure called the mean-squared prediction error (MSPE) is used. Let $\tilde{Y}_{\mathrm{Tmp}}$ be any linear predictor of $Y_{\mathrm{Tmp}}$; then the BLUP, $Y_{\mathrm{Tmp}}^{*}$, minimizes the MSPE since

A general formula for the MSPE of the BLUP is given later in this section.

Returning to the generic MIP, the steps to carry out optimally weighted linear inference follow the same line of development as in the temperature example: (1) specify the first-moment and second-moment parameters in the ANOVA model; (2) derive the variances and covariances of all random quantities in terms of the parameters; (3) weight the outputs (essentially) inversely proportional to their variances to obtain a BLUE/BLUP (i.e., to obtain a best SUPE); (4) quantify uncertainty (variance of the BLUE and MSPE of the BLUP). These four steps require an initialization that might be called a zeroth step, namely the estimation of the parameters $(\sigma_{f}^{(1)})^{2},\ldots,(\sigma_{f}^{(J)})^{2}$ and $\tau^{2}_{f}$; see Section 4. The correlation parameters, $\{\rho_{f}^{(j,j^{\prime})}:j,j^{\prime}=1,\ldots,J\}$, could be known a priori or estimated using a likelihood-based methodology along with the variance parameters; see Appendix B. All together, these steps make up an inference framework for MIPs that we call SUPE-ANOVA.

We now give the optimal estimation and prediction formulas for the generic statistical ANOVA model under the assumption that the correlations between the models, $\{\rho_{f}^{(j,j^{\prime})}\}$, are zero; the general case of non-zero correlations involves matrix algebra and is given in Appendix B. From (2) and (3), it is straightforward to show that $E(Y_{f,i}^{(j)})=\mu_{f}$, $\mathrm{var}(Y_{f,i}^{(j)})=\tau^{2}_{f}+(\sigma_{f}^{(j)})^{2}$, and $\mathrm{cov}(Y_{f,i}^{(j)},Y_{f,i}^{(j^{\prime})})=\tau^{2}_{f}$. From a multivariate version of the Gauss-Markov theorem (e.g., Rencher and Christensen,, 2012), the BLUE of $\mu_{f}$, and its variance, are given by

for $f=1,\ldots,F$, where the optimal weights are

One-sigma and two-sigma confidence intervals for the unknown consensus mean $\mu_{f}$, are given by $\mu_{f}^{*}\pm S_{f}$ and $\mu_{f}^{*}\pm 2S_{f}$, respectively.

From (3), we also want to make inference on $Y_{f,i}$. Recall that $Y_{f,i}=\mu_{f}+\alpha_{f,i}$, where $\{\alpha_{f,i}\}$ are random quantities that capture the departures from the longer-term climatology given by $\mu_{f}$, and they are assumed to have distribution $\mathrm{Dist}(0,\tau^{2}_{f})$. Since the goal of the MIP is to infer a consensus for the underlying process $\{Y_{f,i}\}$, BLUPs of these quantities are sought, along with their MSPEs. Straightforward use of the general Gauss–Markov Theorem (e.g., Rencher and Christensen,, 2012) results in the BLUP $Y_{f,i}^{*}$ and its MSPE, as follows:

where (for the moment) all parameters, including $\mu_{f}$, are assumed known, and

One-sigma and two-sigma prediction intervals for $Y_{f,i}$ are given by $Y_{f,i}^{*}\pm\sqrt{\mathrm{MSPE}(Y_{f,i}^{*})}$ and $Y_{f,i}^{*}\pm 2\sqrt{\mathrm{MSPE}(Y_{f,i}^{*})}$, respectively.

In the next section, we give likelihood-based estimates of the model parameters, which we notate as $\hat{\mu}_{f}$, $\{(\hat{\sigma}_{f}^{(j)})^{2}:j=1,\ldots,J\}$, and $\hat{\tau}_{f}^{2}$. Then the empirical BLUE/BLUP (written as EBLUE/EBLUP) equations are just those in (13)–(18) with the model parameters, $\mu_{f}$, $\{(\sigma_{f}^{(j)})^{2}:j=1,\ldots,J\}$, and $\tau_{f}^{2}$, replaced by their estimated quantities. This yields the EBLUE for $\mu_{f}$, which we denote by $\hat{\mu}_{f}$, and an estimated variance denoted by $\hat{S}_{f}^{2}$. Further, it yields an EBLUP for $Y_{f,i}$, which we denote by $\hat{Y}_{f,i}$, and an estimated mean-squared prediction error, denoted by $\widehat{\mathrm{MSPE}}(\hat{Y}_{f,i})$. These formulas complete the specification of SUPE-ANOVA.

The likelihood-based estimates given in the next section require specification of statistical distributions and, to do this, we assume $\mathrm{Dist}(\cdot,\cdot)$ is $\mathrm{Gau}(\cdot,\cdot)$ for all distributions. When we assume Gaussianity, it should be noted that (16) is best among all possible predictors, not just linear ones (e.g., Cressie and Wikle,, 2011, Sec. 4.1.2).

Estimation of the parameters of the ANOVA model given by (2) and (3) is performed under the assumption that all unknown distributions are Gaussian, which yields a likelihood function of the model’s parameters. The correlation parameters $\{\rho_{f}^{(j,j^{\prime})}\}$ can be modeled and estimated, but they are assumed to be zero in the application in Section 5; here we omit them from the estimation equations, but see Appendix B where they are included. For the remaining parameters, it is known that the maximum likelihood estimates of variance parameters are biased and, to avoid this, we use restricted maximum likelihood (REML) estimation (Patterson and Thompson,, 1971; Harville,, 1977). REML estimation maximizes the restricted likelihood, denoted by $L^{(r)}(\bm{\tau}^{2},\bm{\sigma}^{2}\mid\mathbf{Y})$, where $\mathbf{Y}\equiv(Y_{1,1}^{(1)},Y_{1,1}^{(2)},\ldots,Y_{F,I(F)}^{(J)})^{\prime}$ is the $(J\cdot\sum_{f=1}^{F}I(f))$-dimensional vector of all outputs from the MIP, $\bm{\sigma}^{2}\equiv((\sigma_{1}^{(1)})^{2},(\sigma_{1}^{(2)})^{2},\ldots,(\sigma_{F}^{(J)})^{2})^{\prime}$ is an $(F\cdot J)$-dimensional vector, and $\bm{\tau}^{2}=(\tau_{1}^{2},\ldots,\tau_{F}^{2})^{\prime}$ is an $F$-dimensional vector. The definition of the restricted likelihood and a derivation of its analytical expression are given in Appendix B.

Small sample sizes can lead to variance estimates that need to be regularized in a manner similar to Twomey-Tikhonov regularization. This is achieved through the penalized restricted likelihood:

for

where $\mathbf{b}\equiv(b_{1},\ldots,b_{F})^{\prime}$ and $g_{a}(x;b)=x^{-a-1}\exp(-b/x)$ for fixed $a>0$. The function $g_{a}(x;b)$ is proportional to the density of an inverse gamma distribution with shape $a$ and scale $b$, which is a common choice when making Bayesian inference on variance parameters. The value of the scales $\mathbf{b}$ in the penalty $p(\bm{\sigma}^{2};\mathbf{b})$ is considered unknown and is estimated, whereas the shape $a$ is fixed and pre-specified to guard against one team’s output being over-weighted in the BLUEs/BLUPs. That is, for a given factor $f\in\{1,\ldots,F\}$, the penalty puts a soft limit on the ratio between the largest and the smallest team-specific variances, $\{(\sigma_{f}^{(1)})^{2},\ldots,(\sigma_{f}^{(J)})^{2}\}$, to a degree specified by $a$, while the actual scale of the variances as a whole is left unknown and estimated through $b_{f}$.

Calibration of the shape $a$ can be achieved by interpreting the penalty as a Bayesian prior on the variances. In that case, we specify that the ratio of the prior 97.5th percentile of the variances to the prior 2.5th percentile of the variances is equal to four, which yields the fixed value of $a=8.48$. This choice means that it is unlikely that the largest variance is more than four times larger than the smallest variance.

The expression $L^{(p)}(\bm{\tau}^{2},\bm{\sigma}^{2}\mid\mathbf{Y};\mathbf{b})$ in (19) is maximized with respect to the model parameters $\bm{\tau}^{2}$ and $\bm{\sigma}^{2}$ and with respect to the penalization scale $\mathbf{b}$, to yield the penalized REML estimates $\{(\hat{\sigma}_{f}^{(j)})^{2}:f=1,\ldots,F,\enskip j=1,\ldots,J\}$ and $\{\hat{\tau}_{f}^{2}:f=1,\ldots,F\}$. Appendix B gives the formulas for the general case where correlations are present. For $f=1,\ldots,F$, these estimates are substituted into (13), the BLUE for $\mu_{f}$, to give the EBLUE of $\mu_{f}$, denoted by $\hat{\mu}_{f}$, with empirical variance $\hat{S}^{2}_{f}$. For $f=1,\ldots,F$, $i=1,\ldots,I(f)$, the EBLUP of $Y_{f,i}$, denoted by $\hat{Y}_{f,i}$, with empirical MSPE, denoted by $\widehat{\mathrm{MSPE}}(\hat{Y}_{f,i})$, is obtained similarly.

The statistical methodology given in Sections 2, 3, and 4 is applied in the next section to a MIP comprising nine flux inversion–based estimates of global carbon dioxide fluxes.

The rise in the atmospheric concentration of the greenhouse gas carbon dioxide (CO₂) is the main determinant of global warming (IPCC,, 2013), and estimation of the sources and sinks (called fluxes) of CO₂ is a high priority for the global community (Friedlingstein et al.,, 2014). Fluxes of CO₂ arise from both anthropogenic and natural processes, and the latter are hard to observe directly, as they occur over large spatial and temporal scales. One way to estimate the natural fluxes, called flux inversion, is through observation of their influence on atmospheric concentrations of CO₂. As natural fluxes are so difficult to measure otherwise, flux inversions for CO₂ play an important role in quantifying Earth’s carbon budget.

Flux inversion is an ill-posed inverse problem, in which an atmospheric transport model is used to simulate the transport of the trace gas from its sources to the times and places of its observation. The transport model is combined with a procedure, usually based on a Bayesian state-space model, to work backwards from observations to fluxes. The data used by inversions are typically a combination of surface, aircraft, and ship observations, and data collected by spaceborne instruments. In this section, we shall use the output from a MIP based on CO₂ data from NASA’s Orbiting Carbon Observatory-2 (OCO-2) satellite (Crowell et al.,, 2019) to obtain consensus flux estimates and predictions.

Flux inversion is subject to many forms of model uncertainty, which complicate interpretation of the estimated fluxes. These uncertainties stem from questions about the fidelity of simulated atmospheric transport, about how best to parameterize the unknown fluxes, and about the choice of prior for the flux field, among others. To address these in a rigorous manner, the OCO-2 project established a series of MIPs where an ensemble of flux outputs was assembled from submissions by different teams, who performed global flux inversions under a common protocol. The MIPs involve the use of OCO-2 retrievals of the vertical column–average CO₂ concentration, as well as flask and in situ observations of CO₂ concentrations.

The first of the MIPs concluded in 2019, and its results were reported by Crowell et al., (2019). The satellite data in this MIP were based on version 7r of the OCO-2 data, so henceforth we refer to it as MIPv7. Nine teams participated in MIPv7. The results were released as a set of aggregated monthly non-fossil-fuel fluxes spanning 27 disjoint, large regions derived from TransCom3 regions, shown in Figure A1 (where the unlabeled areas are assumed to have zero flux) and listed in Table A1. The outputs from the MIP covered each month over the two years from January 2015 to December 2016, inclusive. Each team was tasked with producing estimated fluxes from three observation types: from an inversion based on flask and in situ measurements only (IS), from an inversion based on OCO-2 retrievals made in land glint (LG) mode, and from an inversion based on OCO-2 retrievals made in land nadir (LN) mode.

To apply the SUPE-ANOVA framework to MIPv7, the fluxes from each mode were considered separately. Then, the factor of interest is $f=(s,r)$, where $s=1,\ldots,4$ is a season (starting with December-January-February [DJF]), and $r=1,\ldots,27$ is a spatial region. The replicates $i=1,\ldots,6$ are the six months within a season for the two-year study period, so $I(f)=6$ for all factors. From (2) and (3), the flux output of the $j$th team is modeled by,

for $j=1,\ldots,9$, where $\{\eta_{s,r,i}^{(j)}:j=1,\ldots,9\}$ are uncorrelated; we assume they are uncorrelated for simplicity, but, if desired, it is relatively easy to incorporate correlations using the equations in Appendix B. The consensus flux is modeled by

The term $Y_{s,r,i}$ represents the MIP consensus for the underlying geophysical flux for month $i$ in season $s$ and region $r$; $\mu_{s,r}$ is the more-slowly-varying consensus climatological flux for the season and region; $\tau_{s,r}^{2}$ is the variance of the deviations of the monthly fluxes from the climatological flux; and $(\sigma_{s,r}^{(j)})^{2}$ is the variance of the $j$th team’s flux output within each region and season.

The six replicates within each factor $(s,r)$ are too few for reliable estimation of the variance parameters, but this can be addressed straightforwardly by grouping regions with similar variability within a season. Hence, for the four seasons $s=1,\ldots,4$, variance parameters $\tau_{s,r}^{2}$, $\{(\sigma_{s,r}^{(j)})^{2}:j=1,\ldots,J\}$, and the penalization scales $b_{s,r}$ were constrained to have identical values within grouped regions. We present a clustering approach to construct groups in Appendix C; the chosen groups are shown in Figure A2.

The parameters of the model in (20) and (21) were estimated according to the initialization step given in Section 4, and the SUPE-ANOVA framework was used to obtain the consensus flux predictions $\{\hat{Y}_{s,r,i}:s=1,\ldots,4;\enskip r=1,\ldots,27;\enskip i=1,\ldots,6\}$. To check whether the model has good fit to the fluxes, we define the standardized errors, $\hat{\eta}_{s,r,i}^{(j)}=(Y_{s,r,i}^{(j)}-\hat{Y}_{s,r,i})/\hat{\sigma}_{s,r}^{(j)}$, for $s=1,\ldots,4$, $r=1,\ldots,27$, $i=1,\ldots,6$, and $j=1,\ldots,9$. Under the Gaussian assumptions used to fit the model, the standardized errors $\{\hat{\eta}_{s,r,i}^{(j)}\}$ should be approximately $\text{i.i.d.}\,\mathrm{Gau}(0,1)$. We diagnosed this visually using Q-Q plots, by plotting $\{\hat{\eta}_{s,r,i}^{(j)}\}$ for each season and group of regions against the theoretical quantiles of the $\mathrm{Gau}(0,1)$ distribution; under the assumptions, the relationship should be approximately on a 45$\degree$ line. Figure A3 shows these Q-Q plots, which do not indicate any obvious violations of the model assumptions.

As an example of predicting consensus fluxes at the regional level, Figure 2(a) shows the SUPE-ANOVA predictions as a time series for one region, T05a (Temperate Northern Africa), as well as the time series of fluxes from each of the nine MIP participants; the unweighted averages, $\bar{Y}_{s,r,i}=\frac{1}{9}\sum_{j=1}^{9}Y_{s,r,i}^{(j)}$, are also shown. Separate panels give the time series for each observation type (IS, LN, and LG). Figure 2(b) repeats Figure 2(a) but with the participating teams’ fluxes omitted to allow the vertical axis to be re-scaled. It also shows 95% prediction intervals for the predicted consensus fluxes, where the upper and lower values are $\hat{Y}_{s,r,i}\pm 1.96\times\sqrt{\widehat{\mathrm{MSPE}}(\hat{Y}_{s,r,i})}$. For the IS fluxes in the region T05a, the predicted consensus hardly differs from the unweighted average. By contrast, for the LN and LG fluxes, the predicted consensus sometimes differs significantly from the unweighted average.

The SUPE-ANOVA predictions of the consensus fluxes and the unweighted-average fluxes at the regional level were also aggregated spatially and temporally to yield the annual global land and annual global ocean fluxes. These are shown in Figure 3, which also shows the land–ocean aggregates for the nine teams in the MIP. Fluxes are given for all three observation types, IS, LN, and LG, and for both study years, 2015 and 2016. The SUPE-ANOVA consensus fluxes imply a larger land sink and a smaller ocean sink than do the unweighted averages. Also, for the IS observation type, there are more extreme values from the TM5-4DVAR team then the other teams, and these have been downweighted in the consensus estimate, relative to the weights they received in the unweighted average.

The different teams’ weights used in the SUPE-ANOVA framework vary between seasons and regions, although they are identical within the cluster groups. For a given $(s,r)$, define the climatological weights, $\hat{w}_{s,r}^{(j)}=(\hat{\sigma}_{s,r}^{(j)})^{-2}/\left[\sum_{j^{\prime}=1}^{J}(\hat{\sigma}_{s,r}^{(j^{\prime})})^{-2}\right]$, for $j=1,\ldots,9$. From (13), these weights correspond to the relative weighting of each participant in estimating the consensus mean climatological flux, $\mu_{s,r}$. Figure 4 shows the climatological weights for the three observation types (IS, LN, LG), for each season and grouping of regions within each season. From inspection of the weights in Figure 4, the UT team and the UoE team generally have higher weights due to their associated estimated variances being consistently smaller.

A reader of this article might ask whether in the North American temperature example given in Section 3, the simple unweighted average of the MIP outputs provides a reasonable estimate of $\mu_{\mathrm{Tmp}}$. Under the SUPE-ANOVA framework, that estimate $\bar{\mu}_{\mathrm{Tmp}}$ given by Equation (5) is unbiased. However, it is not a statistically efficient estimator under our consensus model. This can have a profound effect on the design of the MIP. Using the optimal weights, namely, $w_{\mathrm{Tmp}}^{(j)}$ inversely proportional to $(\sigma_{\mathrm{Tmp}}^{(j)})^{2}\equiv\mathrm{var}(Y_{\mathrm{Tmp}}^{(j)})$ for $j=1,\ldots,J$, means that the precision of $\mu_{\mathrm{Tmp}}^{*}=\sum_{j=1}^{J}w_{\mathrm{Tmp}}^{(j)}Y_{\mathrm{Tmp}}^{(j)}$ can be matched to the precision of $\bar{\mu}_{\mathrm{Tmp}}$ but with $J$ (much) less than 20. Since $\mu_{\mathrm{Tmp}}^{*}$ is also unbiased, this optimally weighted mean has a distinct inferential advantage over the unweighted mean, $\bar{\mu}_{\mathrm{Tmp}}$, when both are based on the same $J=20$ outputs $\{Y_{\mathrm{Tmp}}^{(1)},\ldots,Y_{\mathrm{Tmp}}^{(20)}\}$. Furthermore, the simple variance, $\sigma^{2}_{\mathrm{Tmp}}/20$, in Equation (5) is an incorrect expression for $\mathrm{var}(\bar{\mu}_{\mathrm{Tmp}})$ when it is no longer true that $(\sigma_{\mathrm{Tmp}}^{(1)})^{2}=\ldots=(\sigma_{\mathrm{Tmp}}^{(20)})^{2}=\sigma_{\mathrm{Tmp}}^{2}$. The correct expression is

which will give valid consensus inferences and correct coverage for the one-sigma and two-sigma intervals presented in Section 3, albeit wider than the intervals based on $\mathrm{var}(\mu_{\mathrm{Tmp}}^{*})$.

A number of articles reviewed in Section 1 have taken a Bayesian viewpoint (Tebaldi et al.,, 2005; Manning et al.,, 2009; Smith et al.,, 2009; Tebaldi and Sansó,, 2009; Chandler,, 2013). A Bayesian spatial ANOVA was used in Sain et al., (2011), Kang et al., (2012), and Kang and Cressie, (2013) to analyze projected temperature in a MIP called the North American Regional Climate Change Assessment Program (NARCCAP; see Mearns et al.,, 2012). This approach adds another layer of complexity to consensus inference in a MIP, where prior distributions on mean and variance parameters need to be specified, sensitivity to the specification assessed, and posterior-distribution sampling algorithms developed. In contrast, our SUPE-ANOVA framework does not depend on priors, is fast to compute, and results in easy-to-interpret weights on the $J$ MIP outputs within the ANOVA factors (e.g., season, region).

We return to the general SUPE-ANOVA framework developed in Sections 2 and 3 to show how the optimal weights could be used in conjunction with verification data. Let $\{Z_{k}:k=1,\ldots,K\}$ denote the verification data, with associated measurement-error variances $\{\mathrm{var}(Z_{k}):k=1,\ldots,K\}$. Each datum $Z_{k}$ can be indexed according to which factor $f$ and replicate $i\in\{1,\ldots,I(f)\}$ it belongs. Hence, rewrite $Z_{k}$ as $Z_{f(k),i(k)}$, which will be compared to the $J$ outputs, $\{Y_{f(k),i(k)}^{(j)}:j=1,\ldots,J\}$. Then a simple and statistically justifiable way to assign weights to each team’s output is to define

where $\sum_{j=1}^{J}\nu^{(j)}=1$. Equation (22) shows how the SUPE-ANOVA weights, $\{w_{f,i}^{(j)}\}$, can be used to define a prior in a Bayesian Model Averaging (BMA) scheme, where $\{\nu^{(j)}:j=1,\ldots,J\}$ can be interpreted as the posterior weights assigned to the $J$ models, given the verification data $\{Z_{k}:k=1,\ldots K\}$. That is, the $J$ models are weighted based on verification data without making the common ‘truth-plus-error’ assumption (Bishop and Abramowitz,, 2013), or even without making the more general assumption of co-exchangeability between the MIP outputs and the data (Rougier et al.,, 2013). Now, since the BMA weights have lost their dependence on the $F$ factors, they will not yield optimal estimates of $\mu_{f}$, or optimal predictions of $\{Y_{f,i}\}$; instead, they can be used to suggest which model is the most reliable overall in the sense of REA.

In summary, basic statistical principles have been followed to make optimal consensus inference in a MIP. Once the variance/covariance parameters in the SUPE-ANOVA framework have been estimated in an initialization step, the weights are immediate; code in the R programming language (R Core Team,, 2020) is available at https://github.com/mbertolacci/supe-anova-paper for the OCO-2 MIP analyzed in Section 5.

This project was largely supported by the Australian Research Council (ARC) Discovery Project (DP) DP190100180 and NASA ROSES grant NNH20ZDA001N-OCOST. Noel Cressie’s research was also supported by ARC DP150104576, and Andrew Zammit-Mangion’s research was also supported by an ARC Discovery Early Career Research Award (DECRA) DE180100203. The authors are grateful to the OCO-2 Flux Group for their feedback and suggestions.

The data for the OCO-2 MIP analyzed in Section 5 are described by Crowell et al., (2019) and may be found on the OCO-2 MIPv7 website hosted by the NOAA Global Monitoring Laboratory (https://gml.noaa.gov/ccgg/OCO2/).