2.1 From bivariate to multivariate measures of causality

Monitoring an ecosystem with continuous observations of net ecosystem exchange (NEE), the underlying gross primary production (GPP), and ecosystem respiration (R_eco) together with the relevant drivers, i.e. global radiation (R_g), surface air temperature (T), and soil moisture (SM), allows us to study the dynamics of the carbon cycle in terrestrial ecosystems. To foster its understanding, a fundamental question is how these variables causally depend on each other. This requires the identification of directional dependencies such as the well-known effects of SM → GPP and GPP → R_eco and their differentiation from physically implausible links such as R_eco → GPP. Graphical causal models (Spirtes et al., 2001) provide a framework to represent and identify causal relations based on conditional independence relations in data streams of this kind. In the case of an ecological monitoring site as described here, we can exploit the temporal information of the observations for identifying a time series graph as a type of graphical model (Runge, 2018). Formally this can be stated as follows. The variable $X_t^i$ is part of a multivariate stochastic process X (where i is the variable index, in the example $i\in \mathit{\{}{R}_{\mathrm{g}},T,\mathrm{GPP},R_{\mathrm{eco}},\mathrm{SM}\mathit{\}}$, and t is the time index). A time series graph 𝒢 visualises how the individual variables Xⁱ∈X depend on each other at specific time lags τ, i.e. $X_{t-\mathit{\tau }}^i$ with $\mathit{\tau }\in \mathit{\{}\mathrm{1},\mathrm{\dots },{\mathit{\tau }}_{max}\mathit{\}}$ (see Runge, 2018, for definitions). In the following, we refer to a variable $X_{t-\mathit{\tau }}^i$ that is causally affecting a variable $X_t^j$ as “parent” or “driver” and the latter as “receiver” or “target”. To come to a causal interpretation, it is important to exclude dependencies between two variables that are due to common drivers ($X_{t-{\mathit{\tau }}_{\mathrm{1}}}^i\leftarrow X_{t-{\mathit{\tau }}_{\mathrm{2}}}^s\to X_t^j$) or indirect paths ($X_{t-{\mathit{\tau }}_{\mathrm{2}}}^i\to X_{t-{\mathit{\tau }}_{\mathrm{1}}}^s\to X_t^j$). For instance, when estimating the effects of GPP on R_eco and R_g on R_eco using a bivariate measure, one likely obtains implausible results like links that are too strong or even unexpected links because T, respectively, as the common driver and mediator (indirect path), is not accounted for. To exclude dependencies due to common drivers or indirect paths, conditional independence tests are used, denoted as $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j|\mathcal{S})$, with some conditioning set 𝒮. If any variable (or their combination) in 𝒮 explains the dependence between $X_{t-\mathit{\tau }}^i$ and $X_t^j$, then the CI statistic is 0.

Two prominent methods that aim for directional dependencies are Granger causality and transfer entropy (Granger, 1969; Schreiber, 2000). Granger causality is typically estimated as a vector auto-regressive model and thus captures only linear links. Transfer entropy, based on information theory, also captures non-linear dependencies. It can be shown that for multivariate Gaussians, transfer entropy is equivalent to Granger causality (Barnett et al., 2009). Both can be phrased as testing for conditional independence (Runge et al., 2019a). In their original bivariate form, neither of these two methods accounts for third variables, but both can also be extended to deal with multivariate time series as required here (Runge et al., 2012; Granger, 1969). There are even non-linear and spectral modifications of Granger causality which have been applied to study biosphere–atmosphere interactions (Papagiannopoulou et al., 2017a; Detto et al., 2012; Claessen et al., 2019). However, the estimation of multivariate transfer entropy is challenging due to the “curse of dimensionality” (Runge et al., 2012), and multivariate Granger causality also exhibits low link detection power for larger number of variables (higher dimensions) and limited sample size, as is the case in our application (Runge et al., 2019b). The strong decrease in detection power happens when using the whole past ${\mathbf{X}}_t^{-}=({\mathbf{X}}_{t-\mathrm{1}},{\mathbf{X}}_{t-\mathrm{2}},\mathrm{\dots })$ of $X_t^j$, truncated at a maximum lag τ_max, as a conditioning set 𝒮. The problem is that this set can contain a high number of conditions that are irrelevant. For example, when assessing the effect of R_g at a specific time lag τ on GPP using multivariate Granger causality one would create a vector auto-regressive model comprised of all variables, i.e. R_g, T, SM, GPP, and R_eco, at each available lag, but R_eco, dominated by heterotrophic respiration, is not expected to affect gross primary productivity and could be removed to decrease the dimensionality. However, manually selecting conditions is not desirable when the underlying dependence structure is unknown, which is why ideally the conditioning set is identified automatically.

