Regularized Fingerprinting in Detection and Attribution of Climate Change with Weight Matrix Optimizing the Efficiency in Scaling Factor Estimation

Detection and attribution analyses of climate change are critical components in establishing a causal relationship from the human emission of greenhouse gases to the warming of planet Earth (e.g., Bindoff et al., 2013). In climate science, detection is the process of demonstrating that a climate variable has changed in some defined statistical sense without providing a reason for that change; attribution is the process of evaluating the relative contributions of multiple causal factors to a change or event with an assignment of statistical confidence (e.g., Hegerl and Zwiers, 2011). Casual factors usually refer to external forcings, which may be anthropogenic (e.g., greenhouse gases, aerosols, ozone precursors, land use) and/or natural (e.g., volcanic eruptions, solar cycle modulations). By comparing simulated results of climate models with observed climate variables, a detection and attribution analysis evaluates the consistency of observed changes with the expected response, also known as fingerprint, of the climate system under each external forcing.

Optimal fingerprinting is the most widely used method for detection and attribution analyses (e.g., Hasselmann, 1997; Allen and Stott, 2003; Hegerl et al., 2010). Fingerprinting is a procedure that regresses the observed climate variable of interest on the fingerprints of external forcings, and checks whether the fingerprints are found in and consistent with the observed data. The central target of statistical inferences here is the regression coefficients, also known as scaling factors, which scale the fingerprints of external forcings to best match the observed climate change. Historically, the method was “optimal” in the context of generalized least squares (GLS) when the precision matrix of the regression error is used as weight, such that the resulting estimator of the scaling factors have the smallest variances. It was later recognized that the fingerprint covariates are not observed but estimated from climate model simulations. This leads to an errors-in-variables (EIV) issue, which has been approached by the total least squares (TLS) (Allen and Stott, 2003) with both the response and the covariates “prewhittened” by the covariance matrix of the error. The covariance matrix represents the internal climate variation. In practice, it is unknown and has to be estimated from climate model simulations (Allen and Tett, 1999; Ribes et al., 2013), which is generally handled preliminarily and independently from the regression inference (e.g., Hannart et al., 2014).

Estimating the error covariance matrix in fingerprinting is challenging because the available runs from climate model simulations are small relative to the dimension of the covariance matrix. The dimension of the covariate matrix is often reduced by considering averages over 5-year or 10-year blocks in time and/or over aggregated grid-boxes in space. For example, decadal average over a 110-year period in 25 grid-boxes would lead to a dimension of $275\times 275$. The sample size from climate model simulations that can be used is, however, at most a few hundreds. The sample covariance matrix is not invertible when the sample size of climate model simulations is less than its dimension. Earlier methods project data onto the leading empirical orthogonal functions (EOF) of the internal climate variability (Hegerl et al., 1996; Allen and Tett, 1999). More recently, the regularized optimal fingerprinting (ROF) method avoids the projection step by a regularized covariance matrix estimation (Ribes et al., 2009), which is based on the linear shrinkage covariance matrix estimator of Ledoit and Wolf (2004). The ROF method has been shown to provide a more accurate implementation of optimal fingerprinting than the EOF projection method (Ribes et al., 2013).

The uncertainty in estimating the error covariance matrix has important implications in optimal fingerprinting. The optimality in inferences about the scaling factor in optimal fingerprinting was historically based on the assumption that the error covariance matrix is known. The properties of the scaling factor estimator obtained by substituting the error covariance matrix with an estimate have not been thoroughly investigated in the literature. For example, it is not until recently that the confidence intervals for the scaling factors constructed from asymptotic normal approximation (Fuller, 1980; Gleser, 1981) or bootstrap (Pesta, 2013) were reported to be overly narrow when the matrix is only known up to a scale (DelSole et al., 2019) or completely unknown (Li et al., 2019). A natural, fundamental question is: when the error covariance matrix is estimated, are the confidence intervals constructed using the optimal approach under the assumption of known error covariance matrix still optimal? Since optimality under unknown error covariance matrix is practically infeasible, to be precise, we term fingerprinting for optimal fingerprinting and regularized fingerprinting (RF) for ROF in the sequel.

