Interactive comment on “ On the data-driven inference of modulatory networks in climate science : an application to West African rainfall ”

May 24, 2014 Review of “On the data-driven inference of modulatory networks in climate science: an application to West African rainfall” by Gonzalez et al. General comments: This contribution employs a portfolio of statistical approaches including LASSO, CHARM and Dynamic Bayesian networks that are applied to climate data sets. The aim of this manuscript is to apply novel statistical techniques that have been developed and utilized on many domain applications and are extended here to issues of climate. These approaches are sequentially applied to a set of fixed climate data sets obtained from various climate data centers. There is an assessment of the pathway significance


Introduction
The climate system is inherently complex, due to the existence of non-linear interactions, or couplings, between its subsystems (e.g., the ocean and the atmosphere), global scale temperature anomalies (e.g., El Niño-Southern Oscillation), and other climate behaviors.Such a system exhibits hierarchical modularity of its organization and function (Havlin et al., 2012): each constituent subsystem performs a simi-lar function and does not act in isolation; instead, they interact or cross-talk.The challenge is to discover the key subsystems and their cross-talk mechanisms; that is, the pos-35 itive and negative feedbacks that collectively modulate the dynamic behavior of the system through a sophisticated network of modulatory pathways that ultimately define the system's functional response.
For example, the rainfall anomaly in the Sahel region of 40 western Africa, which is the focus of this study, represents a "functional response" for the climate system, which is in actuality the predictant of a model (such as Lasso, DBN, etc.).Rainfall in the Sahel is dependent on global Sea Surface Temperature (SST) patterns, as well as on local climate 45 variability.There is a multitude of complex associations between various subsystems that drive the Sahel's climate response mechanisms.Some of these associations have been discovered throughout more than two decades of hypothesisdriven and/or first-principles based research.These associ-50 ations include a diverse range of climate mechanisms.For example, warmer temperatures in the Mediterranean Sea region lead to increased evaporation, and southward moisture advection in the lower troposphere toward the Sahel (Rowell, 2003).On a more global scale, the Atlantic Multidecadal 55 Oscillation (AMO) displaces the Intertropical Convergence zone (ITCZ) further northward, bringing more moisture to the Sahel region (Zhang and Delworth, 2006).The North Atlantic Oscillation (NAO) has been linked to the moisture budget in Northern Africa (Hurrell, 1995), through a direct influ-60 ence on the Sea Level Pressure (SLP), although this mechanism remains underexplored.In the Pacific, a warm ENSO event is associated with enhanced trade winds over the tropical Atlantic and weaker moisture advection over West Africa, consistent with a weaker monsoon system strength (Janicot 65 et al., 2001).Figure 1 illustrates an overview of the climate : modulatory network, which is a collection of modulatory pathways, with some mechanisms driving rainfall in the Sahel known to be directly/indirectly associated, and some not fully understood.Comprehending these mechanisms is particularly important due to the influence of rainfall variability in the region.Severe drought occurred throughout the 1970's and 1980's, leading to severe disruption of agriculture and major food shortages (Mortimore and Adams, 2001).Dry conditions (low rainfall anomaly) also lead to the spread of meningitis, as in wet conditions, higher humidity during both the spring and summer seasons strongly reduce disease risk by decreasing the transmission capacity of the bacteria (Sultan et al., 2005).These issues make the Sahel particularly vulnerable to fluctuations in rainfall, and provide motivation to improve domain scientist's knowledge of the contributing factors (Tetteh, 2012).
For mechanistic understanding of functional responses such as African Sahel rainfall, we posit that a data-driven approach may facilitate the discovery of key players that might cross-talk by identifying candidate modulatory pathways and/or suggesting new factors and relationships with the proper characterization of their inductive or suppressive roles.The goal of our approach is to elucidate the putative modulatory pathways that suggest cross-talking mechanisms controlling a system's functional response.More specifically, given the key climate drivers and their modulatory directions on the response, we must infer (a) the putative pathways of modulatory events (e.g., Pacific ENSO → AMO → Sahel Rainfall in Fig. 1) and (b) the modulatory signs (e.g., induction vs. suppression, such as a positive anomaly sign of EATL, EATL HIGH being related to the negative anomaly sign of Sahel rainfall, Rainfall LOW ) that collectively define the network of modulatory pathways for the response.Furthermore, given that there is a variety of methodologies that can be used to find such modulatory relationships, we must provide a consensus result that accounts for all evidence of a given relationship.To the best of our knowledge, this is a novel proposition in the field of knowledge discovery in the physical science domain, in general, and climate extremes (e.g., droughts), in particular.Moreover, this data-driven approach could contribute, in the long run, to the identification and characterization of more comprehensive and predictive models of the physical phenomenon under study.

