Causal inference from cross-sectional earth system data with geographical convergent cross mapping

Domain: earth_science
Published: 2023-09-21
Editorial outcome: accept

Abstract
Causal inference in complex systems has been largely promoted by the proposal of some advanced temporal causation models. However, temporal models have serious limitations when time series data are not available or present insignificant variations, which causes a common challenge for earth system science. Meanwhile, there are few spatial causation models for fully exploring the rich spatial cross-sectional data in Earth systems. The generalized embedding theorem proves that observations can be combined together to construct the state space of the dynamic system, and if two variables are from the same dynamic system, they are causally linked. Inspired by this, here we show a Geographical Convergent Cross Mapping (GCCM) model for spatial causal inference with spatial cross-sectional data-based cross-mapping prediction in reconstructed state space. Three typical cases, where clearly existing causations cannot be measured through temporal models, demonstrate that GCCM could detect weak-moderate causations when the correlation is not significant. When the coupling between two variables is significant and strong, GCCM is advantageous in identifying the primary causation direction and better revealing the bidirectional asymmetric causation, overcoming the mirroring effect.

Introduction
A sufficient and precise understanding of causal associations between the target variable and influencing variables is an important component for effectively utilizing natural resources and achieving sustainable development. Natural and social experiments are reliable approaches for proving the existence of causation between the cause and effect variable. However, in Earth system sciences, it is usually not feasible to conduct experiments to explore the causal associations between the target and influencing variables at large spatiotemporal scales. In this case, investigation on observational data has become a commonly acceptable alternative to infer reliable causal associations. Traditional statistics is a mainstream approach to analyze observational data, yet it denied a hypothesis of causation until the 1980s. In the past decade, causal inference has become a topic of interest in many disciplines and a diversity of approaches has been proposed, implemented and adapted for better attributing major social, ecological, and economic issues.

Structural causal modeling, which originated from Wright and was developed by Pearl, and the potential outcome frameworks, which originated from Neyman and Fisher and was developed by Rubin, are the two main frameworks to identify the existence of causal relationships and estimate causal effects (the strength of the causal relationship). The former emphasizes the combination of causal graphs and statistical tools, while the latter attempts to imitate the randomized controlled trials (RCT) by adjusting the variables. Causal network learning algorithms developed from the structural causal modeling framework, such as PCMCI, can automatically build the causal networks and estimate casual effects. The potential outcome framework was proposed to embed the instrumental variables to eliminate the confounders and has been widely used in social and economic science. These two frameworks are mainly for stochastic processes. But for earth system sciences, many systems (e.g., climate system or ecosystems) contain deterministic or dynamical interactions, which are not captured well by the assumptions of these frameworks and thus challenge their applicability. For those nonlinear complex systems, causal inference models based on the theory of state space reconstruction were developed. Amongst these models, Convergent Cross-mapping (CCM) has been well accepted and massively employed in a diversity of disciplines, including epidemic diseases, socio-economic issues, and atmospheric pollution. According to Takens’ embedding theory, CCM employs time series data to reconstruct the attractor, based on which the causation can be identified and measured through the cross-mapping prediction. In addition, other models which extract causation following a similar principle were also developed, such as cross-mapping smoothness (CMS), partial cross-mapping (PCM).

However, these models can solely work with time series data, and cannot be applied to spatial cross-sectional data. Spatial cross-sectional data, which often includes rich spatial information but has little information on temporal changes, is more easily available than time series data in Earth System Science. Spatial cross-sectional data records spatial processes and their interactions, and spatial order (or variation) are important indicators for understanding causal associations, which may be ignored by temporal models due to insufficient length of time series or insufficient variation across time series. Our previous research suggested that major temporal causation models, including CCM and Granger Causality, failed to extract NPP (Net Primary Productivity)-temperature (precipitation), a clear existing causation. This was mainly attributed to the slight inter-annual variations of temperature (precipitation). Meanwhile, NPP-climate association was effectively identified by spatial models, thanks to notable spatial variations of temperature and precipitation. Existing models, such as “spatial difference in difference” and “spatial regression discontinuity”, may in some cases measure the average causal effect of binary treatment under strong assumptions. They are based on the stable unit treatment value assumption (SUTVA) (or ignorability), which assumes that (1) the potential outcomes of different spatial units do not interfere with each other, and (2) all treatment levels are included so that the outcome should remain unchanged if the same treatment is assigned. However, due to spatial heterogeneity and spatial spillover, which mean that potential outcomes of one spatial unit can influence its spatial neighbors, these assumptions are easily violated in Earth System Sciences. Another major challenge lies in the identification of the direction of causal associations, which is largely disturbed by the contemporaneous dependence between spatial variables, known as the mirroring effect. Therefore, it remains unclear whether, and more importantly, how reliable causal associations can be inferred from spatial cross-sectional data.

Given the demand for causal inference from spatial cross-sectional data in Earth System Sciences and the limitation of existing temporal and spatial causation models, here we propose a Geographical Convergent Cross Mapping (GCCM) model by adopting dynamical systems theory and generalized embedding theory. GCCM is capable of identifying the causal associations between spatial cross-sectional variables and estimating the corresponding causal effects. By exploring nonlinear associations, GCCM outputs can be more robust than linear correlation, which may easily lead to spurious correlation, and can identify causations neglected by Linear Non-Gaussian Acyclic Model (LiNGAM), a frequently used structural causal model with enhanced assumptions. Based on three typical cases, including clearly existing causation that cannot be measured through temporal causation models, we demonstrate that GCCM can detect weak-moderate causation when the correlation is not significant. On the other hand, when the correlation between two variables is significant and strong, GCCM is disadvantageous of identifying the primary causal direction and better revealing the unidirectional asymmetric causation, overcoming the mirroring effect. GCCM opens up a new vision to infer causal associations from spatial cross-sectional data, by bridging Earth data with the methodologies of complex system, and can be a powerful tool for revealing nonlinear complex relationships in Earth systems.

Results
We illustrate the implementation and interpretation of GCCM with three typical cases, demonstrating the application GCCM to different types of spatial data and causal assoications. The first case is to infer and casual associations between the soil heavy metal and two influencing factors (the density of industry pollutants and residence) in Illinois and Indiana of US. The second case is to verify the casual associations between county-level population density and environmental factors in China. The third case consider the farmland NPP-temperature (precipitation) associations in China, which are clear existing causations but fail to be identified by mainstream temporal causation models. Correlation analysis, which is the most frequently used method for understanding causal associations from spatial cross-sectional data and LiNGAM from the structural causal modeling framework, which is commonly used to identify the direction of causations, are employed for the comparative analysis.

Causal inference with GCCM
Causal associations are important components of inner mechanism, which can be recognized by observing and analyzing the phenomena they present. The spatial distribution is an important phenomenon to extract causal associations in supplement to temporal changes. The corresponding spatial cross-sectional data record the spatial processes and their interactions, and the spatial difference (order) provides a valuable reference for understanding causal associations. There have been some famous examples of inferring causation from the spatial perspective. For instance, Charles Darwin developed the theory of evolution by noticing the spatial difference of animals among Galapagos Islands, while Vasili Vasilievich Dokuchaev proposed the soil formation theory (Pedogenesis) according to the spatial difference of soil along the latitude. Thoses successes are the fruit of human wisdoms applied to spatial cross-sectional data, but only with formalized mathematical methods, can the inferring framework become easily understandable and transferable for researchers from various disciplines, or for artificial intelligence (AI) to infer causal associations from big data.

The ubiquitous complex nonlinear associations make it challenging to infer causation in Earth Systems. Fortunately, theories for dynamic system provide useful tools to reveal the nonlinear and intertwined relationships. Since Lorenz proposed a strategy to address the nonlinear complex problems in the state (phase) space, many prediction and causation detection methods have been developed based on Takens’ theory for time-series data. Takens’ theory proves that, for a dynamic system \phi, if its trajectory converges to an attractor manifold M, which are consisted by a bounded and invariant set of states, then the mapping between \phi and M can be built and time series observations of \phi can be used to construct M.

Earth systems are dynamic and spatial cross-sectional data can be treated as a snapshot of the dynamical system. In contrast to time series data, which are observations from a fixed spatial unit at different times, spatial cross-sectional data are observations from different spatial units at the same time. Both of them record different states of the dynamical system, only from different perspectives. The generalized embedding theorem, which is a theory on the mapping between differentiable manifold M and Euclidean space \mathbb{R} can be adopted to construct the state space with spatial cross-sectional data. Whitney proved that M can be mapped into \mathbb{R} by observation functions, and points in \mathbb{R} can be used to construct M. According to the theorem, for a compact d-dimensional manifold M and a set of observation functions h_1, h_2, ..., h_L, the map \psi_{\phi,h}(x) = (h_1(x), h_2(x), ..., h_L(x)) is an embedding of M with L = 2d + 1. Here embedding means a one-to-one map resolving all singularities of the original manifold. The elements h_i can be lags of observations from single time series observations, lags of observations from multiple time series, or multiple observation functions. The first two constructions are only special cases of the third one. By taking the measured values at one specific unit and its neighbors (named as spatial lags in spatial statistics) as a set of observation functions, the focal unit can then be uniquely determined by its value and its neighbors in a set of order of spatial lags, h_s(x), and h_{s(i)}(x), are their observation functions respectively. Hereinafter, we will use \psi(x,s) to present \psi_{\phi,h}(x,s) for short. In time series data, the lag means a shift from the observation at focal period t to past observation at t-k\tau. Similarly, in spatial data, the lag means a shift to spatial neighbors from the focal spatial unit. For raster data as illustrated in Fig. 1a, the first-order lags are adjacent units in eight directions in yellow color, whose first-order lags (removing those units already included) in turn constitute the second-order lags of the focal units. Similarly, for polygon data in Fig. 1b, the first-order lags are adjacent units sharing common edges or vertexes with the focal unit, and the lags of next order are the first-order lags of those adjacent units (excluding those already included). As the spatial lags in each order contain more than one spatial units, the observation function can be set as the mean of the spatial units or other summary functions considering the spatial direction, to assure the one-to-one mapping of the original manifold M.

For two spatial variables X and Y on the same set of spatial units, organized as regular grids (raster data) or irregular polygons (vector data), their values and spatial lags can be regarded as observation functions reading values from each spatial unit. Following the generalized embedding theorem, their shadow manifold M_x and M_y can be constructed by assembling above-defined \psi(x or y, s). According to the dynamical systems theory, if X and Y are observed from the same dynamical system, they are governed by the same manifold. Consequently, the one-to-one map between the reconstructed M_x and M_y can be deduced because they resolve all trajectories of the same original manifold without crossings. Set one point \psi(x,s) in M_x as a focal state, its corresponding state \psi(y,s) in M_y can be acquired according to the mutual spatial locations, as well as its close states. For example, \psi(x,s_1), \psi(x,s_2), \psi(x,s_3), and \psi(x,s_4) in M_x can be identified as close states of \psi(x,s) according to the distance in the state space. The corresponding states of them in M_y are \psi(y,s_1), \psi(y,s_2), \psi(y,s_3), and \psi(y,s_4) respectively based on the mutual spatial locations. Due to the mutual evolution, in M_y, \psi(y,s_1), \psi(y,s_2), \psi(y,s_3) are also close neighbors of \psi(y,s) as illustrated in Fig. 2a. Therefore, for a given x, the value of y can be predicted according to its close neighbors identified from M_x. This type of prediction based on nearest mutual neighbors is defined as cross-mapping prediction, as Eq. (1).

\hat{Y}_s|M_x = \sum_{i=1}^{L+1} w_{s_i} Y_{s_i} | M_x (1)

where s represents a spatial unit at which the value of Y need to be predicted, \hat{Y}_s is the prediction result, L is the number of dimensions of the embedding, s_i is the spatial unit used in the prediction, Y_{s_i} is the observation value at s_i and simultaneously the first component of a state in M_y, noted as \psi(y,s_i). In further, \psi(y,s_i) is determined by its one-to-one mapping point \psi(x,s_i), which is in turn one of the L+1 nearest neighbors of the focal state \psi(x,s) in M_x. w_{s_i} is the corresponding weight defined in Eq. (2).

w_{s_i}|M_x = \frac{weight(\psi(x,s_i), \psi(x,s))}{\sum_{i=1}^{L+1} weight(\psi(x,s_i), \psi(x,s))} (2)

where weight(*,*) is the weight function between two states in the shadow manifold, defined as Eq. (3).

weight(\psi(x,s_i), \psi(x,s)) = exp(-dis(\psi(x,s_i), \psi(x,s))) (3)

dis(\psi(x,s_i), \psi(x,s)) = \frac{1}{L} |h_{s_i}(x) - h_s(x)| + \sum_{k=1}^{L-1} abs[h_{s_i(k)}(x), h_{s(k)}(x)] (4)

where |*| means the absolute value of a real number, and abs[*,*] represents the distance function between two vectors, as the first element h_{s_i}(x) in \psi(x,s_i) corresponds to the spatial focal units, while \psi(x,s_i) and corresponding \psi(y,s_i) are close neighbors of \psi(x,s) and \psi(y,s) respectively. For example, h_{s(1)}(x) in Fig. 1(a) has eight members in yellow color, h_{s(1)}(x) = {u_{nw}, u_n, u_{ne}, u_e, u_{se}, u_s, u_{sw}, u_s}. The concrete form of abs[*,*] for raster data and polygon data are specified as abs_r and abs_v in Eq. (5) and Eq. (6) respectively.