The contributions of this paper are two-fold. First, we develop a new method to construct the weight matrix in RF by minimizing directly the uncertainty of the resulting scaling factor estimators. The weight matrix is the inverse of a nonlinear shrinkage estimator of the error covariance matrix inspired by Ledoit and Wolf (2017a). We first extend the validity of their nonlinear shrinkage estimator to the context of RF via GLS regression with no EIV. We show that the proposed method is asymptotically optimal as the sample size of climate model simulations and the matrix dimension go to infinity together with a fixed ratio. When there is EIV, as is the case in practice, we show that the proposed weight is more efficient than the existing weight in RF (Ribes et al., 2013) in terms of the asymptotic variance of the scaling factor estimator when RF is conducted with generalized TLS (GTLS). This is why we refer to the current practice by RF instead of ROF. The second contribution is practical recommendations for assumptions about the structure of the error covariance matrix under which the sample covariance is estimated before any shrinkage in RF based on findings of a comparison study under various realistic settings. An implementation of both linear shrinkage and nonlinear shrinkage estimators is publicly available in an R package dacc (Li et al., 2020) for detection and attribution of climate change.

The rest of this article is organized as follows. After a review of RF in Section 2, we develop the proposed weight matrix and the theoretical results to support the asymptotic performance of proposed method in Section 3. A large scale numerical study assessing the performance of the proposed method is reported in Section 4. In Section 5, we demonstrate the proposed method with detection and attribution analysis on global and regional scales. A discussion concludes in Section 6. The technical proofs of the theoretical results are relegated in the Supplementary Materials.

Fingerprinting takes the form of a linear regression with EIV

where $Y$ is a $N\times 1$ vector of the observed climate variable of interest, $X_{i}$ is the true but unobserved $N\times 1$ fingerprint vector the $i$th external forcing with scaling factor $\beta_{i}$, $\epsilon$ is a $N\times 1$ vector of normally distributed regression error with mean zero and covariance matrix $\Sigma$, $\tilde{X}_{i}$ is an estimate of $X_{i}$ based on $n_{i}$ climate model simulations under the $i$th external forcing, and $\nu_{i}$ is a normally distributed measurement error vector with mean zero and covariance matrix $\Sigma/n_{i}$, and $\nu_{i}$’s are mutually independent and independent of $\epsilon$, $i=1,\ldots,p$. The covariance matrices of $\nu_{i}$’s and $\epsilon$ only differ in scales under the assumption that the climate models reflect the real climate variation. No intercept is present in the regression because the response and the covariates are centered by the same reference level. The primary target of inference is the scaling factors $\beta=(\beta_{1},\ldots,\beta_{p})^{\top}$.

The “optimal” in optimal fingerprinting originated from earlier practices under two assumptions: 1) the error covariance matrix $\Sigma$ is known; and 2) $X_{i}$’s are observed. The GLS estimator of $\beta$ with weight matrix $W$ is

where $X=(X_{1},\ldots,X_{p})$. The covariance matrix of the estimator $\hat{\beta}(W)$ is

The optimal weight matrix is $W=\Sigma^{-1}$, in which case, $\hat{\beta}(\Sigma^{-1})$ is the best linear unbiased estimator of $\beta$ with covariance matrix $X^{\top}\Sigma^{-1}X$. Since $\Sigma$ is unknown, a feasible version of GLS uses $W=\hat{\Sigma}^{-1}$, where $\hat{\Sigma}$ is an estimator of $\Sigma$ obtained separately from controlled runs of climate model simulations.

Later on it was recognized that, instead of $X_{i}$’s, only their estimates $\tilde{X}_{i}$’s are observed and that using $\tilde{X}_{i}$’s in place of $X_{i}$’s leads to bias in estimating $\beta$ (Allen and Stott, 2003). If $\Sigma$ is given, the same structure (up to a scale $1/n_{i}$) of the covariance matrices of $\nu_{i}$’s and $\epsilon$ allows precise pre-whitening of both $Y$ and $\tilde{X}_{i}$’s . Then the TLS can be applied to the pre-whitened variables. Inferences about $\beta$ can be based on the asymptotic normal distribution of the TLS estimator of $\beta$ (Gleser, 1981) or nonparametric bootstrap (Pesta, 2013), as recently studied by DelSole et al. (2019). Similar to the GLS setting, a feasible version of the GTLS procedure relies on an estimator of $\Sigma$.

