A New Method of Significance Testing for Correlation-Coefficient Fields and Its Application

Considering a geographical region D covered by a lattice (site) network composed of m grid lattices (sites) ($ s=1,\dots ,m $), the field r composed of the correlation coefficient $ {r}_{s},s=1,\dots ,m $ in the lattice (site) network becomes the distribution field of the correlation coefficient (still noted as r). r can be the distribution field of the correlation coefficient between the time series of a factor or factors $ {x}_{q},q=1,\dots ,n $ (noted as x) and the time series of the previous factor field $ {y}_{q},q=1,\dots ,n $ (noted as $ {\boldsymbol{Y}}_{m\times n} $) in D, which can also be the local correlation-coefficient distribution field of time series $ {\boldsymbol{X}}_{m\times n} $ and $ {\boldsymbol{Y}}_{m\times n} $ of the two factor fields on D. In the analysis of atmospheric circulation and climate anomalies, r is widely used for analyzing the correlation between atmospheric circulation and contemporaneous or lagged climatic anomalies (Feng, 1978; Chen and Song, 2000; Gong et al., 2002; Hua and Ma, 2009; Ma et al., 2009; Wu et al., 2011; Huang and Wang, 2012; Xu and Wu, 2012; Liu et al., 2015).

In order to infer whether there is a correlation between the parent population from r of sample x(X)~Y, r needs to be tested for significance. Livezey and Chen (1983) developed an r significance test method (known as the Livezey method), which assumes that k of m correlation coefficients $ {r}_{s},s=1,\dots ,m $ constituting r pass the significance test of reliability $ \alpha $, and k is the number of significant-correlation lattices (stations). Assuming that the matrix of x(X) and Y are not correlated, the critical value $ {k}_{\alpha } $ of k can be determined by m and $ \alpha $ according to Bernoulli’s theorem. When $ {k\geqslant k}_{\alpha } $, the correlation-coefficient distribution field passes the significance test of reliability $ \alpha $.

However, the Livezey method implies that $ {k}_{\alpha } $ is based on two assumptions: (1) a uniform space of lattice (site) points; and (2) the matrix space of the local sequence X and Y has no correlation. However, for all planetary-scale two-dimensional lattice (site) networks, neither of these assumptions is true. Therefore, the Livezey method for r significance testing has defects.

This paper is based on research results regarding the correction of spatial homogenization of climatic data (Chung and Nigam, 1999; Luo et al., 2011) and introduces a new statistic $ \widehat{D} $, the significant-correlation area of r, which is used to replace the number of total stations k of the significant correlation of r in the Livezey method and eliminate the influence of the lattice (site) network being uneven in reality. Furthermore, the Monte Carlo method to test the significance of $ \widehat{D} $ is replaced with an empirical Monte Carlo method, eliminating the influence of the X(Y) parent space correlation, so as to obtain a reasonable and operable new method of r significance testing.

The correlation-coefficient field r is given in the field area $ \mathrm{\Omega } $, the significant-correlation in the region $ \widehat{\mathrm{\Omega }}(\widehat{\mathrm{\Omega }}\in \mathrm{\Omega }) $ is $ (\left|r\right|\geqslant {r}_{\alpha }) $, and $ D\left(\widehat{D}\right) $ is the area of $ \mathrm{\Omega }\left(\widehat{\mathrm{\Omega }}\right) $. $ \widehat{\mathrm{\Omega }} $ is the significant-correlation region of r, and $ \widehat{D} $ is the significant-correlation area of r.

Fields of r used in atmospheric circulation or climate analysis are obtained from the discretized data given on the lattice (site) network, and the most commonly used lattice network is a rectangular network with uniform spacing on plane $ \lambda -\phi $ ($ \lambda $ is the longitude, ϕ is the latitude). Figure 1a is a uniform grid distance rectangular lattice network in the Northern Hemisphere, and the network area is $ D=2\pi {a}^{2} $(a is the radius). Figure 1b shows a 160-station network in China with an area of about 9.6 million km². Obviously, the distribution of lattice points on the network is uneven. In order to calculate $ \widehat{D} $, it is first necessary to give the measurement method of the lattice network inhomogeneity. Luo et al. (2011) introduced the element area $ {d}_{s} $ of the s-th lattice (site), which is the part of the domain area D represented by the s-th lattice (site). The larger (smaller) $ {d}_{s} $ is, the sparser (denser) the dot distribution is near point s. $ {d}_{s} $ is the basis for calculating the area of the significant-correlation area, and its algorithm should be given first.

