A new statistical approach to climate change detection and attribution

Abstract

We propose here a new statistical approach to climate change detection and attribution that is based on additive decomposition and simple hypothesis testing. Most current statistical methods for detection and attribution rely on linear regression models where the observations are regressed onto expected response patterns to different external forcings. These methods do not use physical information provided by climate models regarding the expected response magnitudes to constrain the estimated responses to the forcings. Climate modelling uncertainty is difficult to take into account with regression based methods and is almost never treated explicitly. As an alternative to this approach, our statistical model is only based on the additivity assumption; the proposed method does not regress observations onto expected response patterns. We introduce estimation and testing procedures based on likelihood maximization, and show that climate modelling uncertainty can easily be accounted for. Some discussion is provided on how to practically estimate the climate modelling uncertainty based on an ensemble of opportunity. Our approach is based on the “models are statistically indistinguishable from the truth” paradigm, where the difference between any given model and the truth has the same distribution as the difference between any pair of models, but other choices might also be considered. The properties of this approach are illustrated and discussed based on synthetic data. Lastly, the method is applied to the linear trend in global mean temperature over the period 1951–2010. Consistent with the last IPCC assessment report, we find that most of the observed warming over this period (+0.65 K) is attributable to anthropogenic forcings (+0.67 $\pm$ 0.12 K, 90 % confidence range), with a very limited contribution from natural forcings ($-0.01\pm 0.02$ K).

Appendix

1.1 Model (1)–(3) as a linear regression Gaussian model

Based on the notation used in (1)–(3), we define

$$\begin{aligned} Z = \left( \begin{array}{c} Y \\ X_1 \\ \vdots \\ X_{n_f}\end{array} \right) , \quad \quad Z^* = \left( \begin{array}{c} X^*_1 \\ \vdots \\ X^*_{n_f}\end{array} \right) , \quad {\text { and }} \quad \varepsilon = \left( \begin{array}{c} \varepsilon _Y \\ \varepsilon _{X_1} \\ \vdots \\ \varepsilon _{X_{n_f}} \end{array} \right) . \end{aligned}$$

(35)

which are vectors of size $({n_f}+1)n$, ${n_f}n$, and $({n_f}+1)n$, respectively. Then, (2) and (3) may be written simultaneously as

$$\begin{aligned} Z = A Z^*+\varepsilon , \end{aligned}$$

(36)

where

$$\begin{aligned} A= \left( \begin{array}{cccc} I_n &{} I_n &{} \dots &{} I_n \\ I_n &{} 0 &{} \dots &{} 0 \\ 0 &{} I_n &{} \ddots &{} \vdots \\ \vdots &{} \ddots &{} \ddots &{} 0 \\ 0 &{} \dots &{} 0 &{} I_n \end{array} \right) , \end{aligned}$$

(37)

is the $({n_f}+1)n \times {n_f}n$ design matrix of the linear Gaussian model (36).

This representation is useful to derive some of the statistical properties of our model. It is used, for instance, in “Appendix 8.2.3”. $A'A$ has also a simple closed form. However, it is difficult to derive the precision matrix $\left ( A'A \right )^{-1}$ in closed form, which would be required to deduce, e.g., the MLEs more directly than in Sect. 3.3.

1.2 Hypothesis testing details

1.2.1 Goodness of fit tests

The tests proposed in Sect. 3.5 are goodness of fit tests constructed as follows.

If $Z \in {\mathbb {R}}^n$ is a random vector with distribution $N(\mu ,{\varSigma })$ under a given model, with $\mu$ and ${\varSigma }$ known, we will call “goodness of fit” test of this model the deviance test, i.e. the likelihood ratio test (LRT) with respect to a saturated alternative hypothesis $H_1: Z \sim N(\theta ,{\varSigma })$, where $\theta \in {\mathbb {R}}^n$ is unknown. The log-likelihood is 0 under this alternative, so the LRT only involves $-2$ log-likelihood under $H_0$, the considered model. Therefore the LRT statistic is

$$\begin{aligned} T = (Z-\mu )' {\varSigma }^{-1} (Z-\mu ), \end{aligned}$$

(38)

which follows a $\chi ^2(n)$ distribution under $H_0$.

1.2.2 Minimized likelihood under $H_0: Y^* = 0$ in model (1)–(3)

Under $H_0: Y^*=0$, each $X_i^*$ has to be estimated under the additional constraint that $\sum _{i=1}^{{n_f}} X_i^* = 0$. Under $H_0$, (6) becomes

$$\begin{aligned} \ell _{H_0}(X^*_i) = Y' {\varSigma }_Y^{-1} Y + \sum _{i=1}^{{n_f}} \left (X_i - X^*_i\right )' {\varSigma }_{X_i}^{-1} \left (X_i - X^*_i\right ). \end{aligned}$$

