1 Introduction
1.1 Our contribution and related work
We propose a very simple yet effective method to infer the ancestor variables in a linear structural equation model from observational data.

Consider a response variable of interest
and covariates in a linear structural equation model. The procedure is as follows. For a nonlinear function , for example , run a least squares regression of versus and all covariates : the pvalue corresponding to the th covariate is measuring the significance that is an ancestor variable of, and it provides type I error control.
We refer to this methods as ancestor regression. Its power (i.e., type II error) depends on the nature of the underlying datagenerating probability distribution. Obviously, the proposed method is extremely simple and easy to be used; yet, it deals with a difficult problem of finding the causal order among random variables. In particular, the proposed method does not need any new software and it is computationally very efficient.
Structure search methods based on observational data for the graphical structure in linear structural equation models have been developed extensively for various settings: for the Markov equivalence class in linear Gaussian structural equation models spirtes2001causation, Chapter 5.4; chickering2002optimal or for the single identifiable directed acyclic graph in nonGaussian linear structural equation models (shimizu2006linear; gnecco2021causal)
or for models with equal error variances
(peters2014identifiability). All these methods do not come with pvalues and type I error control. In addition, for the identifiable cases, the corresponding algorithms assume certain assumptions such as nonGaussian errors, while our procedure does not rely on such conditions but automatically exploits whether the structure is identifiable or not.Our contribution can be summarized in short as follows. We propose ancestor regression which, as ordinary least squares regression, is very easy to be used and equipped with type I error control for inferring ancestor variables in linear structural equation models. The method does not rely on additional assumptions such as nonGaussianity: it automatically avoids to claim causal directions if the problem is nonidentifiable.
1.2 Notation
We present some notation that is used throughout this work. We use upper case letters to denote a random variable, e.g., or . We use lower case letters to denote i.i.d. copies of a random variable, e.g., . If , then . With a slight abuse of notation, can either denote the copies or realizations thereof. We write to denote the jth column of matrix and to denote the element in row and column . With , we emphasize that an equality between random variables is induced by a causal mechanism.
2 Ancestor regression
2.1 Model and method
Let be given by the following linear structural equation model (SEM)
(1) 
where the are independent and centered random variables. We assume that such that the covariance matrix of exists and has full rank. We use the notation , , and for ’s parents, children, ancestors, and descendants. Further, assume that there exists a directed acyclic graph (DAG) representing this structure.
Let with be a variable of interest; it has been denoted as response in Section 1. Consider a nonlinear function . The following result describes the population property of ancestor regression, with general function .
Theorem 1.
Assume that the data follows the model (1). Consider the ordinary least squares regression versus , denote the according OLS parameter by and assume that it exists. Then,
Importantly, itself must be included in the set of predictors as well. The beauty of Theorem 1 lies in the fact that no assumptions on the distribution of the or the size of the
apart from existence of the moments must be taken. This allows to control against false discovery of ancestor variables. The proof of this and the subsequent theoretical results is available in the supplementary material.
Typically, holds for ancestors since a nonlinear function of that ancestor cannot be completely regressed out by the other regressors using only linear terms. For ancestors that are much further upstream, this effect might become vanishingly small. However, this is not such an issue since when fitting a linear model using the detected ancestors, those indirect ancestors are assigned a direct causal effect of anyway.
Based on Theorem 1, we suggest testing for in order to detect some or even all ancestors of . Doing so for all , requires nothing more than fitting a multiple linear model and using its corresponding ztests for individual covariates.
Let be i.i.d. copies from the model 1. Define the following quantities
(2) 
where is meant to be applied elementwise in .
Theorem 2.
Due to this limiting distribution, we suggest testing the null hypothesis
with the pvalue(3) 
where
denotes the cumulative distribution function of the standard normal distribution.
For ancestors, for which , the absolute zstatistics increases as . In typical setups, one can thus detect all ancestors. Having found all ancestors, one can easily infer the parents with an ordinary least squares regression of versus , using the test for assigning the significance of being a parental variable.
2.2 Adversarial setups
There might be cases where does not hold true for some ancestors leading to reduced power of the method. We discuss such cases in the following.
Gaussian .
It is wellknown in causal discovery for linear SEMs that Gaussian error terms lead to nonidentifiability.
Define
i.e., the children of through which a directed path from to begins.
Proposition 1.
Thus, if every directed path from to starts with an edge for which the nodes on either end have Gaussian noise terms, we have no power to detect this ancestor relationship. However, we neither detect the opposite direction as guaranteed by Theorem 1 and thus, control against false positives is guaranteed.
Error terms with the same distribution.
A pathological case occurs, if a child’s error term has the same distribution as the inherited contribution from the parent’s error term. Then, all ancestors to which all directed paths from start with that given child are not detectable.
Proposition 2.
Assume that the data follows the model (1). Let and . Then, it holds
For the variables discussed here, the limiting Gaussian distribution as stated in Theorem
2 is not guaranteed even though ; see also the proof in the supplemental material.2.3 Simulation example
We study ancestor regression in a small simulation example. We generate data from a SEM that can be represented by the directed acyclic graph in Fig. 1 and aim to infer the ancestors of . We chose the precise setup such that neither of the adversarial cases discussed above is present. We let the sample size vary from to executing 200 simulation runs for each. As a nonlinear function, we use . By zstatistics, we mean as defined in Theorem 2. We calculate pvalues according to (3) and apply a BonferroniHolm correction (without cutting of at for the sake of visualization) to them.
In Fig. 1, we see the desired growth of the absolute zstatistics for the ancestors, while for and their sample averages are close to their theoretical mean. The indirect ancestor is slightly harder to detect than the two parents. Although the null distribution is only asymptotically achieved, the type I familywise error rate is controlled for every sample size, supporting our method’s main benefit, i.e., robustness against false causal discovery.
3 Ancestor detection in graphs: nodewise and recursive
3.1 Algorithm
In the previous section, we assumed that there is a (response) variable that is of special interest. However, this is not always the case, instead, one might be interested in inferring the full graph of causal connections between the variables. Naturally, our ancestor detection technique can be extended to graphs by applying it nodewise. We suggest the procedure described in Algorithm 1.
In BuildRecursive the graph is recursively constructed by adding the estimated ancestors of every estimated ancestor. This recursive construction facilitates the detection of all ancestors. This procedure cannot increase the type I familywise error rate compared to just using the significant pvalues because a false causal discovery can only be propagated if it existed in the first place.
Since there is no guarantee that the recursive construction does not create directed cycles, we need to address this. In FindStructure, the significance level is gradually reduced until no more directed cycles are outputted. This means, that the output becomes somewhat independent of the significance level, e.g., in a case with two variables and and , we would never claim no matter how large is chosen.
In practice, when there is no a priori knowledge on the structure of the data, one can also consider the largest such that no loops are created as a pvalue for the null hypothesis that the modeling assumption (1) holds true if the multiplicity correction is designed to control the familywise error rate, e.g., BonferroniHolm. We summarize the properties of our algorithm.
3.2 Simulation example
We extend the simulation from Section 2.3 to estimating every variable’s ancestors using Algorithm 1. We compare our method to LINGAM using the implementation provided in the Rpackage pcalg (kalisch2012causal). Based on the graph given in Fig. 1, we consider two slightly different linear SEMs. In the first, which uses the same data as the simulation in Section 2.3, only follows a Gaussian distribution while as in the second, both and do such that Proposition 1 applies. For our method, we apply a BonferroniHolm correction (without cutting off at for the sake of visualization). As LINGAM provides an estimated set of parents, we additionally apply BuildRecursive to the output to get an estimated set of ancestors to enable comparison with our method.
For the model with only one Gaussian error variable, we can reliably detect all ancestors without any false causal claims for . Not all curves reach a power of even when letting the significance level become arbitrarily large. This can be explained by the possible insensitivity to the significance level as sketched in Section 3.1.
We are able to control the familywise error rate even for low sample size, which is not the case for LINGAM. The power of LINGAM approaches faster than ancestor regression and one might consider it to be favorable for or in the model with one Gaussian noise term. The picture changes when looking at the model with multiple Gaussian noise terms. LINGAM is still more powerful for some sample sizes but does not control the error at all. Ancestor regression reaches a power of almost without any false causal discovery. The lacking are due to the fact that out of ancestor relationships depicted in Fig. 1 fulfills the conditions of Proposition 1 for the nonidentifiable Gaussian case.
4 Real data example
Causal effect  Suggested  Passing HOLS  Significant in linear model  minimum pvalue 