PCMCI addresses this issue by reducing the set of conditions 𝒮 prior to quantifying the dependence between two variables. The two-step approach utilises a variant of the PC algorithm (Spirtes and Glymour, 1991) and the momentary conditional independence measure (MCI) (Runge et al., 2019b). More detailed descriptions are given in Sect. 2.4 and 2.5, respectively (a full description of PCMCI including proofs and quantitative comparisons with other methods are provided in Runge et al., 2019b). A schematic of the PCMCI approach is given in Scheme A1 in Appendix A. PCMCI belongs to the family of causal graphical models (Spirtes et al., 2001; Pearl, 2009) and follows the assumptions listed in Sect. 2.2. In the limit of infinite time series length, PCMCI converges to the true graph of dependencies, which is why we use the term “causal”. As we deal with finite sample length and partially unfulfilled assumptions, spurious links can still appear (beyond the expected false positive rate) and therefore each detected link has to be interpreted with caution.

2.2 Assumptions

PCMCI assesses the causal structure of a multivariate dataset or process X by estimating its time series graph. To draw causal conclusions from observational data, any causal method must adopt a number of assumptions (Pearl, 2009; Spirtes et al., 2001). For the time series case, here we assume time order, the causal Markov condition, faithfulness, causal sufficiency, causal stationarity, and no contemporaneous causal effects (Runge et al., 2018). PCMCI is applied in combination with the ParCorr linear independence test based on partial correlations (see Sect. 2.3). This application additionally requires stationarity in the mean and variance and linear dependencies. In the following, we briefly discuss these assumptions (further details in Runge, 2018; Runge et al., 2019b).

The time order within the time series allows us to orient directed links which are only pointing forward in time. This accounts for causal information propagating forward in time only, i.e. the cause shall precede the effect. Therefore, a directed causal link $X_{t_{\mathrm{1}}}^i\to X_{t_{\mathrm{2}}}^j$ can only exist between two nodes $X_{t_{\mathrm{1}}}^i,X_{t_{\mathrm{2}}}^j$ if t₁<t₂. When a contemporaneous link is found, i.e. t₁=t₂, it is considered to be undirected. In ecological language it means that in order to claim that R_g is driving GPP, any change in R_g that is affecting GPP must be measured at a time step before the change in GPP occurs. The causal Markov and faithfulness assumptions relate the underlying physical causal mechanisms to statistical relationships manifest in the data. The causal Markov condition states that if two processes are not directly connected by some physical mechanism, then they should be statistically independent conditional on their direct drivers, like R_g and R_eco conditional on T. The faithfulness assumption concerns the other direction: if two processes are statistically independent, then there cannot be a direct physical mechanism. The causal sufficiency assumption implies that every common cause of two or more variables Xⁱ∈X is included in X. If this is not the case, detected links may be indirect or due to an unobserved common driver. However, the absence of a link in the detected graph still implies that no direct link is present (as this only requires the assumption of faithfulness). For example, R_g is expected to influence R_eco via T, the indirect path. However, a link between R_g and R_eco might be detected if T is not included in the analysis. However, a missing link between T and R_eco might indicate conditions inhibiting respiratory processes, i.e. very cold temperatures with frozen surfaces or very dry conditions with dead vegetation coverage. Causal stationarity refers to the existence of links over time. In a deciduous forest, for example, the ecosystem's CO₂ exchange is not causally stationary as the link R_g → GPP is given in summer but not in winter. Formally, a process X with graph 𝒢 is called causally stationary over a time index 𝒯 if and only if for all links $X_{t-\mathit{\tau }}^i\to X_t^j$ in the graph the condition $X_{t-\mathit{\tau }}^i\textcolor{red}{}{X}_t^j|{\mathbf{X}}_t^{-}\backslash \mathit{\{}{X}_{t-\mathit{\tau }}^i\mathit{\}}$ holds for all t∈𝒯.

2.3 Independence test

At the core of PCMCI there are conditional independence tests $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j,\mathcal{S})$ to evaluate whether $X_{t-\mathit{\tau }}^i\textcolor{red}{}{X}_t^j|\mathcal{S}$ given a conditioning set 𝒮. Within the PCMCI software package TIGRAMITE (Runge, 2017), several independence tests are implemented. Here, we focus on the linear independence test called ParCorr. The ParCorr conditional independence test is based on partial correlations and a t test. This assumes the following model:

with coefficients β and Gaussian noise ϵ. This leads to the following residuals:

with estimated $\widehat{\mathit{\beta }}$. ParCorr removes the influence of 𝒮 on Xⁱ and X^j via ordinary least-squares regression and tests for independence of the residuals using the Pearson correlation with a t test. The independence test returns a p value and test statistic value I, i.e. the correlation coefficient in case of ParCorr. Thus, to identify the effect of GPP on R_eco that does account for their common driver T∈𝒮, ParCorr will perform a linear regression of T on both GPP and R_eco accounting for time lags. The p value of the residuals' partial correlation test can be used to assess whether the two variables are dependent.

2.4 PC algorithm

