Ancestral causal learning in high dimensions with a human genome-wide application

TL;DR

This paper introduces a supervised learning method for inferring ancestral causal relationships in high-dimensional data, demonstrated on human genome-scale data with over 19,000 variables, showing robustness and scalability.

Contribution

It presents a scalable, supervised approach to learn ancestral causal structures in high dimensions, leveraging known relationships and applying it to human genome data.

Findings

Method is effective on human genome data with 19,000 variables.

Results are robust to perturbations in prior knowledge.

Approach scales to large, high-dimensional biological data.

Abstract

We consider learning ancestral causal relationships in high dimensions. Our approach is driven by a supervised learning perspective, with discrete indicators of causal relationships treated as labels to be learned from available data. We focus on the setting in which some causal (ancestral) relationships are known (via background knowledge or experimental data) and put forward a general approach that scales to large problems. This is motivated by problems in human biology which are characterized by high dimensionality and potentially many latent variables. We present a case study involving interventional data from human cells with total dimension $p \! \sim \! 19{,}000$. Performance is assessed empirically by testing model output against previously unseen interventional data. The proposed approach is highly effective and demonstrably scalable to the human genome-wide setting. We consider sensitivity to background knowledge and find that results are robust to nontrivial perturbations of the input information. We consider also the case, relevant to some applications, where the only prior information available concerns a small number of known ancestral relationships.

Equations