(39)

The gradient of $\ell _{H_0}$ with respect to $X_i^*$ is

$$\begin{aligned} {\varSigma }_{X_i}^{-1} (\widehat{X}_i^* - X_i). \end{aligned}$$

(40)

The gradient of the constraint is ${\mathbbm {1}}_n$, i.e. the “all-one” vector of size n. The theory of Lagrange multipliers imposes that at the constrained minimum, these two gradients are proportional, so that

$$\begin{aligned} {\varSigma }_{X_i}^{-1} (\widehat{X}_i^* - X_i) = \lambda {\mathbbm {1}}_n, \end{aligned}$$

(41)

where $\lambda$ is some constant. From (41),

$$\begin{aligned} \widehat{X}_i^* = X_i + \lambda {\varSigma }_{X_i} {\mathbbm {1}}_n, \end{aligned}$$

(42)

and thus the constraint $\sum _{i=1}^{{n_f}} \widehat{X}_i^* =0$ gives

$$\begin{aligned} \lambda {\mathbbm {1}}_n = - {\varSigma }_X^{-1} X, \end{aligned}$$

(43)

where $X = \sum _{i=1}^{{n_f}} X_i$ and ${\varSigma }_X= \sum _{i=1}^{{n_f}} {\varSigma }_{X_i}$. The maximum likelihood estimate of $X_i^*$ under $H_0$ is therefore

$$\begin{aligned} \widehat{X}_i^* = X_i - {\varSigma }_{X_i} {\varSigma }_X^{-1} X. \end{aligned}$$

(44)

Finally, the minimized value of $-2$ log-likelihood under $H_0$ is

$$\begin{aligned} \ell _{H_0}(\widehat{X}^*_i) = Y' {\varSigma }_Y^{-1} Y + X' {\varSigma }_X^{-1} X. \end{aligned}$$

(45)

This result is used in “Appendix 8.2.3”, Eq.  (46), to derive the LRT of $H_0: Y^*=0$ within the statistical model defined by (1)–(3).

1.2.3 Choice of the null and the alternative hypotheses in the detection test

The “detection” test deals with the null hypothesis of no change, i.e. $H_0: Y^*=0$. The detection test we propose in (18) is a goodness of fit test based on (2) only. In particular, it doesn’t consider (1) and (3). Here, we present two alternatives that might be considered to test the same null hypothesis. Both are LRTs between two nested well-defined hypotheses on the data (Y, X).

• In the statistical model (1)–(3), consider the null-hypothesis $H_0: Y^*=0$ versus the alternative hypothesis $H_1:$(1)–(3), ie the same model with an unspecified $Y^*$. Following (45) and (14), this test would be based on the statistic

  $$\begin{aligned} \ell _{H_0} - \ell _{H_1} = Y' {\varSigma }_Y^{-1} Y + X' {\varSigma }_X^{-1} X - (Y-X)' ({\varSigma }_Y+{\varSigma }_X)^{-1} (Y-X) \sim _{H_0} \chi ^2(n). \end{aligned}$$

  (46)

  Note that the distribution under $H_0$ is known to be $\chi ^2(n)$ since $H_0$ is a linear sub-hypothesis of $H_1$, and they differ by a dimension of n (see 8.1).

• In the statistical model (1)–(3), consider the null-hypothesis $H_0: Y^*=0$ versus the saturated alternative hypothesis. This test is a goodness of fit test as defined above. Such a test would be based on the statistic (see (45))

  $$\begin{aligned} \ell _{H_0} = Y' {\varSigma }_Y^{-1} Y + X' {\varSigma }_X^{-1} X \sim _{H_0} \chi ^2(2n), \end{aligned}$$

  (47)

  where the null-distribution is deduced directly from (2) and (3).

Both of these tests would treat information from Y and X very symmetrically because, given (1), $Y^*=0$ implies not only that Y is small, but also that X is small (as $X^*=0$). Therefore, rejection of these tests may happen because either Y or X are “large”. In the case where X is large while Y is not, detection would not be a consequence of an abnormal observation, but rather the consequence of too large a response simulated in climate models. This, of course, is not consistent with the definition of detection. In our opinion, it is then more appropriate to discuss the first term in the right hand side of (45) separately.

We further argue that there is a fundamental distinction between detection and attribution in this respect. As detection only assesses whether observations are consistent with internal variability, historical simulations by climate models (X) are not required for detection. The only requirement is to quantify internal variability—which is usually done based on other simulations by climate models. Attribution, however, definitively relies on historical simulations to disentangle contributions from different forcings based on some physical knowledge. This fundamental difference is the main reason why we propose to base detection on (2) only, while (1)–(3) are considered as a whole to perform attribution, and in particular to estimate the contributions of individual forcings.