Methods
In our previous work, we proposed an approach for the aforementioned data-driven, semi-automatic inference of phenomenological physical models based on Lasso multivariate regression (Pendse et al., 2012).This approach was applied to quantify the influence of key factors on the Sahel rainfall anomaly.The results obtained enabled the formulation of the North Atlantic Oscillation (NAO)-driven hypothesis, among others, which theorizes that the NAO modulates the drivers of West African climate, the Atlantic Dipole and the EATL, via the low-level westerly (LLW) jet.

120
We extended this work by developing coupled heterogeneous association rule mining (CHARM), which allowed us to mine higher-order couplings of climate relationships and to capture the anomaly phases with which each climate factor is related to each other (e.g., a negative anomaly of LLW 125 may be related to a positive anomaly of EATL, and the presence of both factors may be associated with a negative Sahel rainfall anomaly) (Gonzalez et al., 2013) (Sect. 2.1).Such relationships are not typically captured from modulatory inference frameworks, let alone traditional association rule min-130 ing (ARM) methodologies.
Here we propose to extend CHARM by incorporating other existing methodologies, namely Lasso multivariate regression (Tibshirani, 1994) (Sect.2.2) and Dynamic Bayesian Networks (Murphy, 2002) (Sect.2.3), as comple-135 mentary approaches to increase the confidence of the inferred modulatory relationships.Moreover, in order to obtain a consensus as to which of the relationships identified have the most evidence of being present, we treat the results of each methodology as individual pieces of evidence in an informa-140 tion fusion approach, and combine them into a unified, coherent result.This unified result can provide us a means with which to increase the confidence of the relationships identified throughout the different methodologies.This should allow us to contrast the methodologies by studying how each of 145 their results differ, and to correlate these results with known relationships found in literature.Furthermore, the application of this unified result to the climate network may allow the identification of previously-undiscovered relationships, which can then be analyzed from a traditional climate per-150 spective.In the Section 3 (specifically Table 3), we present the application of such a method to the climate indices affecting Sahel rainfall.CHARM is an extension of ARM that enables the discovery of climatologically-relevant modulatory pathways from spatio-temporal climate data.The traditional ARM methodology CHARM is based on is presented in Sect.2.1.1,and the limitations of ARM that CHARM aims to address are described in Sect.2.1.2.

Traditional Association Rule Mining (ARM)
Traditional ARM was pioneered by Agrawal et al. (1993) as a methodology for capturing the frequency with which two items are present within transactions in market basket data.
For instance, Fig. 2 presents a set of transactions that indicate whether or not an item was purchased.ARM takes this information and organizes each transaction as a combinatorial set of items  3.
The aforementioned relationship can be equally represented as {Bread → Milk} or {Milk → Bread}, utilizing an arrow to capture that the presence of the antecedent (e.g., the items on the left side) implies the consequent (e.g., the items on the right side) will be present with the given support.However, if we use a metric such as confidence, which captures the conditional probability of the consequent being present given the presence of the antecedent, the direction of the rule carries more weight.For example, {Bread → Milk} has a confidence of 0.75 (of the four times bread is present, milk is only present thrice) while {Milk → Bread} has a confidence 1 (bread is present every time milk is present) (Tan et al., 2006).As such, the metric used to measure the interestingness of mined rules affects their overall interpretation and should be selected carefully (Sect.2.1.4)(Tan et al., 2001).
ARM is an increasing area of interest for domain sciences, because of the growing need to mine data to identify the cooccurrence of important events (Agrawal and Srikant, 1994;Tan et al., 2001).ARM, unlike other methodologies for inference of phenomenological models, takes into account the latent but vital signals embedded in the intermediary pathways associated with the system's functional response.However, the application of ARM to spatio-temporal climate data puts forth a series of challenges.In Sect.2.1.2,we outline these challenges and describe CHARM as a means to address 210 them.