To efficiently estimate $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j|\mathcal{S})$ the conditioning set 𝒮 should be as small as possible, which means that it should only contain relevant conditions, which allows us to isolate the unique influence of $X_{t-\mathit{\tau }}^i$ on $X_t^j$. For an estimation of $\mathrm{CI}({R_{\mathrm{g}}}_{t-\mathit{\tau }},{\mathrm{GPP}}_t|\mathcal{S})$, for example, 𝒮 should contain T and SM (at certain lags), as they influence the ability of an ecosystem to perform photosynthesis. Likewise, when estimating $\mathrm{CI}(T_{t-\mathit{\tau }},{\mathrm{GPP}}_t|\mathcal{S})$, 𝒮 should include R_g and SM for the same reasons. A sufficient set of relevant conditions includes the drivers and parents of the variable $X_t^j$. Consequently, the aim of the PC step is to identify an as small as possible superset of the parents of each variable included in the process. The algorithm uses a variant of the PC algorithm (Spirtes et al., 2001); a comprehensive pseudo-code of this procedure is given in the Supplement of Runge et al. (2019b). In the limit of infinite sample size, the relevant conditions indeed converge to the true causal parents, although, practically, an estimate that contains a few irrelevant conditions, like R_eco, is sufficient as well.

The PC step starts by initialising the whole past of a process: $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)={\mathbf{X}}_t^{-}=\mathit{\{}{X}_{t-\mathit{\tau }}^i:i=\mathrm{1},\mathrm{\dots },N,\mathit{\tau }=\mathrm{1},\mathrm{\dots },{\mathit{\tau }}_{max}\mathit{\}}$. Next, by evaluating $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j,\mathcal{S})$, conditions $X_{t-\mathit{\tau }}^i$ are removed from $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$ that are independent of $X_t^j$ conditionally on a subset $\mathcal{S}\in \stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)\backslash \left\{X_{t-\mathit{\tau }}^i\right\}$. 𝒮 starts as the empty set ∅ and is iteratively increased. For instance, let $X_{t-\mathit{\tau }}^i$ be R_eco (at a specific lag) and $X_t^j$ be GPP. The conditional independence between GPP and R_eco will be estimated first by using no conditions. If GPP and R_eco appear related, one variable will be included in the conditioning set. If the residuals are still dependent, a second variable is included and so on. When T is part of the conditioning set, the residuals of GPP and R_eco might not be dependent anymore and R_eco is removed from the estimated set of parents of GPP. The PC algorithm adopted in PCMCI efficiently selects those conditioning sets to limit the number of tests conducted.

Every conditional independence test is evaluated at a significance threshold α_pc, which is usually set to a liberal value between 0.1 and 0.4. Alternatively, in TIGRAMITE one can let α_pc be unspecified. PCMCI then evaluates the best choice of α_pc $\in \mathit{\{}\mathrm{0.1},\mathrm{0.2},\mathrm{0.3},\mathrm{0.4}\mathit{\}}$ based on the Akaike information criterion (AIC), which is further explained in Runge et al. (2018).

2.5 MCI tests

MCI is the actual causal discovery step that ascribes a p value and strength to each possible link. MCI iterates through all pairs $(X_{t-\mathit{\tau }}^i,X_t^j),i=\mathrm{1},\mathrm{\dots },N,\mathit{\tau }=\mathrm{0},\mathrm{\dots },{\mathit{\tau }}_{max}$, and calculates $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j,\mathcal{S})$, where 𝒮 consists of the two (super-)sets of parents $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$ and $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_{t-\mathit{\tau }}^i\right)$ obtained in the PC step. $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_{t-\mathit{\tau }}^i\right)$ is constructed by shifting the time series of $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^i\right)$ by τ. In cases where $X_{t-\mathit{\tau }}^i\in \stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$, $X_{t-\mathit{\tau }}^i$ has to be removed from $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$. If τ=0, conditional dependence is estimated for contemporaneous nodes $X_t^j$ and $X_t^i$. Due to missing time order, a dependence would be left undirected. Further, as the parents $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^i\right)$ and $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$ used in each conditional dependence test are defined to lie in the past of $X_t^i$ and $X_t^j$, links, both contemporaneous and lagged, can be spurious due to contemporaneous common drivers or contemporaneous indirect paths. The absence of a link, though, means that a physical (contemporaneous) link is unlikely (assuming faithfulness; see Runge et al., 2018). For simplicity, the previously given examples were omitting the time lag. Thus, if R_eco responds instantaneously (considering the sampling temporal resolution) to changes in T but T responds with a time lag to R_g, both variables will likely appear contemporaneously coupled to R_g.

The link strength in the PCMCI framework can be given by the effect size of the conditional independence test statistic measure CI used in combination with MCI. In case of ParCorr, the effect size is given by the partial correlation value, which is between −1 and 1. Assuming a linear Gaussian model the partial correlation value is shown to directly depend on the receiver's and driver's variance as well as the coupling coefficient (Runge et al., 2019a):