Figure 1.  Example of lattice (site) network. (a) Uniform grid distance rectangular lattice network in the Northern Hemisphere（$ \mathrm{\Delta }\lambda ,\Delta \phi =\pi /18 $）；(b) 160-station network in China (Wang et al., 2011)

(1) Calculation of element area $ {d}_{s} $

Chung and Nigam (1999) pointed out that the element area of the lattice point located at latitude $ \varphi $ on the lattice network in Fig. 1a was proportional to $\mathrm{c}\mathrm{o}\mathrm{s}\varphi$. Figure 1a shows the uniform spacing on the $ \lambda -\varphi $ plane of the Northern Hemisphere ($\Delta \lambda \;\mathrm{a}\mathrm{n}\mathrm{d}\;\Delta \varphi$ are constant). The zonal lattice number is $ u(\Delta \lambda =2\pi /u) $, and the meridional lattice number is $(v+1\left)\right(\Delta \varphi =\pi /2v)$. The two-dimensional lattice point ordinal number $ \left(i,j\right)=\left(\mathrm{0,0}\right) $ is taken on the intersection point $ \left(\lambda ,\varphi \right)=\left(\mathrm{0,0}\right) $ between the prime meridian and the equator, and by using Fig. 2, the corresponding element area can be easily written (i is the longitude coordinates, j is the latitude coordinates):

Figure 2.  A schematic diagram of the calculation of the area of s lattice element. The shaded area of the Arctic is only $ {d}_{s}/4 $.

where the units of $ \Delta \lambda $ and $ \Delta \varphi $ are radians (rad), and the units of d(i,j) and ${d}\left(i,j\right)$ are the radius of the Earth squared. Because $ \Delta \lambda $ is constant, $ d\left(i,j\right) $ only changes with j. Using the relation S~(i, j),

where $ d\left(i,j\right) $ can be easily converted to $ {d}_{s} $. Note that the North Pole is processed into $ u $ lattice points, ${d}_{s}=\overline {1,m}$, $ m=u\left(v+1\right) $, $\Delta \lambda =\Delta \varphi =\pi /$18, $ u $=36, ${v}$=9, and $m=u (v+1)$=360.

Luo et al. (2011) presented the calculation scheme of S station area $ {d}_{s} $ on a Chinese 160-station network (Fig. 1b), which is briefly introduced here with the help of Fig. 3. In the figure, S represents the test station located at $ {(\lambda }_{s},{\phi }_{s}) $. The shaded area in Fig. 3 is the spherical part of the spherical crown with S as the pole, and the area $ {S}_{0} $ is called the area of the crown [(${S}_{0}=2\pi (1-\cos\theta ){a}^{2}$, where $ a $ is the radius of the Earth and $ \theta ={\angle SOS}^{*} $ is the center angle corresponding to the pole of the spherical cap to the edge point)]. For the selected $ {S}_{0} $, S station crown area $ {\mathrm{\Omega }}_{s} $ is determined, and its area is $ {S}_{0} $. The number of stations in the statistics is determined to be $ {\widehat{m}}_{s} $, and then the area $ {D}_{s} $ that belongs to the coverage area of station network $ {\mathrm{\Omega }}_{s} $ is calculated. Thus, the area $ {{d}}_{s} $ of station element S and the density of the station network are obtained, which are defined as:

Figure 3.  A schematic diagram of the calculation of the area of s station. The shaded area is $ {\mathrm{\Omega }}_{s} $.

Figure 4 is the diagram representing the calculation of $ {{d}}_{s} $ and ${{m}}_{s}$ for the Lanzhou (inland), Altay (land border), and Shanghai (ocean boundary) stations. For the inland station, ${{D}}_{s}$=1, and for the boundary station, $ {{D}}_{s} $<1. Table 1 shows ${{d}}_{s}$, ${{m}}_{s}$, and the intermediate calculation results of $ {\widehat{m}}_{s} $ and ${{D}}_{s}$ for the three stations.