Coupling of climate indices
Due to the complexity of the climate system, building comprehensive models over climate data is not trivial, in part because of the interactions between its subsystems, the dimen-215 sionality and structure of its underlying data, and the quality of such data.
The key drivers of a climate system are spatially distributed and active at different temporal phases in the modulatory network of the system's functional response.State-of-220 the-art mining methodologies are not well-equipped to handle such diversity of spatio-temporal alignment between the system's features.For example, due to the transactional nature of the ARM methods, each spatial grid point (specified by latitude, longitude, and/or altitude) at a given point in time 225 (e.g., month and year) defines a transaction ID, or a row in the transaction matrix.This requires that all features be aligned with respect to their transaction IDs, complicating the use of multi-resolution, multi-variate, spatio-temporal climate data by these methods.For this reason, we leverage climate indices, known to be a valid abstraction of the underlying subsystem's zonal climate behavior (Hallett et al., 2004), thus significantly reducing the number of features needed to capture spatial data.However, these climate indices that capture data for different subsystems are located in different parts of 235 the globe (see Fig. 3).Some climate variables from observations or simulations (e.g., SST) are defined only over the ocean; yet others (e.g., rainfall) are defined over land.Hence, considering both features as columns in a transactional matrix is impossible, 240 given that they have no common grid points.Even if they share some spatial region, they are often still not perfectly aligned due to variation of their grid resolutions.While math-: ematical methods (e.g., interpolation or extrapolation) exist to facilitate data alignment, they introduce uncertainty and instability, affecting the interpretability of the results (Gonc ¸alves, 2002).Subsequently, when a new feature is integrated into the study, realignment and the aforementioned mathematical operations must be performed again.
An ARM-based approach for the discovery of relationships among climate variables was proposed by Tan et al. (2001).This approach studies only spatially-aligned datasets, and affixes climate indices alongside them.However, due to its inherently grid-based nature, this approach assigns each climate index's locally observed anomalies to all grid points.
That is, it assumes that the anomaly equally affects the entire globe.While such an assumption may hold for some climate drivers, it can increase the number of false positives due to its inherent amplification bias.For example, in the representative case shown in Fig. 4, the high anomaly of the SOP index (SOI-HI) and the low anomaly of the NP index (NP-LO) would occur in so many transactions (see Sect. 2.1.1)that it would be present in every itemset for time t 1 , which complicates understanding the information gained from any resulting rules that includes them.
Climate scientists often study climate factors in a coupled manner and relate certain variables to others.Traditional lagged climate techniques employed for coupled pattern analysis include singular value decomposition, grid point correlations, among others (Polo et al., 2008).For example, Principal Component Analysis (PCA) has been used to determine the relationship between Indian Ocean Dipole and East African rainfall (Schreck and Semazzi, 2004;Manatsa et al., 2012).Hence, we adopt a similar approach by coupling climate indices.We take a quotient of these relationships be-275 fore identifying any anomaly, to capture the anomaly in the relationships between these variables.
For each climate index λ to be used as a coupling listener, we iterate through all other climate indices δ as coupling inciters, and calculate their ratio, δ/λ, as a data coupling that intends to capture the behaviors of the logical sentence "how abundant is δ given the presence of λ?"An issue in calculating these ratios is the potential emergence of large values due to the denominator possibly being orders of magnitude smaller than the numerator.To handle this, we normalize the 285 resulting data couplings such that they range between −1 and 1, allowing us to avoid wide-ranging quotients that could affect the abstraction of anomalous events (Sect.2.1.3).
We note that each tuple (row in the database) now represents a specific coupling inciter λ and time, while each col-290 umn represents a particular coupling listener δ, and each cell contains the relevant data coupling value.For ARM, this data must be binned, after which the resultant dataset cells indicate the presence or absence of anomalies for each previously calculated coupling, described further in Sect.2.1.3.We ad-295 dress the increase in dimensions this leads to in Sect.2.2.1, by using only the most prominent temporal phases in order to reduce the search space.

Identifying anomalous events
As suggested by NOAA (2014) and based on our interac-300 tions with climate scientists, we identify anomalies as any set of values below the 33.33 rd percentile or above the 66.67 th percentile for any given variable.Given that the data was normalized before calculating ratios, we identify the anomalies using the aforementioned norm, based on the phase-wise 305 groupings.We take each tuple corresponding to each unique combination of δ and λ, and identify high anomalies as being those ratios in the upper 66.67 th percentile, and low anomalies as being ratios in the lower 33.33 rd percentile.
Since we are trying to identify the presence or absence 310 of anomalous events, we divide each column into two separate high and low cases, and assign a binary 1 when either anomaly occurs, and a 0 otherwise.This results in a very sparse matrix, as no particular year can fit in both high and low categories, and it is likely the majority of years have most 315 variables falling into a non-anomalous category.
Figure 5 represents a particular (i, j)th iteration.In this example, for the coupling of λ i = EATL 8 and δ j = NAO 3 , we identify high and low anomalies and assign transaction IDs that indicate that the cells pertain to the coupling of that 320 year's data for the listener (shown in the column header)  against the stated inciter.This transaction ID shows that each row in the matrix consists of the anomaly of the ratios calculated over each possible coupling δ x /λ i , where x indicates that all values in this row were divided by the same λ i .