$$y_k=\{\begin{array}{ll}1& \text{if }\text{i(k) is a causal ancestor of j(k) according to Π}\\ 0& \text{otherwise.}\end{array}$$

$$\hat{\theta }(X,\Pi )=\underset{\theta \in \Theta }{\mathrm{argmin}}\sum _{k\in \mathcal{T}(\Pi )}\mathrm{RL}(f({\varphi }_k(X);\theta ),y_k(\Pi )),$$

$${\hat{y}}_{k^{\prime }}=f({\varphi }_{k^{\prime }};\hat{\theta }(X,\Pi ))$$

$${\hat{G}}_{ij}=\{\begin{array}{ll}f({\varphi }_{k(i,j)};\hat{\theta }(X,\Pi ))& \text{if }k(i,j)\notin \mathcal{T}(\Pi )\\ y_{k(i,j)}(\Pi )& \text{otherwise}\end{array}$$

$${\hat{G}}_{ij}(X,\Pi )=f({\varphi }_{k(i,j)};\hat{\theta }(X,\Pi )),\forall k\in \mathcal{K}.$$

$${\tilde{y}}_k=\{\begin{array}{ll}0& \text{if }k\in {\tilde{\mathcal{T}}}^{(1)}\\ 1& \text{if }k\in {\tilde{\mathcal{T}}}^{(0)}\\ y_k& \text{otherwise}\end{array}$$

Full text

Keywords:   causal learning, ancestral causality, human gene expression data, interventional data, supervised learning

1 Introduction

We consider the learning of ancestral causal relationships in high dimensions. Suppose $V=\{1\ldots p\}$ indexes a set of variables of interest, with $p$ potentially large. Our goal is to learn a directed graph $\hat{G}$ comprising vertex set $V$ and edge set $E$, with edges $(i,j)\!\in\!E$ having the interpretation that node $i$ has an ancestral causal influence on node $j$. Here, an ancestral causal relationship refers to the existence of a directed causal path from the ancestor to the descendant (possibly via latent variables). Such relationships are of particular relevance to experimental studies due to their amenability to empirical verification against experimental data (Zhang, 2008). In particular, a change in the variable with index $j$ induced by an intervention on $i$ implies the existence of such an ancestral path from $i$ to $j$.

Many causal learning methods are based on graphical models, with models based on directed acyclic graphs (DAGs) playing a key role (Spirtes et al., 2001; Pearl, 2009). The PC algorithm is a prominent example of such a method (Spirtes et al., 2001). It estimates an equivalence class of DAGs – encoded as a completed partially directed acyclic graph or CPDAG – via a series of conditional independence tests. The PC output can in turn be used to estimate bounds on quantitative total causal effects between nodes via an algorithm known as IDA (Intervention calculus when the DAG is Absent; Maathuis et al., 2009). Greedy Interventional Equivalence Search (GIES; Hauser and Bühlmann, 2011) is a score-based approach that allows for the inclusion of interventional data. FCI (Fast Causal Inference; Spirtes et al., 2001) and RFCI (Really Fast Causal Inference; Colombo et al., 2012) allow for latent variables and learn equivalence classes of ancestral graphs (encoded as partial ancestral graphs or PAGs, see Richardson and Spirtes, 2002). Recently, Malinsky and Spirtes (2017) showed how quantitative total causal effects can be obtained from PAGs using the LV-IDA algorithm.

In contrast to these approaches, which are rooted in data-generating models of the causal system, there has been recent work with an emphasis on learning to tell apart causal and non-causal relationships in a data-driven fashion. Work in this “discriminative” direction has included Lopez-Paz et al. (2015), Mooij et al. (2016) and Hill et al. (2019). Our work follows in this line. Our goals are to put forward a scalable approach for causal learning and investigate whether it can be effective in the challenging context of high-dimensional human genomic data.

In a nutshell, our approach works as follows. Available information on some ancestral causal relationships (via background knowledge or experimental data) are treated as “labels” that are combined with a featurization of the data to train a classifier. The fitted classifier is then used to obtain labels (or probabilistic scores) across the entire problem which are intended to encode ancestral causal relationships. The output is a directed graph (or associated probabilistic scores) where the presence of an edge $(i,j)$ means that $i$ is inferred to be a causal ancestor of $j$. Our work builds on aspects of the framework presented in Hill et al. (2019). The key methodological difference is that the present paper focuses on high-dimensional settings to which Hill et al. (2019) does not scale and on a supervised framework. Furthermore, we present a biological case study using new, genome-wide data. This represents a relevant use-case that goes beyond the scale to which their methods can be applied. We include a wide range of experiments investigating behavior in high dimensions and in terms of practically relevant factors, including sensitivity to incorrect prior information.

Our work is motivated by, and applied to, data from molecular biological experiments on human cells. Such experiments involve measuring molecular quantities across many thousands of variables (in our example the variables represent gene expression levels and $p\!\sim\!19{,}000$), including under intervention. The question of interest is to understand causal relationships between these variables. Ancestral relationships are of particular relevance in the context of data of this type (see e.g. Zhang, 2008). Human gene regulation is complex and many of the underlying processes have been investigated in considerable mechanistic detail, including via detailed biochemical and biophysical studies. The known set of processes and factors is enormously rich, with many types of variables beyond those included in the available data that are known to influence gene expression levels. These cover, among others, diverse non-coding RNAs (RNAs that are not translated into proteins) as well as proteins and their post-translational modifications (Mercer et al., 2009; Vaquerizas et al., 2009; Rinn and Chang, 2012). From a causal perspective, these mechanisms represent a high dimensional set of latent variables that may influence the observables in complex ways. High dimensionality and the richness of regulatory influences are characteristic of human biology and present key challenges for large-scale causal learning in biomedical research.

Our approach differs from graphical models-based methods in several key aspects. First, our method copes relatively easily with large numbers of variables, as demonstrated (in terms of both performance and computational efficiency) in our case study. Graphical models-based approaches are theoretically principled and benefit from well-developed causal semantics but the high-dimensional case poses well-known challenges for estimation. That said, it should be noted that while graphical models-based approaches usually aim to model the underlying data-generating process, our approach learns only ancestral edges, which is a more limited form of output. Second, our approach is not hurt by increasing dimensionality, since the number of objects being classified grows with the number of variables and the computations turn out to be very tractable, as we discuss below. We show empirically using the high-dimensional human data that the performance of our approach is stable as dimension increases up to $p\!\sim\!19{,}000$. Third, our approach allows for cycles in the graph. This stands in contrast to PAG-based ancestral approaches like FCI and RFCI. In biological settings, cycles are common, for example due to feedback regulation in the underlying systems (Alon, 2006). Cyclic causation has been discussed in, for example, Spirtes (1995); Richardson (1996); Hyttinen et al. (2001).

The remainder of the paper is organized as follows. We first introduce notation and describe the methods. We then present empirical results and conclude with a discussion.

2 Methods

2.1 Notation

Variables and graphs. The set of observed variables of interest is indexed by $V\!=\!\{1\ldots p\}$ and the variables themselves are denoted $(X_{1},\ldots,X_{p})$. The variables will often be identified with vertices in a directed graph. Where useful, the vertex and edge sets of a directed graph $G$ are denoted $V(G)$ and $E(G)$, respectively. Occasionally we overload $G$ to refer also to the corresponding binary adjacency matrix, as will be clear from context. Then, for a graph $G$, $G_{ij}$ refers to the entry $(i,j)$ of the corresponding adjacency matrix.

Linear indexing. To facilitate mapping to the machine learning problem, we use linear indexing of variable pairs. Specifically, an ordered pair $(i,j)\!\in\!V\!\times\!V$ has an associated linear index $k\in\mathcal{K}=\{1\ldots K\}$, where $K$ is the total number of variable pairs of interest. Where useful to make the mapping explicit, the linear index corresponding to a pair $(i,j)$ is denoted as $k(i,j)$ and the variable pair corresponding to a linear index $k$ as $(i(k),j(k))$. The linear indices of variable pairs whose ancestral causal relationships are unknown and of interest are denoted $\mathcal{Q}\subset\mathcal{K}$. That is, for all pairs $\{(i,j):k(i,j)\in\mathcal{Q}\}$ we want to learn whether or not $X_{i}$ has an ancestral causal influence on $X_{j}$.

Prior or background information. Some ancestral causal relationships $\mathcal{T}(\Pi)\subset\mathcal{K}$ are known in advance via background knowledge $\Pi$. In all empirical experiments the sets $\mathcal{T}(\Pi)$ and $\mathcal{Q}$ are disjoint, i.e. no prior causal information is available on the pairs $\mathcal{Q}$ of interest.

2.2 Ancestral causal learning via a supervised formulation

Our approach is based on framing the task of learning ancestral causal relationships as a supervised learning problem. We first formulate the problem in a generic manner without specifying implementation details. This makes it clear how specific combinations of featurization and supervised learning could be used to perform ancestral causal learning. We then describe the particular formulations used in experiments below. Our goal is to learn ancestral causal relationships between the variables $V$. These relationships are amenable to experimental verification (Zhang, 2008; Kocaoglu et al., 2017) and are of particular interest for follow-up experiments.

2.2.1 Problem setting

The inputs for the given problem setting are: (i) an $n\times p$ data matrix $X$ containing data on $p$ variables indexed $V\!=\!\{1\ldots p\}$; (ii) knowledge $\Pi$ concerning ancestral causal relationships on a subset $\mathcal{T}(\Pi)\!\subset\!V{\times}V$ of variable pairs; and (iii) a set $\mathcal{Q}\!\subset\!V{\times}V$ of variable pairs whose ancestral causal relationships are unknown and of interest. It is assumed that $\mathcal{T}$ and $\mathcal{Q}$ are disjoint. This corresponds to the assumption that no prior information is available on the causal relationships of interest, i.e. that the causal relationships in $\mathcal{Q}$ are entirely unseen. Similarly, while no specific assumption is made on the data $X$ (in particular the data need not be i.i.d.), it is assumed that it does not contain interventional data corresponding to the pairs in $\mathcal{Q}$.

2.2.2 Featurization

In order to formulate the problem as a machine learning task, we require a featurization that can be applied to both the pairs $\mathcal{T}$ whose causal status is known from background information $\Pi$ and the unknown pairs $\mathcal{Q}$. Let $\phi_{k}=\phi_{k}(X)\in\mathbb{R}^{d}$ generically denote such a featurization of the $p$-dimensional data (for generality this is written as a function of the entire data matrix). The subscript $k$ indicates that the features are specific to a variable pair $k\in\mathcal{K}$; this is important to allow linking to labels encoding causal status (see below). We consider a specific, very simple featurization below via vectorized bivariate histograms following Hill et al. (2019). Collecting the features over all variable pairs $\mathcal{K}$ gives rise to a $K\times d$ feature matrix $\Phi=[\phi_{1}(X)\ldots\phi_{K}(X)]^{\mathrm{T}}$.

2.2.3 Learning

Labels. The background information $\Pi$ provides information on ancestral causal relationships for the variable pairs $\mathcal{T}(\Pi)$. In the specific case of available interventional data, this means that prior interventional experiments reveal, for each pair $k\in\mathcal{T}$, whether or not $i(k)$ is a causal ancestor of $j(k)$ (recall that $i(k),j(k)$ is the ordered pair corresponding to the linear index $k$). This information is used to form (training) labels $(y_{k})_{k\in\mathcal{T}}$. Specifically,

[TABLE]

Collecting together these labels gives a label vector $Y(\Pi)=(y_{k})_{k\in\mathcal{T}}$, where the notation emphasizes that the labels are derived from $\Pi$.

The prediction function. Thus, labels are available for all vertex pairs $k\!\in\!\mathcal{T}$. For one such pair, the corresponding features are $\phi_{k}\!\in\!\mathbb{R}^{d}$ (the $k^{\mathrm{th}}$ row of the feature matrix $\Phi$). That is, $(\phi_{k},y_{k})_{k\in\mathcal{T}}$ form (feature, label) couples in the usual supervised learning sense. Let $f:\mathbb{R}^{d}\rightarrow\{0,1\}$ denote a prediction function that maps features to a binary label (or probabilistic score, in which case the range of $f$ is $[0,1]$). Assume $f$ can be fully specified by an unknown parameter $\theta$, giving the classification function as $f(\,\cdot\,;\theta)$ (here we use $\theta$ in a general sense to encode all information needed to specify $f$, potentially including parameters, hyper-parameters, model architecture, etc.). Estimation of $\theta$ – and thereby of the prediction function – is done in a supervised manner via

[TABLE]

where $\mathrm{RL}$ denotes generically a loss function with regularization (specific choices correspond to various supervised learning approaches) and $\Theta$ denotes the parameter space. Note that the learning step (2) requires only the variable pairs $\mathcal{T}$ for which the labels are available.

Now, using $f$, the estimated label for a new pair $k^{\prime}$ is

[TABLE]

with the notation emphasizing that estimation (of the parameters $\theta$ that in turn specify the function $f$) is based on the data and background information.

The estimate $\hat{G}$. Finally, the fitted classifier $f(\,\cdot\,;\hat{\theta}(X,\Pi))$ is used to obtain labels for all unknown pairs $\mathcal{Q}$ via (3) using the features $\{\phi_{k}\}_{k\in\mathcal{Q}}$. This gives an estimate $\hat{G}(X,\Pi)$ with edges specified by

[TABLE]

where $(i,j)$ are ordered variable pairs. Note that the overall estimate depends only on the data $X$ and background information $\Pi$, via the parameters of $f$, as estimated via (2). No change is made for pairs $\mathcal{T}$ whose status was known at the outset.

2.2.4 Correcting prior information

In some settings, available information $\Pi$ may not be entirely trustworthy. In this case, it may be desirable to allow the learner to attempt to correct the input labels. We consider doing so by simply applying the fitted classifier to all points $\mathcal{K}$ (i.e. not only those in $\mathcal{Q}$). This gives an alternative to (4) with edges estimated as

[TABLE]

Note that in contrast to (4) this means that for a pair $k\in\mathcal{T}$ whose causal status is (thought to be) known from background knowledge it may be that the output $\hat{G}_{ij}\neq y_{k(i,j)}$ (“error correction”).

2.3 Computational scaling

Several computational aspects make the proposed approach highly efficient for large-scale applications. In particular:

• •

The learning step (2) is restricted to the set $\mathcal{T}$ which in most applications will be much smaller than the complete set of pairs $\mathcal{K}$. If very large amounts of interventional data were available for the training step, (2) could be handled using standard techniques for supervised learning for large data, such as batch-wise, stochastic gradient type methods (these are already used in the neural network examples below).

• •

For problems with large $p$, the estimation of the graph structure (4), which spans to the complete set of pairs $\mathcal{K}$, involves applying a fitted classifier to a potentially large number $|\mathcal{Q}|$ of points. This can typically be done very efficiently computationally.

• •

For featurizations that involve dealing with variable pairs, construction of the feature matrix $\Phi$ is $O(p^{2})$, however this can usually be parallelized (this is the case for the featurization considered below).

Taken together, the foregoing points mean that the proposed approach can be effectively scaled to very high dimensional problems. To give an idea of practical compute times, the largest experiments in the case study below, involving $p\sim 19{,}000$ variables, took roughly 10 minutes to solve on a standard workstation (24 CPU cores and 256GB of RAM). Figure 8 in the Appendix shows run-times for our method and various standard algorithms as a function of $p$: the favourable scaling due to the reasons mentioned above can be clearly seen.

2.4 Specific model formulations used in the experiments

The overall approach we propose is modular in the sense that, in principle, any featurization $\phi$ could be combined with any predictor $f$ to provide a learner $\hat{G}$. In this Section, we provide some specific formulations that we use in experiments below. For featurization, we use a simple bivariate histogram approach following Hill et al. (2019). In brief, this involves constructing a bivariate histogram for each pair $k$, with $\phi_{k}$ being obtained as a PCA reduction of vectorized histogram bin counts down to $d=100$ dimensions.

We consider the following specific classifiers, which we refer to as Scalable Causal Learning (SCL):

• •

SCL-L1. Featurization as above, with prediction done using a standard $\ell_{1}$-regularized generalized linear model, as implemented in glmnet (Friedman et al., 2010).

• •

SCL-NN. As in SCL-L1, but with prediction performed using a feedforward neural network (implemented via tensorflow (Abadi et al., 2016), containing a total of 5 dense layers and 132,865 parameters).

3 Results

3.1 Overall set-up

Data. We consider a case study involving human gene expression data. After standard pre-processing and trimming, the data consisted of sequencing-based measurements of gene expression levels of $p=24{,}052$ genes in human cells (the neuroblastoma cell line BE(2)-M17) under intervention by short hairpin RNAs. The main interventional data consisted of 70 interventions, under each of which all $p$ gene expression levels were measured, and an additional panel (comprising 35 interventions) that was used solely to define gene-specific thresholds (see below). Specific genes affected by too many interventions could make the learning problem too easy, we therefore excluded genes affected by half or more of the interventions, leading to a total of approximately $19{,}000$ genes.

We sampled the data to obtain specific problem instances in the following way. A subset of 35 interventions was used to train and test the models (as described above, in terms of the set notation this means $|\mathcal{K}|=35\!\times\!p$). The remaining 35 interventional data were used to populate the data matrix $X$ along with the observational data. The intention was to mimic the realistic setting in which the available data contains variation of unspecified source. Note that the interventions included in $X$ are neither used to provide background information $\Pi$ nor to test the models.

The general set-up is that we use some interventional data to train the learners and then test the output against unseen interventional data. A gene $i$ was said to have an ancestral causal effect on gene $j$ if $j$’s expression level under intervention on $i$ was larger or smaller than any measurement (for the same gene $j$) in the panel of additional interventional data. This emphasizes salient changes whilst ensuring an appropriate scale for each gene $j$. The data were strictly split in the sense that (i) no data used to define the true ancestral relationships against which the model output was tested appear in the data matrices $X$ and (ii) in all experiments the pairs $\mathcal{Q}$ on which the model output is tested are disjoint from those pairs $\mathcal{T}$ whose ancestral relationships are provided as auxiliary background information.

Comparisons. We compare the proposed method against a panel of existing causal learning methods, including PC, IDA, RFCI and GIES (implemented in the R package pcalg (Kalisch et al., 2012)). We obtained total causal effects from the output of RFCI using the LV-IDA method described in Malinsky and Spirtes (2017). These methods differ in the nature of their output and, as noted at the outset, should not be regarded as direct competitors to our approach. Rather, we include them as important examples of causal learning methods whose output can reasonably be tested with respect to interventional experiments of the kind we consider here. While IDA and RFCI/LV-IDA provide output that can reasonably be directly compared with the interventional data (since they estimate total causal effects), this is not the case for PC and GIES. In the case of PC and GIES we therefore considered an additional transitive closure step to extract ancestral statements from the CPDAG or essential graph output.

Additionally, we considered some non-causal statistical methods including Gaussian graphical models (GGMs; estimated using a shrinkage estimator as implemented in the R package corpcor (Schäfer and Strimmer, 2005)), and Pearson and Kendall correlation coefficients. Although non-causal, these approaches are included to provide a set of simple data-driven baselines. We note that non-causal methods of this kind remain widely used in many applied settings with scientific goals that are arguably causal (including in bioinformatics), hence it is interesting to consider their empirical behavior with respect to human interventional data.

In the experiments below, we consider in turn varying a number of key factors, including total problem dimension and the amount of available background information. We go on to assess sensitivity to incorrect background information, where we intentionally induce errors in the input information and examine whether such incorrect background information can be corrected using the “error correction” approach described above. Finally, we look at a set-up, relevant in several applications, where background information is very sparse in the sense of comprising information on only a small number of known ancestral causal relationships (rather than the entire set of changes under a given intervention).

3.2 Varying dimension $p$

Random sampling. Figure 1(a) shows performance as a function of problem dimension $p$, with the fraction of background information $\rho$ set to 0.5 throughout (in terms of the set notation above, $\rho$ equals $|\mathcal{T}|/|\mathcal{K}|$, i.e. the fraction of all pairs whose causal status is known to the learner as background information; we consider varying $\rho$ below). Figures 1(b-d) show the ROC curves for three values of $p$. Results from PC (GIES is only appropriate for an “intervention-wise” problem set-up, see below) are shown as points on the ROC plane. For PC and RFCI/LV-IDA we show results for two different significance levels $\alpha\in\{0.01,0.5\}$; additional results for different $\alpha$’s appear in the Appendix (Figure 9).

The proposed methods perform well in absolute terms and relative to other approaches. Note that the performance initially improves with increasing dimension $p$, in line with the arguments presented above with respect to the relationship between dimension and the learning task. Computational cost is moderate even for the largest $p$ (see Figure 8 in the Appendix for wall clock times).

Intervention-wise sampling. The results above are based on random sampling of pairs (in the sense that $\mathcal{T}$ is a random subset of $\mathcal{K}$). Figure 2(a) shows the corresponding result for intervention-wise sampling. The latter corresponds to the case where the learner is asked to generalize to an unseen intervention (and not just an unseen causal relationship or pair) and is in a sense a harder problem. Suppose $I$ is the initial set of available interventions. Intervention-wise sampling is done by sampling a random subset $I^{\prime}\subset I$ of the interventions and setting $\mathcal{T}=\{k:i(k)\in I^{\prime}\}$. This ensures that for any test pair $q\in\mathcal{Q}$, it always holds that $i(q)\notin I^{\prime}$ or equivalently $q\in\mathcal{Q}\implies\nexists j:k(i(q),j)\in\mathcal{T}$. This amounts to training on some interventions and testing on interventions that are entirely unseen in the sense that for any test pair $(i,j)$, the training pairs $\mathcal{T}$ do not include the effect of intervening on $i$ on any other variable.

We note that GIES requires the interventional data itself, whereas our approach requires only information on known ancestral relationships (encoded in $\Pi$). For these experiments, to enable a reasonable comparison, we provided GIES the interventional data upon which the information in $\Pi$ is based, i.e. precisely the interventions $I^{\prime}$. Hence, in this example, GIES receives additional data (the interventional data for interventions $I^{\prime}$) that the other methods do not. We note also that since GIES requires the interventional data on all variables, it cannot be obviously applied to the random sampling case above.

3.3 Varying amount of available background information

The previous results are based on a fixed fraction $\rho$ of pairs on which background information is available. Next, we considered varying the fraction $\rho$. This is practically relevant since in many settings information may only be available on a relatively small fraction of pairs. Figure 3 shows performance as a function of $\rho$ for the random sampling case. Figure 4(a)(a) shows the corresponding result for the intervention-wise case. In both cases, we find that performance improves rapidly with the fraction $\rho$. Since only GIES from among the other methods varies with $\rho$, we show only results for the proposed approaches and GIES (in the intervention-wise case). We were not able to run GIES for the very high dimensional examples due to computational demands. For the other methods, performance would be as shown in the previous figures, since for these methods $\rho$ plays no role.

3.4 Sensitivity to background information

The learner relies on knowledge $\Pi$ of some ancestral causal relationships. In this section we consider sensitivity to perturbation of such information. The intention is to investigate the effect of including in $\Pi$ information that is in fact incorrect. This was done by intentionally perturbing the input labels in such a way as to control the fraction of labels altered whilst keeping the overall sparsity (i.e. the fraction of ancestral causal edges in the full graph) fixed.

To clarify how this was done, we introduce some additional notation. In what follows let $\tilde{\Pi}$ denote a perturbed version of the background knowledge. For the experiments all learning was done as described above, but with the perturbed $\tilde{\Pi}$ in place of $\Pi$. Denote by $\mathcal{T}^{(1)}=\{k:k\in\mathcal{T}\wedge y_{k}(\Pi)=1\}$ the set of all causal pairs known from background knowledge $\Pi$ and by $\mathcal{T}^{(0)}=\mathcal{T}\setminus\mathcal{T}^{(1)}$ the complement (i.e. pairs known to be non-causal from background knowledge). Let $\tilde{\mathcal{T}}^{(1)}\subset\mathcal{T}^{(1)}$ be a random subset of the positive cases and similarly $\tilde{\mathcal{T}}^{(0)}\subset\mathcal{T}^{(0)}$ for the negatives (with $|\tilde{\mathcal{T}}^{(0)}|=|\tilde{\mathcal{T}}^{(1)}|$) and $\tilde{\mathcal{T}}=\tilde{\mathcal{T}}^{(0)}\cup\tilde{\mathcal{T}}^{(1)}$ be the complete set of perturbed pairs. Now, we perturb the labels as

[TABLE]

This means that exactly $|\tilde{\mathcal{T}}^{(1)}|$ of the known causal pairs are perturbed whilst leaving the total number of “known” casual pairs as input to the learner unchanged. The above procedure guarantees that (i) causal sparsity is unchanged in the sense that $\frac{1}{|\mathcal{T}|}\sum_{k\in\mathcal{T}}y_{k}=\frac{1}{|\mathcal{T}|}\sum_{k\in\mathcal{T}}\tilde{y}_{k}$, and (ii) input information on the perturbed pairs is indeed incorrect, i.e. $\forall k\in\tilde{\mathcal{T}}:\tilde{y}_{k}\neq y_{k}$.

Figure 5 shows performance as a function of the fraction $|\tilde{\mathcal{T}}^{(1)}|/|\mathcal{T}^{(1)}|$ of the input causal information that was corrupted as described above. Performance is highly robust for the $\ell_{1}$ approach but the neural network is dramatically more fragile. We note that this is probably not a fundamental difference between neural networks and sparse regression and could likely be rescued with a different architecture and/or regularization strategy. Rather, the point is that a highly flexible model can easily fit the input information, hence for corrupted input information the simplicity of the $\ell_{1}$ model is an advantage and amounts to stronger regularization.

3.5 Correcting background information

As described in Section 2.2.4, the learner could in principle be used to correct prior information. Here, we investigate the behavior of this approach in the perturbed input setting described in the previous Section, taking advantage of the fact that we have intentionally induced errors. Figure 6 shows performance with respect to error correction, as a function of the fraction $|\tilde{\mathcal{T}}^{(1)}|/|\mathcal{T}^{(1)}|$ of the input causal information that was corrupted. Note that here performance is assessed with respect to the set of perturbed pairs $\tilde{\mathcal{T}}$. Although these are a subset of the training data, the inputs are by design incorrect. Hence, in the present case an AUC of unity implies that the learner corrects the input labels to their true (i.e. pre-perturbation) values. We find that the $\ell_{1}$ learner is highly effective at error correction in the high dimensional settings: this is a blessing of dimensionality and is due to the fact that in larger problems there is more information and therefore more opportunity to detect regularities that can then be used to cope with incorrect inputs and ultimately, as shown here, rectify such inputs.

3.6 Knowledge of a small number of ancestral relationships

In fields such as economics, epidemiology or the social sciences, although large scale interventional experiments of the kind now performed in molecular biology are not possible, some ancestral relationships are typically known from the science. However, these tend not be the total effects of a certain intervention on a large number of variables, but rather isolated facts (corresponding in our notation to specific pairs $(i,j)$). To take one example: an expert epidemiologist might not be able to provide prior information on the effect of any intervention on all measured variables in a high-dimensional study, but might be able to specify some pairs (e.g. connecting specific exposures to specific phenotypes). Here, we take advantage of the molecular data in our case study to ask whether learning is possible if the background information $\Pi$ contains only a small number of such ancestral relationships. That is, we consider the case where $\Pi$ provides information on only a small number of pairs $(i,j)$ where it is known that $i$ is a causal ancestor of $j$ (i.e. only positive information).

To study this case, rather than giving the learners access to background information as in the examples above, we allowed access only to information on a specified fraction of (positive) ancestral relationships. Specifically, let $\mathcal{T}^{+}\!\subset\!\mathcal{T}^{(1)}$ be a (possibly small) random subset of the causal pairs known from background knowledge (notation as in Section 3.4). We applied the proposed methods using this information to construct the input labels $y_{k}$ by simply assuming all other pairs as non-causal, i.e. setting $\tilde{y}_{k}=1$ if $k\in\mathcal{T}^{+}$ and to zero otherwise. Figure 7 shows results for the random sampling and intervention-wise cases. Performance improves rapidly with the number $|\mathcal{T}^{+}|$ of positive examples available and the gain is more rapid for higher dimensional examples. The $\ell_{1}$ approach is more effective (likely due to model complexity as discussed above). As a point of comparison, we considered the case of locating the “known” ancestral examples entirely at random (that is, randomly setting $\tilde{y}_{k}=1$). Learning failed in all cases, with AUCs approximately 0.5 (see Figure 12 in the Appendix).

4 Discussion

In this paper, we presented a simple, discriminative approach to causal learning that we showed can be effectively applied at human genome-wide scale. The general approach we propose is modular in the sense that a variety of specific featurizations and learners could be used. The specific approach of SCL that we focused on in our empirical work scaled demonstrably well, both in terms of computational demand and statistical performance.

Although effective at the specific tasks that we considered, in contrast to graphical models-based methods, SCL is less general in the sense that it only learns ancestral relationships, but cannot provide a full range of probabilistic output (e.g. post-intervention distributions). Nevertheless, its good performance and favourable scaling suggest that it may be fruitful to investigate combining SCL with graphical models-based approaches in future.

In molecular biology, the question of how to learn causally meaningful, systems-level networks remains open. If this were possible at scale, the implications for study design would be considerable. Our results suggest that such scalable causal learning may be possible, but further empirical work will be needed in future to fully understand the scope and limitations of our approach in diverse biological contexts.

Appendix A Additional results