The current practice of fingerprinting consists of two separate steps. First, estimate $\Sigma$ from controlled runs of climate model simulations under the assumption that the climate models capture the internal variability of the real climate system. Second, use this estimated matrix to pre-whiten both the outcome and covariates in the regression model (1)–(2), and obtain the GTLS estimator of $\beta$ on the pre-whitened data. Nonetheless, estimation of $\Sigma$ in the first step is not an easy task. The dimension of $\Sigma$ is $N\times N$, with $N(N+1)/2$ parameters if no structure is imposed, which is too large for the sample size $n$ of available climate model simulations (usually in a few hundreds at most). The sample covariance matrix based on the $n$ runs is a start, but it is of too much variation; when $N>n$ it is not even invertible. The linear shrinkage method of Ledoit and Wolf (2004) regularizes the sample covariance matrix $\hat{\Sigma}_{n}$ to in the form of $\lambda\hat{\Sigma}_{n}+\rho I$, where $\lambda$ and $\rho$ are scalar tuning parameters and $I$ is the identity matrix. This class of shrinkage estimators has the effect of shrinking the set of sample eigenvalues by reducing its dispersion around the mean, pushing up the smaller ones and pulling down the larger ones. This estimator has been used in the current RF practice (Ribes et al., 2009, 2013).

Substituting $\Sigma$ with an estimator introduces an additional uncertainty. The impact of this uncertainty on the properties of resulting ROF estimator has not been investigated when the whole structure of $\Sigma$ is unknown (Li et al., 2019). The optimality of the optimal fingerprinting in its original sense is unlikely to still hold. Now that the properties of the resulting estimator of $\beta$ depends on an estimated weight matrix, can we choose this weight matrix estimator to minimize the variance of the estimator of $\beta$? The recently proposed nonlinear shrinkage estimator (Ledoit and Wolf, 2017a, 2018) has high potential to outperform the linear shrinkage estimators.

We consider constructing the weight matrix by inverting a nonlinear shrinkage estimator of $\Sigma$ (Ledoit and Wolf, 2017b) in the fingerprinting context. New theoretical results are developed to justify the adaptation of this nonlinear shrinkage estimator of $\Sigma$ to minimize the uncertainty of the resulting estimator $\hat{\beta}$ of $\beta$. Assume that there are $n$ replicates from climate model simulations (usually pre-industrial control runs) that are independent of $Y$ and $\tilde{X}_{i}$’s. Let $Z_{i},Z_{2},\ldots,Z_{n}\in\mathbb{R}^{N}$ be the centered replicates so that the sample covariance matrix is computed as $\hat{\Sigma}_{n}=n^{-1}\sum_{i=1}^{n}Z_{i}Z^{\top}_{i}$. Our strategy is to revisit the GLS setting with no EIV first and then apply the result of the GTLS setting to the case under EIV, the same order as the historical development.

The finite sample performance of the proposed method in comparison with the existing practice in RF needs to be assessed to make realistic recommendations for detection and attribution of climate change. We conducted a simulation study similar to the setting of a study in Ribes et al. (2013). The observed climate variable of interest is 11 decadal mean temperatures over 25 grid-boxes, a vector of dimension $N=275$. Two $N\times 1$ fingerprints were considered, corresponding to the anthropogenic (ANT) and natural forcings (NAT), denoted by $X_{1}$ and $X_{2}$, respectively. They were set to the average of all runs from the CNRM-CM5 model as in Ribes et al. (2013). To vary the strength of the signals, we also considered halving $X_{1}$ and $X_{2}$. That is, there were two levels of signal-to-noise ratio corresponding to the cases of multiplying each $X_{i}$, $i\in\{1,2\}$, controlled by a scale $\lambda\in\{1,0.5\}$. The true scaling factors were $\beta_{1}=\beta_{2}=1$. The distribution of the error vector $\epsilon$ was multivariate normal $\mathrm{MVN}(0,\Sigma)$. The distribution of the measurement error vector $\nu_{i}$, $i\in\{1,2\}$, was $\mathrm{MVN}(0,\Sigma/n_{i})$, with $(n_{1},n_{2})=(35,46)$ which are the values in the detection and attribution analysis of annual mean temperature conducted in Section 5.