CHARM pathway significance assessment
As mentioned before, rule interestingness in ARM is estimated using metrics that quantify the importance of each rule.Selecting which metric to use depends on the information to be obtained (Tan et al., 2001).Support and confidence are commonly used to measure rule quality, and although there are other possible metrics to measure interestingness, none is regarded a "catch-all" for high quality rules (Tan et al., 2001).Such metrics also require predetermined thresholds to cut off rules deemed "uninteresting", which reduces the accuracy in retention of significant rules.Hence, we first prune based on bare-minimum thresholds for support and confidence (the rule appearing in a single transaction), and if any given rule needs to be further analyzed (based on domain knowledge or otherwise), we perform a Monte Carlo simulation to test against the null hypothesis of observing the rule at random.Thus, we define a rule to be significant and interesting if it meets the following criteria: p value ≤ 0.01, support ≥ 6%, and confidence ≥ 75 %.
This criteria constrains the search space, and trims the result space, pruning unimportant rules.Once a set of possibly interesting rules is identified, the computationally more demanding, but embarrasingly parallel, statistical significance test is applied to further prune insignificant rules.On average, this removed 20-30 % of the generated rules from the result set.

Coupling heterogeneity
Data coupling creates a large set of transactions, covering each year studied for each possible coupling inciter.Rules that only have sufficient support when counted over multiple inciters would be difficult to interpret, thus we must heterogeneously generate rules for each coupling inciter separately, as shown in Fig. 6.This allows us to identify preferential bias towards a particular coupling inciter, and preserves information relating to each data coupling individually.

CHARM computational complexity
Finding all frequent itemsets for ARM is an NP-complete problem that when bounding transaction length, becomes linear with complexity O(r • n • 2 l ), where n is the transaction count, l is the maximum itemset length, and r is the number 365 of maximal frequent itemsets (Zaki, 2000).Rules are generated such that a user-specified minimum confidence and/or a minimum support is satisfied.Thus, for an itemset of length k, there are 2 k − 2 potentially confident rules, making the complexity O(c•2 q ), where c is the number of frequent item-370 sets, and q is the length of the longest frequent itemset (Tan et al., 2001;Zaki, 2000).CHARM leverages the Apriori algorithm, ensuring only maximal frequent itemsets are considered for rule generation (Agrawal et al., 1993), while leveraging sequential ARM 375 to identify such relationships across different temporal instances (Huang et al., 2008).As mentioned in Sect.2.1.5,each inciter is studied heterogeneously.Thus, the method operates in smaller parallel executions with low overhead.

Lasso multivariate regression 380
Least absolute shrinkage and selection operator (Lasso) multivariate regression is an approach pioneered by Tibshirani (1994) that takes a set of inputs and an outcome measurement and fits a linear model, seeking to shrink the regression and sparsify the predictor feature space.This is achieved by 385 constraining the L 1 norm of the β parameter vector B = {β 1 , β 2 , . . ., β n }, calculated as in Eq. ( 1), such that it is no greater than a given s value to be minimized (Tibshirani, 1994).
In the context of this study, this process highlights the prominent phases of the features (Sect.2.2.1).It derives the temporal phases of predictors lagged behind a response of interest, generating predictor coefficients indicating the mag-395 nitude and type of the modulatory relationships with said response (Pendse et al., 2012).Recent work on inference of modulatory relationships based on Lasso multivariate regression of temporal and spatio-temporal data includes means to improve upon the 400 Lasso methodology.We apply the method proposed by Pendse et al. (2012), given that it incorporates prominent phase detection and significance assessment.Pendse et al. (2012) presents an approach toward a data-driven, semiautomatic inference of phenomenological physical models 405 based on the Lasso multivariate regression model and quantifies the influence of key "players" on the response of interest (e.g., Sahel rainfall anomaly) through use of the Expected Causality Impact (ECI) score.The work presented in Pendse et al. (2012) also proposes methods for search space pruning, significance estimation and impact analysis that provide quantifiable metrics in terms of predictors' contributions to the rainfall variability and their probability of detections (PODs).

Prominent phase detection 415
We employ the methodology suggested by Pendse et al. (2012) to identify the most prominent phases (i.e., seasons) in the data.For the benefit of reproducibility, we utilized the supplemental material provided therein (Pendse et al., 2012).The results obtained by Pendse et al. (2012) were consistent with many well-known modulatory relationships from prior climate knowledge (Chang et al., 2006;Marshall et al., 2001;Sutton et al., 2000).These results complement the existing physical models and may help climate scientists categorize the correct season for the response of interest (e.g., Sahel rainfall variability).Hence, by leveraging these prominent phases (shown in Table 1), we can focus on features that should have a stronger influence over the response.