For raster data, as the number and location of neighbors in a certain order are fixed, abs_r[*,*] can be defined as the averaged absolute difference of each spatial unit with the consideration of anisotropy. For polygon data, since the number and direction of neighbors in a fixed order of two spatial focal units differ frequently, abs_v[*,*] is defined as the absolute difference of spatial lags. When the anisotropy can be neglected, Eq. (6) can also be used as the absolute value function of raster data.

abs_r(h_{s_i(k)}(x), h_{s(k)}(x)) = \frac{1}{D} \sum_{d} |u_{s_i(k,d)}(x) - u_{s(k,d)}(x)| (5)

abs_v(h_{s_i(k)}(x), h_{s(k)}(x)) = |\frac{1}{D_1} \sum_{d} u_{s_i(k,d)}(x) - \frac{1}{D_2} \sum_{d} u_{s(k,d)}(x)| (6)

where u_{s_i(k,d)}(x) is the spatial unit of the kth-order spatial lags of s_i in the direction d. D is the number of spatial units (or directions) in the kth-order.

The skill of cross-mapping prediction is measured by the Pearson correlation coefficient between the true observations and corresponding predictions, defined in Eq. (7)

\rho = \frac{Cov(Y, \hat{Y})}{\sqrt{Var(Y) Var(\hat{Y})}} (7)

The prediction skill \rho varies by setting different sizes of libraries, which means the quantity of observations used in reconstruction of the shadow manifold. For raster data, since spatial units are regularly arranged with an equal area, the window size is used to represent the size of library. For polygon data, the number of spatial units is used to represent the size of library due to the irregular arrangement and unequal areas. Sugihara et al. suggested to use the convergence of \rho to infer the causal associations. For GCCM, the convergence means that \rho increases with the size of libraries and is statistically significant when the library becomes largest. By plotting a line graph of \rho versus libraries’ size (as Figs. 3–5), the increasing trend can be determined from it. In further, the null hypothesis for the significance test of \rho is that H_0: \rho = 0, and the statistics t in Eq. (8) can be used with the Student t distribution to get the p-value. Meanwhile, the confidence interval of \rho can be estimated based the statistics z in Eq. (9) with the normal distribution. In following case studies, the significance level is set to 0.05 and the confidence interval is 95%.

t = \rho \sqrt{\frac{n-2}{1-\rho^2}} (8)

where n is the number of observations to be predicted.

z = \frac{1}{2} ln \frac{1+\rho}{1-\rho} (9)

If X and Y are from the same dynamic system, they are causally linked and can be predicted with the cross-mapping method. To assure the causal associations between X and Y are bidirectional or they are both the effects of a shared cause, it is a required condition that both \psi(x,s_i) and corresponding \psi(y,s_i) are close neighbors of \psi(x,s) and \psi(y,s) respectively. In this case, with the increase of libraries, the points in the local area containing \psi(x,s) (or \psi(y,s)) become denser and the cross-mapping prediction becomes more reliable, as demonstrated in Fig. 2a. If the causal association is unidirectional and X is an external cause (noted as X \rightarrow Y), the original manifold contains dynamic information of both X and Y, while Y also contains information of X and itself. Meanwhile, X does not contain the dynamic information that solely belongs to Y. Consequently, M_x is a lower dimensional sub-manifold of the complete system, while M_y is a one-to-one map of M. It is not guaranteed that \psi(y,s_i), corresponding to \psi(x,s_i), is close to \psi(y,s), as illustrated in Fig. 2b. So X can be predicted reliably with Y, but not vice versa, i.e. the accuracy of the cross-mapping prediction of X with Y converges faster and notably stronger than the inverse effect of Y with X, as in Figs. 3–5. If Y cross-mapping predicting X (i.e. noted as Y xmap X) is higher, then we can interpreted as X causes Y (noted as X \rightarrow Y). However, if the causal influence of X is strong enough that Y becomes subordinated to X, they become synchronized (the well-known phenomenon of “synchronization”). Consequently, the capability of X cross-mapping predicting Y (noted as X xmap Y) is also strong. Therefore, when the cross-mapping prediction skills are strong in both directions, it either means that the causal associations are bidirectional, or the effect variable is enslaved by the cause variable. The interpretation of GCCM outputs under different scenarios are illustrated through three following case studies.

Extracted causations between soil pollution and multiple influencing factors
Soil heavy metal pollution exerts a strong influence on soil quality and public health. The contents of soil heavy metals are generally stable for a long period and large-scale surveys are very expensive and time-consuming. Thus, the datasets of heavy metal may be updated every 10 or 20 years, making the attribution of soil pollution based on the time series data and temporal statistical models not feasible. Meanwhile, due to the wide variety of pollution sources and complex factors that influence their accumulation, the correlations between the causal factors and soil pollution are always weak and insignificant. In this study, we select four heavy metals (the soil concentrations of Cu, Cd, Mg and Pb) to present soil pollution and two influencing factors, the density of industrial pollutants and residential pollutants. The residential pollutants is usually measured by the residential density, which could be reflected by the intensity of night light. Due to the lack of data concerning residential density data, the intensity of night light has been widely as the proxy variable of residence density. Hereinafter in this manuscript, we use nightlight intensity to delegate the residential density. As is well known, industry and residential wastes are two major causes of soil pollution. Since the ultimate challenge for verifying causation models is the lack of observable evidence, the clear existence of causation between the concentration of major heavy metals and industrial (or residential) density provides a valuable reference for assessing the reliability of causation models.

The output of causal inference using GCCM is displayed as Fig. 3. As shown in Fig. 3a–c, the spatial distribution of soil Cu concentrations, the density of industrial pollutants and nightlight, present little similarity and shared patterns. Meanwhile, the reconstructed manifold of soil Cu concentrations, the density of industrial pollutants and nightlight, as shown in Fig. 3d–f, present a similar pattern. It demonstrates that the reconstruction process of GCCM makes some unclear spatial associations between variables clearer in the phase space. Based on the reconstructed manifold, GCCM reveals the unidirectional causation between industrial and residential density and Cu concentrations. The \rho (cross-mapping prediction skill) of Cu xmap industrial density, which measures the causal effect of industrial pollution density on Cu (noted as industry \rightarrow Cu), in Fig. 3g, and Cu xmap nightlight, which measures the causal effect of residential pollution on Cu (residence \rightarrow Cu), in Fig. 3h are 0.37 (p=0.00) and 0.25 (p=0.00). Meanwhile, the \rho of industrial pollution density xmap Cu (Cu \rightarrow industry) and nightlight xmap Cu (Cu \rightarrow residence) are weak and insignificant, as 0.05 (p=0.09) and 0.00 (p=0.43) respectively. The larger and significant \rho of Cu xmap industrial density (industry \rightarrow Cu) and Cu xmap nightlight (residence \rightarrow Cu) suggest that industrial and residential pollution exert a strong influence on Cu concentrations. Conversely, the lower and insignificant \rho of industrial pollution density xmap Cu (Cu \rightarrow industry) and nightlight xmap Cu (Cu \rightarrow residence) suggest that Cu concentrations are not the cause of the density of industrial and residential pollutants.

Due to the lack of similarity in Fig. 3(a), (b) and (c), the commonly employed Pearson correlation analysis fail to extract significant correlations between Cu concentrations and two influencing factors, with correlation efficient as -0.08 (p=0.43) and -0.03 (p=0.08) for industrial pollution density and nightlight density respectively. The linear coefficients matrix (B pruned) outputted by LiNGAM is a zero matrix, indicating that none causal associations were identified between Cu concentrations and the two influence factors. The detailed coefficients and uncertainties of GCCM, correlation analysis and LiNGAM are shown in Supplemental Table S2. The causation inference outputs of other heavy metals (Cd, Mg and Pb) are similar to Cu and can also be found in Supplemental Table S2 and Fig. S2.

Extracted causations between population density and multiple influencing factors
In this case, we examine the county-level population density organized as polygon data, to evaluate the reliability of GCCM. Different from soil pollution, we select four environmental factors closely correlated with population density, including precipitation, temperature, elevation and slope, as shown in Supplemental Fig. S3. These factors are all significantly correlated with population density (Supplemental Table S3). It is worth mentioning that similar to CCM, the running of GCCM also requires a non-linear relationship between variables. To address this issue, we conduct a linearity-removal pre-processing for these variables before the running of GCCM. The detailed methodology for linearity-removal can be found in the Supplemental Section S3. The causation between population density and environmental variables revealed by GCCM are displayed in Fig. 4. The \rho of population density xmap precipitation, temperature, elevation and slope (precipitation \rightarrow population density, temperature \rightarrow population density, elevation \rightarrow population density and slope \rightarrow population density respectively) are large and significant, as 0.64 (p=0.00), 0.63 (p=0.00), 0.64 (p=0.00) and 0.48 (p=0.00). As a comparison, the \rho of precipitation, temperature, elevation and slope xmap population density (population density \rightarrow precipitation, population density \rightarrow temperature, population density \rightarrow elevation and population density \rightarrow slope) are also significant yet much smaller, as 0.27 (p=0.00), 0.31 (p=0.00), 0.36 (p=0.00) and 0.26 (p=0.00). It is worth mentioning that the much smaller \rho of environmental factors xmap population density (population density \rightarrow environmental factors) do not necessarily prove the existence of a reverse-direction causation. As explained by Suighara et al., when there is a strong causation from X to Y (X \rightarrow Y), then even if there is no feedback causal influence from Y to X, there can also be a smaller \rho of X xmap Y (Y \rightarrow X), resulting from the subordinating effect under the strong-association scenario. Therefore, the smaller \rho of X xmap Y can either indicate a weaker causal influence of Y on X (feedback causation), or the spatial variation of Y can have a small reflection on the spatial variation of X (reflection, not causation). The much smaller \rho of environmental factors xmap population density (population density \rightarrow environmental factors) did not necessarily prove the existence of a reverse-direction causation. When there is a strong causation from X to Y (X \rightarrow Y), then even if there is no feedback causal influence from Y to X, there can also be a smaller \rho of X xmap Y (Y \rightarrow X), resulting from the subordinating effect. Therefore, the smaller \rho of X xmap Y can either indicate a weaker causal influence of Y on X (feedback causation), or the spatial variation of Y can have a small reflection on the spatial variation of X (reflection, not causation). For this research, since population density can influence the temperature through intensive human activities, and the precipitation by releasing particulate matters, which are important nuclei for rain, the small \rho of precipitation xmap population density and temperature xmap population density can be regarded as a weaker feedback causation. On the other hand, population density cannot directly alter county-level elevation. Instead, population density can partially reflect the county-level elevation. For instance, when noticing a city with a population density notably lower than other cities in China, we can presume that the elevation of this city is very likely large and located in plateau regions. So the small \rho of elevation xmap population density can be regarded as a weak reflection. Although GCCM itself cannot distinguish the feedback causation or reflection in strong coupling cases, it can identify the leading direction of the casual associations.

As a comparison, the linear correlation coefficients between population densities and precipitation, temperature, elevation, and slope are 0.06 (p=0.00), 0.15 (p=0.00), -0.18 (p=0.00) and -0.23 (p=0.00) respectively. Although the linear correlations are also significant, the coefficients of GCCM are much higher and effectively reflect the asymmetric association between two variables. As a comparison, the fixed value of the correlation coefficient cannot help understand the asymmetric direction and strength of causation between two variables. The linear coefficients of LiNGAM between precipitation, temperature, elevation and slope, and population density are 0.00, 0.00, -1.08, and -364.03. It demonstrates that LiNGAM, with linear relationship and independent non-Gaussian distributed noise assumption, fails to identify the causal associations between population density and precipitation, and temperature. The detailed results of three method can be in found in Supplemental Table S3.

Extracted causations between farmland NPP and climate factors
In addition to the above cases which illustrate the performance of GCCM when processing raster data with weak correlations and polygon data with strong correlations respectively, we also check whether GCCM could effectively reveal the well-known NPP-Climate causation, a case explored in our previous study. As is widely acknowledged, precipitation and temperature are the main driving forces of farmland NPP. However, Gao et al. proved that mainstream temporal causation models such as CCM and Granger Causality Test (GCT) failed to infer causal associations between farmland NPP and precipitation (or temperature) due to insignificant temporal variations of these variables from 2000 to 2015. In this case, we employ the same datasets, including the average farmland NPP, precipitation, and temperature across China, to infer NPP-Climate associations using GCCM. Considering the strong linear correlation between NPP and precipitation (temperature), we conducted the same linearity-removal pre-processing (as in Supplemental Section S4) before the running of GCCM. The causal inference outputs from GCCM are displayed in Fig. 5, and detailed outputs of correlation analysis and LiNGAM can be found in Supplemental Table S4.

For GCCM, the prediction skill \rho of farmland NPP xmap precipitation (precipitation \rightarrow farmland NPP) and temperature (temperature \rightarrow farmland NPP) are 0.48 (p=0.00) and 0.50 (p=0.00). Conversely, the \rho of precipitation and temperature xmap farmland NPP (farmland NPP \rightarrow precipitation and farmland NPP \rightarrow temperature) are 0.00 (p=0.50) and 0.17 (p=0.00) respectively. These outputs prove a strong casual influence of precipitation and temperature on farmland NPP, which is consistent with the well-known facts, yet cannot be extracted by major temporal causation models. The zero \rho of precipitation xmap farmland NPP indicate that farmland NPP is not a cause of precipitation. And the much-smaller \rho of temperature xmap farmland NPP majorly result from the above-introduced enslaved effect from the strong causal influence of temperature on farmland NPP. In other words, farmland NPP can partially reflect temperature.