Two settings of true $\Sigma$ were considered. In the first setting, $\Sigma$ was an unstructured matrix $\Sigma_{\mathrm{UN}}$, which was obtained by manipulating the eigenvalues but keeping the eigenvectors of the proposed minimum variance estimate from the same set of climate model simulations as in Ribes et al. (2013). Specifically, we first obtained the eigen decomposition of the minimum variance estimate, and then restricted the eigenvalues to be equal over each of the 25 grid-boxes (i.e., only 25 unique values for the $N=25\times 11$ eigenvalues) by taking averages over the decades at each grid-box. The pattern of the resulting eigenvalues is similar to the pattern of the eigenvalues of a spatial-temporal covariance matrix with variance stationarity and weak dependence over the time dimension. Finally, the eigenvalues were scaled independently by a uniformly distributed variable on $[0.5,1.5]$, which results in a more unstructured covariance matrix setting similar to the simulation settings in Hannart (2016). In the second setting, $\Sigma$ was set to be $\Sigma_{\mathrm{ST}}$ whose diagonals were set to be the sample variances from the climate model simulations without imposing temporal stationarity; the corresponding correlation matrix was set to be the Kronecker product of a spatial correlation matrix and a temporal correlation matrix, both with autoregressive of order 1 and coefficient 0.1.

The observed mean temperature vector $Y$ and the estimated fingerprints $(\tilde{X}_{1},\tilde{X}_{2})$ were generated from Models (1)–(2). The control runs used to estimate $\Sigma$ were generated from $\mathrm{MVN}(0,\Sigma)$ with sample size $n\in\{50,100,200,400\}$. For each replicate, the two GTLS estimators of $\beta$ in Theorem 3 were obtained. For each configuration, 1000 replicates were run.

Figure 1 displays the boxplots of the estimates of the ANT scaling factor $\beta_{1}$ from the two TGLS approaches based on pre-whitening with $\hat{\Sigma}_{\mathrm{LS}}$ (denoted as M1) and $\hat{\Sigma}_{\mathrm{MV}}$ (denoted as M2) and, respectively. Both estimators appear to recover the true parameter values well on average. The variations of both estimators are lower for larger $n$, higher $\lambda$, and more structured $\Sigma$ (the case of $\Sigma_{\mathrm{ST}}$). These observations are expected. A larger $n$ means more accurate estimation $\Sigma$; a higher $\lambda$ means stronger signal; a more structured $\Sigma$ means an easier task to estimate $\Sigma$. In the case of $\Sigma=\Sigma_{\mathrm{UN}}$, the M2 estimates have smaller variations than the M1 estimate, since the eigenvalues were less smooth and, hence, favored the nonlinear shrinkage function. For the case of $\Sigma=\Sigma_{ST}$ where the covariance matrix is more structured, both methods estimate the true covariance matrix much more accurately, and the differences between methods are less obvious. More detailed results are summarized in Table D.1 and Table D.2, the latter of which had smaller ensembles in estimating the fingerprints with $(n_{1},n_{2})=(10,6)$. The standard deviations of the M2 estimates are over 10% smaller than those of the M1 estimates for both cases.

Confidence intervals are an important tool for detection and attribution analyses. It would be desirable if the asymptotic variance in Theorem 2 can be used to construct confidence intervals for the scaling factors. Unfortunately, it has been reported that the confidence intervals for the scaling factors based on $\hat{\Sigma}_{\mathrm{LS}}$ have coverage rates lower than, sometimes much lower than, their nominal levels (Li et al., 2019). The under-coverage issue remains for the estimator based on $\hat{\Sigma}_{\mathrm{MV}}$. To give confidence intervals with correct coverage rates, Li et al. (2019) proposed a calibration procedure which enlarges the confidence intervals based on the asymptotic normality of the estimators by an appropriate scale tuned by a parametric bootstrap to achieve the desired coverage rate. We applied this calibration to both M1 and M2 estimators. Figure 2 shows the empirical coverage rates of the 95% confidence intervals after the calibration. The coverage rate of a naive confidence interval could be as low as 70% (not shown). After the calibration, the coverage rates are much closer to the nominal levels. The agreement is better for larger $n$ and more structured $\Sigma$. The calibrated confidence intervals from M2 are about 10% shorter to those from M1 overall in both settings of true $\Sigma$, except for the case of $\Sigma=\Sigma_{\mathrm{UN}}$ and sample size $n=50$ where the confidence intervals suffer from under-coverage issue..