where ${\mathit{\sigma }}_{X_{t-\mathit{\tau }}^i,X_t^j}$ are the variances of the noise terms driving $X_{t-\mathit{\tau }}^i$ and $X_t^j$, respectively, and c is their coupling coefficient. In practice, non-linear links can often also be well detected with ParCorr in so far as they are linearly approximated. In case the linear part is even “stronger” than the non-linear part, ParCorr might also have a better detection power than a non-linear independence test (Runge et al., 2018).

2.6 Data

3.1 Test model

As an example, in Fig. 1 we show PCMCI and lagged correlation networks in the form of process graphs. It is clearly visible that many more spurious links pass the significance threshold of 0.01 using lagged correlation as compared to using PCMCI. Those spurious links can complicate the analysis or lead to false conclusions and misleading hypotheses. We examined four cases of different time series lengths: 91 [$\mathrm{183}\le \mathrm{doy}\le \mathrm{274}$] and 120 [$\mathrm{153}\le \mathrm{doy}\le \mathrm{274}$] d and 1 and 5 years for daily data (doy: day of the year). For each time series length and each parameter set, the causal network structure was estimated for 100 realisations of the model (each based on a realisation of intrinsic noise), which allowed the estimation of false positive (FPR) and true positive (TPR) detection rates. The detection rates are calculated for each tower, FPR in general and TPR link-wise. The TPR for each link is its sum of detections among 100 realisations divided by 100. The FPR is the number of falsely detected links divided by the number of all possible false links and 100. The summary of the experiments, i.e. the overall false positive rate (FPR) and the distributions of the link's true positive rate (TPR) across sites are given in Figs. 2 and 3, respectively. The blue violin plots always report the case of normal distributed (non-heteroscedastic) intrinsic noise and the corresponding orange violin plots summarise the case of heteroscedastic noise. The effect of heteroscedasticity and seasonality is then assessed by comparing the distributions obtained from the baseline dataset to the results of the seasonality dataset.

Figure 1The artificial datasets are generated with a prescribed interaction structure (true network), which is obtained by fitting the test model to the FLUXNET sites. Here we show for four time series lengths the process graphs estimated via both lagged correlation and PCMCI. The data used stems from the homoscedastic realisation of the seasonality dataset of the Hainich site. The significance level was set to 0.01. The number of time lag labels were limited to five in the correlation networks. However, for the longest time series typically the whole range of lags (0–25) was significant.

Download

Figure 2The distribution of false positive detection rates estimated for the baseline dataset, the seasonality dataset and the anomalised seasonality dataset (mean seasonal cycle subtracted). The distributions are given for different time series length (number of datapoints). Additionally, the distributions are split to show the impact of heteroscedastic noise (orange) compared to normal distributed noise (blue). The significance level of 0.01 is given by a blue horizontal line.

Download

Figure 3The distribution of the true positive detection rates for each link in the test model estimated for the baseline dataset, the seasonality dataset, and the anomalised seasonality dataset. The distributions are given for different time series lengths (number of datapoints). Additionally, the distributions are split to show the impact of heteroscedastic noise (orange) compared to normally distributed noise (blue).

Download

The FPR of homoscedastic time series in the baseline dataset is in the expected range of 0.01, the chosen significance level, indicating a well-calibrated test due to fulfilled assumptions. The assumption of stationarity is violated as soon as heteroscedasticity or seasonality is present. The effect on the FPR is an increase above 0.01 for time series length of 1 and 5 years, with a much stronger increase due to seasonality (factors of 4 and 8, respectively) than for heteroscedasticity (factor of 2).

The effect of non-stationarities on the TPR differs among the links. The detection of linear links (ℛ_g→𝒯 and 𝒯→𝒯) is not affected by seasonality and slightly improves for heteroscedastic dynamical noise. The detection of non-linear links is improved by seasonality, with the strongest effects in the link 𝒯→𝒢𝒫𝒫. The link 𝒯→ℛ_eco has a stronger non-linearity, and therefore the detection rate shows a weaker effect on seasonality. Furthermore, the coupling coefficient c₅, the base of T, can be close to or be exactly 1 (see Fig. A3). This would actually cause, on the one hand, 𝒯→ℛ_eco to vanish, rendering a detection impossible and, on the other hand, result in a linear dependence of 𝒢𝒫𝒫 on ℛ_eco, which improves its detection. Heteroscedasticity seems to have a slight negative effect on non-linear links. In general, the TPRs in the seasonality dataset are quite high, even for non-linear links, and predominantly above 80 % and often reaching 100 %.