Lasso pathway significance assessment
To assess pathway significance, we follow the method de-430 scribed in Pendse et al. (2012).That is, we apply the Monte Carlo method to estimate the statistical significance of the relationships found between the input features and the response in terms of the null hypothesis, by iteratively permuting the response and performing Lasso multivariate regression for 435 this permuted data.This method allows us to prune insignificant edges in the Lasso network, represented by higher p values.

Lasso computational complexity
Given Lasso's iterative nature in finding appropriate λ and 440 β values, the computational complexity of Lasso is reliant on the q parameters and n observations provided by the source data, as it would need to attain solutions for all subsets M k , k ∈ 1, . . ., m (Meinshausen, 2007).Hence, the computational complexity of this methodology 445 is O (n • q • min {n, q}) (Meinshausen, 2007).Furthermore, since all variables must be at some point evaluated as the Lasso response r ∈ q, this is multiplied by a factor q.However, the q value that affects the actual Lasso execution would also grow smaller, as considerations are made to remove q 450 parameters that temporally cannot modulate r.

Dynamic Bayesian networks
DBNs expand upon Hidden Markov Models (HMMs) and Kalman Filter Models (KFMs), indexing instances of arbitrary variables.DBNs are represented as a structure similar 455 to that of Bayesian Networks with the added benefit of incorporating the temporal space (Dean and Kanazawa, 1989;Murphy, 2002).DBNs are a very popular means with which to mine and represent modulatory relationships in spatial and temporal data, given that the conditional probability distri-460 bution of each node can be estimated independently (Friedman et al., 1998;Murphy, 2002;de Kock et al., 2008).The model's dynamicity is obtained by combining a traditional Bayesian network with a temporal Bayesian network that allows for capturing behaviors of the Bayesian network over 465 the temporal space, and is not to be confused with the idea that the model changes over time (Murphy, 2002).

DBN pathway significance assessment
To assess pathway significance, we again apply the Monte Carlo method to estimate the statistical significance of edges 470 representing modulatory relationships.This affects the computational complexity of the methodology, as each random combination must be mined individually.Hence, to somewhat alleviate this matter, we verify only columns for which relationships were found by the base method, and omit the 475 features for which no relationships were found.: 7

DBN computational complexity
Several different implementations of DBN inference exist, each with varying degrees of complexity.To mine the DBNs for our problem, we leverage the toolkit provided by Zou and Feng (2009), built upon the design proposed in Murphy (2002).This toolkit allows us to infer the network structure of the DBNs in O(T ), where T is the length of the sequence to be mined, which could be exponentially large depending on the number of possible feature combinations a se-485 quence could contain (Murphy, 2002).Given our utilization of the toolkit, we abide by this complexity for our estimation, only restricting the execution by disallowing temporallyinfeasible edges (i.e.edges are only allowed between two nodes if the originating node occurs temporally before the destination node).In doing so, we ensure that the directionality of the network is temporally sound and fits proper modulatory relationships.

495
in a different manner, affecting the interpretability of the information they provide.Hence, the resulting relationships between climate factors should be structured such that all possible modulatory pathways are captured in a comparable context, while preserving the information given by each 500 method.
The results provided by DBN capture a network of relationships between climate factors as a directed, acyclic graph (DAG).Such a graph includes a set of vertices and directed edges, which in the context of this study represent the climate factors and the relationships between them, respectively.Furthermore, there are no cycles in the graph (i.e.following a path originating at a given node will never lead back to that node).This structure provides an intuitive visualization of the behaviors in the system, as each edge represents the existence of a relationship between the climate factors it connects.Therefore, we adapt the results provided by CHARM to a similar structure, by building a network where the edges represent each possible combination of high/low anomalies, directed from antecedent to consequent.However, given that CHARM uses coupled climate indices, we must ensure the the networks generated for the three methods can be equally interpreted.Hence, each Lasso and DBN experiment will also use such coupled data and will also be executed heterogeneously, as described in Sect.2.1.5.This allows us to directly use the results provided by DBN, given that it already adheres to the proper network structure.As for the results from the Lasso experiments, we generate the network of modulatory pathways by drawing a directed edge from vertex A to B when a β coefficient was found for 525 an execution where B was the response and A was a predictand.Given the temporal window constraints set upon this problem, we can follow the graph backwards from our de-sired response to study all relationships, both direct and indirect.530