We apply the proposed approach to the detection and attribution analyses of annual mean temperature of 1951–2010 at the global (GL), continental, and subcontinental scales. The continental scale regions are Northern Hemisphere (NH), NH midlatitude between $30^{\circ}$ and $70^{\circ}$ (NHM), Eurasia (EA), and North America (NA), which were studied in (Zhang et al., 2006). The subcontinental scale regions are Western North American (WNA), Central North American (CNA), Eastern North American (ENA), southern Canada (SCA), and southern Europe (SEU), where the spatio-temporal correlation structure is more likely to hold.

In each regional analysis, we first constructed observation vector $Y$ from the HadCRUT4 dataset (Morice et al., 2012). The raw data were monthly anomalies of near-surface air temperature on $5^{\circ}\times 5^{\circ}$ grid-boxes. At each grid-box, each annual mean temperature was computed from monthly temperatures if at least 9 months were available in that year; otherwise, it was considered missing. Then, 5-year averages were computed if no more than 2 annual averages were missing. To reduce the spatial dimension in the analyses, the $5^{\circ}\times 5^{\circ}$ grid-boxes were aggregated into larger grid-boxes. In particular, the grid-box sizes were $40^{\circ}\times 30^{\circ}$ for GL and NH, $40^{\circ}\times 10^{\circ}$ for NH 30-70, $10^{\circ}\times 20^{\circ}$ for EA, and $10^{\circ}\times 5^{\circ}$ for NA. For the subcontinent regions, no aggregation was done except for SCA, in which case $10^{\circ}\times 10^{\circ}$ grid-boxes were used. Details on the longitude, latitude and spatio-temporal steps of each regions after processing can be found in Table 1.

Two external forcings, ANT and NAT, were considered. Their fingerprints $X_{1}$ and $X_{2}$ were not observed, but their estimates $\tilde{X}_{1}$ and $\tilde{X}_{2}$ were averages over $n_{1}=35$ and $n_{2}=46$ runs from CIMP5 climate model simulations. The missing pattern in $Y$ was used to mask the simulated runs. The same procedure used to aggregate the grid-boxes and obtain the 5-year averages in preparing $Y$ was applied to each masked run of each climate model under each forcing. The final estimates $\tilde{X}_{1}$ and $\tilde{X}_{2}$ at each grid-box were averages over all available runs under the ANT and the NAT forcings, respectively, centered by the average of the observed annual temperatures over 1961–1990, the same center used by the HadCRUT4 data to obtain the observed anomalies.

Estimation of $\Sigma$ was based on $n=223$ runs of 60 years constructed from 47 pre-industrial control simulations of various length. The long-term linear trend was removed separately from the control simulations at each grid-box. As the internal climate variation is assumed to be stationary over time, each control run was first split into non-overlapping blocks of 60 years, and then each 60-year block was masked by the same missing pattern as the HadCRUT4 data to create up to 12 5-year averages at each grid-box. The temporal stationarity of variance at each grid implies equal variance over time steps at each observing grid-box, which is commonly incorporated in detection and attribution analyses of climate change (e.g., Hannart, 2016). Both M1 and M2 estimates based on linear and nonlinear shrinkage, respectively, were obtained for comparison. Pooled estimation of the variance at each grid-box was considered in each of the shrinkage estimation to enforce the stationary, grid-box specific variance.