Comparing the TPRs of the non-linear links shows some disparity. The links 𝒯→𝒢𝒫𝒫 and 𝒯→ℛ_eco experience zero detection in the baseline dataset but partially considerable rates in the seasonality dataset with a strong dependence on the time series length. On the contrary, the median of the TPR of the links ℛ_g→𝒢𝒫𝒫 and 𝒢𝒫𝒫→ℛ_eco is above 95 % in the seasonality dataset, even for time series length as short as 91 d, but remains high in the baseline dataset. The removal of the seasonal cycle keeps the TPRs largely unaffected but reduces the FPR. Nevertheless, it still remains above the significance level by factors of 4 and 2 for 5-year and 1-year time series lengths, respectively.

In summary, the seasonality dataset exhibits high TPR even for non-linear links. Compared to stationary time series, the detection of non-linear links actually benefits from seasonality. The high detection, though, comes at the cost of a high false positive rate for time series length of and above 1 year. To a certain degree, the increase in FPR can be counteracted by anomalisation.

3.2 Majadas de Tiètar dataset

First, we look at the link consistency by comparing networks that were obtained for each tower within a month. The comparison is done for 2 months with strongly differing climate conditions: April and August. In Fig. 4 we compare the estimated link strengths (effect size estimated via partial correlation) as long as the corresponding links are significant in at least one network. The confidence intervals are overlapping for the majority of links, suggesting that the uncertainty of the fluxes is much smaller than the observed effects (El-Madany et al., 2018). Exceptions are found for only a few links ($R_{\mathrm{g}}\stackrel{\mathrm{0}}{⟶}T$, $R_{\mathrm{g}}\stackrel{\mathrm{1}}{⟶}LE$, $\mathrm{VPD}\stackrel{\mathrm{1}}{⟶}\mathrm{VPD}$, $H\stackrel{\mathrm{2}}{⟶}H$, $\mathrm{NEE}\stackrel{\mathrm{0}}{⟶}LE$, $H\stackrel{\mathrm{2}}{⟶}\mathrm{NEE}$, $\mathrm{NEE}\stackrel{\mathrm{4}}{⟶}\mathrm{NEE}$; the number above the arrow indicates the lag) where the detection rates do not or barely overlap. Cross-links (a link from one variable to another) with two or more significant appearances are predominantly at zero lag. Approximately half of the links with lag 1 are auto-dependent links (a link from the past of a variable to its present). Comparing the links between the months of April and August, distinct differences can be noticed. First, August has slightly fewer significant links compared to April. Second, the only links remaining that are significant in two or three towers are between atmospheric variables. Third, the remaining link strengths tend to be weaker in August than in April.

Figure 4Comparison of the networks of three eddy covariance measurement stations (LMa, LM1, LM2) located in Majadas de Tiètar (Spain). Links that are found to be significant in one of the three networks are included. For each link, the calculated strength of all three networks is plotted with its 90 % confidence interval. The colours blue, orange, and green correspond to the towers LMa, LM1, and LM2, respectively. The significance threshold is 0.01. If a link does not pass the significance, it is marked by a black dot. The links are grouped into lag 0 (top), lag 1 (middle), and all lags greater than 1 (bottom). Negative NEE is associated with carbon uptake by the ecosystem. Links at lag 0 are left undirected (–), yet as R_g is set as main driver, links incorporating R_g at lag 0 are directed (→).

Download

The difference among the seasons is further investigated in Fig. 5, which shows process graphs for each month of the year. We combine the networks of the three towers to one process graph by plotting every link that is significant in at least one tower. The process graphs in Fig. 5 visualise clearly gradual changes within the interaction structure of the biosphere–atmosphere system during the course of a year. During the main growing season from February to May, NEE is coupled strongly to the energy fluxes latent (LE) and sensible heat (H). These connections weaken, disappear, or even switch sign during the dry season. Less regularly, NEE also shows connections to radiation (R_g) and temperature (T). Between the atmospheric variables, a basic network between VPD, T, R_g, and H remains intact and relatively constant in strength. The dominance of contemporaneous links is found as well, as already seen in Fig. 4. Aside from the decoupling of NEE from any variable in the dry period, there are additional interesting patterns. For example, the positive reappearance of the link between NEE and LE in September. Here, the onset of precipitation events (see Fig. A10) occurred, which led to strong respiration peaks (Ma et al., 2012). Creating such a network via lagged correlation would result in much more significant links (causing the network to be not interpretable as opposed to PCMCI) and NEE does not decouple from the atmosphere in August (see Fig. A7).

Figure 5To visualise the gradual changes in interaction structure the networks of the three towers are combined for each month. The number of significant occurrences of a link is given by its width. The link strength, given by the link colour, is calculated by averaging the significant links of the towers. The link's lag is shown in the centre of each arrow, sorted in descending order of link strength. The resulting graphs are shown for April 2014 until March 2015. The significance threshold is 0.01. The networks of April and August, illustrated in Fig. 4, are highlighted by a box.

Download