As a last remark, the test defined in (18) may also be considered as a test of the null hypothesis $H_0: ``Y^*=0$ and (1) does not hold”. Indeed, removing (1) implies that the second term in 45 is zero, and then our test (18) is actually a LRT against a saturated alternative hypothesis.

1.3 Estimation of ${\varSigma }_{\text {m}}$ and ${\varSigma }_X$ with unbalanced data

This section deals with the realistic case where models have run ensembles of historical simulations of various sizes. For instance, in the CMIP5 archive, the number of historical simulations with all external forcings varies from 1 to 10 depending on the model considered. Here, we mainly discuss the estimation of ${\varSigma }_{\text {m}}$, which is not explicitly addressed in the main text. We also mention what ${\varSigma }_X$ should be considered in this unbalanced case.

Consistent with Sect. 4, we assume that simulation k of model j can be decomposed as

$$\begin{aligned} w_{jk} = \mu + m_j + \epsilon _{jk}, \quad j=1,\dots ,n_m, \quad k=1,\dots ,n_j, \end{aligned}$$

(48)

where $m_j \sim N(0,{\varSigma }_{\text {m}})$, and $\epsilon _{jk} \sim N(0,{\varSigma }_{\text {v}})$, leading to

$$\begin{aligned} w_{jk} \sim N \left (\mu ,{\varSigma }_{\text {m}}+{\varSigma }_{\text {v}}\right ). \end{aligned}$$

(49)

We then introduce

$$\begin{aligned} w_{j.} = \frac{1}{n_j} \sum _{i=1}^{n_j} w_{jk} \sim N \left (0, {\varSigma }_{\text {m}}+\frac{{\varSigma }_{\text {v}}}{n_j} \right ). \end{aligned}$$

(50)

This kind of framework is actually a multivariate linear mixed model, and useful references could be found in Rao and Kleffe (1988). Within such models, the main challenge comes from the estimation of variance components (here ${\varSigma }_{\text {m}}$ and ${\varSigma }_{\text {v}}$), and no optimal estimator is known in the general unbalanced case (i.e. if the ensemble sizes $n_j$ are not all equal). We propose to use a method of moments approach very similar to that proposed by Henderson (1953), but with a couple of specific features:

• Consistent with common practice in climate science, we estimate the fixed effect $\mu$ as the mean of the ensemble means from each model. In this way each model is given equal weight, disregarding the number of simulations performed. Consequently, we consider $\widehat{\mu } = \overline{w}= \frac{1}{n_m} \sum _{j=1}^{n_m} w_{j.}$, whereas the common approach in statistics uses $\widehat{\mu } = w_{..} = \frac{1}{n} \sum _{j,k} w_{jk}$.

• As a large number of unforced simulations are available to estimate ${\varSigma }_{\text {v}}$ (Table 1), we assume that this term is already known, and only focus on the estimation of ${\varSigma }_{\text {m}}$. Note that within-ensemble differences might also be used to estimate ${\varSigma }_{\text {v}}$.

One has:

$$\begin{aligned} \overline{w}= \frac{1}{n_m} \sum _{j=1}^{n_m} w_{j.} \quad \sim \quad N \left ( \mu ,\frac{1}{n_m} {\varSigma }_{\text {m}}+ \frac{1}{n_m^2} \sum _{j=1}^{n_m} \frac{1}{n_j} {\varSigma }_{\text {v}}\right ). \end{aligned}$$

(51)

$$\begin{aligned} {\text {Var}} \left ( w_{j.} - \overline{w}\right ) = \left (1-\frac{1}{n_m} \right ) {\varSigma }_{\text {m}}+ \left ( \frac{1}{n_j}-\frac{2}{n_m n_j}+\frac{1}{n_m^2} \sum _{j=1}^{n_m} \frac{1}{n_j} \right ) {\varSigma }_{\text {v}},\end{aligned}$$

(52)

$$\begin{aligned} {\text {E}}\ \left ( \sum _{j=1}^{n_m} (w_{j.} - \overline{w})^2 \right ) = (n_m-1) {\varSigma }_{\text {m}}+ \frac{n_m-1}{n_m} \sum _{j=1}^{n_m} \frac{1}{n_j} {\varSigma }_{\text {v}}. \end{aligned}$$

(53)

Using the method of moments, we estimate this quantity with

$$\begin{aligned} SSM = \sum _{j=1}^{n_m} (w_{j.} - \overline{w})^2. \end{aligned}$$

(54)

