Abstract

Summary To generate a climate data set, temperature data collected at the Earth's surface must be adjusted to remove non-climatic effects such as urbanization and measurement discontinuities. Some studies have shown that the post-1980 spatial pattern of temperature trends over land in prominent climate data sets is strongly correlated with the spatial pattern of socioeconomic development, implying that the adjustments are inadequate, leaving a residual warm bias. This evidence has been disputed on three grounds: spatial autocorrelation of the temperature field undermines significance of test results; counterfactual experiments using model-generated data suggest such correlations have an innocuous interpretation; and different satellite covariates yield unstable results. Somewhat surprisingly, these claims have not been put into a coherent framework for the purpose of statistical testing. We combine economic and climatological data sets from various teams with trend estimates from global climate models and we use spatial regressions to test the competing hypotheses. Overall we find that the evidence for contamination of climatic data is robust across numerous data sets, it is not undermined by controlling for spatial autocorrelation, and the patterns are not explained by climate models. Consequently we conclude that important data products used for the analysis of climate change over global land surfaces may be contaminated with socioeconomic patterns related to urbanization and other socioeconomic processes. Research Supported by Social Sciences and Humanities Research Council of Canada Grant Number 430002 2 Key words: global warming, data quality, spatial autocorrelation, economic activity DISCLAIMER: The authors declare that this is an original work and is not under consideration for publication elsewhere. 3 Socioeconomic Patterns in Climate Data 1 Introduction Background Temperatures are monitored at thousands of sites around the world, for a variety of purposes: weather analysis and forecasting, tracking aviation conditions at airports, assessing the risk of smog formation in urban areas, etc. In most cases, the actual temperature recorded at a specific site, such as an airport, is the measure sought. But for the purpose of measuring long term climatic trends, such as for use in studies of global warming, raw temperature data are not appropriate. This is because, over time, local mean temperatures can change for reasons unrelated to the climate. If we consider a monitoring site in a region which is known to have experienced no climatic warming, there might still be a trend in the temperature record if a nearby city expanded and encompassed the site, or the type of thermometer changed and the equipment was not carefully calibrated, or the time of day at which readings are taken was changed, etc. These local changes (called "inhomogeneities") must be filtered out to reveal, if possible, the underlying trends attributable solely to long-term climatic changes. The ideal monitoring site would be one which has undergone no changes to the local land surface since the start of modern temperature monitoring (roughly, since the late 1800s), and which has not experienced any local air pollution, which can cause small but measurable local cooling or warming, and which has had the same equipment and recording methodology across time. Because there are few such sites in the world, long temperature records are typically taken from sites that fail one or more of these conditions. Hence the records must be modified in some way to reveal the climatic trend. The dangers of using raw temperature data for climate analysis has long been noted. In 1953, commenting on the growing use of weather records for measuring climate change, J. Mitchell Jr. 1.2 Climate Data Adjustments The Climatic Research Unit (CRU) at the University of East Anglia disseminates the world's most widely-used surface temperature data sets, including those used in reports of the influential 4 Intergovernmental Panel on Climate Change (IPCC, Brohan et al. [3] assumes any inhomogeneity uncertainties are symmetric around zero (p. 6). [3] Section 2.3.3 states that to properly adjust the data for urbanization bias would require a global comparison of urban versus rural records, but classifying records in this way is not possible since "no such complete meta-data are available" (p. 11). The authors instead invoke the assumption that the bias is no larger than 0.006 degrees per century. [8] likewise offers little information about the data adjustments. They discuss combining multiple site records into a single series, but do not discuss removing nonclimatic contamination. Moreover, the article points out (page 208) that it is difficult to say what homogeneity adjustments have been applied to the raw data since the original sources do not always include this information. The two reports cited ("Jones et al. 1985("Jones et al. , 1986c are technical reports submitted to the US Department of Energy some 25 years ago. They only cover data sets ending in the early 1980s, whereas the data now under dispute is the post-1979 interval. Even if the adjustments were adequate in the pre-1980 interval it is likely impossible to have estimated empirical adjustments in the early 1980s that would apply to changes in socioeconomic patterns that did not occur until the 1990s and after. In sum, the CRU cautions that its unadjusted temperature data products (TS) are inappropriate for climatic analysis, and refers users to the CRUTEM products. Yet the accompanying documentation does not appear to explain the adjustments made or the grounds for claiming the CRUTEM products are reliable for climate research purposes. Neverthetheless, the assumption that the adjustments are adequate is widely held. For example, Jun et al. 5 Inhomogeneities in the data arise mainly due to changes in instruments, exposure, station location (elevation, position), ship height, observation time, urbanization effects, and the method used to calculate averages. However, these effects are all well understood and taken into account in the construction of the data set. In its 4 th Assessment Report, as in the previous three, the IPCC also claimed that their data are adjusted to remove non-climatic contamination. This forms an essential assumption behind all its key conclusions. Global temperature trends were presented in "Urban heat island effects are real but local, and have a negligible influence (less than 0.006°C per decade over land and zero over the oceans) on these values." The 0.006°C is referenced back to Systematic instrumental errors, such as changes in measurement practices or urbanisation, could be more important, especially earlier in the record (Chapter 3), although these errors are calculated to be relatively small at large spatial scales. Urbanisation effects appear to have negligible effects on continental and hemispheric average temperatures (Chapter 3). The citation to IPCC [6] Chapter 3 is uninformative. That chapter does not describe the data adjustments and only briefly mentions two studies (McKitrick and Michaels [15] and de Laat and Maurellis 1.3 Evidence of Data Problems Schmidt Regarding the first point, it is true that numerous different data sets show an upward trend after 1980. But the surface data show a relatively large trend compared to data collected from satellites On the second claim, we do find that use of RSS data rather than UAH data weakens the MM07 coefficients, although removal of a small number of outliers from the data set largely eliminates this 7 contrast. S09 did not present joint significance tests on which the core conclusions were based, and using RSS data these still uphold the MM07 findings, albeit at reduced significance. On the third point, S09 reported significant socioeconomic coefficients in a regression using GISS-E data. However, we show that the significance of individual coefficients disappears when the residuals are treated for SAC, something not done in the S09 analysis. In addition, the coefficients estimated on GISS-E data, as well as those estimated on the ensemble means of a much larger suite of climate models, are of opposite signs and different magnitudes compared to those estimated on observations. This provides further evidence against the view that the socioeconomic correlations are spurious, since the coefficient pattern on observed data is significantly different from that on data generated by climate models operating on the assumption that local socioeconomic process do not influence surface trends. An additional piece of evidence comes from applying the filtering methodology of MM07 to the GISS-E data. The methodology uses the regression coefficients from the socioeconomic variables to estimate the trend distribution after removing the estimated non-climatic biases in the temperature data. On observational data this reduces the mean warming trend by between one-third and one-half, but it does not affect the mean surface trend in the model-generated data. Again this is consistent with the view that the observations contain a spatial contamination pattern not present in, or predicted by, the climate models. Finally, we look at the differences between observed surface trends and the predicted values from an ensemble mean of a large suite of GCM's. If the models explain the observations, and if the observations have been filtered to remove socioeconomic influences, these trend differences should be independent of the socioeconomic variables. But we find that the differences are highly correlated with the socioeconomic indicators, and the coefficients are very close to those estimated on the observed trends themselves. This strengthens the argument that the socioeconomic pattern in the data is not accounted for by the processes in the climate models. Taken together we find significant evidence against the view that surface climate data are free of biases due to socioeconomic development and other inhomogeneities. Instead, measures of socioeconomic influence appear to be an essential component of a well-specified model of the spatial trend pattern in climate data over land. The coefficient pattern on observational data differs in both sign and magnitude from that predicted by climate models as a response to natural oscillations and anthropogenic (greenhouse) forcing. Hence we consider the standard interpretation of climatic data to be untenable. In the next section we explain the data sets used throughout this paper. In Section 2 we model spatial autocorrelation and give detailed results for the data configurations of interest. In Section 3 we explore the mismatch between the regression results from model-generated and observed data. Section 4 presents further specification tests and Section 5 concludes. 1.3 Data sets Most data sets used herein are taken from MM07 and S09. 1 Readers should consult both these papers for detailed explanations; only a brief summary will be presented herein. For ease of notation we will drop the gridcell subscript i when doing so does not create ambiguity. Equation (1) explains the spatial pattern of temperature trends using three main variable groups: temperature trends in the lower tropospheric layer about 5 km above the surface, fixed geographical factors and socioeconomic variables. The geographical variables include latitude, coastal proximity, mean air pressure, etc. The socioeconomic variables measure factors that influence data quality, land use change, etc. The standard interpretation of climate data is that their effects have been filtered out of climatic data products like CRUTEM. Summary statistics for the data are in The S09 data set comprises surface and tropospheric grid cell trends like those in MM07, except the surface trends are from later CRU compilations and the tropospheric trends are from RSS The tropospheric data used in MM07 and S09 were at a 2.5x2.5 degree level, one-fourth of the 5x5 CRU surface grid size, so the top-right tropospheric cell was used. For some of our calculations herein we retain the 2.5 degree scale aloft where our intent is to replicate earlier results. Otherwise, in order to reconcile the spatial scales between surface and tropospheric gridcells we develop matched 5x5 9 grid cells. We denote the data series in which four tropospheric cells have been combined to yield a 5x5 grid cell as UAH4 and RSS4. S09 also provided synthetic trends from GISS-E. For a description of this model see CCSP ([4] Sct. 2.5.3) and The average GISS-E land surface trend is 0. With a vector of trend terms on both the left-(surface) and right-hand (troposphere) side there are 24 possible data combinations: CRU, CRU2v, CRU3v at the surface, UAH, UAH4, RSS and RSS4 aloft, and GISS and the all-GCM averages. Additionally results can be run with no spatial autocorrelation terms, or with corrections on either the error or lagged dependent variables, making 72 possible model configurations. Since there are many common results across different specifications we will only report those central to our argument, but other results are available on request. For instance, since CRU2v was not used in MM07 and has been superseded by CRU3v we will not report CRU2v results, and we will generally use UAH4 and RSS4 rather than UAH and RSS. 2 Robustness across multiple data sets: the non-SAC case 2.1 Observational data sets We begin this section by looking at whether the MM07 results were unique to the particular data configuration used therein. Estimations throughout this paper were done using STATA version 8.2, running on a Dell Studio laptop with an Intel i7 quad core 64-bit processor. All data and code are available in the supplementary information. The socioeconomic coefficient estimates are quite similar across the observational columns (2 though 9). Use of CRU3v data does not yield much of a change in socioeconomic coefficients compared to CRU2v and CRU. Use of RSS data rather than UAH data yields smaller and less significant coefficients, though for some reason leaves a greater component to be explained by the latitude variable. Use of reconciled gridcell sizes also yields smaller and less significant coefficients compared to the 2.5x2.5 tropospheric grids. Across all these specifications the coefficient sizes and signs remain comparable and the socioeconomic effects taken as a group P(g-c=0) remain jointly significant. The largest drop in the coefficient magnitudes is associated with using the CRU3v-RSS4 pair, yet the effect is 10 apparently due to a relatively small number of outliers. The C3/RSS4x column repeats the CRU3v/RSS4 results with outliers removed. The Ordinary Least Squares "hat matrix" was evaluated, and an observation was flagged as an outlier if it exceeded twice the mean diagonal element of the hat matrix (see 2.2 Model-generated data sets S09 hypothesized that these effects arise from a lucky match between the spatial pattern of socioeconomic activity and the spatial pattern of enhanced natural and greenhouse forcing of the climate. Since the GCM does not contain a socioeconomic component, if, upon using GISSES and GISSET in place of observations in the MM07 regression model, significant coefficients of the same approximate size and sign emerge on the socioeconomic variables, then correlations such as those in MM07 could be dismissed as coincidental. It is worth quoting the argument in S09 directly to make this point clear. "There is a relatively easy way to assess whether there is any true significance to these correlations. We can take fully consistent model simulations for the same period and calculate the distribution of the analogous correlations. Those simulations contain no unaccounted-for processes (by definition!) but plenty of internal variability, locally important forcings and spatial correlation. If the distribution encompasses the observed correlations, then the null hypothesis (that there is no contamination) cannot be rejected. (S09, p. 2, emphasis added) However, as is shown in columns 10 and 11 of Table 2, the regression coefficients estimated on the data generated by the GISS-E ensemble and the all-GCM ensemble are quite different from those estimated on observational data. Coefficients in With regard to the quoted paragraph, the distributions of the coefficients estimated on GCM data do not encompass the coefficients from either the MM07 data set or any other observational grouping in 3 Spatial autocorrelation of the trend field 3.1 Testing framework The analysis in the previous section ignores the issue of spatial dependence, so we must now reconsider the models taking it into account. Both S09 and Benestadt [2] point out, correctly, that the surface temperature field is spatially autocorrelated, and argue that this can, in principle, bias the inferences from regressions on the spatial trend field. They both concluded on this basis that the results in MM07 and MM04 were unreliable. However neither one formulated the argument as a testable hypothesis, though S09 presented variograms of the dependent variable and some independent variables from MM07. It is insufficient to observe autocorrelation in a dependent variable and conclude that the inferences from a regression model are therefore biased. An additional step in the argument is required, namely a test showing that the regression residuals also exhibit SAC. As we will show, they do when model-generated data are used, as in S09, but they do not when observational data are used, as in MM07. The contrast is important. Inferences concerning the coefficients in a regression model are based on the statistical properties of the residuals, not the dependent variable. Thus, the absence of SAC in the residuals of a regression model in which the dependent variable is spatially autocorrelated is evidence in support of the specification, i.e. that the right hand side variables do have explanatory power. We test for residual spatial dependence as follows. The regression model (1) can be rewritten in matrix notation as where T is a 440x1 vector of temperature trends in each of 440 surface grid cells, X is a 440xk matrix of climatic and socioeconomic covariates, b is a kx1 vector of least-squares slope coefficients and u is a 440x1 residual vector. Spatial autocorrelation in the residual vector can be modeled using where λ is the autocorrelation coefficient, W is a symmetric n n × matrix of weights that measure the influence of each location on the other, and e is a vector of homoskedastic Gaussian disturbances, (Pisati [20]). The rows of W are standardized to sum to one. n equals 440 except in some regressions where grid cells are missing, as noted below. A test of 0 : measures whether the error term in (1) is spatially independent. As argued in S09, it is likely the dependent variable is spatially autocorrelated. Anselin et al. Autocorrelation can affect either the dependent variable or the residuals or both. Florax et al. [6] discuss a sequential testing and estimation regime for models of the form combining equations Here φ is the spatial lag and λ is the spatial error term. If neither one is significant then a least-squares model with no spatial autocorrelation term can be used. If only one of them is significant, then (5) should be used retaining only the significant lag term. If both are significant then the one that has a higher significance level should be chosen, implying either the spatial lag model (4) or the spatial error model (2, 3). Hypothesis tests, and any subsequent parameter estimations, are conditional on the assumed form of the spatial weights matrix W in (3). Denote the great circle distance between the grid cell centers from which observation i and observation j are drawn as ij g . The weighting function is µ − ij g where µ determines the rate at which the relative influence of one cell on adjacent cells declines. The function was estimated using a grid search to find the maximum-likelihood value. For observational data groupings the likelihood function was maximized at values of µ between 2.5 and 2.7, whereas for the GISS-E and all-GCM model-generated data, the likelihood function was maximized at µ values of 3.2 or 3.0 respectively. The optimal exponent values are listed in Spatial Autocorrelation Testing Results The Stata command 'spatwmat' was used to generate the row-standardized weights and eigenvalues, then the command 'spatdiag' was used to generate the test scores. The second column (exponent) reports the maximum likelihood value of µ . The next column reports the robust LM score for a null hypothesis of no spatial dependence on the dependent variable, while the final column reports the robust LM score for the test on the residuals. In each column the corresponding p value is shown in parentheses. The first row of results refers to the original configuration in MM07: the CRU gridded trends regressed on the UAH tropospheric trends and the rest of the MM07 model variables in Equation (1). In 13 this case the robust LM score is significant for both the dependent variable and the residuals, but much more so for the dependent variable, indicating that a spatial lag model is appropriate. The second and third rows show the test scores using CRU3v surface data and either UAH4 or RSS4 tropospheric series. In both cases the dependent variable lag is significant while the residual lag term is insignificant. Again this indicates that the spatial lag model is appropriate, and also indicates that the regression model is well-specified in the sense that the SAC is removed from the error terms. The next two rows report the results after substituting in model-generated data on both sides of Equation (1). The results change in an interesting way. The no-SAC test in the residuals is now strongly rejected, and the ranking of the significance also changes, so that the appropriate estimation model is now the spatial error form rather than the spatial lag. This is suggestive of a deficiency in the specification of (1), whereby the terms on the right hand side fail to explain the spatial lag pattern for model-generated data the way they do when observational data are used. The final two rows refer to a test applied in the next section and will be discussed below. To summarize, in regressions using observational temperature data, the right-hand side terms in equation 3.3 Regression results with SAC controls The pattern observed in Note that the coefficients in Another interesting contrast concerns the tropospheric coefficient (first row). The usual expectation from climate modeling studies is that the trend aloft will be slightly stronger than that at the surface, but as shown in 3.4 Do GCM's predict the observed temperature-industrialization correlation pattern? Returning to Schmidt's [22] argument that the distribution around coefficients estimated on modelgenerated needs to encompass those estimated on observed data, The spatial autocorrelation results in the previous section, and the mismatch between coefficients estimated on observed and model-generated data, point to the likely existence of a non-climatic contamination pattern in the observed surface trend data. Further evidence on this point is obtained by repeating the filtering experiment of MM07 on the GISS-E and all-GCM data. The MM07 method sets the contaminating influences to zero and assumes each country has equivalent measurement resources as the United States. Since the climate model does not contain any contaminating processes or qualitycontrol variations we should not expect much difference between the raw data from the models and that obtained by applying the MM07 method for removing socioeconomic effects. 15 Finally, to investigate the ability of the all-GCM ensemble to explain the surface trend pattern we took the differences between the observed surface trends (CRU3v) and the all-GCM mean surface trends and examined if the differences can be explained by the socioeconomic variables in equation 4 Further specification tests MM07 presented a series of specification tests using the UAH data. We repeated all of them using CRU3v and RSS4. The results closely follow those reported in MM07. Details are available on request, and can be summarized as follows. • The RESET test does not reject a null hypothesis of no un-modeled residual nonlinearity (P = 0.241). • The Hausman test does not reject a null hypothesis of no endogeneity bias (P = 0.999). • The outlier test (described in Section 3) flags 26 observations as influential. When these are removed the individual and joint socioeconomic coefficient tests become more significant, yet we do not reject a null hypothesis that the coefficient vectors with and without the outliers are equivalent (P = 0.278). • Coefficient results are individually and jointly significant in rich countries but not poor countries, and in economies with growing but not declining incomes. • After removing a randomly-selected third of the data set and re-estimating the model, the prediction of the withheld sample scatters along a 45-degree line with the observed values. In 500 repetitions, a regression of the predicted and observed values has a constant of 0.011 and a slope of 0.961, and a test of a perfect fit (constant = 0, slope = 1) obtains an average P value of 0.407, i.e. does not reject on average. The SAC tests reveal that one of the specification tests in MM07 was done incorrectly. In Section 4.6 of MM07, an alternative estimation is presented in which the surface trends were replaced by the UAH-derived lower tropospheric trends. Had the socioeconomic coefficients retained their size and significance it would provide evidence that the surface results might be spurious. Some variables, such as growth in coal consumption (c) and population growth (p) can yield regional changes that affect the 16 lower troposphere. One mechanism is the changing atmospheric aerosol load (Li et al. [13]); another is induced changes in regional precipitation downwind from major cities (Shepherd et al. [24]). However, the tropospheric trends vector is spatially autocorrelated, so the regression equation needs to be augmented with an SAC correction, which was not done in MM07. As in MM07, x changes sign and in one case (RSS4 with a spatial lag model) acquires significance, indicating that in this regression it is likely acting as a proxy for some unrelated spatial pattern. In the full regression x is never significant (see The results in 5 Conclusions We have examined the question of whether spatial trend patterns in surface temperature data can be explained in part by non-climatic, socioeconomic processes of the kind that are supposed to have been filtered out of the gridded data products. We have shown that a coefficient pattern connecting temperature trends to indicators of industrialization is robust across a wide range of data configurations in the surface and lower troposphere, but is absent in climate mode-generated data. The failure to 17 reproduce this pattern in models indicates that it is not a natural feature of the climate system nor a response to greenhouse gas-induced forcing. One strand of argument against earlier findings on this issue was that spatial autocorrelation of the temperature field reduces the effective number of degrees of freedom, biasing significance calculations. We have estimated robust SAC test statistics and have shown that while the trend field is spatially autocorrelated, SAC is not found in the residuals when updated observations are used, but is strongly present when model-generated data are used. This again points to a qualitative difference between observations and models that may be explained by the patterns of socioeconomic development. After re-estimating our results with the appropriate corrections for spatial dependence we find the socioeconomic coefficients remain significant in observations, but, again, disappear in model-generated data. Also, we find no overlap between the implied distribution of coefficients from estimations on model-generated data and coefficients estimated on observed data, refuting a key conjecture in Schmidt The data set presented in MM07 includes trends up to the end of 2002, and includes coarse resolution of some socioeconomic variables at the national level. Further investigation of the potential surface climatic data problems we have identified herein could involve a reconstruction of the MM07 data base using updated socioeconomic and climatic variables, use of cross-sectional time series (panel) regression rather than trend fields, and use of regional, rather than national, socioeconomic data where available. Acknowledgments We thank Chad Herman for assistance in obtaining the PCMDI archive trends. Financial support from the Social Sciences and Humanities Research Council of Canada is gratefully acknowledged.