Consensus modulatory network inference via information fusion
To infer a consensus modulatory network for a functional system response, we must combine the modulatory networks inferred by CHARM, Lasso multivariate regression,

535
and DBN into a single unified network that captures the consensus of the results.The field of evolutionary biology has leveraged methods related to information fusion to combine evidences found for specific gene classifications in collected field data (Bailey and Gribskov, 1998;Li et al., 2008).Of 540 the methods in this field, we chose to combine the resulting p values of each edge for each modulatory inference methodology by overlaying the resulting graphs from each methodology upon one another, and performing Fisher's combined probability test, shown in Eq. ( 2) where p i represents the p value for the ith independent test.This score presents a large χ 2 test statistic when p i values are smaller, suggesting the null hypotheses are not true for every 550 test.In contrast, when all the null hypotheses are true, and the p i are independent, χ 2 has a chi-squared distribution with 2k degrees of freedom, where k is the number of tests being combined.This can then be used to determine the p value for χ 2 (Fisher, 1932).

555
After obtaining the combined p value, we compute an ARM-inspired support count to quantify the amount of methods providing evidence for this result.With this, we determine which edges are worthwhile of inclusion, opening the realm for climate scientists to determine which amount of ev-560 idence constitutes a satisfiable minimum for which an edge is acceptable, and additional information can be obtained from the underlying individual results.For example, if ARM found some A HIGH → B LOW relationship between features A and B, Lasso or DBN also found evidence of some A → B 565 relationship, and we obtain a significant Fisher statistic, we can state that this relationship is founded, since 3 out of the 6 possible method results have evidence of such relationship.Furthermore, given the information provided by ARM's result highlighting specific phases, domain scientists can inves-570 tigate the A HIGH → B LOW relationship in further detail.
We use these statistics to determine a consensus in the relationships found by the methods employed and build a network to capture said consensus.Note that the number of relationships in the consensus result will not be restricted by 575 the methodology that identifies the fewest relationships.Instead, each methodology serves as evidence for the consensus result and affects the strength of the evidence provided for a particular relationship.Hence, each method contributes to the consensus result, with no specific methodology act-580 ing as a determining factor that would bias the result towards that specific methodology.Determining a bound at which to remove rules from the consensus result based on the amount of evidence provided is a topic of future work.

585
Setting the Sahel rainfall anomaly as the system's response presents an ideal model for assessing the climatological relevance of modulatory pathways identified by CHARM.We evaluate the computational validity of CHARM, using the criteria defined in Sect.2.1.4,by studying the mined rules 590 for the connections it identifies, and compare the results to those of the other approaches and study the effect of combining the rules as described in Sect.2.5.

Data
Table 1 presents the data used for this study, which was obtained from the NCEP/NCAR (2014), along with rainfall data obtained from the UDel AirT Precip dataset provided by NOAA/OAR/ESRL (2014).Along with these, our study uses new indices created by climate scientists (items 8-9) using empirical orthogonal function (EOF) techniques (Wilks, 2006) to isolate the dominant mode(s) in reanalysis data (PSD, 2014).The inclusion of the new indices is based on the fundamental knowledge that Sahel climate is modulated by different climatic drivers (Hurrell, 1995;Rowell, 2003).

605
These drivers originate from the ocean, atmosphere, land surface, and vegetation, where they interact intricately, and ultimately exert a strong influence over the Sahel region.However, the tropical Pacific, the Atlantic, and the Indian Oceans, as well as the Mediterranean Sea and the overlying atmosphere are key drivers of Sahel climate, so the creation of such indices ensures they are given equal chance to participate in the experiment, represented in Table 2. Hence, where climate literature suggested a teleconnection between a given climate variable and Sahel rainfall, but a representative cli-615 mate index for it was not readily available from NCAR, an EOF analysis of the 850mb height field was created instead, using reanalysis data.Each mode is represented as feature # (i.e., the first mode of variance over the Mediterranean Sea is referred to as MSEA1).
We select the most prominent seasons for these indices, as described in Sect.2.2.1, and utilize eastern Sahel rainfall over the season July-August-September (JAS) as the desired response.These variable-phase combinations are denoted as feature phase , where the subscript corresponds to temporal 625 phases (i.e., 1 = Jan-Feb-Mar, 2 = Feb-Mar-Apr. . ., 12 = Dec-Jan-Feb).Thus, we can contrast these to provide a description of the sub-region's climate variability in association to this response for the period of 1950-2008.

Network interpretations 630
We will discuss specific use cases of our experiment herein.Images for the generated networks for each method, their associated DAG matrices and combination metrics can be obtained via supplemental material1 .