The above results demonstrate that PCMCI is sensitive enough to capture seasonal differences and certain physiological reasonable biosphere behaviour. Moreover, PCMCI yields a better interpretable network structure than pure correlation approaches.

3.3 Gridded global dataset

The significant lags and MCI values of each climatic variable on NDVI were subject to inspection. In Fig. 6 the maximal MCI value and the corresponding lag are plotted for the links $X_{t-\mathit{\tau }}\to \mathrm{NDVI}:\mathit{\tau }\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}$, with X being one of the climatic drivers: radiation, temperature, or precipitation. The chosen significance threshold was set to 0.05. Figure 7 shows the climatic driver with the largest MCI per grid point. PCMCI detects a regionally varying influence of climatic drivers. As expected, the boreal regions are strongly driven by temperature instantaneously, while (semi-)arid regions, which correspond predominantly to grass- or prairie-dominated areas, respond strongest to precipitation at a time lag of 1 month. Radiation is found to have a comparatively low spatial effect, with hotspots in southern and eastern China, central Russia, and eastern Canada.

Figure 6Influence of climatic drivers on NDVI as calculated by PCMCI. The first column shows the estimated causal influence given as maximal absolute MCI value of climatic drivers on NDVI. The second column gives the time lag at which the maximal absolute MCI value occurs (in months).

Figure 7Map of the strongest climatic driver (largest absolute MCI value) per grid point.

The dominant lags are found to be 0 and 1. Only a very small fraction of the total area shows a maximal MCI value at a higher lag of two or three months. The lags are also not equally distributed among the climatic drivers. Radiation and temperature are predominantly strongest at lag 0, while precipitation has a much larger fraction of area showing the strongest response at lag 1. Regions where the impact of R_g on NDVI is strongest at lag 1 tend to respond negatively to R_g but positively to precipitation at lag 1. On the other hand, a large part of the regions with the strongest impact of precipitation at lag 0 respond negatively to it but positively to radiation.

In summary, PCMCI estimates coherent interaction patterns that match well with anticipated behaviour based on vegetation type and prevailing climatic conditions.

4.1 Lessons learned from the test model

The probability of detecting a link with PCMCI depends strongly on a link's MCI effect size, which is larger for strong variance in the driver and a low variance in the receiver (see Sect. 2.5). Several results can be explained by this observation. First, the variance of three out of five drivers of cross dependencies in the test model are either directly or indirectly (via 𝒢𝒫𝒫) influenced by ℛ_g, which has the highest variance of all variables. Consequently, the detection power of the three links is large, almost 100 %. In comparison, the other variables' variances are weaker, since they are influenced by 𝒯, which results in a lower detection power. This is the origin of the disparity in detection rates of the non-linear links. Second, the partially strong increase in TPR of non-linear links (influenced in a multiplicative way by 𝒯) from the baseline dataset to the seasonality dataset can also be explained by this increase in variance. A multiplicative link is actually not generally expected to be found by ParCorr (Runge, 2018), but the value of the multiplicative factor is dominated by the seasonal value, and not the dynamical noise, which might cause a scaling of the dynamical noise terms rather than a random distortion. Third, the dependence on the variance ratio can also explain the difference in TPR between homoscedastic (equal error variance) and heteroscedastic (error variance changing over time) time series, i.e. the variance of ℛ_g and 𝒢𝒫𝒫 exhibits a strong seasonality, with its peak in summer, while the variance of 𝒯 is rather constant. This explains, for example, the strong decrease in TPR for the link 𝒯→𝒢𝒫𝒫 at 91 d time series length when comparing homoscedasticity to heteroscedasticity. The decrease in TPR is less pronounced when another season, implying a different variance, is chosen for this comparison. As links with weak driver variance and strong response variance are more likely to be missed, one may ask what effect this will have on the detection of feedback loops where one variable has low and the other high variance. Here lies a limitation of the test model where no feedback loops were implemented.