Figure 4.  A schematic diagram of the structure of three types of stations. The elliptical line is the $ {\mathrm{\Omega }}_{s} $ boundary, the shallow shadow area is the ocean, the deep shadow area is the land, and the white curve is the land boundary/sea land demarcation line. (a) Lanzhou (inland) station, (b) Aletai (land border) station, and (c) Shanghai (ocean boundary) station.

(2) Calculate of the significant-correlation area D

Figure 1 shows the correlation-coefficient field r in the lattice (site) network composed of m phase relationship values $ {{r}}_{s} $. The significance test of reliability $ \alpha $ is carried out for each ${{r}}_{s}$, where k-th $ {r}_{s} $ are significant. We introduce 0 and 1 functions,

in which we take 1 through the test lattice (site) point $ {p}_{s} $; otherwise, we take 0. Thus, the formula for calculating the area of the significant-correlation area in the correlation-coefficient field is obtained as

where $ \widehat{D} $ is the new statistic used in this paper for the significance test of correlation-coefficient field.

The empirical Monte Carlo (EMC) method for significance testing a correlation-coefficient field r on a lattice network (Fig. 1a) and site network (Fig. 1b) is the same as using statistic $ \widehat{D} $ (area of significant-correlation area) to replace k (number of significant-correlation points) in the Livezey method. Here, the EMC method is given by taking the significance test of the correlation-coefficient field r of a 160-station network in China (Fig. 1b) as an example. Meanwhile, the Livezey method is introduced in order to compare it (Livezey and Chen, 1983; Wang et al., 2007, 2011) with the test results.

r is the correlation-coefficient distribution fields of the time series x and field series Y or the local correlation-coefficient distribution fields of the sequence of fields X and Y in the same grid (site) network at the same time (or delayed) in time series $ q=1,\dots ,n $. x, X, and Y are expressed as

where $ {x}_{q} $ is the t-th time value of sequence x, and $ {x}_{q} $ and $ {y}_{q} $ are the $ t $-th time fields of field orders X and Y. The test steps of reliability $ \alpha $ of r using the EMC method are as follows:

(1) According to the method described in section 2, the significant-correlation area $ \widehat{D} $ of sample correlation-coefficient field r is found.

(2) The sample sequence $ {x}_{q},q=1,\dots ,n $ or $ {\mathit{x}}_{q}, $ $q=1,\dots ,n $ is randomly selected L times by a random function, where L = 1000 or more. Then, the correlation-coefficient field of $ {\mathit{x}}_{l} $ (or $ {\mathit{X}}_{l} $), $ l=1,\dots ,L $, and Y are found, and the following is obtained:

(3) The significant area $ {\widehat{D}}_{l},l=1,\dots ,L $ of $ {r}_{l},l=1,\dots ,L $ is found, and it is arranged in non-ascending order to get $ {\widehat{D}}_{l},l=1,\dots ,L $. Then, making ${l'}_{\alpha }=\alpha \times L$, ${\widehat{D}}_{{{l}^{{'}}}_{\alpha }}$ is the critical value of $ \alpha $ for the reliability of $ {\widehat{D}}_{l} $. For $ \alpha $=0.05, L=1000, ${l'}_{0.05}$=50, $ {\widehat{D}}_{50} $ is a value that the reliability of $ \widehat{D} $ is the critical value of 0.05.

(4) The significance of r is tested with statistic $ \widehat{D} $: If $ \widehat{D}\geqslant {\widehat{D}}_{{l}_{\alpha }} $, then r is significant with respect to $ \alpha $ under the reliability, and x(X) and Y are from the relevant parent; otherwise, r is not significant, and x(X) and Y come from an uncorrelated parent.

The key of the EMC method is to use the random sequence of the sample sequence x(X) to generate its simulated sequence $ {\mathit{x}}_{l}\left({\mathit{X}}_{l}\right) $, which is the basis of the word “empirical” in its name. This process makes the statistical properties of $ {\mathit{x}}_{l}\left({\mathit{X}}_{l}\right) $ closest to the sample and the statistical properties of the sample closest to its parent. For example, the $ {\mathit{r}}_{l} $ obtained by the EMC method in step (2) should be well organized. Therefore, the EMC significance test of r considers the spatial correlation of X (or Y).

For the 160-station network of China, when the reliability $ \alpha $=0.05 and L=1000, the critical value of k for the number of significantly correlated stations of in the correlogram is $ {k}_{0.05} $=13 (Wang et al., 2007). When $ {k > k}_{0.05} $, r passes the significance test, and x(X) and Y come from the correlation matrix.

According to the method of Sun et al. (2015), we obtained the intensity sequence $ {P}_{t},t=1,\dots ,70 $ (P) of five atmospheric activity centers (ACAs) during a 70-yr period (1951–2020) in the Northern Hemisphere in winter, which is x in Eq. (5). Then, the temperature and precipitation field sequences (T and R) were obtained at 160 stations in China during the same period, which is Y in Eq. (5). Figures 5 and 6 show the simultaneous correlation-coefficient field r of x-T and x-R, respectively.

Figure 5.  Correlations between ACA intensity index in the Northern Hemisphere and temperature at 160 stations in China in winter from 1951 to 2020. The shaded area is the significant-correlation area that passes the t test at the $ \alpha =0.05 $ level. (a) Iceland low. (b) North Atlantic High. (c) Mongolia High. (d) Aleutian Low. (e) North Pacific High.

Figure 6.  Correlations between ACA intensity index in the Northern Hemisphere and precipitation at 160 stations in China in winter from 1951 to 2020. The shaded area is the significant-correlation area that passes the t-test at the $ \alpha =0.05 $ level. (a) Iceland low. (b) North Atlantic High. (c) Mongolia High. (d) Aleutian Low. (e) North Pacific High.

Table 2 shows the significance test results of the EMC method with the area $ \widehat{D} $ and $ \alpha $ = 0.05 of the significance test area in Figs. 4 and 5. Two qualitative ACAs in the North Atlantic (Iceland low and North Atlantic High) and the Mongolia High intensity in eastern Eurasia are significantly correlated with the corresponding winter temperature in China. Quantitatively, the correlation between the Mongolian High and temperature is the strongest, followed by the Icelandic Low and the North Atlantic High. The correlation between intensity and China’s winter precipitation in the same period is much weaker, and only the Icelandic Low intensity is significantly correlated with China’s precipitation in the same period. The above conclusions are consistent with the size of the significant-correlation area in Figs. 5 and 6.

Table 3 shows the Livezey significance test results of the number of significantly correlated stations k and $ \alpha $=0.05 in Figs. 5 and 6. It can be seen that the test results of P-T are consistent. However, the test results of P-R vary greatly. Regarding P-R, in addition to the Icelandic Low, the Mongolian High, Aleutian Low, and North Pacific High were also significantly correlated. As can be seen from Figs. 6c–e, although the area of the significant-correlation area above them is small, they are all located in the dense area of the site in Fig. 1b. The inhomogeneity of the station network seriously affects the significance test of r.

By comparing Tables 2 and 3, it is reasonable and effective to test the significance of the correlation-coefficient field r of an example by using the EMC method and through the significant-correlation area $ \widehat{D} $.

In summary, this paper analyzes the defects of the Livezey method for significance testing of a correlation-coefficient field r, constructs a new statistic $ \widehat{D} $ (the significant-correlation area of r), and presents an EMC method for significance testing of r using $ \widehat{D} $. The Livezey method implies an assumption of the spatial independence of the correlation-coefficient field. In fact, the spatial distribution of the lattice or site network of the circulation and climate field is generally not uniform, and the field quantity on the adjacent lattice or site must have a certain correlation; otherwise, the system in the circulation and climate field cannot be given. The significance test by the EMC method can eliminate the influence of parent-space correlation. The significance test results from applying the method in this paper to the correlation chart r of the intensity index of ACAs in the Northern Hemisphere winter and the temperature and precipitation in China in the same period show that the test method in this paper is more reasonable than the Livezey method.

Acknowledgements. This work was supported by the National Key Research and Development Program of China (Grant No. 2018YFC1505602), the National Natural Science Foundation of China (Grant Nos. 41705055 and 41505088), the Project of Scientific Creation of Post-Graduates of Jiangsu (Grant No. CXZZ12_0485), the Creative Teams of Jiangsu Qinglan Project, and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). The authors are very grateful for the constructive comments of two reviewers.