Figure 3 summarizes the GTLS estimates of the scaling factors $\hat{\beta}_{1}$ and $\hat{\beta}_{2}$ for the ANT and NAT forcings, respectively. The estimates from pre-whitening weight matrix $\Sigma_{\mathrm{LS}}$ and $\Sigma_{\mathrm{MV}}$ are denoted again as M1 and M2, respectively. The 95% confidence intervals were obtained with the calibration approach of Li et al. (2019). The point estimates from M1 and M2 are similar in all the analysis. The confidence intervals from the M2 method are generally shorter than those from the M1 method in the analyses both at continental and subcontinental scale. More obvious reduction in the confidence interval lengths is observed at the subcontinental scales, e.g., the ANT scaling factor in EA/NA/SCA and the NAT scaling factor in NA/WNA/SCA. This may be explained by that signals at subcontinental scale are weaker and that the error covariance matrix has non-smooth eigenvalues that form some clustering patterns due to weak temporal dependence, as suggested by the simulation study. Although the detection and attribution conclusions based on the confidence intervals remain the same in most cases, the shortened confidence intervals means reduced uncertainty in the estimate of the attributable warming (Jones et al., 2013) and other quantities based on detection and attribution analysis, such as future climate projection and transient climate sensitivity (Li et al., 2019).

Optimal fingerprinting as the most commonly used method for detection and attribution analyses of climate change has great impact in climate research. Such analyses are the basis for projecting observationally constrained future climate (e.g., Jones et al., 2016) and estimating important properties of the climate system such as climate sensitivity (e.g., Schurer et al., 2018). The original optimality of optimal fingerprinting, which minimizes the uncertainty in the resulting scaling factor estimator, is no longer valid under realistic settings where $\Sigma$ is not known but estimated. Our method constructs a weight matrix by inverting a nonlinear shrinkage estimator of $\Sigma$, which directly minimizes the variation of the resulting scaling factor estimator. This method is more efficient than the current RF practice (Ribes et al., 2013) as evident from the simulation study. Therefore, the lost optimality in fingerprinting is restored to a good extent for practical purposes, which helps to reduce the uncertainty in important quantities such as attributable warming and climate sensitivity.

There are open questions that we have not addressed. It is of interest to further investigate how the asymptotic results under $N,n\to\infty$ and $N/n\to c$ can guide the RF practice. The temporal and spatial resolution that controls $N$ can be tuned in RF practice, which may lead to different efficiencies in inferences and, hence, different results in detection and attribution. Is there an optimal temporal/spatial resolution to obtain the most reliable result? Goodness-of-fit check is an important element in detection and attribution analyses. The classic approach to check the weighted sum of squared residuals against a chi-squared distribution under the assumption of known $\Sigma$ is not valid when $\Sigma$ has to be estimated. Can a test be designed, possibly based on parametric bootstrap, to take into account of the uncertainty in regularized estimation of $\Sigma$? These questions merit future research.

JY’s research was partially supported by NSF grant DMS 1521730. KC’s research was partially supported by NSF grant IIS 1718798.

Supplementary Materials

Consider the singular value decomposition of $N\times p$ design matrix $X$:

where $U$ is a $N\times N$ orthogonal matrix, $D$ is a $N\times p$ matrix with $p$ non-negative singular values on the diagonal, and $V$ is a $p\times p$ orthogonal matrix.

Let $X^{*}=UD=XV=(X_{1}^{*},\ldots,X_{p}^{*})$ and $\beta^{*}=V^{\top}\beta=(\beta^{*}_{1},\ldots,\beta^{*}_{p})$. The columns of $X^{*}$ are orthogonal. The linear regression can be expressed as

The estimator of $\beta^{*}$ is

and the corresponding covariance matrix is

The orthogonality of $V$ ensures that $\operatorname{Tr}(\operatorname{\mathbb{V}}(\hat{\beta}))=\operatorname{Tr}(\operatorname{\mathbb{V}}(\hat{\beta}^{*})$ and $\det(\operatorname{\mathbb{V}}(\hat{\beta}))=\det(\operatorname{\mathbb{V}}(\hat{\beta}^{*})$. Therefore, we only need to consider the orthogonal design in the regression model. In other words, we may assume that the columns of the design matrix have columns such that $X_{i}^{\top}X_{j}=0$ for any $i\neq j$.

It is critically important to estimate the covariance matrix precisely as the covariance matrix of $\hat{\beta}$ depends on $\hat{\Sigma}$. In the estimation of the covariance matrix, an important question is how to quantify the quality of the covariance matrix estimator. We use loss functions to measure the quality of the covariance matrix estimator. For example, one loss function is the Frobenius norm of the bias of the estimator

and another loss function is Stein’s loss

which is the Kullback–Leibler divergence under the normal assumption.