Seasonality and heteroscedasticity constitute violations of the stationarity assumption underlying the independence test ParCorr. Seasonality constitutes a common driver in this model. In general, such common drivers increase the dependence among the variables and hence lead to a higher detection rate for true links (TPR) as well as a higher false positive rate (FPR) for absent links if this driver is not conditioned out properly. This additionally causes the TPR and the FPR rate to increase in the seasonality model. As shown in Runge (2018), including the cause of the non-stationarity as an exogenous driver in the analysis allows PCMCI to regress out its influence on the other variables. However, for ParCorr this is only valid if the dependence on the non-stationary driver is linear. Therefore, the regression on ℛ_g fails for 𝒢𝒫𝒫 and ℛ_eco in the test model. With this ill-posed setting, the probability to detect false links increases with increasing time series length or when more periods are included. Stationarity in the mean is obviously also not fully guaranteed when subtracting the seasonal mean. Here we observe that the FPR stays above the significance level for the anomalised seasonality dataset. One can ask whether the FPR stays above the significance threshold because subtracting the seasonal mean does not remove the heteroscedasticity. However, we attribute this high FPR to a not fully removed seasonality since the FPR of both homoscedastic and heteroscedastic time series decreases by roughly the same amount in the anomalised seasonality dataset and the effect of heteroscedasticity is rather weak in the baseline dataset. The increasing FPR with increasing time series length can further raise doubts regarding the analysis of long time series. For such an analysis, though, the assumption of causal stationarity should first be assessed. For example, the link from radiation to GPP vanishes in winter as there is mostly no active plant material left. To account for causal stationarity, the analysis should be limited to time series sections where the causal structure is expected to be similar. This is typically done by limiting the analysis to a specific time period (i.e. “masking”), e.g. a specific season, month, or time of the day. Such masking additionally further reduces influences of remaining seasonality or heteroscedasticity. One can argue, as in Peters et al. (2017), that the causality of a system is invariant even between seasons because the physical mechanism is the same in all seasons. Yet, while the physical, i.e. functional relationship might be constant over time, its imprint in the time series might vary. For example, a functional dependence f(x) might be “flat” for small values of x and linearly increasing for larger values. If only small values occur in the winter season, then the link is absent, while it “appears” only in the summer season. Across all seasons, this can be considered a non-linear functional dependence f(x). In practice, restricting an analysis to different seasons can help in interpreting the mechanism and is performed here using a linear framework.

Summarising the results of the test model, the different detection rates, disparity among non-linear links, and the detection of multiplicative links are largely explainable via the effect of the variance on the link detection. Yet, the discussion revealed the need for further research in several aspects. On the one hand, feedback loops are not included in the test model yet are an important aspect in natural systems. On the other hand, removing non-stationarities is essential to keep the false positive rate in the expected range, but standard procedures of subtracting the mean seasonal cycle are not sufficient. Further, the effect of non-stationarity on the causal network structure needs to be investigated.

4.2 Causal interpretation of estimated networks from observational data

In both the half-hourly time-resolved eddy covariance data and the monthly global dataset, the predominant type of dependence found is contemporaneous. PCMCI leaves these undirected since no time order indicating the flow of causal information is available. Further, as discussed in Sect. 2.2, contemporaneous common drivers or mediators are not accounted for. The consequence is that both spurious contemporaneous and spurious lagged links can appear if they are due to contemporaneous variables. For interactions that are contemporaneous in nature, since they occur on considerably shorter timescales than the time resolution, PCMCI is not the optimal choice regarding a causal interpretation and other methods should be applied in conjunction (Runge et al., 2019a). Further, we faced a trade-off between fulfilling causal assumptions and detection power. In practice, accounting for causal stationarity (by limiting the analysis to certain periods of the dataset) means decreasing the number of available data points but accounting for causal sufficiency leads to an increase in dimensionality by adding variables and increasing the maximal lag. Both will lead to a decrease in detection power, which can affect the network structure. PCMCI alleviates the curse of dimensionality by applying a condition selection step, but one still cannot indefinitely add more variables. Another important factor that affects detection power and dimensionality is the time resolution. There are several points in favour and against increasing time resolution. On the one hand, increasing time resolution can resolve contemporaneous links and potentially increases the detection power due to an increased number of datapoints. On the other hand, the dimensionality increases if the maximal lag is adapted. Further, causal information might be split apart and distributed over more lags, rendering the links at each individual lag less detectable. This can cause links to disappear, but links can also appear if new processes are resolved at a higher timescale. Finally, observational noise (measurement errors) might be larger in higher-resolution data than in aggregated data, as it is averaged out in the latter and thus affects link detection less. Consequently, when comparing network structures based on different settings, i.e. maximal lag, included variables, time resolution, and considered time period, the (dis)appearance of single links among specific variables can stem from several factors, i.e. a change in detection power, a changed (conditional) dependency, or due to a common driver. These factors together with a non-zero false positive detection rate are challenging for a causal interpretation. Therefore, detected links should be interpreted with care and can give rise to new hypotheses and analyses involving further variables. Generally, a causal interpretation is more robust regarding the absence of links (see Sect. 2.2). Specifically, it does not require that all common drivers are observed.