Based on the same dataset, the Pearson correlation coefficient between NPP and temperature and precipitation are 0.76 (p=0.00) and 0.66 (p=0.00), both significant. However, this coefficient cannot help understand the asymmetric direction and strength of causation between two variables and can be largely biased due to other confounding factors. As a comparison, LiNGAM fail to identify the existence of causal associations between farmland NPP and precipitation (temperature), with all coefficients as 0.00. The incapacity of LiNGAM in extracting NPP-climate causation, is mainly caused by the fact that its strong assumption on the linear relationship and noise distribution is easily violated in Earth System Science. The detailed outputs of three method can be in found in Supplemental Table S4.

Case 2 and Case 3 suggest when variables are strongly and significantly correlated, GCCM can also extract reliable causation between them. The leading cause-effect relationship between two variables (the notably larger \rho) can be well identified and measured through GCCM based on spatial cross-sectional data. Meanwhile, the feedback causal influence or the reflection capability from the target variable to influencing variables can as well be extracted as a smaller \rho, forming the unidirectional or bidirectional asymmetric associations between two variables. Since the correlation between two variables can be largely biased by other confounding variables, the correlation coefficient can be notably underestimated or overestimated, especially in highly complex ecosystems. In this case, the \rho value generated by GCCM, which is drawn from the nonlinear associations and thus more robust, can effectively reduce the influencing of other variables and provide a more reliable and complete reference for better understanding and comparing the strength of causation between the target variable and multiple influencing variables.

Methods
Data sources
We employ three typical cases to illustrate the effects, reliability, and interpretation of GCCM in the real world. The first is the causation between soil pollution and industrial pollution and residential density based on raster data. The second is the causation between population density and a series of environmental factors based on vector data. The third is the causation between farmland NPP and climate factors. The soil pollution data are obtained from the National Geochemical Survey database of the U.S. covering Illinois and Indiana. The population density data of China are obtained from the National Bureau of Statistics of China, and the climate and topographical data are obtained from published dataset. The NPP data are MOD17A3V055 masked with farmland of China, which are obtained from the China Multi-period Land Use Land Cover Remote Sensing Monitoring DataSet (CNLUCC). The dominant linear trend is removed from the population density and farmland NPP with the Ordinary Least Squares regression method. The detailed setting and complete result of the experiment are presented in the supplementary documents.

GCCM is implemented with R language, with one implementation for raster data and another for polygon data. All the cases are run in R 4.2.0, the stats and pcalg package are adopted to conduct the linear correlation and LiNGAM method respectively.

Brief steps for conducting GCCM
GCCM aims to identify and measure the causation from spatial cross-sectional data. It reconstructs the manifold in state space with spatial observations and their spatial lags, and infer causation according to the cross-mapping prediction in the state space. With the major principles and parameters of GCCM already introduced in details, the implementation of GCCM is realized in the following major steps:

1. Construct the embeddings: Firstly, set the dimension of the embeddings as L; Then go through each spatial unit (pixel or polygon) s and find its spatial lags of different orders (as Fig. 1a, b), each of which is marked as s(1), s(2), ..., s(L-1), and contains a set of neighboring units; Following this, get values of X and Y at the corresponding units and organize them as vectors, {h_s(x), h_{s(1)}(x), ..., h_{s(L-1)}(x)} and {h_s(y), h_{s(1)}(y), ..., h_{s(L-1)}(y)}; Finally, assemble those vectors together as a matrix according to their spatial orders.

2. Predict Y based on X: Set a sequence of library sizes (i.e. the widow size for raster data and the number of spatial units for polygon data); For each library size, predict Y as Eq. (1) by searching near points in the state space and confined by the library size; Calculate prediction skill \rho as Eq. (7), and finally get a series of \rho of Y with different library sizes.

3. Predict X based on Y: predict X using the similar approach in the above step and get a series of \rho of X with different library sizes.

4. Extract and interpret causations: Plot the prediction outputs from above two step as line graph using \rho on the vertical axis and library sizes on the lateral axis (as demonstrated in Figs. 3–6); Identify the existence, direction and strength of bidirectional (or unidirectional) causation based on the convergence, significance and confidence interval of \rho with the largest library size. It should be noted, which may be misunderstood, that the line of Y xmap X is the basis to determine whether X is a cause of Y (noted as X \rightarrow Y). If the line of Y xmap X is much higher with the increase of the library size, then X causes Y (noted X \rightarrow Y) is the leading direction in a bidirectional causation, or the sole direction in a unidirectional causation.