Rainfall interconnections 635
Table 3 captures the relationships known from reference material in contrast to the findings of the evaluated methods for the EATL 8 coupling inciter, which implies that for the given data couplings, the coupling inciter used was Atlantic ENSO in the temporal phase of August-September-October..This 640 table serves to present length of pathway from the variable in question to the expected rainfall response, so as to verify the findings of this experiment in terms of known relationships gathered over two decades of research (Tetteh, 2012).The dynamical substance in the processes involved and tele-645 connections in these mining techniques is highlighted.
Lasso reveals that four oceanic modes, the Pacific (represented by MEI, Nino 3.4 and Nino 4), IOD, MSEA3, and AMO influence the Eastern Sahel rainfall (Rowell, 2003).The dynamical processes inferred from warm ocean sur-650 face anomalies associated with the IOD (Lu, 2009), MSEA3 (Rowell, 2003) and AMO (Zhang and Delworth, 2006) are related to an increase in the magnitude of the main rainfall season in the Sahel.The IOD and MSEA3 specifically facilitate positive moisture advection whereas the AMO displaces 655 the Intertropical Convergence Zone (ITCZ) to its climatological position over the Sahel.These mechanisms are tied to moisture transport from the tropical Atlantic by LLW2 and LLW3.On the contrary, a warming of the Pacific is generally associated with rainfall diminution over the Sahel Jani-660 cot et al. (1996).
role or impact is modulated directly or indirectly by distinct phases of Nino3, MSEA1 and MSEA2.While warm (cold) phase of Nino3 suppresses (enhances) moisture flux over the Sahel, MSEA1 and MSEA2 have a competing effect, with positive and negative moisture transport over the Mediterranean Sea respectively, and are involved in negative and positive moisture transport over the Mediterranean Sea.The model also reveals high (low) phase of LLW1 over the Atlantic is associated with strengthening (weakening) of westerly moisture flow.This co-occurs with GHT 1, 2, and 3, which determine troughs and ridges that govern high and low rainfall anomalies.

680
Lasso and DBN coincide in capturing extra-tropical NAO forcing.Although the NAO is known to impact Sahel rainfall (Hurrell, 1995), the mechanism by which this occurs is unclear.A link to the tropical Atlantic, particularly through the LLW's, is suggested by the results here.It is possible that the 685 moisture flux from the tropical Atlantic is dependent on the phase of the NAO.On a finer scale, the model also predicts a direct link to the NAO.The association between the NAO and Sahel rainfall may be multifaceted, and our results are being further investigated by the authors based on an NAO-690 driven hypothesis over the entire West African Sahel (Tetteh, 2012).
Figure 7 captures the aggregate number of relationships found by each method, averaged across all coupling inciters studied.We find that per each coupling inciter, Lasso 695 is the more sensitive methodology, finding edges for most possible feature combinations.Given the design of our CHARM experiment capturing phase-specific relationships (i.e., A HIGH → B LOW ), as described in Sect.2.1.3,we group all possible high/low combinations to merely visualize how 700 many relationships CHARM found as a whole.We note that these highly coincide with the findings of Lasso, but find their share of unique relationships that contribute to the final result.Lastly, we find that DBN produces very few relationships, but the majority of these contribute to the Fisher 705 statistic for the three methods, as 97 % of the found relationships coincide with either CHARM or LASSO, while 60 % coincides with both.The central area of Fig. 7 would lead to further study, as it indicates all three methods provided evidence of relationships in this area.

710
Given the intent is to find drivers for the rainfall feature, Fig. 10 captures the average number of direct relationships to the rainfall response found by each method.This again highlights the sensitivity of Lasso, as of a maximum 20 features, an average 17.67 were selected, whilst CHARM and DBN 715 find less such relationships.When evaluating the coupling inciter EATL 8 (see Fig. 9), we see this in further detail, as while Lasso captures 17 direct relationships to the Rainfall response, CHARM and DBN capture 8 and 1, respectively.Hence, Lasso appears to detract from discovering indirect re-720 lationships, unless β values are inspected directly.This especially affects the fused network, as most features are marked as directly associated with the response (see Fig. 8).  4 highlight that AMM has the highest betweenness centrality in the case of CHARM and the fused network, meaning it is found in all shortest paths in the network, marking AMM's importance in the network, as 10 : Fig. 8. Resulting combined network for coupling inciter EAT L8 stated in Sect.3.2.1.The clustering coefficient for CHARM and DBN have standard deviations σ = 0.02 and σ = 0.11, respectively, but σ = 0.01 when fused together, indicating Lasso's sensitivity plays a key role in the fused result, as in the individual Lasso graph, most nodes were hubs for many inbound connections.Furthermore, we find that the fused network performed best at creating a cluster with the desired Rainfall response, seemingly best influenced by the CHARM network, although given Fig. 9, we know the influence of Lasso again came into play.Hence, exploring a means to limit Lasso's influence may be a beneficial next step for future work.