Finally, given an estimate of ${\varSigma }_{\text {v}}$, we estimate ${\varSigma }_{\text {m}}$ with

$$\begin{aligned} \widehat{{\varSigma }}_{\text {m}} = \frac{1}{n_m-1} \left ( SSM - \frac{n_m-1}{n_m} \sum _{j=1}^{n_m} \frac{1}{n_j} {\varSigma }_{\text {v}}\right )_+, \end{aligned}$$

(55)

where $A_+$ means that negative eigenvalues of A are set to 0. $\widehat{{\varSigma }}_{\text {m}}$ is a truncation of a quadratic unbiased estimator, very similar to Henderson (1953) or MIVQUE estimators (Rao and Kleffe 1988). A nice property of this approach is that it can be used even if the $w_{j.}$ only are known (e.g. the $w_{jk}$ are not observed). This happens, for instance, if one ensemble of simulations has not been performed, and the response of each model is estimated by subtraction, e.g. $w^{ANT} = w^{ALL} - w^{NAT}$. In such a case, we could assume, for model j:

$$\begin{aligned} w^{ANT}_{j.} = w^{ALL}_{j.} - w^{NAT}_{j.} \sim N \left ( 0,{\varSigma }_{\text {m}}^{ANT}+\frac{{\varSigma }_{\text {v}}}{n_j^{ANT}} \right ), \end{aligned}$$

where $n_j^{ANT}$ is defined by

$$\begin{aligned} \frac{1}{n_j^{ANT}} = \frac{1}{n_j^{ALL}} + \frac{1}{n_j^{NAT}}. \end{aligned}$$

(56)

Finally, under this unbalanced assumption, putting (51 52 53) together with $(\mu - w^*) \sim N(0,{\varSigma }_{\text {m}})$ suggests considering

$$\begin{aligned} {\varSigma }_X = \left ( 1 + \frac{1}{n_m} \right ) \widehat{{\varSigma }}_{\text {m}} + \frac{1}{n_m^2} \sum _{j=1}^{n_m} \frac{1}{n_j} {\varSigma }_{\text {v}}. \end{aligned}$$

(57)

Fig. 5

Same as Fig. 5 based on the “models are centered on the truth” paradigm

1.4 Application to global mean temperature using the “models are centered on the truth” paradigm

In this section, we present and briefly discuss the results obtained in the analysis of global mean temperature, when the “models are centered on the truth” paradigm is used instead of the “models are statistically indistinguishable from the truth” paradigm, which was used in Sect. 6.2.

Under this paradigm, outputs from each climate model, $w_j$, may be regarded as being sampled independently from the same distribution, which is centered on the truth, i.e. $w_j \sim N(w,{\varSigma }_m)$, where $w$ is the true value of the simulated parameter, and ${\varSigma }_m$ describes the climate modelling uncertainty on this parameter. The distribution of the multimodel mean is then given by $\overline{w}\sim N(w, {\varSigma }_m/n_m)$. Finally, ${\varSigma }_X = {\varSigma }_m /n_m$ has to be considered under this paradigm (if internal variability is neglected).

The primary consequence of considering this alternative paradigm is to narrow the climate modelling uncertainty. Consequently, the assessment of consistency between models and observations is usually more demanding, while more weight is given to models in the estimation of individual forcing contributions.

Under this revised assumption, panel a) is unchanged.

In panel b), observations are found to be barely consistent with models (p value 0.09). The expected ALL warming, based on climate models only, is $[+0.74$K$, +0.86$K] (90 % confidence interval). This is a much narrower interval than reported in Sect.  6.2 ($[+0.44$K$, +1.16$K]). The weak consistency with observations mentioned above is then related to internal variability, which has to be added to these numbers. After inference, the estimated past warming lies between $[+0.72$K$, +0.83$K], which is still substantially greater than the observed value of $+0.65K$. This can be understood as follows: because there is little uncertainty in the estimate provided by climate models, the method considers that internal variability is partly responsible for the low observed value. The overall forced change is then estimated to be higher than that found in raw observations.

In panel c), results are similarly impacted. The expected (climate models only) NAT response is unchanged, and the ANT response is expected to lies within $[+0.67$K$, +0.93$K], which is again quite a narrow interval. After the inference is performed, the ANT response is estimated to be within $[+0.64$K$, +0.82$K]. If compared to the results given in Sect. 6.2, the impact of changing the paradigm is limited. The main impact of changing the paradigm is to discard the lowest values, from $+0.55$K to $+0.64$K.

Overall, these results suggest that the “models are statistically indistinguishable from the truth” paradigm, which was used in the main text, is more appropriate to ensure consistency between models and observations, and avoids over emphasizing the climate models outputs.