RAF MEK  2  1  1  0 
Akt Erk  7  2  2  5.1e207 
p38 PKC  4  1  1  1.1e143 
PKA Akt  3  1  1  5e119 
MEK RAF  5  2  2  7.1e115 
Akt PKA  4  2  2  2.1e44 
PIP3 PIP2  7  1  1  3.7e41 
PIP3 PLCg  4  2  1  9.6e35 
JNK p38  3  2  2  1.6e17 
PKC p38  4  1  1  3.6e16 
PKA Erk  3  1  1  6.6e12 
PKC JNK  4  1  1  3.1e09 
PIP2 PLCg  1  1  0  1 
PLCg PIP2  5  1  0  1 
level in the respective linear model fit, i.e., a direct (parental) causal edge, correcting for the total number of suggested connections in all environments. The last column reports the minimum of the pvalues from linear regression among environments, where the edge connection is suggested and passes the HOLS check.
We analyze the flow cytometry dataset presented by sachs2005causal. It contains cytometry measurements of 11 phosphorylated proteins and phospholipids. Data is available from various experimental conditions, some of which are interventional environments. The dataset has been further analyzed in various projects, see, e.g., mooij2013cyclic, meinshausen2016methods and taeb2021perturbations. Following these works, we consider data from 8 different environments, 7 of which are interventional. The sample size per environment ranges from to .
For each environment, we fit a graph using Algorithm 1 with nonlinear function and BoferroniHolm correction. For the claimed set of ancestors, we subsequently fit an ordinary linear regression to see which ancestors might be parents.
The outputted ranges from to . This happens due to the occurrence of estimated directed cycles and thus, there is strong indication that the data does not follow the model (1). The deviation can be in terms of hidden variables, nonlinear effects, or noise that is not additive. Therefore, we want to additionally incorporate the HOLS check introduced in schultheiss2022assessing that aims to find locally unconfounded linear effects in regression models (here, when regressing on the estimated ancestor variables). We find some edge connections that pass this subsequent HOLS check at the level in certain environments. We omit the multiplicity correction for this check to lower the tendency to falsely claim causal detection. Thus, we find that the assumption of locally unconfounded linear effects is not unrealistic. We summarize our findings in Table 1 which reports suggested edge connections that pass the HOLS check in at least one environment.
Acknowledgment
C. Schultheiss and P. Bühlmann have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme.
References
Appendix A Proofs
a.1 Additional notation
We introduce additional notation that is used for these proofs.
Subindex , e.g., denotes a matrix with all columns but the th. is the
dimensional identity matrix.
denotes the orthogonal projection onto and denotes the orthogonal projection onto its complement. is the orthogonal projection onto all .For some random vector
, we have the moment matrix . This equals the covariance matrix for centered . We assume this matrix to be invertible. Then, the principal submatrix is also invertible. We denote statistical independence by .a.2 Previous work
We adapt some definitions from and results proven in schultheiss2022assessing.
where  (4)  
where  
Using these definitions, we have from partial regression. We cite a Lemma fundamental to our results.
a.3 Proof of Theorem 1
Let . Then, we can write for a suitable with if . We can now find using this representation.
The third equality follows from the independence of the . Naturally, for all we have such that . Further, if . To see this, note that would be lower triangular, if denoted a causal order. Then, its inverse would be lower triangular as well. Naturally, the same principle applies for every other permutation. Thus,
such that if , i.e., if is not an ancestor of .
Alternatively, we could invoke Lemma 1 to see that for a nonancestor . Then,
a.4 Proof of Theorem 2
Define
Since we assume the covariance matrix to be bounded, we find
where we use invertibility and the continuous mapping theorem. This then implies
and analogously
We get a better rate for than for since we assume existence of the fourth moments. Then,
(5) 
The second line is restricted to covariates for which this variance exists, which includes all nonancestors as . Using Slutsky’s theorem, we have
(6) 
which proves the first part of the theorem.
It remains to consider the variance estimate. Similar to above, we have
Further,
Similar to before
Combined, we find
proving the second part of the theorem.
a.5 Proof of Proposition 1
Using the representation from above, we have
We find
The last remaining term corresponds to the only set such that , i.e., the only set that could potentially induce dependence. However, in the setup of Proposition 1, both sums are Gaussian random variables such that uncorrelatedness implies independence. Therefore, such that and .
By the same argument for all mediators between and . Generally, . As is independent of all mediators but also of all nondescendants of itself, it is independent of all ancestors of except itself. Therefore, it is also independent from
such that the variance as in (5) is consistently estimated.
a.6 Proof of Proposition 2
Following the derivation from the previous section, we see that dependence between and can only be induced through with . Since we exclude mediated paths in the assumptions of Proposition 2, dependence between and can only be induced by and . It follows
For the last equivalence, we use the assumption on the distribution.
Furthermore,
since contributions of in must be mediated through with the according weights. We get
Since and the two terms appear with the same weight in , the two expectations cancel each other such that . However, in general such that is not generally true. Then, the limiting distribution of the estimator is not the same as for nonancestors, see also (5).
a.7 Proof of Corollary 1
Let be such that . This means that there is at least one set
such that , where . If at least one of these must correspond to a false causal discovery, i.e., there is an such that but . We conclude
Let denote the number of true null hypotheses. By the construction of BonferroniHolm
Let as used in Theorem 2. We find
which proves the first part of the corollary. The second to last equality uses Theorem 2 and the continuous mapping theorem.
As the model (1) excludes the possibility of directed cycles, any output of BuildRecursive that contains cycles must include at least one false causal detection. If , it corresponds to the maximal pvalue such that including the corresponding ancestor relationship creates cycles. Therefore, it must hold such that
using similar arguments as above.