Nevertheless, robust patterns were identified in our studies that are also consistent with other studies. Furthermore, a causal analysis has the advantage of an enhanced interpretability compared to correlative approaches. First of all, we could show that the networks' estimated link strengths are consistent for observational data, even though measurement error affects the data. The dataset used was suitable for this analysis, as the measurement stations are located in a reasonably homogeneous ecosystem that shows only little spatial variation (El-Madany et al., 2018). Thus, also the interaction between biosphere and atmosphere is expected to change only marginally across space within this ecosystem. Second, the gradual changes in plant activity that are taking place in the ecosystem of Majadas de Tiètar throughout the year do emerge very well in the coupling strength of daytime NEE to the atmospheric variables. The observed decoupling during the dry season is in accordance with the one of a soybean field during drought conditions observed by Ruddell and Kumar (2009). The gradual changes in ecosystem activity are not visible in a pure (lagged) correlation analysis or are only visible in colour or density changes but the large number of significant links prevents any detailed interpretation on the physical mechanisms and changes thereof. The large number of significant links compared to the PCMCI networks stems solely from the absence of conditioning in common drivers or mediating variables, which often leads to an overestimation of the link strength in correlation networks. As a result, processes, such as the decoupling of NEE during the dry period, stay hidden. To reduce the effect of confounding, analyses often utilise partial correlation (see, e.g. Buermann et al., 2018). However, a partial correlation can introduce new dependencies (as opposed to removing them) if one conditions on causal effects of the variables under consideration (the “marrying parents” effect). This issue is avoided in PCMCI by only conditioning on past variables. Additionally, PCMCI chooses only relevant variables as conditions by applying the PC condition selection step, which is especially valuable in high-dimensional study cases and improves detection power and computation time (Runge, 2018).

The global study of climatic drivers of vegetation shows a general pattern of lags and dependence strengths of vegetation on climatic variables that is easily interpretable. The boreal regions appear energy limited and especially driven by temperature (see Fig. 6c), while the strongest dependence of (semi-)arid regions on precipitation reflects their limitation in water supply. Two recent studies performed a similar analysis. Both Wu et al. (2015) and Papagiannopoulou et al. (2017b) investigated lagged effects and dependence strengths of NDVI on precipitation, temperature, and radiation. Wu et al. (2015) estimated the lags of the strongest effects via an univariate regression of the climatic drivers on NDVI and subsequently used those lags to fit a multivariate regression model of the climatic drivers on NDVI and determined their relative effects. Papagiannopoulou et al. (2017b) applied a non-linear Granger causality framework utilising a random forest predictive model; the method was presented by Papagiannopoulou et al. (2017a). We recognise that similar patterns are observed in Wu et al. (2015), but the lags at the maximal MCI value are usually lower than the one found in Wu et al. (2015), which stems from the methodical differences. Besides having used anomaly values, PCMCI regresses both NDVI and the climatic drivers on their parents before calculating the MCI value (see Sect. 2.5). This especially removes the influence of autocorrelation. Runge et al. (2014) shows how autocorrelation affects the correlative lag causing it to be larger for stronger autocorrelation; thereby, the correlative lag may become larger than the causal lag. Therefore, according to Fig. 6b, the causal information embedded in monthly resolution is predominantly received within 1 month. Finding the strongest causal links at a time lag up to 1 month appear in agreement with Papagiannopoulou et al. (2017b). Additionally, the spatial distribution of the strongest climatic influences compares well, but there are certain noteworthy differences that do not necessarily stem from masking differences, i.e. that we took only values belonging to the growing season while Papagiannopoulou et al. (2017b) took the whole time series. First, there is little significant Granger causality of water availability found in boreal regions while there are significant negative causal dependencies detected via PCMCI. Second, NDVI in arid regions is not (or is only partially) Granger caused by radiation and temperature but in part shows a negative PCMCI value for those variables. There might be physiological reasons that can explain the PCMCI patterns, i.e. waterlogging or temperatures that are too high. To explain the differences though, we could identify two possible reasons. First, Papagiannopoulou et al. (2017b) masked out negative influences of radiation, arguing that radiation is not negatively affecting NDVI. They found that negative influences of R_g are usually a consequence of poor conditioning on other variables. Second, a precipitation event in boreal regions coincides with a reduction in radiation and temperature. Boreal regions usually do not suffer from water shortages. Thus, they respond more strongly to the reduction of radiation and temperature than precipitation. As precipitation is coupled negatively to radiation and temperature at lag 0, the effect of precipitation on NDVI is found to be negative. Thus, the link $P\stackrel{-}{\to }\mathrm{NDVI}$ might be an effect of the contemporaneous common driver scheme $P\stackrel{-}{\to }{R}_{\mathrm{g}}\stackrel{+}{\to }\mathrm{NDVI}$ and therefore would not be causal. In fact, a similar argumentation can be given for the negative impact of temperature and radiation on NDVI in arid regions.

In summary, we pointed out the need for careful interpretations in applying causal discovery methods and especially highlighted the challenges linked to the study of biosphere–atmosphere interaction via PCMCI. We demonstrated that the network structures estimated from observational data are explainable with respect to plant physiology and climatic effects. Finally, our study shows that causal methods can deliver better interpretability and a much improved process understanding in comparison to correlation and bivariate Granger causality analyses that are ambiguous to interpret since they do not account for common drivers.

4.3 Outlook

The preceding discussion has shed light on the merits of PCMCI, as well as the challenges of applying causal discovery methods. Runge et al. (2019a) discuss further challenges and methods and give an outlook of how multiple methods can be combined to alleviate limitations.