The dimension parameter L can be adjusted automatically to find the largest \rho. The library sizes can be commonly set as Arithmetic Sequence to reveal the trend of prediction skill.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Summary:\nIn the article under review, the authors adapt convergent cross mapping (CCM), a dynamical systems-\nbased method for inferring cause-effect relationships between time series, to spatial data. As the\nauthors correctly point out in the introduction, although CCM has been applied in several different\nfields of research (e.g. the Earth sciences), it can often be difficult to obtain time series data of\nsufficient length, for instance if measurements only started recently compared to the time scale one is\ninterested in. Spatial cross-sectional data on the other hand might be available more abundantly or at\nleast can be obtained more easily. Therefore, the method presented in this work, GCCM, can be a\nvaluable and potentially impactful contribution to the literature on causal inference for dynamical\nsystems.\nHowever, while the idea behind GCCM is original and the real data examples to which the method is\napplied are interesting, I also have serious criticism concerning the interpretation of the method’s\nresult, the discussion of the method’s applicability and the presentation of the material.\nRegarding the former, as I explain in more detail below, I believe that the results of the first real-data\nexample should be interpreted more carefully. In addition, to be of real value to domain experts, any\ncausality method should be accompanied by a thorough discussion of its underlying assumptions, its\nlimitations, and its failure modes. Unfortunately, such a discussion is mostly lacking here. For\nexample, CCM-based methods have been found to be vulnerable to dynamic noise, see e.g. (Bartsev\net al., 2021; Cui, Moore, 2021; Krakovska et al, 2018; Cobey, Baskerville, 2016), and this is likely to\nbe a weakness of GCCM as well. Moreover, the criteria that the method invokes (convergence and rho\n> 0) are somewhat vague (for instance, when is the method considered to converge exactly?) and it\ndoes not come with a measure of uncertainty that conveys to a user how trustworthy the output is.\nLastly, the article contains many typos and other minor issues, and the methodology section is not\nparticularly easy to follow.\nIn conclusion, I unfortunately must vote to reject the publication of this work. Due to the novelty of\nthe idea, a resubmission might be justified if the issues above are addressed properly.\n* Bartsev, S., et al. \"Imperfection of the convergent cross-mapping method.\" IOP Conference Series:\nMaterials Science and Engineering. Vol. 1047. No. 1. IOP Publishing, 2021.\n* Cui, Lipeng, and Jack Murdoch Moore. \"Causal network reconstruction from nonlinear time series: A\ncomparative study.\" International Journal of Modern Physics C 32.04 (2021): 2150049.\n* Krakovska, Anna, et al. \"Comparison of six methods for the detection of causality in a bivariate time\nseries.\" Physical Review E 97.4 (2018): 042207.\n* Cobey S, Baskerville EB (2016) Limits to Causal Inference with State-Space Reconstruction for\nInfectious Disease. PLoS ONE 11(12): e0169050.\nDetailed comments to the authors:\n• As a major motivating example, you mention the causal link between Net Primary Productivity and\nprecipitation which can’t be identified by time series-based causality methods. Yet you don’t apply\nyour method to this example and only refer to other spatial methods in the introduction. Is GCCM able\nto find this causal link?\n• At this point, the baseline to which you compare the performance of GCCM is correlation analysis. It\nwould be good to add another baseline method to your analysis that aims to infer causal links. This\ncould be a naïve method, e.g. treating the spatial observations as i.i.d. observations of one random\nvariable and then using a bivariate causality method such as LiNGaM. I expect GCCM to clearly\noutperform this baseline, but it would be nice to a have an empirical demonstration.\n• The weaknesses and limitations of GCCM should be discussed more extensively. How vulnerable is\nthe method to noisy data? Did you encounter real world or simulated data examples where GCCM\ncould not identify the correct causal direction? It would be helpful if an example that is a failure mode\nof the method would be included to illustrate the methods limitations more clearly.\n• The methods section ‘The establishment of GCCM’ feels a bit rushed. It would improve readability if\nyou took more time to explain your notations and disentangle some of the sentences such as in lines\n397-400. Computing different quantities explicitly for a toy example would also help.\nComments on real data examples:\nI like the choice of the two real data examples. However, there are additional aspects that should be\ndiscussed.\n• In the first example, the role of proxies should be discussed more clearly. What I mean is that you\nargue that there is a unidirectional causal link from nightlights to soil pollution whereas clearly, the\nfact that there is light at night does not pollute the soil. Therefore, in the interventional meaning of\ncausation, the variable nightlights is not a cause of soil pollution. Rather, nightlights are a good proxy\nfor an actual underlying cause, the number of residents in the area. In a way, you don’t detect a\ncausal link but a confounded connection between soil pollution and nightlights, the confounder being\nthe number of residents. However, the variable nightlights is so informative about the confounder (i.e.\nsuch a close proxy) that it inherits the causal asymmetry between residents and soil pollution. This is\na reason why the asymmetry that CCM-based methods such as yours detect should be interpreted\nmore carefully.\n• In the second example, there seems a mismatch between the labels of the plot and the causal\ndirection that they identify (the labels should be switched). Am I right in assuming that this is not a\nconceptual issue but rather a mix up when that happened when generating the labels of the causal\ndirections in the plot?\n• Is there a way to assess how significant the results of your method are and how much confidence\none should attribute to them? I agree that the plots in your examples look convincing but those were\nrelatively clear-cut cases from the start. In less clear-cut cases, it would be important to express the\nuncertainty about the results in order to interpret them properly or integrate them into other\napproaches as you propose in the discussion section.\nMinor observations:\n• P.1,l.35: typo, sate -> state\n• P.2, l.51: typo, casual -> causal\n• P.2, l.53 : denies -> denied,\n• P.2, l.53 : the 1980s\n• P.2. l 57ff: Crucially, counterfactual SCMs in the sense of Pearl and potential outcomes are two sides\nof the same coin as they can be shown to be logically equivalent, see (Peters et al.)\n• P.2, l.69: you should elaborate more clearly here why the models that are mentioned here are not\napplicable.\n• P.3, l.97, typo, spillover7\n• P.3, l. 104, proposed -> propose, generalized -> generalize\n• P.3, l. 123, typo: bias15\n• P.4, l. 148, typo: as well known -> as is well known\n• P.5, l. 153, the sentence suddenly ends\n• P.5, l. 164, the figure references need to be switched (g)<-> (h)\n• P.6, l.186: why does this demonstrate to effectively reveal bidirectional associations? This example\nis very clearly unidirectional as you also point on a bit further above.\n• P.9, l.299: research6\n• P.9, l. 309: research45\n• P. 9, l. 338 and l 350: casual", "reviewer_2": "Review of: Inferring causal associations from cross-sectional Earth System data\nGao et al\nFebruary 2023\nThis paper extends convergent cross mapping (CCM) to a spatial dimension, to identify causal\nassociations within “cross-sectional” or spatial datasets. This method is based on constructing\nmanifolds and considering lags as geographical distances between points. The topic is highly relevant\nin that there are many cases where we have high resolution datasets of spatial information, but\nlimited time-series data points, such that we need methods to draw causal conclusions from these\ndatasets. This was an interesting paper that presented two case studies based on pollution and\npopulation to illustrate the utility of the proposed GCCM. However, I found much of the paper and\nmethods difficult to understand, and many details are lacking. With this, I am not convinced that the\npresented GCCM results make sense, or constitute a robust way to define causal associations from a\nspatial perspective. These weaknesses are further described below. Unfortunately I cannot\nrecommend this paper for publication at this point in Nature Communications.\nDescription of GCCM is somewhat lacking, particularly for the meaning and significance of a mapping\nat a certain “spatial lag”, which is the corollary to a temporal lag in time-series analysis. After reading\nthe paper I still cannot understand why values of rho tend to increase with spatial lags for most\nrelationships in the case studies, since a spatial lag should imply some closeness within an attractor\nfor pixels that are a certain distance away from each other. More generally, the extension of CCM to a\nspatial context could be interesting, but the lack of details in this paper and weaknesses in the case\nstudies makes it difficult to tell.\nThis paper would benefit from more definitions that do not rely on the reader being already specialized\nin CCM.\nThe case studies based on pollution and population are not very clear, and it is not convincing that the\nresults are novel. Details are missing regarding the actual datasets used (e.g. in the first one, it seems\nlike it is data from China but it is actually from Indiana). Additionally, the conclusions are not\nparticularly interesting, since a causal association is similarly found in all cases in the expected\ndirection, and meanwhile the aspect of associations at different spatial lags are not well explained.\nThese weaknesses are highlighted below in line-by-line comments.\nLine by line comments:\nLine 28: This is a very long sentence, split up. In general recommend to go through the paper and\nrephrase several sentences like this.\nLine 80: Here is where you first bring up “cross-sectional data” and I am not sure this is the best\nterm, or one that is very familiar. It seems like all your examples are just “spatial data” or something\nsimilar to that. I think cross-sectional is specifically referring to different aspects (e.g. population,\nindustry, and pollution) at a snapshot in time, but I think the key aspect in this study is really the\nspatial one, hence the “geographic” term in GCCM.\nLine 91: “time series presenting” sentence does not make sense\nLine 92-101: There are a lot of references to models with no explanations, e.g. Graphical Detector,\nSTUVA, spatial spillover. These made this paragraph difficult to interpret.\nLine 104: “resorting to” is an odd choice of wording, would say “used” or similar\nLine 123: Here and several places, seems like a referencing typo (“bias15”)\nLine 126: Here you mention China, which led me to infer your pollution study was based there,\nhowever according to the maps in Fig 1 and data statement at the end, it is Indiana and Illinois. This\ncontext needs to be more clear.\nLine 142: first mention of spatial lags, which are not intuitive. Can a spatial lag be defined as a\ndistance (or radius) away from the point of interest? And what do we expect to see with these relative\nto the provided result plots? See main comment above on this. Potentially this could use some sort of\ntoy or illustrative example with generated spatial data with a certain dependency embedded.\nLine 153: Figure 1 and.\nFigure 1: I am not sure how to interpret the x-axis of panels G and H, in that they are the spatial lags.\nWhy is the CCM value increasing with spatial lag for nightlights xmap Cu and industry xmap Cu?\nParticularly if L can be translated from number of pixels to distances (km), what do these values\nmean? Similar, if D, E, and F panels are showing the manifolds for just 1 and 2 spatial lags, what\nshould we be reading from these?\nFigure 2: These results seem backwards relative to what is described in the paper. Isn’t “population\ndensity xmap Elevation” the influence of population on elevation? And if so, why is it actually larger\nthan in the opposite direction? These could simply be labeled incorrectly. Additionally, the notation “X\nxmap Y” is awkward and there could be a better way to write this directional relationship.\nI have a similar question as for Figure 1 on how to interpret L on the x-axis, and the increasing\n“prediction outputs” with L. Finally, in this study, it seems like the different “input” variables behave\nextremely similarly which makes me suspicious – but I also see that the y-axis are not the same for all\npanels (panel A has a greater range).\nLine 263: The word “complicated ecosystems” is often used, but better to use complex?\nLine 364: redundant sentence, states the meaning of “dynamic”\nLine 442: “will be illustrated” is strange, since this comes after the main results section.\nMethods section: the value “rho” is never defined here, but it is the main vertical axis value of all the\nresults plots.\nLine 463: The website with the data (https://www.resdc.cn) is all in Chinese, which makes sense, but\nit seems like for an English-speaking audience there should be a way to obtain the underlying data.\nAs it is, it would be difficult for most readers to reproduce this study. Here, I was trying to see the\nspatial resolution of the population etc data (details that are missing in the paper) and realized it would not\nbe possible to do that.\nSupplementary Information: After reading the paper, I realized some of the information I was looking\nfor was in the SI, but the SI was not referenced at all the appropriate points. It would be more useful\nto put the SI into sections and refer to these within the paper, or put critical aspects regarding the\ndatasets into the main paper. The SI also has a very large table that lacks any references."}, "author_response": "To Reviewer 1:\nThanks so much for your very professional and detailed comments on our\nmanuscript. Clearly, you are an authority in causation inference, as many\ncomments proved you know the principle, advantage and limitations of\nmainstream causation models well. Against this background, we feel very\nencouraged that you gave us positive comment on the innovation, robustness and\npotential implementations of our new GCCM and we fully believe that your\nconstructive comments can largely improve this manuscript. In the past two\nmonths, According to your comment, we have fully revised this manuscript by\nadding the NPP-climate case and the uncertainty test. Meanwhile, you mentioned\nthat we should consider the influence of random noise in the data sources on the\nperformance of GCCM, which we believe is a very important point. Therefore, we\nconducted a specific experiment to test the effects of gradually increased noise in\nthe data sources on GCCM outputs. The results proved what you said, the\nperformance of GCCM decreased with the increase of noise. But this experiment\nalso proved GCCM remains valid until the noise was very high. Furthermore,\nsince you suggested the use of such causation models as LiNGaM to compare with\nGCCM, we also employed LiNGaM for the three cases and LiNGaM failed to\nextract those causations, which further proved the major advantage of GCCM.\nFurthermore, since you suggested some details of this methodology should be\nbetter explained, we not only added much details to the revised manuscript, but\nalso largely revised the structure of the manuscript to highlight the major details\nfor an easier and clearer understanding of GCCM.\nIn general, according to these comments, substantial revisions have been made\nand the manuscript has been improved significantly. Thank so much again for\nyour constructive comments, which are the essential for making this manuscript\nqualified for the esteemed NC journal.\nPlease feel free to contact us if additional revisions are required and we are more\nwilling to conduct any revision accordingly.\nComment 1: In the article under review, the authors adapt convergent cross mapping\n(CCM), a dynamical systems-based method for inferring cause-effect relationships\nbetween time series, to spatial data. As the authors correctly point out in the introduction,\nalthough CCM has been applied in several different fields of research (e.g. the Earth\nsciences), it can often be difficult to obtain time series data of sufficient length, for\ninstance if measurements only started recently compared to the time scale one is\ninterested in. Spatial cross-sectional data on the other hand might be available more\nabundantly or at least can be obtained more easily. Therefore, the method presented in\nthis work, GCCM, can be a valuable and potentially impactful contribution to the\nliterature on causal inference for dynamical systems.\nHowever, while the idea behind GCCM is original and the real data examples to which\nthe method is applied are interesting, I also have serious criticism concerning the\ninterpretation of the method’s result, the discussion of the method’s applicability and\nthe presentation of the material.\nRegarding the former, as I explain in more detail below, I believe that the results of the\nfirst real-data example should be interpreted more carefully. In addition, to be of real\nvalue to domain experts, any causality method should be accompanied by a thorough\ndiscussion of its underlying assumptions, its limitations, and its failure modes.\nUnfortunately, such a discussion is mostly lacking here. For example, CCM-based\nmethods have been found to be vulnerable to dynamic noise, see e.g. (Bartsev et al.,\n2021; Cui, Moore, 2021; Krakovska et al, 2018; Cobey, Baskerville, 2016), and this is\nlikely to be a weakness of GCCM as well. Moreover, the criteria that the method\ninvokes (convergence and rho > 0) are somewhat vague (for instance, when is the\nmethod considered to converge exactly?) and it does not come with a measure of\nuncertainty that conveys to a user how trustworthy the output is. Lastly, the article\ncontains many typos and other minor issues, and the methodology section is not\nparticularly easy to follow.\nResponse: Thanks so much for your comments. According, we have interpreted\nthe casual relationship in the first real-data case (Line 295-304). For this case, the\nnightlight was a proxy variable, which reflected residential density and residential\nwaste, well-known as a major cause of soil pollution.\nWhat you mentioned here is very important and we are sorry that we did not make\nthe limitations of GCCM clearer in the previous manuscript. In the revised version,\nwe fully revised the structure and put much detailed introduction of GCCM in the\nfirst part of the Results(Line 132-285). In this case, readers can better understand\nGCCM. Meanwhile, we have added a specific paragraph to the discussion part\nconcerning the limitations of GCCM (Line 522-542). As revealed by other scholars,\nthe limitations of CCM-based methods proposed mainly include: being sensitive\nto periodic fluctuations and vulnerable to dynamic noise. The periodic variation is\nnot a main issue for GCCM, as periodic spatial phenomena are very rare in\npractice. According to your comment, to investigate the influence of noise on\nGCCM, we added different levels of noise to the heavy metal data (the first case)\nand ran the GCCM on the synthetic data. By increasing the added random noise\nfrom 10% to 30% , 60% and finally 90% of the original observed value, the\nprediction skill of heavy metals based on industry pollution and residential density\ndecreased. However, GCCM remained robust in identifying the causation with low\nor medium noise. Only when the added noise reached 60%, GCCM failed to\nextract causation of residential density on Cu.\nThe uncertainty test was really a good idea. In the revised manuscript, according\nto your suggestions, we employed two criteria (the significance and the\nconfidential interval) to demonstrate the uncertainty of the results of GCCM. We\nhave introduced the testing method (Line 251-252) and added it to the\ncorresponding R codes.\nIn the revised manuscript, according to your comment, we added the definition of\nconvergence(Line 248-251). Since GCCM was the spatial extension of CCM, we\nhad used the original definition of convergence from CCM (Sugihara et al., 2012)\nas “ρ increases with the size of libraries and is statistically significant when the\nlibrary becomes largest”.\nWe believe your comment that the principle and interpretation of GCCM can be\nbetter explained are very important and useful. To address this issue, we have fully\nre-organized the structure of the manuscript according to similar Method papers\nin NC, re-rewritten the method section by adding many more details, and invited\nall authors, including native-speaker scholars to carefully polish the manuscript.\nThanks again for all your constructive general and detailed comments.\nComment 2: In conclusion, I unfortunately must vote to reject the publication of this\nwork. Due to the novelty of the idea, a resubmission might be justified if the issues\nabove are addressed properly.\n* Bartsev, S., et al. \"Imperfection of the convergent cross-mapping method.\" IOP\nConference Series: Materials Science and Engineering. Vol. 1047. No. 1. IOP\nPublishing, 2021.\n* Cui, Lipeng, and Jack Murdoch Moore. \"Causal network reconstruction from\nnonlinear time series: A comparative study.\" International Journal of Modern Physics\nC 32.04 (2021): 2150049.\n* Krakovska, Anna, et al. \"Comparison of six methods for the detection of causality in\na bivariate time series.\" Physical Review E 97.4 (2018): 042207.\n* Cobey S, Baskerville EB (2016) Limits to Causal Inference with State-Space\nReconstruction for Infectious Disease. PLoS ONE 11(12): e0169050.\nResponse: Thanks so much for your encouragement on the innovation of this\nmanuscript. And we have fully revised the manuscript according to all your\ngeneral and detailed comments and believe it has been largely improved. Thanks\nagain for your great help on our research. These references have all been added to\nthe revised manuscript.\nComment 3: As a major motivating example, you mention the causal link between Net\nPrimary Productivity and precipitation which can’t be identified by time series-based\ncausality methods. Yet you don’t apply your method to this example and only refer to\nother spatial methods in the introduction. Is GCCM able to find this causal link?\nResponse: This is a good idea and thanks so much for pointing it out. In the revised\nmanuscript, we added NPP-climate causation as the third case(Line 422-482), and\nthe result proved GCCM successfully identified the casual links. Meanwhile,\nLiNGAM, which could be applied to the spatial cross-sectional data, also failed to\nidentify the NPP-climate causation. Thanks again for this constructive comment,\nwhich further proved the advantage of GCCM.\nComment 4: At this point, the baseline to which you compare the performance of\nGCCM is correlation analysis. It would be good to add another baseline method to your\nanalysis that aims to infer causal links. This could be a naïve method, e.g. treating the\nspatial observations as i.i.d. observations of one random variable and then using a\nbivariate causality method such as LiNGaM. I expect GCCM to clearly outperform this\nbaseline, but it would be nice to a have an empirical demonstration.\nResponse: Again, this is a good suggestion. In the revised manuscript, we employed\nLiNGaM as a comparison with GCCM. For the three cases, LiNGaM failed to\nextract any causations in the first and third case, and simply identified part of\ncausations between population and influencing factors(Line 332-334, 408-412,455-\n459). Thanks so much for this valuable comment, which further proved that few,\nif not none, existing causation models were effective for causation inference in\ncross-sectional data.\nComment 5: The weaknesses and limitations of GCCM should be discussed more\nextensively. How vulnerable is the method to noisy data? Did you encounter real world\nor simulated data examples where GCCM could not identify the correct causal direction?\nIt would be helpful if an example that is a failure mode of the method would be included\nto illustrate the methods limitations more clearly.\nResponse: Again, thanks so much for this constructive comment. What you\nmentioned here is very important and we are sorry that we did not make the\nlimitations of GCCM clearer in the previous manuscript. In the revised version,\nwe fully revised the structure and put many detailed introduction of GCCM in the\nfirst part of the Results. In this case, readers can better understand GCCM.\nMeanwhile, we have added a specific paragraph to the discussion part concerning\nthe limitations of GCCM(Line 522-542). As revealed by other scholars,\nthe limitations of CCM-based methods proposed mainly include: being sensitive to\nperiodic fluctuations and vulnerable to dynamic noise. The periodic variation is\nnot a main issue for GCCM, as periodic spatial phenomena are very rare in\npractice. According to your comment, to investigate the influence of noise on\nGCCM, we added different levels of noise to the heavy metal data (the first case)\nand ran the GCCM on the synthetic data. By increasing the added random noise\nfrom 10% to 30% , 60% and finally 90% of the original observed value, the\nprediction skill of heavy metals based on industry pollution and residential density\ndecreased. However, GCCM remained robust in identifying the causation with low\nor medium noise. Only when the added noise reached 60%, GCCM failed to\nextract some causations between heavy metals and industry (residence), which\nindirectly proved the robustness of GCCM.\nComment 6: The methods section ‘The establishment of GCCM’ feels a bit rushed. It\nwould improve readability if you took more time to explain your notations and\ndisentangle some of the sentences such as in lines 397-400. Computing different\nquantities explicitly for a toy example would also help.\nResponse: Thanks so much for this comment. In the revised manuscript, we fully\nrevised the structure of this manuscript and put the description of GCCM earlier\nin the manuscript (similar to other method papers in NC), so it is easier for readers\nto understand the methodology of GCCM. In Fig.1 and Fig.2, toy data was\ndisplayed to explain GCCM. Meanwhile, many details have been added to the\nrevised manuscript (Line 132-285), thus we believe the readability of the\nmethodology has been largely improved. Thanks again for your help.\nComment 7: Comments on real data examples: I like the choice of the two real data\nexamples. However, there are additional aspects that should be discussed.\nIn the first example, the role of proxies should be discussed more clearly. What I mean\nis that you argue that there is a unidirectional causal link from nightlights to soil\npollution whereas clearly, the fact that there is light at night does not pollute the soil.\nTherefore, in the interventional meaning of causation, the variable nightlights is not a\ncause of soil pollution. Rather, nightlights are a good proxy for an actual underlying\ncause, the number of residents in the area. In a way, you don’t detect a causal link but\na confounded connection between soil pollution and nightlights, the confounder being\nthe number of residents. However, the variable nightlights is so informative about the\nconfounder (i.e. such a close proxy) that it inherits the causal asymmetry between\nresidents and soil pollution. This is a reason why the asymmetry that CCM-based\nmethods such as yours detect should be interpreted more carefully.\nResponse: Thanks so much for your comment. Yes, your understanding is perfect,\nand what we meant was exactly the same. We are sorry we did not make this\nclearer in the previous version. In the revised manuscript, we had made it very\nclear that nightlight was a proxy variable of the residential density, which was one\nmajor influencing factor for heavy metal pollution (Line 295-300). Thanks again\nfor your comment, which helped us address the potential confusion.\nComment 8: In the second example, there seems a mismatch between the labels of the\nplot and the causal direction that they identify (the labels should be switched). Am I\nright in assuming that this is not a conceptual issue but rather a mix up when that\nhappened when generating the labels of the causal directions in the plot?\nResponse: Thank you very much for pointing out this problem. Yes, this is a mixed\nlabeling and we have checked it all through the text and figures, and make sure\nthere was no any potential confusion or mix up in the revised manuscript. Yes, as\nexplained in CCM, X xmap Y means the influence of Y on X (YX), which\nsuggests Y is the cause and X is the effect. In the revised manuscript, we have\nimproved the expression of the causation labeling for a better and clear\nunderstanding. Thanks so much for your valuable comment, which help us reduce\npotential confusion.\nComment 9: Is there a way to assess how significant the results of your method are\nand how much confidence one should attribute to them? I agree that the plots in your\nexamples look convincing but those were relatively clear-cut cases from the start. In\nless clear-cut cases, it would be important to express the uncertainty about the results\nin order to interpret them properly or integrate them into other approaches as you\npropose in the discussion section.\nResponse: The uncertainty test was really a good idea. In the revised manuscript,\naccording to your suggestions, we employed two criteria (the significance and the\nconfidential interval) to demonstrate the uncertainty of the results of GCCM. We\nhave introduced the testing method (Line 251-252) and added it to the\ncorresponding R codes.\nComment 10:\nMinor observations:\n• P.1,l.35: typo, sate -> state\n• P.2, l.51: typo, casual -> causal\n• P.2, l.53 : denies -> denied,\n• P.2, l.53 : the 1980s\n• P.2. l 57ff: Crucially, counterfactual SCMs in the sense of Pearl and potential outcomes\nare two sides of the same coin as they can be shown to be logically equivalent, see\n(Peters et al.)\nResponse: Thanks so much for your comment. We have checked and revised all\nthese typos and checked through the manuscript. We quite agree with you that\nthey were logically equivalent. The main differences between SCMs and potential\noutputs is the emphasis on casual graph. For the causal graph, Imbens and Rubin\n(2015) said that “we have not found this approach to aid drawing of causal\ninferences,”. And Pearl and Mackenzie wrote that “Rubin has steadfastly\nmaintained over the years that diagrams serve no useful purpose”. Therefore,\nconsidering that most readers do not have such an incisive understanding of\ncausation models, we divided them into two categories. Again, we agreed with you,\nand we think that the SCMs and potential outcomes could be combined together,\nas well as GCCM. We discussed this in the discussion part.\nReference:\nImbens GW, Rubin DB. Causal Inference for Statistics, Social, and Biomedical\nSciences: An Introduction. Cambridge University Press (2015).\nPearl J, Mackenzie D. The Book of Why: The New Science of Cause and Effect. Basic\nBooks (2018).\nComment 11: P.2, l.69: you should elaborate more clearly here why the models that\nare mentioned here are not applicable.\nResponse: Thanks so much for your comment. We have revised the statement to\nbetter express it(Line 65-67). What we wanted to explain here is that two\nframeworks are mainly for stochastic processes, while earth system sciences often\ncontain deterministic or dynamical interactions, which breaks the assumption of\nstochastic processes (e.g. information separability and probability distribution).\nComment 12:\n• P.3, l.97, typo, spillover7\n• P.3, l. 104, proposed -> propose, generalized -> generalize\n• P.3, l. 123, typo: bias15\n• P.4, l. 148, typo: as well known -> as is well known\n• P.5, l. 153, the sentence suddenly ends\n• P.5, l. 164, the figure references need to be switched (g)<-> (h)\nResponse: Thanks so much for your comment. We have revise the typos and\nchecked through the manuscript.\nComment 13: P.6, l.186: why does this demonstrate to effectively reveal bidirectional\nassociations? This example is very clearly unidirectional as you also point on a bit\nfurther above.\nResponse: Thanks so much for your comment. You are right that it is\nunidirectional. It was a typo and we have revised it to “unidirectional associations”\nin the revised manuscript.\nComment 13:\n• P.9, l.299: research6\n• P.9, l. 309: research45\n• P. 9, l. 338 and l 350: casual\nResponse: Thanks so much for your comment. We have revised the typos and\nchecked through the manuscript.\nTo Reviewer 2:\nThanks so much for your very professional and detailed comments on our\nmanuscript. And many thanks for the encouragement and acknowledgement from\nyou and another reviewer, who are clearly both experts in spatial analysis and\nmodels, on the innovation and robustness of GCCM. Based on a careful analysis\nof your comments, we found they are all very constructive for use to improve the\nmanuscript. Some your concerns (such as the data availability links) are caused\nby our unclear expression that easily lead to misunderstanding. In the revised\nmanuscript, we rewrote them to clarify it according to your suggestion.\nThe most important concern you mentioned is that due to the lack of details, the\nprinciple, parameter setting and methodology of GCCM were not clearly\nexplained. Yes, due to limited spaces and a lot of details concerning GCCM setting,\nwhich is similar to CCM, we did not explain the setting of GCCM clearly enough.\nBut we fully agreed with you that even if some principles of GCCM have been\nintroduced before in previous paper concerning CCM (Sugihara et al. 2012), it\nremains necessary for us to explicitly explain our models in this manuscript.\nTherefore, according to your comment, we not only added many details to the\nrevised manuscript, but also largely revised the structure of the manuscript to\nhighlight the major details for an easier and clearer understanding of GCCM. All\nyour raised issues have been addressed or explained in the revised manuscript.\nMeanwhile, according to the suggestions from the other reviewer, we added\ncontent concerning another case study, uncertainty test, noise test and another\ncausation model LiNGaM f. Thanks so much to your professional comments. In\nthe past two months, we have made substantial efforts to revise the manuscript\naccordingly and we believe the manuscript has been improved significantly.\nPlease feel free to contact us if additional revisions are required and we are more\nthan willing to conduct revisions according to your constructive suggestions.\nComment 1: This paper extends convergent cross mapping (CCM) to a spatial\ndimension, to identify causal associations within “cross-sectional” or spatial datasets.\nThis method is based on constructing manifolds and considering lags as geographical\ndistances between points. The topic is highly relevant in that there are many cases where\nwe have high resolution datasets of spatial information, but limited time-series data\npoints, such that we need methods to draw causal conclusions from these datasets. This\nwas an interesting paper that presented two case studies based on pollution and\npopulation to illustrate the utility of the proposed GCCM. However, I found much of\nthe paper and methods difficult to understand, and many details are lacking. With this,\nI am not convinced that the presented GCCM results make sense, or constitute a robust\nway to define causal associations from a spatial perspective. These weaknesses are\nfurther described below. Unfortunately I cannot recommend this paper for publication\nat this point in Nature Communications.\nDescription of GCCM is somewhat lacking, particularly for the meaning and\nsignificance of a mapping at a certain “spatial lag”, which is the corollary to a temporal\nlag in time-series analysis. After reading the paper I still cannot understand why values\nof rho tend to increase with spatial lags for most relationships in the case studies, since\na spatial lag should imply some closeness within an attractor for pixels that are a certain\ndistance away from each other. More generally, the extension of CCM to a spatial\ncontext could be interesting, but the lack of details in this paper and weaknesses in the\ncase studies makes it difficult to tell. This paper would benefit from more definitions\nthat do not rely on the reader being already specialized in CCM.\nResponse: Thanks so much for your comment. Yes, due to limited spaces and a lot\nof similar details to CCM, we did not explain the setting of GCCM clearly enough.\nBut we fully agreed with you that even if some principles of GCCM have been\nintroduced before in previous paper of CCM (Sugihara et al. 2012), it remains\nnecessary for us to explicitly explain our models in this manuscript. Therefore,\naccording to your comment, we not only added many details to the revised\nmanuscript, but also largely revised the structure of the manuscript to highlight\nthe major details for an easier and clearer understanding of GCCM. All your\nraised issues have been addressed or explained in the revised manuscript.\nAs you suggested, the spatial lag is indeed the corollary to a temporal lag in time-\nseries. So we have added very detailed description of spatial lag in the revised\nmanuscript and specifically explained it here. In time series data, the lag k means\na shift from the observation at focal period t to past observation at t-k . Similarly,\nin spatial data, the lag means a shift to spatial neighbors from the focal spatial unit.\n𝝉\nThe spatial lags are used to compose a vector with the focal point. And the vector\nare then used to reconstruct the manifold in the state space. As illustrated in Figure\n3, due to the weak correlation, the spatial distribution of soil Cu concentrations,\nthe density of industrial pollutants, and the nightlight images present little\nsimilarity. But in the reconstructed manifold in the state space, they presented a\nsimilar pattern. GCCM is based on the one-to-one map in the reconstructed\nmanifold in the state space. For example, to predict the density of industrial\npollutants at spatial unit s, noted as Ys, GCCM employs the one-to-one mapping\nrelationship of corresponding points and in the reconstructed\nmanifold, which contains s as the first com𝑴po𝒙,n𝒔ents, an𝑴d t𝒚h,𝒔e neighbors of and\n. In fact, the lateral axis of the GCCM results is the size of librari𝑴es,𝒙 ,w𝒔 hich\n𝑴me𝒚a,𝒔ns the quantity of observations used in reconstruction of the shadow manifold.\nWe did not clarify this point in previous version, and we introduced this point in\ndetail in this revised version according to your suggestion. For raster data, since\nspatial units are regularly arranged with equal areas, the window size is used to\nrepresent the size of library. For polygon data, the number of spatial units is used\nto represent the size of library due to the irregular arrangement and unequal areas.\nWith the increase of the libraries, the points in the local area containing (or\n) of state space become denser and the cross-mapping prediction 𝑴be𝒙c,o𝒔mes\n𝑴mo𝒚r,𝒔e precise, as demonstrated in Fig. 1(c). Therefore, rho tends to increase with\nlibrary size when causal associations exist. The details and definitions were added\nto introduce GCCM more clearly.\nThanks again for your comment, according to which the manuscript has been\nimproved significantly.\nComment 2: The case studies based on pollution and population are not very clear, and\nit is not convincing that the results are novel. Details are missing regarding the actual\ndatasets used (e.g. in the first one, it seems like it is data from China but it is actually\nfrom Indiana). Additionally, the conclusions are not particularly interesting, since a\ncausal association is similarly found in all cases in the expected direction, and\nmeanwhile the aspect of associations at different spatial lags are not well explained.\nThese weaknesses are highlighted below in line-by-line comments.\nResponse: Thanks so much for your comment.\nFirstly, according to your suggestions, we added more details concerning the\ndataset and deleted the misleading words.\nSecondly, we believe there is a major misunderstanding here. Yes, you are right,\nthe extracted metal pollution-industry (residential) causations and population-\nenvironment causations were not novel to the society. However, the novel and\nundiscovered causations were not the focus of this research. Instead, all these cases\nwere employed to test the suitability of GCCM in different situations. As\nintroduced in our manuscript, it is a common challenge that there is no reliable\ntrue data to verify causation models. So the well-accepted prior-knowledge, such\nas industry-pollution causation, can be used as reference to verify the specific\nsuitability and interpreting capability of GCCM. Therefore, these not-new and\nwell-known causations, as you pointed out, worked effectively to reveal the\nadvantage of GCCM. As the other reviewer, who seems also an authority in\ncausation inference, mentioned he liked the choice of the two cases and he even\nsuggested us to a third case (which is also a well-known causation) to the revised\nmanuscript. According to his/her comments, we added another case concerning\nfarmland NPP-climate (Line 422-482). For the three case studies, GCCM\nsuccessfully extracted all these well-known causations and effectively identified\nthe direction and strength of these causations based on the cross-sectional data.\nMeanwhile, these clearly-existing causations cannot be identified by other\ntemporal causation models, such as the commonly employed model LiNGAM (a\nmodel suggested by the other reviewer). Generally, since no actual reference data\nare available for us to properly verify GCCM, we can only consider the strategy\nto test whether GCCM can identify those clearly existing causations. So, in the\nrevised manuscript, we further emphasized that the aim of these clearly existing,\nnot-new causations was to verify the suitability and robustness of GCCM, which\nis the major novelty and contribution of this research.\nThanks again for pointing out your concern and we believe that with clearer\nexplanation, readers can better understand our choice of case studies.\nThirdly, according to previous comments, we believe the “different spatial lags”\nyou mentioned here means “different library size”. As explained in details in the\nabove comment, the rho increase with library size when causal associations exist,\nbecause with the increase of the libraries, the points in the local area containing\n(or ) of state space become denser. To avoid potential confusions, we\n𝑴ex𝒙p,𝒔lained 𝑴th𝒚i,s𝒔 more clearly in the revised manuscript (Line 244-275).\nThanks again for all these comments. By addressing these confusions and\nmisunderstanding, this manuscript has been largely improved.\nComment 3: Line 28: This is a very long sentence, split up. In general recommend to\ngo through the paper and rephrase several sentences like this.\nResponse: Thanks so much for pointing this out. We have rephrased this sentence\nand checked throughout the manuscript to polish the language. Thanks again for\nthis valuable comment.\nComment 4: Line 80: Here is where you first bring up “cross-sectional data” and I am\nnot sure this is the best term, or one that is very familiar. It seems like all your examples\nare just “spatial data” or something similar to that. I think cross-sectional is specifically\nreferring to different aspects (e.g. population, industry, and pollution) at a snapshot in\ntime, but I think the key aspect in this study is really the spatial one, hence the\n“geographic” term in GCCM.\nResponse: Thank you so much for your suggestion. The key aspect in this study is\nthe spatial data. Therefore, we revised the cross-sectional data to spatial cross-\nsectional data in the revised manuscript.\nComment 5: Line 91: “time series presenting” sentence does not make sense\nResponse: Thank you so much for your suggestion. We have deleted this sentence\nin the revised manuscript.\nComment 6: Line 92-101: There are a lot of references to models with no explanations,\ne.g. Graphical Detector, STUVA, spatial spillover. These made this paragraph difficult\nto interpret.\nResponse: Thank you so much for this constructive comment. We have rewritten\nthis part and added explanation to all these key concepts(Line 93-100) .\nComment 7: Line 104: “resorting to” is an odd choice of wording, would say “used”\nor similar\nResponse: Thanks for pointing this out. We have revised it accordingly.\nComment 8: Line 123: Here and several places, seems like a referencing typo\n(“bias15”)\nResponse: Thanks for pointing this out. We have checked and corrected these.\nComment 9: Line 126: Here you mention China, which led me to infer your pollution\nstudy was based there, however according to the maps in Fig 1 and data statement at\nthe end, it is Indiana and Illinois. This context needs to be more clear.\nResponse: Thanks so much for pointing this out. Yes, the description in the\nprevious version may be misleading and thus we fully revised this paragraph to\navoid unnecessary confusions.\nComment 10: Line 142: first mention of spatial lags, which are not intuitive. Can a\nspatial lag be defined as a distance (or radius) away from the point of interest? And\nwhat do we expect to see with these relative to the provided result plots? See main\ncomment above on this. Potentially this could use some sort of toy or illustrative\nexample with generated spatial data with a certain dependency embedded.\nResponse: Thank you so much for your comment. Despite a limited space, we\nshould still try to make the details of the model setting clear. So as you suggested,\nmany details have been added to the revised manuscript (Line 167-181).\nYes, the spatial lag is defined as a distance (or radius) away from the point of\ninterest. In time series data, the lag k means a shift from the observation at focal\nperiod t to past observation at t-k . Similarly, in spatial data, the lag means a shift\nto spatial neighbors from the focal spatial unit. The spatial lags are used to\n𝝉\ncompose a vector with the focal point. To explain the spatial lags, we defined it in\ndetails, and illustrated in Fig.1a and b for raster data and polygon data based on\nsome toy data (not the data from three cases). The spatial lags do not join the result\nplots directly. It is the number of the vectors they composed were used as the\nlateral axis, i.e. the size of libraries. With the increase of the libraries, the points\nin the local area containing (or ) of state space become denser and the\ncross-mapping prediction be𝑴co𝒙m,𝒔 es m𝑴or𝒚e, 𝒔precise, as demonstrated in Fig. 1(c).\nTherefore, rho tend to increase with library size when causal associations exist.\nThanks so much again for your comment.\nComment 11: Line 153: Figure 1 and.\nResponse: Thanks for pointing this out. We have checked and corrected these.\nComment 12: Figure 1: I am not sure how to interpret the x-axis of panels G and H, in\nthat they are the spatial lags. Why is the CCM value increasing with spatial lag for\nnightlights xmap Cu and industry xmap Cu? Particularly if L can be translated from\nnumber of pixels to distances (km), what do these values mean? Similar, if D, E, and F\npanels are showing the manifolds for just 1 and 2 spatial lags, what should we be\nreading from these?\nResponse: Thanks for your point. I think there is some unclear parts here and we\nare giving more details to explain it here and in the revised manuscript(Line253-\n275 and Line 599-629).\nThe x-axis of panels G and H is the size of libraries, which means the quantity of\nobservations used in reconstruction of the shadow manifold. The spatial lags are\nused to compose a vector with the focal units. To simplify this problem, here we\nthe treat the observation function as the average of spatial lags with the same order.\nFor example , where the first component is the X\n〈𝒉𝒔(𝒙),𝒉𝒔(𝟏)(𝒙),….,𝒉𝒔(𝑳−𝟏)(𝒙)〉\nvalue of the focal unit s, is the average of X values of first order\n𝒉𝒔(𝟏)(𝒙)\nneighbors of s, is the average of X values of second order neighbors of s.\nThe spatial lags𝒉 a𝒔r(e𝟐 )u(s𝒙e)d to form the elements of the vectors, i.e. the dimension of\npoints in Fig.3 d-f (The Fig.1 in previous version), while the size of library means\nto the number of points in Fig.3 d-f. Fig.3 d-f plotted all observations, i.e. the\nlargest size of library. From Fig.3 d-f, it could be found that, due to the weak\ncorrelation, the spatial distribution of soil Cu concentrations, the density of\nindustrial pollutants, and the nightlight images present little similarity. But in the\nreconstructed manifold in the state space, they presented a similar pattern.\nTherefore we can predict them based on the one-to-one map in the reconstructed\nmanifold in the state space. For example, to predict the density of industrial\npollutants at spatial unit s, noted as Ys, GCCM employs the one-to-one mapping\nrelationship of corresponding points and in the reconstructed\nmanifold, which contains s as the first com𝑴po𝒙,n𝒔ents, an𝑴d t𝒚h,𝒔e neighbors of and\n. With the increase of the libraries, the points in the local area contain𝑴in𝒙g,𝒔\n𝑴𝒚,𝒔 𝑴𝒙,𝒔\n(or ) of state space become denser and the cross-mapping prediction becomes\nmor𝑴e p𝒚r,𝒔ecise, as demonstrated in Fig. 1(c). Therefore, rho tends to increase with\nlibrary size when causal associations exist.\nThanks again for your comment. The details and definitions were added to\nintroduce GCCM more clearly (Line 132-285 and 599-629).\nComment 13: Figure 2: These results seem backwards relative to what is described in\nthe paper. Isn’t “population density xmap Elevation” the influence of population on\nelevation? And if so, why is it actually larger than in the opposite direction? These could\nsimply be labeled incorrectly. Additionally, the notation “X xmap Y” is awkward and\nthere could be a better way to write this directional relationship. I have a similar\nquestion as for Figure 1 on how to interpret L on the x-axis, and the increasing\n“prediction outputs” with L. Finally, in this study, it seems like the different “input”\nvariables behave extremely similarly which makes me suspicious – but I also see that\nthe y-axis are not the same for all panels (panel A has a greater range).\nResponse: Yes, thanks so much for your careful checking and we acknowledged\nthis was a typo (or a mixed labelling). In the revised manuscript, we have checked\ncarefully and make sure no such mixed use happened again. Meanwhile, to make\nthe expression clearer and avoid potential confusions, we further improved the\nexpression in the revised manuscript according to your good advice.\nAs defined in CCM, X xmap Y is a specific defined term (Sugihara et al, 2012),\nmeans the influence of Y on X (YX), which indicates that Y is the cause and X is\nthe effect.\nSo Yes, “Population density xmap Elevation” means the casual influence of\nElevation on population density (Elevationpopulation density).\nBy revising the expression of causations, especially causation direction in the\nrevised manuscript, readers can better understand the direction of causations\nbetween variables.\nThanks so much again for your constructive comment.\nComment 14: Line 263: The word “complicated ecosystems” is often used, but better\nto use complex?\nResponse: Thank you so much for your suggestion. We have corrected accordingly.\nComment 15: Line 364: redundant sentence, states the meaning of “dynamic”\nResponse: Thank you so much for your suggestion. We have corrected accordingly.\nComment 16: Line 442: “will be illustrated” is strange, since this comes after the main\nresults section.\nResponse: Thank you so much for your suggestion. We have corrected accordingly.\nComment 16: Methods section: the value “rho” is never defined here, but it is the main\nvertical axis value of all the results plots.\nResponse: Thank you so much for your suggestion. We have corrected accordingly.\nWe defined the “rho” in the newly added Equation (7) and explained it in the\nrevised manuscript.\nComment 17: Line 463: The website with the data (https://www.resdc.cn) is all in\nChinese, which makes sense, but it seems like for an English-speaking audience there\nshould be a way to obtain the underlying data. As it is, it would be difficult for most\nreaders to reproduce this study. Here, I was trying to see the spatial resolution of the\npopulation etc data (details that are missing in the paper) and realized it would not\nbe possible to do that.\nResponse: Thanks for pointing it out. I think there is a major misunderstanding\nhere. In the previous manuscript, we had already uploaded a specific dataset,\nwhich included all the data sources used for this research and English data\ndescription, and introduced it in the data availability part as “The Code for\nGCCM and the data employed for this research are available at\nhttps://doi.org/10.6084/m9.figshare.21782201”. The website with the data\n(https://www.resdc.cn) is the original source we got the climate data and we listed\nit just to demonstrate its authority. The default language for this website is Chinese,\nwhich may cause such difficulties. To avoid such misunderstanding, we listed the\nwebsite with English instead (Line 586-591).\nSo generally, all readers can easily get the required data from the dataset we\nuploaded or directly download it from the website.\nThanks again for raise your concern.\nComment 18: Supplementary Information: After reading the paper, I realized some of\nthe information I was looking for was in the SI, but the SI was not referenced at all the\nappropriate points. It would be more useful to put the SI into sections and refer to these\nwithin the paper, or put critical aspects regarding the datasets into the main paper. The\nSI also has a very large table that lacks any references\nResponse: Thanks so much for your comment. According to your comment, for a\nclearer understanding, we not only moved many contents to the main text from SI,\nbut also largely restructured the manuscript and added many details. We also\nadded many references to SI tables as you suggested.\nThanks again for advising this, which improved the readability of this paper\nsignificantly.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nAnswers to rebuttal/ comments on first revised version:\nIn the revised version of their submission, the authors have incorporated several changes that\nsignificantly improve their work compared to the original version. Nevertheless, some issues remain,\nand these should be addressed before I can recommend this paper for publication. Given that the\nauthors have been able to improve the previous version in several places, I would encourage them to\nprovide a second revised version that addresses the remaining problems.\nBelow, when referring to earlier comments, I will employ the numbering that you used in your rebuttal\nto the first review.\nMajor comments to the authors:\n• Compared to the earlier draft, language and grammar have clearly improved in the main document.\nThere are still several grammatical mistakes and imprecise formulations, some of which I point out\nbelow, but these affect readability much less than before. Unfortunately, the supplement (especially\nthe first section) is still in a poor state language-wise and this does seriously affect intelligibility. This\nsection needs to be thoroughly edited for language structure and clarity, or with the help of software\ntools. In addition, the section feels rushed, and you use a lot of terminology without explaining it. I\nwould suggest that you either explain different concepts and assumptions in more detail, or you only\nmention their existence very shortly and refer to existing work for detailed explanations. Currently,\nyou use an in-between approach which is often hard to follow.\n• Comments 3,4 & 5: I like that you included the NPP – precipitation example and that you evaluated\nthe soil pollution example once more with different levels of added noise. Similarly, adding LinGaM as\nan additional baseline strengthens your argument. These changes improve the paper substantially in\nmy opinion.\n• Comment 6: The section on methodology that explains the theory behind GCCM has improved and it\nwas a good idea to move it to an earlier section. However, it is still hard to digest and sometimes too\nhand-wavy on the mathematics behind the method. I think you should explain the mathematical\ndetails in much larger detail and give more precise formal definitions of the objects you employ. It\nwould also help if some of the underlying results, such as the generalization of Takens’ theorem were\nincluded, at least in the supplement. From a mathematician’s point of view, some choices of notations\nare also suboptimal. This is the most fundamental section of your article, and it warrants a more\ncareful and precise discussion.\n• Comment 7: I think there is a bit of a misunderstanding here. The fact that your method finds a\ncausal link nightlights -> soil pollution is a challenge to the interpretation of your method’s output.\nNightlights are not a cause of soil pollution in the interventional sense of Pearl: If we would intervene\non nightlights by switching them off, this would not reduce soil pollution by itself. Rather the correct\ncausal graph is something like nightlights <- residential density -> soil pollution and the asymmetry\nyou detect is a consequence of the fact that the association between nightlights and residential density\nis very strong. If you applied your method on new data (X,Y) without much background knowledge\nand you find a link X->Y, you cannot claim that this is a causal link, but only that X is sufficiently close\nto a cause X’ of Y. This is still useful information and does not invalidate your method, but it needs to\nbe discussed much more thoroughly than what you currently do.\n• Comment 8: fixed.\n• Comment 9: I am glad that you added a significance test to your method. However, this warrants\nmore than two lines of text. The explanation on how you test is too hand-wavy in my opinion as this is\ncrucial information. Please be more specific on the following: what do you mean by ‘rho is statistically\nsignificant when the library becomes largest’? Do you mean ‘rho is significantly different from 0’? What\nis your test statistic for this and what is your choice of significance level? How do you measure\nconcretely whether there is an increase of rho with increasing size of L?\n• Comment 10 – 13: fixed.\nMinor comments:\n• For some reason, the abstract is written in past tense. I would suggest switching to present tense.\nThere are also some sentences in the introduction where you should change to present tense, e.g. line\n151.\n• Abstract: Causation inference in complicated systems -> Causal inference in complex systems.\nGenerally, the common term in the field is causal inference, not causation inference.\n• Abstract: split last sentence in two.\n• Line 44: the prerequisite -> a prerequisite\n• l. 58: ‘the strength of the cause’ -> the strength of the causal relationship\n• l. 60: ‘the randomized control trial’ -> randomized control trials\n• l. 67: ‘… which render these two frameworks unapplicable’. I find this too harsh, I would rather say\n‘… which are not captured well by the assumptions of these frameworks and therefore challenge their\napplicability.’\n• l. 69: ‘the spatial heterogeneity’ -> lose ‘the’\n• l. 107: again, lose ‘the’ before the theories\n• l. 133: this sentence is strange, what do you mean by ‘internal essences’?\n• l.158: consists -> is\n• l.161: dynamic system -> dynamical system; different from -> in contrast to\n• l.178: lags of current orders is strange terminology, maybe explain this sentence more clearly.\n• L.224: missing bracket )\n• P.6: I still find your notation conventions quite confusing, especially the notation\nM_{x,si}. You explain a lot in words, but I would like to see cleaner mathematical definitions of the objects that you\nuse. Why don’t you work more with the embedding map \\Psi_h for instance? As far as I understand,\nthe only thing you do here is to define a metric d on the manifold M_x by d(w,z) = d’(\\Psi_h(w),\n(\\Psi_h(z)) for arbitrary points w,z on M_x, where d’ is an averaged sum 1/L \\sum d_s. over the\nEuclidean metrics d_s on R^{|lag s neighbors|}. This is well-defined as \\Psi_h is an embedding.\n• Equation (6): with this definition of h, I don’t see why ||…|| defines a norm in a mathematical sense,\nand consequently why dis() defines a metric on M_x. More precisely, can you explain why ||…||or\ndis() satisfy the property of positivity of metrics here? Or is being a metric not essential? In the latter\ncase, use a different notation than ||…||as this strongly suggests that this is a norm in a mathematical\nsense. Again, it would help to have clearer notation in this section.\n• L.254: ‘it is a required condition’. Required for what exactly?\n• L. 448: in-slaved, do you mean enslaved?\n• L. 526: ‘the periodic time series’ -> ‘periodic time series’\nReviewer #2 (Remarks to the Author):\nSecond review of: Inferring causal associations from cross-sectional Earth System Data\nI think the authors did a very nice job making improvements to this paper, and largely addressed\ncomments and suggestions from my own and another review. The paper reads a lot more clearly as\nthe method is more understandable, and the case studies are better presented as “tests” rather than\nnovel research findings using the GCCM method. There is also a clearer motivation for spatial causality\ntechniques like this. With that, I think this paper would be acceptable for publication after some\nminor/moderate revisions noted below. Mainly I feel like the figures could be improved – some figures\ncould be combined into more complex figures, and figures could also show the comparison of GCCM\nresults with the more traditional methods.\nGeneral comment: Throughout the paper, I found that past and present tense were combined, and\nwould recommend to use present tense for everything the authors did (e.g. We found, proposed,\ndiscovered, this was, these were --> present tense). However this is mainly stylistic.\nLine 89: Here a study on NPP is alluded to but not referenced (I think it gets referenced later in the\npaper)\nLine 122: Should this little paragraph be labeled “Results”? Maybe this is a journal specific section that\nrequires a short review of the results that are explained later.\nEquation 4: missing a parenthesis\nLine 266: Here you mention that X xmap Y is the basis to determine whether Y causes X (so Y --> X,\nor “Y causes X”). To me this is counter-intuitive, but useful to state explicitly like this. Later (page 18)\nthis aspect is repeated and it is noted that it may easily be mistaken – it might be helpful to re-iterate\nat different points in the paper such as in a figure caption (e.g. “Note that the higher Cu xmap\nIndustry indicates a causal effect from industry to soil pollution”)\nFigure 2: It seems like panel c in Figure 1 would go well with this figure – maybe Figure 1c and Fig 2a\nbelong together rather than separate figures?\nFigure 3: in panel (g), should Cd be Cu?\nLine 311: demonstrate\nLine 368: Here is a case where the value does not fit into the confidence interval? “0.48(p=0.00; 95%\nconfidence interval= [0.61,0.66])”.\nIn general, I think the p-values are useful, but the confidence intervals make the sentences difficult to\nread and don’t contribute to understanding the results. I would shift to presenting these results in\ntable form if you want to include all p-values and CI for each relationship. Alternately, it would be\nmost useful to have those values, and the values of linear correlation and LinGam metrics listed\ndirectly inside the white space in each figure panel. This would help illustrate the differences between\nthe metrics and whether the baseline models do or do not find a correlation for each case so the\nreader can refer to the figures and not have to sift through paragraphs for the statistical significance\nand comparison.\nLine 392: “Particulate Matters” should probably be “particulate matter”\nFigure 4: formatting issue with vertical line around panel (d). Also, it seems like this figure could be\ncombined into one potentially more informative figure, with 8 lines (e.g. 4 colors for the different\nvariables, and 2 line styles, or something similar). Finally, suggest to change “slop” to “slope” in figure\ntitle.\nFigure 4 discussion: Explanations are discussed for different possible mechanisms behind the finding\nof significant but smaller Y xmap Population (where Y = precip, temp, elev, or slope), but it should be\nnoted that GCCM doesn’t actually distinguish between any of these possible reasons (for example the\nfact that population could have a physical feedback on air temperature, but not so much with\nelevation).\nReference 53: some typo in this reference\nReviewer #1 (Remarks to the Author):\n• Compared to the earlier draft, language and grammar have clearly improved in the\nmain document. There are still several grammatical mistakes and imprecise\nformulations, some of which I point out below, but these affect readability much less\nthan before. Unfortunately, the supplement (especially the first section) is still in a poor\nstate language-wise and this does seriously affect intelligibility. This section needs to\nbe thoroughly edited for language structure and clarity, or with the help of software\ntools. In addition, the section feels rushed, and you use a lot of terminology without\nexplaining it. I would suggest that you either explain different concepts and\nassumptions in more detail, or you only mention their existence very shortly and refer\nto existing work for detailed explanations. Currently, you use an in-between approach\nwhich is often hard to follow.\nR: Thanks so much for your kind encouragement. We found your comments\nduring the first and second-round review are most professional and constructive.\nBy carefully revising the manuscript fully according to all your detailed comments,\nour manuscript has been largely improved. We are very glad to have you, clearly\na top mathematician who are familiar with causation models and CCM, as the\nmost suitable reviewer for this manuscript.\nSpecifically, thanks so much for suggesting us for further improving the\nsupplemental materials. Yes, we realized that some terminology should be clearly\nexplained to let readers have a better understanding of the background of casual\ninference methods. So according to your comments, we have added many details\nto the supplemental materials and carefully polished the English. Therefore, the\nrevised SI has been largely enhanced.\nThanks again for helping us improving this manuscript. Please feel free to contact\nus if additional revisions are required.\n• Comments 3,4 & 5: I like that you included the NPP – precipitation example and that\nyou evaluated the soil pollution example once more with different levels of added noise.\nSimilarly, adding LinGaM as an additional baseline strengthens your argument. These\nchanges improve the paper substantially in my opinion.\nR: Thanks so much for your encouragement. Yes, your suggestion on adding the\nNPP case and the LiNGam were great ideas. They improved readers’\nunderstanding of the wide suitability and advantage of GCCM and largely\nenhanced the original manuscript. Thanks again for all these constructive\ncomments.\n• Comment 6: The section on methodology that explains the theory behind GCCM has\nimproved and it was a good idea to move it to an earlier section. However, it is still hard\nto digest and sometimes too hand-wavy on the mathematics behind the method. I think\nyou should explain the mathematical details in much larger detail and give more precise\nformal definitions of the objects you employ. It would also help if some of the\nunderlying results, such as the generalization of Takens’ theorem were included, at least\nin the supplement. From a mathematician’s point of view, some choices of notations\nare also suboptimal. This is the most fundamental section of your article, and it warrants a more\ncareful and precise discussion.\nR: Thanks so much for this comment. We have followed your suggestions and\ndouble-checked all the notations in the manuscript. We found that you were indeed\nan excellent mathematician and your comments are precise and very valuable.\nAgain, we are very grateful to have you as the reviewer, who helped us improve\nthis manuscript significantly. We have rewritten the methodology section by\nadding more details of the method. We have carefully read the original works\nexplaining Takens’ theorem with time series and the generalized embedding\ntheorem to present a concise, yet precise and informative introduction of the two\ntheory in the methodology part. In addition, we have checked all relevant\nreferences and realized that you are quite right. We have carefully revised all\nnotations according to your comments and believe these equations are largely\nimproved.\nThanks again for these wonderful and professional comments.\n• Comment 7: I think there is a bit of a misunderstanding here. The fact that your method\nfinds a causal link nightlights -> soil pollution is a challenge to the interpretation of\nyour method’s output. Nightlights are not a cause of soil pollution in the interventional\nsense of Pearl: If we would intervene on nightlights by switching them off, this would\nnot reduce soil pollution by itself. Rather the correct causal graph is something like\nnightlights <- residential density -> soil pollution and the asymmetry you detect is a\nconsequence of the fact that the association between nightlights and residential density\nis very strong. If you applied your method on new data (X,Y) without much background\nknowledge and you find a link X->Y, you cannot claim that this is a causal link, but\nonly that X is sufficiently close to a cause X’ of Y. This is still useful information and\ndoes not invalidate your method, but it needs to be discussed much more thoroughly than what you currently do.\nR: Thanks so much for pointing this out. This is a very good point. We attempted\nto establish the casual links between residence and soil pollution. But due to the\nlack of residential density data, we chose the strategy of employing nightlight\nintensity data as the proxy variable. Actually, nightlight has been widely accepted\nas an effective reflection and employed as a proxy variable of residence density\n(Wang et al. ,2018; Tan et al.,2018; etc).Yes, you are right, the reason and validity\nof using the nightlight as the proxy value should be clearly explained. So in the\nrevised manuscript, we have added relevant references and revised the expression\nto avoid unnecessary misunderstanding. Meanwhile, we have referred nightlight\nas residence throughout the manuscript to better explain the causation.\nThanks again for this constructive comment, which made the methodology part\nmuch clearer.\nWang L, et al. Mapping population density in China between 1990 and 2010 using\nremote sensing. Remote Sensing of Environment 210, 269-281 (2018).\nTan M, et al. Modeling population density based on nighttime light images and\nland use data in China. Applied Geography 90, 239-247 (2018).\n• Comment 8: fixed.\n• Comment 9: I am glad that you added a significance test to your method. However,\nthis warrants more than two lines of text. The explanation on how you test is too hand-\nwavy in my opinion as this is crucial information. Please be more specific on the\nfollowing: what do you mean by ‘rho is statistically significant when the library\nbecomes largest’? Do you mean ‘rho is significantly different from 0’? What is your\ntest statistic for this and what is your choice of significance level? How do you measure\nconcretely whether there is an increase of rho with increasing size of L?\nR: Thanks so much for this comment. We added the details of the significance test\nand confidence interval estimation. The statistical significance means that rho is\nsignificantly different from zero when the library becomes largest. The\nsignificance level in the case studies is set to 0.05 and the confidence interval is\n95%. Following CCM, the increase of rho with increasing size of L is determined\nby plotting a line graph of rho vs L.\n• Comment 10 – 13: fixed.\nMinor comments:\n• For some reason, the abstract is written in past tense. I would suggest switching to\npresent tense. There are also some sentences in the introduction where you should\nchange to present tense, e.g. line 151.\nR: Thanks so much for this comment. Corrected.\n• Abstract: Causation inference in complicated systems -> Causal inference in complex\nsystems. Generally, the common term in the field is causal inference, not causation\ninference.\nR: Thanks so much for this comment. Revised, and we replaced all “Causation\ninference ” with “causal inference”\n• Abstract: split last sentence in two.\nR: Thanks so much for this comment. Corrected.\n• Line 44: the prerequisite -> a prerequisite\nR: Thanks so much for this comment. Corrected.\n• l. 58: ‘the strength of the cause’ -> the strength of the causal relationship\nR: Thanks so much for this comment. Corrected.\n• l. 60: ‘the randomized control trial’ -> randomized control trials\nR: Thanks so much for this comment. Corrected.\n• l. 67: ‘… which render these two frameworks unapplicable’. I find this too harsh, I\nwould rather say ‘… which are not captured well by the assumptions of these\nframeworks and therefore challenge their applicability.’\nR: Thanks so much for this comment and advising us how to improve it. We have\ncorrected it accordingly.\n• l. 69: ‘the spatial heterogeneity’ -> lose ‘the’\nR: Thanks so much for this comment. Corrected.\n• l. 107: again, lose ‘the’ before the theories\nR: Thanks so much for this comment. Corrected.\n• l. 133: this sentence is strange, what do you mean by ‘internal essences’?\nR: Thanks so much for pointing this out. We have corrected according to your\nsuggestion.\n• l.158: consists -> is\nR: Thanks so much for this comment. Corrected.\n• l.161: dynamic system -> dynamical system; different from -> in contrast to\nR: Thanks so much for this comment. Corrected.\n• l.178: lags of current orders is strange terminology, maybe explain this sentence more\nclearly.\nR: Thanks so much for this comment. We have revised it into “the first-order\nlags are adjacent units sharing common edges or vertexes with the focal unit, and the\nlags of next order are the first-order lags of those adjacent units (excluding those\nalready included)”.\n• L.224: missing bracket )\nR: Thanks so much for this comment. Corrected.\n• P.6: I still find your notation conventions quite confusing, especially the notation\nM_{x,si}. You explain a lot in words, but I would like to see cleaner mathematical\ndefinitions of the objects that you use. Why don’t you work more with the embedding\nmap \\Psi_h for instance? As far as I understand, the only thing you do here is to define\na metric d on the manifold M_x by d(w,z) = d’(\\Psi_h(w), (\\Psi_h(z)) for arbitrary points\nw,z on M_x, where d’ is an averaged sum 1/L \\sum d_s. over the Euclidean metrics d_s\non R^{|lag s neighbors|}. This is well-defined as \\Psi_h is an embedding.\nR: Thanks so much for this comment. As explained above, we are very grateful\nthat you gave us these very professional and helpful comments. We have re-\nchecked all required documents and revised all related notations according to your\nsuggestions. Again, thanks so much for this great comments.\n• Equation (6): with this definition of h, I don’t see why ||…|| defines a norm in a\nmathematical sense, and consequently why dis() defines a metric on M_x. More\nprecisely, can you explain why ||…||or dis() satisfy the property of positivity of metrics\nhere? Or is being a metric not essential? In the latter case, use a different notation than\n||…||as this strongly suggests that this is a norm in a mathematical sense. Again, it would\nhelp to have clearer notation in this section.\nR: Thanks so much for this comment. Again, we are very impressed and grateful\nfor giving us this constructive comments. Yes, we checked this and realized that\nyou are quite right. Accordingly, in the revised manuscript, we replaced ||…||\ninto an function called abs. Thank you very much for your reminder.\n• L.254: ‘it is a required condition’. Required for what exactly?\nR: Thanks so much for this comment. It is required for that the causal associations\nbetween X and Y are bidirectional or they are both the effects of a shared cause.\nThat it is a necessary condition for the fact after the “when” (the causal\nassociations between X and Y are bidirectional or they are both the effects of a\nshared cause). We have rephrased this sentence.\n• L. 448: in-slaved, do you mean enslaved?\nR: Thanks so much for this comment. Yes, it should be enslaved.\n• L. 526: ‘the periodic time series’ -> ‘periodic time series’\nR: Thanks so much for this comment. Corrected.\nTo Reviewer #2 (Remarks to the Author):\nR: Thanks so much for acknowledging our efforts on revising this manuscript\nand your very kind encouragement. We found your comments during the first and\nsecond-round review are most professional and constructive. By carefully revising\nthe manuscript fully according to all your detailed comments, our manuscript has\nbeen largely improved. We are very glad to have you as a very suitable reviewer\nfor this manuscript. Thanks again for helping us improving this manuscript.\nPlease feel free to contact us if additional revisions are required.\nGeneral comment: Throughout the paper, I found that past and present tense were\ncombined, and would recommend to use present tense for everything the authors did\n(e.g. We found, proposed, discovered, this was, these were --> present tense). However\nthis is mainly stylistic.\nR: Thanks so much for this comment. This is a very good suggestion. We have\nchecked throughout the whole paper and revised the tense fully according to your\ncomment. Thanks again for pointing this out, which improved the manuscript a\nlot.\nLine 89: Here a study on NPP is alluded to but not referenced (I think it gets referenced\nlater in the paper)\nR: Thanks so much for this comment. We have added the reference to the study\non NPP.\nLine 122: Should this little paragraph be labeled “Results”? Maybe this is a journal\nspecific section that requires a short review of the results that are explained later.\nR: Thanks so much for this comment. We labeled “Results” by referring to\npublished papers on Nature Communications. Thanks again for pointing this out.\nEquation 4: missing a parenthesis\nR: Thanks so much for this comment. We have corrected it accordingly.\nLine 266: Here you mention that X xmap Y is the basis to determine whether Y causes\nX (so Y --> X, or “Y causes X”). To me this is counter-intuitive, but useful to state\nexplicitly like this. Later (page 18) this aspect is repeated and it is noted that it may\neasily be mistaken – it might be helpful to re-iterate at different points in the paper such\nas in a figure caption (e.g. “Note that the higher Cu xmap Industry indicates a causal\neffect from industry to soil pollution”)\nR: Thanks so much for this comment. Yes, as you suggested in this-round and first-\nround, this part may cause confusion to readers. Therefore, according to your\ncomment, we tried to repeat this point at different parts throughout the revised\nmanuscript to remind readers the counter-intuitive results. Thanks again for this\ngood advice.\nFigure 2: It seems like panel c in Figure 1 would go well with this figure – maybe Figure\n1c and Fig 2a belong together rather than separate figures?\nR: Thanks so much for this comment. We removed Figure 1c in to Fig 2, and put\nthem together to illustrate the mutual neighbors of bidirectional and\nunidirectional associations.Thanks again for this good advice.\nFigure 3: in panel (g), should Cd be Cu?\nR: Thanks again for pointing it out. Yes, It is Cu, we have corrected it.\nLine 311: demonstrate\nR: Thanks so much for this comment. Corrected.\nLine 368: Here is a case where the value does not fit into the confidence interval?\n“0.48(p=0.00; 95% confidence interval= [0.61,0.66])”.\nR: Thanks so much for this comment. We have corrected it accordingly.\nIn general, I think the p-values are useful, but the confidence intervals make the\nsentences difficult to read and don’t contribute to understanding the results. I would\nshift to presenting these results in table form if you want to include all p-values and CI\nfor each relationship. Alternately, it would be most useful to have those values, and the\nvalues of linear correlation and LinGam metrics listed directly inside the white space\nin each figure panel. This would help illustrate the differences between the metrics and\nwhether the baseline models do or do not find a correlation for each case so the\nreader can refer to the figures and not have to sift through paragraphs for the statistical\nsignificance and comparison.\nR: Thanks so much for your good comment. It really helps. We have deleted all\nconfidence intervals with in the sentence in the main text, while readers can also\nfind them from the supplementary tables. And your suggestions on adding result\nof correlation and LinGAM to the figure is very useful, as it can give readers a\nvery clear picture of the comparison between GCCM, correlation analysis and\nLinGAM. Since correlation and LinGAM only provided one single value, instead\nof a figure, which may not be properly added to the GCCM figures with a good\nvisual effects. And we cannot directly list the value of correlation and LinGAM on\nthe figure of GCCM, otherwise it may cause confusion for readers when reading\nGCCM outputs. In this case, we added the value of correlation and LinGAM to\nthe detailed note in figures. Consequently, readers can immediately get the\ncomparison between GCCM, correlation analysis and LinGAM without checking\nSupporting information. Thanks again for this valuable comments.\nLine 392: “Particulate Matters” should probably be “particulate matter”\nR: Thanks so much for this comment. Corrected.\nFigure 4: formatting issue with vertical line around panel (d). Also, it seems like this\nfigure could be combined into one potentially more informative figure, with 8 lines (e.g.\n4 colors for the different variables, and 2 line styles, or something similar). Finally,\nsuggest to change “slop” to “slope” in figure title.\nR: Thanks so much for this comment. Yes, your suggestions can help us achieve a\nbetter explanation effect. Accordingly, we have revised to plot and the title. We\ncombined four plots into one, and foud it is not easy to distinguish them. As in the\nwhole paper, we keep the same style that the orange color is the main casual\ndirection, and the blue is not, we retained the four plots to make fast interpretation\nfeasible. Thanks again for you good advice.\nFigure 4 discussion: Explanations are discussed for different possible mechanisms\nbehind the finding of significant but smaller Y xmap Population (where Y = precip,\ntemp, elev, or slope), but it should be noted that GCCM doesn’t actually distinguish\nbetween any of these possible reasons (for example the fact that population could have\na physical feedback on air temperature, but not so much with elevation).\nR: Thanks so much for this comment. In the strong coupling cases, GCCM itself\ncannot distinguish the feedback causation or reflection. GCCM can only identify\nthe leading direction of the casual associations. To decide whether the weaker\ncoupling was a feedback causation or a reflection caused by the strong enslaved\neffects, sufficient prior-knowledge or further investigation are required. Yes, this\nis a very good point when employing GCCM. Accordingly, we have specifically\nadded this to the discussion part to give instructions to readers.\nThanks again for this constructive comments.\nReference 53: some typo in this reference\nR: Thanks so much for this comment. We have fixed typos in this reference.\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThank you for incorporating my suggestions into your work. Overall, I believe that the manuscript has\nimproved a lot since its original submission and may be suited for publication in Nature\nCommunications.\nI will address the three major open issues that I mentioned in my previous comments:\nMathematical notation: I am still not fully convinced by some of your notational choices, but I think\nthat the formal description of your methodology has improved enough to be acceptable. As a minor\nsuggestion, maybe give your abs() functions in equations (5) and (6) subscripts to make clear that\nthese are distinct from each other.\nThe role of proxies: I am still under the impression that you are not fully understanding the exact\nnature of my criticism here. I am <i> not </i> questioning that nightlights are a good proxy for\nresidential density. Instead, I am making the point that the use of proxies is a challenge to\ninterpretating the asymmetry that your method finds as <i>causation</i>. At least that is the case, if\nwe interpret causation as intervening on X induces a change in Y. In fact, we know for sure that\nintervening on the proxy in your example (X = nightlights) does not induce a change in Y (=\npollution). In this particular example, it might be clear what the actual underlying cause is (many\nresidents) but in other examples this might be much less clear. And these other examples may be\nexactly what people would like to use your method for! After all, we want to use it in cases where we\ndon’t know the causal relationship a priori already, as in your toy example!\nI do not expect you to write a philosophical text on the role of proxies in causality but at least it would\nhave been nice if this issue was acknowledged. However, if you prefer not to do this,<i> I will accept\nthat</i> since this point is not significant enough to keep your work from being published. I would ask\nyou though to reverse the change from nightlights to residence as this change is simply hiding the\nproxy issue, rather than addressing it.\nSupplement: The first section of the supplement has clearly improved, good work on that!\nReviewer #1 (Remarks to the Author):\nThank you for incorporating my suggestions into your work. Overall, I believe that the manuscript\nhas improved a lot since its original submission and may be suited for publication in Nature\nCommunications.\nR: Thanks so much for your encouragement and great help during three-round review, without\nwhich the quality of this manuscript cannot meet the high-standard of NC.\nI will address the three major open issues that I mentioned in my previous comments:\nMathematical notation: I am still not fully convinced by some of your notational choices, but I\nthink that the formal description of your methodology has improved enough to be acceptable. As\na minor suggestion, maybe give your abs() functions in equations (5) and (6) subscripts to make\nclear that these are distinct from each other.\nR: We have revised the notations in equations (5) and (6) according to your suggestions. Thanks\nso much for your comment in this-round and last-round review, which helped improve this\nmanuscript so much. .\nThe role of proxies: I am still under the impression that you are not fully understanding the exact\nnature of my criticism here. I am not questioning that nightlights are a good proxy for residential\ndensity. Instead, I am making the point that the use of proxies is a challenge to interpretating the\nasymmetry that your method finds as causation. At least that is the case, if we interpret causation\nas intervening on X induces a change in Y. In fact, we know for sure that intervening on the proxy\nin your example (X = nightlights) does not induce a change in Y (= pollution). In this particular\nexample, it might be clear what the actual underlying cause is (many residents) but in other\nexamples this might be much less clear. And these other examples may be exactly what people\nwould like to use your method for! After all, we want to use it in cases where we don’t know the\ncausal relationship a priori already, as in your toy example!\nI do not expect you to write a philosophical text on the role of proxies in causality but at least it\nwould have been nice if this issue was acknowledged. However, if you prefer not to do this, I will\naccept that since this point is not significant enough to keep your work from being published. I\nwould ask you though to reverse the change from nightlights to residence as this change is simply\nhiding the proxy issue, rather than addressing it.\nR: We are sorry that we did not fully get you in the last-round revisions. We have reverse the use\nof residence to nightlight in the revised manuscript according to your comment. Meanwhile, in\nthe revised manuscript (Discussion part), we have admittedly that we used a proxy variable in\nthe first case to reveal an already known causation, which should be carefully conducted for\nother scholars. We also suggested that there are some criteria to follow for properly selecting\nagent variables. Thanks so much again for your clarification.\nSupplement: The first section of the supplement has clearly improved, good work on that!\nR: Thanks so much for your encouragement."}]

Decision letter