Conclusions
We evaluated three different methods for finding modulatory relationships in spatio-temporal climate data and validated 750 the results obtained against known relationships modulating rainfall in the Sahel region of Western Africa.These results show that each method has its benefits and drawbacks.Noting that significant changes had to take place for utilizing CHARM for this purpose given spatial alignment issues, we  devised data coupling as a means with which to study the relationships in the underlying data.These changes served to make CHARM an efficient methodology for addressing the data-driven discovery of predictive, climatologically relevant, and statistically significant modulatory pathways in the 760 physical model of the Sahel rainfall anomaly.We also evaluated the consensus network obtained after combining the results of these methods via information fusion.In any case, this study served to validate these methods against known relationships from over two decades 765 of hypothesis-driven and first-principles research.The IOD, ENSO, MSEA and AMO were confirmed as important SST anomalies modulating rainfall in the region, as previously discussed in the climate literature.The relationship with the NAO is found to have both direct and indirect components, and is particularly related to equatorial westerlies (LLWs) in the Atlantic, known to influence the region (Nicholson, 2009).It is hypothesized that the NAO modulates the position and strength of the equatorial westerlies, impacting the Tropical Easterly Jet and therefore Sahel rainfall.This hy-775 pothesis is currently under investigation by climate domain scientists (i.e.(Tetteh, 2012)), based on the results of this study, and serves as an example of a relationship which is not fully understood being highlighted by the framework presented here.Furthermore, the presented CHARM approach studies rules heterogeneously, which is for the benefit of understanding rules with significant support within particular coupling inciters.However, future work would include addressing rules across multiple inciters, given the geophysical nature of the underlying data and understanding that climate relationships are in reality affected by multiple inciters.
Additionally, given our findings regarding Lasso's sensitivity to finding relationships at varying beta magnitudes, future work will be directed towards limiting such sensitivity and/or its influence in the fused network.
Finally, some of the modulatory relationships identified by these methods may represent underlying causal pathways in the climate system.Future work will also focus on inferring these causal pathways by leveraging causal modeling frameworks, such as Causal Bayesian Networks.Under this framework, inferring causal relationships becomes a problem of network structure learning.Several score-based and constraint-based algorithms have been proposed to this end (Spirtes, 2010).However, due to the inherent complexity of the climate system, learning this causal network structure is not a simple task.Future work should include identifying an appropriate causal inference algorithm for the problem at hand, by determining which underlying assumptions must hold to infer causal models for the climate domain.This causal inference algorithm should not assume causal sufficiency nor acyclicity of the causal structure (Hyttinen et al., 2013), since latent variables (i.e., confounders) and feedback loops are ubiquitous in the climate system.This algorithm should also be able to handle the high-dimensionality and small sample size of climate data (Bühlmann, 2013).Furthermore, an algorithm that allows to incorporate prior knowledge (i.e., known causal relationships from domain knowledge) would also be desirable (Borboudakis et al., 2011).

Fig. 1 .
Fig. 1.Complex relationships between climate indices and Sahelian rainfall, with some direct and indirect relationships well defined in literature (light arrows) and others not fully understood (dark arrows)

Fig. 2 .
Fig. 2. Simplistic traditional representation of market basket data in the form of transactions.The rule Bread → Milk has a support of 0.6, meaning it appears in 60% of transactions.

Fig. 4 .
Fig. 4. Spatially-defined variables' anomaly presence or absence is affixed to climate index anomalies expanded to represent a global effect.

Fig. 6 .
Fig. 6.Heterogeneous rule sets are generated for each coupling inciter individually, preserving the independence of their anomalies.

Figure 8
Figure 8 captures the resulting network after the different 725

Fig. 9 .
Fig. 9. Relationships directly associated with rainfall for λ = EAT L8 770 WorkIn this work, we used the same set of climate indices for all the individual methodologies employed, to facilitate the : 11 comparison and fusion of results.It is possible that using different datasets or different time/spatial series for each methodology can improve their individual results, and in turn the overall outcome of the models.However, additional considerations would be needed for interpreting which models provide evidence of particular relationships between individual climate indices.

Table 4 .
Network vertex statistics for coupling inciter EAT L8