Conditional Independence Testing with Heteroskedastic Data and Applications to Causal Discovery
Authors: Wiebke Günther
Urmi Ninad
Jonas Wahl
Jakob Runge
Venue: NeurIPS 2022 Accept (2022)
Keywords: conditional independence test
heteroskedasticity
causal discovery
causal model

Abstract: Conditional independence (CI) testing is frequently used in data analysis and machine learning for various scientific fields and it forms the basis of constraint-based causal discovery. Oftentimes, CI testing relies on strong, rather unrealistic assumptions. One of these assumptions is homoskedasticity, in other words, a constant conditional variance is assumed. We frame heteroskedasticity in a structural causal model framework and present an adaptation of the partial correlation CI test that works well in the presence of heteroskedastic noise, given that expert knowledge about the heteroskedastic relationships is available. Further, we provide theoretical consistency results for the proposed CI test which carry over to causal discovery under certain assumptions. Numerical causal discovery experiments demonstrate that the adapted partial correlation CI test outperforms the standard test in the presence of  heteroskedasticity and is on par for the homoskedastic case. Finally, we discuss the general challenges and limits as to how expert knowledge about heteroskedasticity can be accounted for in causal discovery.

Conditional Independence Testing with Heteroskedastic Data and Applications to Causal Discovery
Wiebke Günther German Aerospace Center Institute of Data Science 07745 Jena, Germany wiebke. guenther@dlr.c
Urmi Ninad Technische Universität Berlin 10623 Berlin, Germany urmi rimii.inad@tu-brr//
Jonas Wahl Technische Universität Berlin 10623 Berlin, Germany
wahl@tu-berlin.e
Jakob Runge
German Aerospace Center Institute of Data Science 07745 Jena, Germany and Technische Universität Berlin 10623 Berlin, Germany j jakob. reaebbrungg@dlr.com
Abstract
Conditional independence (CI) testing frequently used data analysis and ma- chine learning for various scientific fields and it forms the basis of constraint- based causal discovery. Oftentimes, CI testing relies on strong, rather unrealistic assump- tions. One of these assumptions is homoskedasticity. in other words, constant conditional variance assumed. We frame heteroskedasticity in structural causal model framework and present adaptation of the partial correlation CI test that works well in the presence of heteroskedastic noise, given that expert knowledge about the heteroskedastic relationships is available. Further, provide theoretical consistency results for the proposed CI test which carry over to causal discovery under certain assumptions. Numerical causal discovery experiments demonstrate that the adapted partial correlation CI test outperforms the standard test the presence of heteroskedasticity and is par for the homoskedastic case. Finally, we discuss the general challenges and limits as to how expert knowledge about heteroskedasticity can be accounted for in causal discovery.
1
Introduction
Conditional independence (CI) testing is frequently used step across wide range of machine learning tasks for various scientific fields. is also very challenging. Discovering causal relationships from purely observational data is even more challenging task and an important topic in sciences where real experiments are infeasible, e.g. in climate research [Ebert-Uphoff' and Deng, 2012, Runge al., 2019a]. One can distinguish several frameworks that address this problem: score-based approaches [Chickering, 2002], restricted structural causal models [Peters et al., 2017], and constraint- based methods [Spirtes et al., 2000] that rely on CI testing. A typical representative is the PC algorithm [Spirtes and Glymour, 1991] which can be combined with any conditional independence
"qqull contribution
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

test and is thus adaptable to a wide range of data distributions. It utilizes the Faithfulness assumption to conclude that no causal link can exist between two variables X and if CI test suggests that they are independent given set Z.
When linear additive noise model may be assumed, popular CI test is the partial correlation test [Lawrance, 1976] which has the advantages of fast computation time and that the null distribution known analytically, As a disadvantage rather strict and often unrealistic assumptions have to be satisfied. One of these assumptions is homoskedasticily, meaning that the variance of the error term constant. When this assumption is violated, i.e., in the case of heteroskedasticity, the variance can, for instance, depend on the sampling index or the value of one or multiple influencing variables.
In regression analysis one distinguishes between impure and pure heteroskedasticity. Impure het- croskedasticity stems from an insufficient model specification, c.g. if unobserved variables or con- founders are present. Another reason might be that the model fails capture the full relationships of the variables, e.g. the underlying model might be linear with multiplicative noise rather than additive. On the other hand, pure heteroskedasticity is non-constant noise variance that is present despite correct model.
The sources for heteroskedasticity in real data are manifold. For example, in environmental sciences precipitation in different areas might exhibit different variances that are unaccounted for by other variables in the system, i.e. location-scalle
noise. Such problem could be introduced by aggregating data of different catchments and not adding a variable that is well enough correlated with catchment location [Merz et al., 2021]. An example for sampling index-dependent heteroskedasticity in the time series case are seasonal effects that are present in many climate variables [Proietti, 2004]. Finally, an example where the noise variance of a variable is dependent on an observed cause is the distance to sea influencing the variability of temperature. This is a special case of state-dependent noise.
One assumption of the PC algorithm is causal sufficiency. meaning that there are no unobserved confounders, formally excluding impure heteroskedasticity. However, in practice this assumption is often violated which can lead to heteroskedastic noise.
If unaccounted heteroskedasticity is present in the data, the estimator of the ordinary least squares (OLS) regression slope parameter is still unbiased. However, the estimator of the covariance matrix of the parameter estimates can be biased and inconsistent under heteroskedasticity [Long and Ervin, 2000]. This can skew subsequent partial correlation significance tests and affect the link detection rate of causal discovery method. Furthermore, heteroskedasticity might even lead to the detection of wrong links (false positives). Moreover, the Gauss-Markon theorem assumes homoskedasticity and, hence, with heteroskedastic data the OLS slope estimator is no longer guaranteed to be the most efficient linear unbiased estimator, which can further harm power of the CI test.
There are multiple ways to treat heteroskedasticity, either already during the modeling step by addressing possible model misspecification, by pre-processing the data, e.g. by applying a log- transform, or post-hoc using robust statistics for CI testing. Adapting the model, for example, by using CI tests allowing for multiplicative dependencies [Runge, 2018] might not always be feasible for limited sample sizes. Pre-processing in real-world problems also comes with drawbacks: The transformed variables are difficult to interpret and can introduce dependencies or spurious links.
In this work we propose an adapted weighted least-squares (WLS) partial correlation variant as CI test for the PC algorithm that is able to deal with particular forms of heteroskedasticity. Our contributions are theoretical consistency results as well numerical experiments demonstrating that this approach yields well-calibrated CI tests leading to controlled false positive rates and also improves upon detection power as compared to the standard partial correlation CI test. Our approach requires expert knowledge in that it needs to be known which of the variables the heteroskedasticity depends on, or if it depends on the sampling index.
2 Related Work
The effect of heteroskedasticity on the standard Pearson correlation test has been investigated, for example, by Wilcox and Muska [2001]. Remedies for heteroskedasticity also have been extensively studied. Hayes [2007] discuss and evaluate heteroskedasticity consssttttt standard errors. Practical approaches to choose weights for WLS, consistency and asymptotic results have been obtained, e.g. by Neumann [1994], Fan and Yao [1998], Carroll [1982], Robinson [1987], Brown and Levine [2007].

---

Romano and Wolf [2017] propose a method that combines WLS with beieroskedastiittycooiisseer standard errors.
Recently, within the causal discovery framework of restricted structural causal models [Peters et al., 2017] several authors [Xu al., 2022, Tagasovska et al., 2020] have relaxed the assumption of homoskedasticity. In these works the authors focus on identifying cause from effect only in the bivariate setting. Specifically, Xu al. [2022] base their inference score on the log- ligliiklioood of regression residuals and apply a binning-scheme on regions where variance is approximately constant.
More closely related to our work is the robust-PC method of Kalisch and Bühlmann [2008] who use an estimator in a recursive partial correlation formula to robustify the PC algorithm against outliers and otherwise contaminated data.
Non-constant conditional variance can also be regarded as a specific kind of distributional shift or context change. The effects of distributional shifts on causal discovery have been investigated, for instance, in Huang et al. [2020], Mooij et al. [2020] where the authors propose framework that includes the environment into the structural causal model formulation using context variables.
3 Problem setting
3.1 Heteroskedasticity in causal models
In this work, we consider discovering causal relationships in linear models in the presence of non- constant error variance which is potentially dependent on the parents or on the sampling index. To translate this into structural causal model (SCM), we represent heteroskedasticity function of the noise variable. In this way, it can also be viewed as state-dependent or multiplicative noise.
scaling
Consider finitely many random variables (X)... ......... x²) with joint distribution Px over a domain X' = X1 .... x Na. Then we are interested in samples from the following SCM with assignments
x = f((Pa(X))) + (+(())))) N.
i d
(1)
where fi are linear functions, E T stands for the sample index, and have the heteroskedasticity functions h; X' x T Rgo- The noise variables N; are assumed independent standard Gaussian distributions. The parent set of the variable Xi is denoted by Pu(X;), and H (X;) C Pa(X;) U {1} which can also be the empty set. Furthermore, are stable over time, i.c. the parent sets (x();) as well as the functions f. are not time-dependent.
make the restriction that the causal relationships
For h3 = const. the homoskedastic version of this SCM is obtained, which simply is linear SCM with Gaussian noise. To learn such a causal model from data, a well-known and easy to implement method is to use constraint- based causal discovery method with CI tests based on partial correlation.
To recap the partial correlation test and illustrate our adaptations, consider the following simple example model for a time-series, or otherwise indexed, SCM (1). The general multivariate case is discussed in section 4.2.
X. -0= cE, + hx (Z,, Nx al,
Z. Nz E, = Ne
Y, bZ, cE, + by Aby(Zz, (Ze, t) Ny
for some constants and and standard normal independent noise terms Nx, .N, N2, and Ng. The constant e set to zero for (conditionally) independent X and Y. Note that here we are only interested in discovering the relationship between X and Y and not in learning the whole causal graph.
An intuitive way to define and understand partial correlation is in terms of the correlation between residuals. The partial correlation between X and Y given a controlling variable Z. or a set thereof, is the correlation between the residuals rx and ry resulting from the linear regression of X on Z and of on Z, respectively. The linear regression that used in the standard variant of the partial correlation test is ordinary least squares (OLS) regression, thus we refer to this CI test as Corr-OLS. The Cov(rxiry)
Pearson correlation coefficient p is then estimated by p(X, Z
We use the
WVrr(r)War(rv))
hat operator to indicate estimators. For testing the null hypothesis Ho p 0 versus the alternative
3

---

-dstribution PLL
0.4
ParCor--USS
distributian
0.4
+distributtion PerCorr-QLS
0.3
danslk
0.3
0.3
ParCor-WiLS
0.
0.2
0.2
0.1
0.1
0.1
0.0
5.0
2.5
0.0
2.5
5.0
7.5
0.0
0.0
Figure 1: Data generated with SCM (2). (Left) Null distribution when X and Y are affected by the same kind of heteroskedasticity. (Middle) Null distribution when only X is affected by heteroskedasticity. (Right) Alternative distribution (assuming dependence with c o 0) when only X affected by heteroskedasticity In all plots the heteroskedasticity a linear function of the value of Z and the solid black line depicts the t-distribution under the null hypothesis.
Hi p O, we use the studentized version of the partial correlation as test statistic T(ộ) This statistic is {(n and n is the sample size.
V1-p2
diitiirrbb where & is the number of variables are conditioning on,
However, for the partial correlation test to be correct, the assumptions of OLS, in particular ho- moskedasticity, have to be fulfilled. In the next section, we investigate what happens if this assumption is violated.
3.2 Effects heteroskedasticity on partial correlation
Under heteroskedasticity, the estimator of slope regression parameters in OLS regression is still unbiased. However, the estimator of the covariance matrix of the parameter estimates can be biased and inconsistent under heteroskedasticity [Wooldridge, 2009, p.268J. This also affects the subsequent residual- based correlation test in ways that we summarize below. Proofs for and further discussions of the following statements are provided in the Supplement A.1.
Effect 1. Under the null hypothesis, the studentized Pearson correlation coefficient is t(n 2 k:) distributed if X is independent of the variables inducing heteroskedasticity in Y.
This means that if only one node is compromised by heteroskedasticity, or, more generally, the heteroskedasticity X is independent of that in Y, the type-I error rate of the t-test will not be affected. Please refer to the middle plot in figure for visualization of the associated null distribution comparison to the analytical one used by the partial correlation test. On the other hand, in the left plot the case where both X and Y are affected by the same kind of heteroskedasticity is shown, leading to a different null distribution.
Effeet 2. If X and Y are dependent and at least one is affected by heteroskedasttcctt, the detection power of the t-test might be degraded.
In particular, if X and Y are dependent, X is affected by linear heteroskedasticity, the mean of the t-distribution is closer zero compared to the distribution of weighted least squares based studentized partial correlation coefficient, which is introduced in the next section and essentially transforms the data be homoskedastic, for fixed sample size. See the right plot figure This effect is an immediate consequence of the reduced efficiency of the OLS slope estimate. Intuitively, heteroskedasticity masks the relationship between X and Y
4 Weighted least squares partial correlation test and causal discovery
We have seen that the standard partial correlation test sensitive to heteroskedastic noise since is based on an OLS regression step. Therefore, we propose to replace the OLS regression by the weighted least squares (WLS) approach which known to be able to handle non constant crror variance. We will refer to the resulting CI test ParCorr-WLS.

---

4.1 Weighted least squares partial correlation test
The idea of WLS is perform re-weighting of each data point depending on how far it is from the true regression line. It is reasonable to assume data points where the enror has low variance to be more informative than those with high error variance. Therefore. ideally the weights are chosen as the inverse variance of the associated error.
To formalize this idea, consider the linear model y X + with Ele X 0, Cov'e X] diag(a (a ):=...., for n observations. The variance the error term 6; is denoted by of- Note, as opposed the assumptions of OLS, the entries of the conditional variance matrix are allowed to differ from each other. Denote the weight matrix by W == diag( )im The WLS method
estimates 3 by solving the adjusted optimization problem
B argmin(y- Xb) Xb).
This objective is quadratic, thus we can write down the solution in closed form
3 = (xTwx)" -! xTyy.
If the true weights are known, WLS is equivalent to applying OLS to a linearly transformed, ho- moskedastic version of the data the weights have the effect of standardizing the scale of the errors. Thus, the following lemma holds. Proofs for these statement can be found, for instance, in Greene [2003].
Lemma 1. If the weights are chosen the reciprocal of the ervor variance per sample, the WLS estimator is consistent, efficient, and asymptotically normal, as well as BLUE.
Moreover, note that the weighted residuals W(y x are homoskedastic.
The goal now is approximate the conditional variance function of both X and Y in SCM (2). Since it works analogously for both cases, we illustrate the approach for X; = aZ + s(, t) Nx, i.c. we approximate o² (z,t)
Var I 2).
For that, we use a residual- -based con-paraettric estimator for the conditional variance, similar to the approach of Robinson [1987]. Motivated by the identity V Var(X) Z) = E[(X - E[X Z Z]. and noting that this is the regression of (X aZ)² on Z, the first step using OLS regression to obtain the squared residuals (X - az)². Afterwards, we use a non-parametric regression method to regress these residuals on Z and thereby predict the conditional mean by using a linear combination of the k residuals closest in Z value. For sampling index-dependes heteroskedasticity this turns into a windowing approach, which essentially smoothes the squared residuals.
Algorithm 1 details the proposed partial correlation CI test based on the feasible WLS approach that employs our weight approximation method. For the test to perform well, is crucial to know the type of heteroskedasticity. more precisely, which of the predictors the variance depends on. In practice, this kind of expert knowledge could be obtained by performing a test for heteroskedasticity, e.g. as suggested in Wooldridge [2009, p.277] or by investigating plots of the residuals. Here we make the limiting assumption that the heteroskedasticity only depends on one of the predictors or on the sampling index. Further extensions are considered in the section 6.
Now, we formulate our main assumption under which it is possible to obtain consistency result for the WLS method that uses our weight approximation method.
Assumption 1 (Heteroskedastic relationships). For each node i d in SCM (1), the skedasticity noise scaling function fú only depends one of the predictors H or the sampling index t, i.e. its domain is one dimensional, or is constant, Furthermore, it is known what depends
We also have to impose a rather technical assumption on the functions h4- Assumption (Weight approximation). Assumptions (3.3) f3om)) Robinson [1987].
The proof of the following lemma can be found in Robinson [1987].
Lemma 2. Under assumptions I and 2, the WLS estimator that uses the reciprocal of the approxi- mated variance weights is consistent and attains the correct covariance matrix asymptotically.
5

---

Assumption 3 (Technical assumptions), See Supplement A.2.
We now state the first main theorem.
Theorem Under the assumptions 3 ParCorr-WI3, CI test with estimated weights (Algorithm 1) consistent for testing the conditional independence between two potentially heteroskedastic variables X and Y conditioned on a set of variables Z.
Algorithm 1: ParCorr-WLS
Data: Expert knowledge as map with keys and and values in [false, 'sampling index', 'heteroskedastic parent H' ). where H E V. observational data with sample size for nodes X, Y, Zi,..., Zz. H, window length X, significance level
Result: boolean indicating whether there is (conditional) dependence between X and Y given Zi...., Zk resid- 0::
for node in /X. Y/ do k then
F = node
else
obtain residuals F = (fi,.... F.) by regressing node on Zi,..., Zz using OLS;
end
if Enoode] == "falsee then
append F to resid:
else if E(node] = 'sampling index' then min((,
compute weights w, = (* 2 j=j=mxx11,,
(3) for 1,..., n obtain residuals r by regressing node on Zi...., append resid:
,....,ZZ, Zk using WLS with the weights w;
else
sort such that their corresponding values of increase;
compute weights = ( E min(n,,5 {)
_ }) (3)-1 for = I,...,n;
revert the sorting in the indices of w; obtain residuals by regressing node on Z ...... Zx using WLS with the weights w;
append w- r to resid;
end
calculate studentized Pearson correlation between rx resid[0] and ry residd [1]: perform 1-test with (two-sided) significance level a, i.e. reject if [[ >> I(1 2. - k)
4.2 Extension the PC algorithm
well-known and widely used algorithm for discovering causal relationships in terms of the com- pleted partially directed acyclic graph (CPDAG) from observational data the PC algorithm as introduced in Spirtes and Glymour [1991]. It consists of two phases: The first one is concerned with learning the skeleton adjacencies based on iterative CI testing. Subsequently, a of rules is applied to determine the orientation of the found links.
To ensure consistency of this method, the following assumptions have be fulfilled. Details be found in Spirtes et al. [2000]. Let G = (V, EE) be graph consisting of a set of vertices V and a set of edges E V V. Let P denote the probability distribution of V. Assumption 4 (PC algorithm). The Causal Markov condition, Causal Faithfulness, and Causal Sufficiency are fulfilled for the SCM (1) with graph G.
Under these assumptions and if the utilized conditional independence test is consistent, it can be shown that the PC algorithm converges in probability to the correct causal structure, i.e.
lim P(G,, * G) = O,
where denotes the ground truth CPDAG and Gn the finite sample output of the PC algorithm. See Kalisch and Bühlmann [2007] for proof.

---

In the following, we need to further restrict Assumption 1 to prove consistency of the PC algorithm under heteroskedassicity.
Assumption 5 (Heteroskedastic relationships regarding PC algorithm), Assumption is further lim- ited to the case where there is only sampling index dependent heteroskedastictty or homoskedasticity,
Given Assumption 5, we can apply proposed method ParCorr- WL. in every CI test of the PC algorithm. Using Lemma I, we thus can establish consistency of the PC algorithm with ParCorr-Wi WLS if the true weights are known. Lemma yields the same result for ParCorr- WLL with estimated weights.
Theorem 2.
Under Assumptions 2,3,4,5 5 the oulput of the PC algorithm, with ParCorr-WIS CI test, converges in probability to the correct causal graph.
Note that we prove consistency of the PC algorithm with ParCorr-W CLL only in the case of sampling index dependent heteroskedasticity. We discuss the challenges to extending this more general types of heteroskedasticity in section 6.
Experiments
In the following, we conduct experiments evaluating our proposed CI test separately and in conjunc- tion with the PC algorithm. Throughout the experiments, heteroskedasticity strength refers the parameter s in the scaling functions h of linear and periodic type given by
l(x)) = + seTo x 1,20
(linear)
(xx) = + se" (ir( + (periodic).
In other words, in case of linear heteroskedasticity strength is the slope of the variance function, and for periodic heteroskedasticity refers to the amplitude of the variance function.
5.1 Conditional independence testing
We generate the data from the SCM (2) where we consider various types of heteroskedasticity, i.e. functions namely linear and periodic as given in Eq. (3). Plots of the simulated data are provided in the Supplement. Each of the types can either be Z- or sampling index dnpex-dependent. A visualization of the considered heteroskedasticctyyyyyp can be found in the Supplement A.5.
We use the Kolmogorov-Sminng Smiinov (KS) statistic to quantify how uniform the distribution of p-values is, and therefore as a metric for type-I errors, in Runge [2018]. Type-II errors measured by the area under the power curve (AUPC). The metrics were evaluated with sample size of 500 from 100 realizations of the SCM. Error bars indicate the bootstrapped standard errors.
Figure shows that ParCorr-We WL. is well calibrated in the presence of heteroskedasticity, regardless if it affects only one both of the variables X and Y. On the other hand, the ParCorr-OLe test becomes ill-calibrated in the case heteroskedastic and as expected by Effeet 1. In particular, this means that if we can choose the weights reasonably well, we are able to overcome multiplicative confounding with our proposed CI test.
Regarding power measured by AUPC, we observe rather rapid decrease for ParCorr-OLS as the heteroskedasticity increases (compare to Effect 2). Our proposed method ParCorr- WLS has higher power in the heteroskedastic scenario, and even for homoskedastic noise (heteroskedasticity strength equal zero) the power comparable to that of ParCorr-OLS. The drop power for ParCorr- WLS can be explained by an overall increase in noisiness of the data as heteroskedasticity strength is growing. Refer to the Supplement A.5 for plot of AUPC of ParCorr-OLS and ParCorr-WIs on data with homoskedastic but increasing noise.
Similar effects are present for all considered types of heteroskedasticity (see Supplement). The results remain very similar if we do not use the ground truth weights but estimate them using the window approach with reasonable window length as detailed in Section 4.1.
5.2 Causal discovery
To test our proposed CI test in a more realistic setting, we apply it within the PC algorithm to recover the causal graph from simulated observational data. We build upon the PC-stable algorithm
7

---

Perent-dependent heteroskedasticnny affecting only
Parent- detent-dependettt heteroskedastehty affecting and Y
Insar
puriodi
Iinerr
perlodic
ParCom-LLS
-..-Par CaCm--OLLS
ParCom-MUL
ParCom-WLL
PerCam-AS
ParCarm-00.5
0.4
FarCardAyS
ParCorn M1.5
FarCor-41.5 (nt true weights)
(wi true weigrtt)
(wl true weiggtes)
Owr true wigghts)
0.2
0.2
0.2
2
44
0.0
0.0
0.0
hoter caddrcadddaniccyy strength
badettccedeellily cit strength
Mataroseeaariic Cily strengh
Ister cocorsssiiciii stroingh
Inear
puriodic
Iinerr
parlodic
-e DLS
1.0
1.0
ParCorm-M
ParCor-4iLs
S0.8
(Ai true weights)
Par-arrr OLS
Fa--aarC DLS
ParCor--LS QLS
ParCort WLS
FarCori WUS
ParCon WLS
0.6
Per Com--WLL (e) thue wcighh52)
0.6
ParCues- 48.55
"Wf true weights)
ParCortMLS
(N frum weights)
hate bycodeatichy abrengg
butorracedenttetyy strength
hetaroszeoasicii strengts
hete cileesscddssiikp strength
Figure 2: Performance of partial correlation CI tests for dependence (c 0.5) and conditional independence (c = 0) between X and Y given Z. In all plots the heteroskedasticity is a function of the confounder Z and only affects X (left two columns) or both X and Y (right two columns). Shown are KS (top row) and AUPC (bottom row) for different strengths of heteroskedasticity for linear (left), periodic (right) noise scaling functions. The ground truth weights are used for WLS, or the weights are estimated using the window approach with window length 10 as detailed in section 4.1. sample size of 500 is used and the experiments are repeated 100 times.
implementation within the Tigramite software package [Runge et al., 2019b] which is published under the GNU General Public License.
The data for these experiments was generated with SCM (1) in the following way. Given a random affected ground truth graph, we fix percentage of heteroskedassicityy afffeee nodes which are then selected uniformly at random. Throughout the experiments this percentage set to 0.3 reduce the chance of parent and child being affected by the same kind of heteroskedasticity. For these affected nodes, we choose as a heteroskedasticity type linear or periodic with equal probability (Eq. (3)) and let the noise variance h either depend on one randomly selected parent or the sampling index. We also set fixed strength per experiment and investigated the effect of increasing the strength on the performance of our method compared to the PC algorithm with PurCor--OLS All linear dependencies have coefficient e = 0.5. We tested for various small to medium sized causal graphs (see also Supplement).
In these experiments, we estimate the weights based on the expert knowledge as required by As- sumption 1, namely which node is affected by heteroskedasticity, and whether the noise variance depends on the sampling index or another node. In the case of the heteroskedasticity depending on another node, we also need to know which node is. We also compare with the PC algorithm which uses Corr-WI.S based on the ground truth weights for all direct heteroskedastic relations. Note, however, that this does not take care of indirect heteroskedasticity due to heteroskedastic parents that are not part of the conditioning set. Note that this challenging setting is outside of the stricter Assumption for which the PC algorithm is consistent.
Figure 3 shows, similar to the experiments of the previous section, that the PC algorithm with ParCorr- WLS continues to have rather small false positive rate (FPR) even though the heteroskedasticity strength increases. In contrast, ParCorr-LLe shows an increase in FPR as the heteroskedasticity strength increases. Additionally, even in this rather complicated setting, we see that the true positive rate (TPR) can be improved by using our method. The PC algorithm with ParCorr-WISS has an average runtime of 0.29 seconds on homoskedastic data compared to 0. seconds for ParCorr-OL.S evaluated on AMD 7763.
Discussion and Outlook
In this work we relaxed the common assumption of a constant variance by explicitly allowing for hetcroskedasticity as a multiplicative scaling of noise in an otherwise lincar SCM. Our proposed partial correlation test based on weighted least squares regression is a linear, computationally fast and
8

---

0.90
ParCorr-OLS
0.013
ParCornWLSS
0.960
PLS weights)
adi precision
0.955
0.012
0.950
0.80
ParCor-OLS
0.011
ParCor-OLS
ParCor-WLLS
0.945
ParCorr-WL
ParCorr-WLS
ParCor--WLS
(w/ true weights)
0.010
0.940
(w/ true weights)
2. 4. Heteroskedasticity strength
5.
0.
4.
5.
Heteroskedasticity strength
Heteroskedasticcty strength
Figure 3: Results for the PC algorithm with the standard ParCorr-OLS CI compared to that with the proposed ParCorrWLLS test. Shown are adjacency TPR (left) and FPR (middle), as well as the adjacency precision (right) for increasing strengths of heteroskedasticity. The graph has 10 nodes and 10 edges, sample size of 500 is used. The significance level a is set to 0.05. The experiment is repeated 500 times. Errorbars show standard errors. Estimated weights with a window length of or ground truth weights are used for WLS.
easy to implement method and constitutes useful standalone method for various data science tasks, but our focus here is its use in causal discovery. The main strengths of our approach are that is able to produce more reliable results both as a CI test and in causal discovery in the presence of heteroskedasticity than the standard partial correlation variant. More reliable results here refers to controlled false positives well higher detection power. Furthermore, the suggested adaptations do not compromise calibratedness and power of the CI test on homoskedastic data.
The main weakness of our method is that the CI test requires substantial expert knowledge in As- sumption and for the consistency of the PC algorithm the even stricter Assumption 5. Assumption 1 requires that the heteroskedasticity only depends on one of the predictors or on the sampling index, i.c. its domain one dimensional. Furthermore, one needs to know which of these types is the case. Assumption 5 only allows for sampling index heteroskedasticity The reason for this is that there can be indirect heteroskedasticity the node X Y that is not induced by a parent of X or Y, but by some other ancestor. This heteroskedasticity then propagates through the causal graph and essentially makes X or Y heteroskedastic whenever this path the expert knowledge does not tell us about this heteroskedasticity and we would not be able to apply ParCorr-V CLL to remove it. See also Section A.4 the Supplement a detailed treatment cases in which the CI tests within the PC algorithm are consistent under more general heteroskedasticity forms. In the following, we discuss alternatives and further avenues of research.
blocked by the conditioning set. In this case,
In future work one may alter the required expert knowledge and weight approximation scheme to overcome the issues induced by indirect heteroskedasticity. Here, a possible remedy could be an iter- ative approach using information about causal relationships from earlier steps. An alternative would be use the CMIknn CI test [Runge, 2018] that is able fully treat heteroskedastic multiplicative confounding. However, this increased generality comes at the price of reduced detection power, also see Figure 10 in the Supplement.
Another open question is how to further improve the weight approximation method. For instance, by iteratively repeating the regression and smoothing steps. Another important consideration the choice of the window length. Linear properties of the variance function allow us to use a larger window length. However, if the variance function shows high variability, crucial information might be lost if the chosen window length is too large. Model selection criteria might be employed to alleviate this problem.
Furthermore, one can consider extensions of ParCorr- WLS to multiple heteroskedastic influencing factors, e.g.. Spokoiny [2002] extend the residual-basee conditional variance function estimation to dimensions larger than one. Another multidimensionall method based on differences is discussed in Cai et al. [2009].
would also be interesting to explore using generalized least squares to able to account for additional correlation of the residuals. Potentially, the weighted least squares partial correlation

---

coefficient could also be combined with permutation- based test to extend the method to problems with non- Con-Gaussian noise.
Concluding, our proposed ParCor--WLS CI test makes constraint-based causal discovery methods better applicable to real world problems where the assumption of homoskedasticity is violated, such as in climate research or neuroscience. Ethically, we believe that our rather fundamental work has low potential for misuse.
Acknowledgments and Disclosure of o Funding
WG was supported by the Helmholtz Al project CausalFlood. UN and JW supported by grant no. 9481 12 Causal Earth of the European Research Council (ERC). This work used resources the Deutsches Klimarechenzentrum (DKRZ) granted its Scientific Steering Committee (WLA) under project ID bd1083. We thank the anonymous reviewers for their helpful comments.
References
L. D. Brown and M. Levine. Variance estimation in nonparametric regression via the difference sequence method. The Annals of Statistics, 35(5):2219-2232, 2007.
Cai, M. Levine, and Wang. Variance function estimation in multivariate nonparametric regression with fixed design. Journal of Multivariate Analysis, 100(1):12--13,, 2009. ISSN 0047-259X.
R. Carroll. Adapting for heteroscedasticity in linear models. The Annals of Statistics, 10(4):1224 12322222211122--1111121112322311 1982.
D. M. Chickering. Learning Equivalence Classes of Bayesian-Network Structures. J. Mach. Learn. Res., 2: 445-498, 2002. ISSN 1532-4435.
1. AbeettUUpoof and Y. Deng. Causal discovery for climate research using graphical models. Journal of Climate, 25:5648-5660 2012.
J. Fan and Q Yao. Efficient estimation of conditional variance functions in stochastic regression. Biometrika, 85 (3):645 16065-660, 1998.
W. Greene, Econometric Analysis. Pearson Education, fifth edition, 2003. ISBN 0-13-066189-9
Hall and R. J. Carroll. Variance function estimation regression: effect estimating the mean, Journal of the Royal Statistical Society: Series B (Methodological), 51(1):3-14, 1989.
A. Hayes. Using beteroskedaasiiityycossse standard error estimators in ols regression: An introduction and software implementation. Behavior reseurch methods, 39:709-22, 12 2007.
B. Huang, K. Zhang. J. Zhang. J. D. Ramsey, R. Sanchez-Romero, C. Glymour, and B. Schölkopf. Causal discovery from helerogeneouv/noostttiionnn data. Mach. Learn. Res., (1(89):1-53, 2020.
M. Kalisch and P. Bühlmann. Estimating high-dimensional directed acyelic graphs with the pe-algorithm. Journal of Machine Learning Research, 8(22))613 636, 2007.
M. Kalisch and P. Bühlmann. Robustification of the pe-algorithm for directed acyelic graphs. Journal of Computational and Graphical Statistics, 17(4))773-789, 2008
A. Lawrance. On conditional and partial correlation. The American Statistician, 146-)1446149, 1976.
J. S. Long and L. H. Ervin. Using heterosecdasticity consistent standard errors in the lincar regression model. The American Statistician, 54(3):217-223 2000.
B. Merz, G. Blöschl, S. Vorogushyn, F. Dottori, J. C. Aerts, Bates, M. Bertola, M. Kemter, H. Kreibich, U. Lall, Causes, impacts and patterns disastrous river floods. Nature Reviews Earth Environment. 2(9):592-600 2021.
J. M. Mooij, S. Magliacane, and T. Claassen. Joint causal inference from multiple contexts. Journal of Muchine Learning Research, 21(99):1-108 2020.
M. H. Neumann. Fully data- daiven nonparametric variance estimators. Statistics: A Journal of Theoretical and Applied Statistics, 25(3):189-211 1994.
10

---

J. Peters. D. Janzing, and B. Schölkopf. Elements of causal inference: foundations and learning algorithums. MIT Press, Cambridge, MA, 2017. ISBN 9780262037310.
T. Proietti. Seasonal specific structural time series. Studies in Nonlinear Dynamics & Econometrics, 8(2). 2004.
P. M. Robinson. Asymptotically efficient estimation the presence of heteroskedasticity of unknown form. Econometrica: Journal of the Econometric Society, pages 875-891, 1987.
J.P. P. Romano and M. Wolf. Resurrecting weighted least squares. Journal of Econometrics, 197(1):1-1, 2017.
J. Runge. Conditional independence lesting based nearest-neighhor estimator of conditional mutual information.
Storkey and F. Pérez-Cruz editors, International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain, volume 84 of Proceedings of Machine Learning Research, pages 938 947. PMLR, 2018.
J. Runge, Bathiany, E. Bollt, G. Camps- Valls,Valls, Coumou, E. Doyle, C. Glymour, M. Kretschmer, M. D. Mahecha, J. Muñoz-Mari, et al. Inferring causation from time series in earth system sciences. Nature communications, 10(1):1-13, 2019a.
J. Runge, P. Nowack, M. Kretschmer, S. Flaxman, and D. Sejdinovic. Detecting and quantifying causal associations in large nonlinear time series datasets. Science advances, s(eauu44, 2019b.
P. Spirtes and C. Glymour. An Algorithm for Fast Recovery of Sparse Causal Graphs. Soc. Sei. Comput. Rev., 9 (2):62-72, 1991. ISSN 0894-4393.
Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT press, 2nd edition, 2000.
Spokoiny. Variance estimation for high-dimensional regression models. Journal Multivariate Analysis, 82 11)3111-133, 2002. ISSN 0047-259X.
N. Tagasovska, V. Chavez-Demoulin, and T. Vatter. Distinguishing cause from effect using quantiles: Bivariate quantile causal discovery. Proceddiggs of the 37th International Conference on Machine Learning, ICMI. 2020, 13-18 July 2020, Virtual Event, volume 119 of Proceedings of Machine Learning Research, pages 9311 93233. PMLR, 2020.
R. Wilcox and Muska. Inferences about correlations when there is heteroscedasticity The British journal of mathematical and statistical psychology, 54:39 47, 06 2001.
J. M. Wooldridge. Introductory Econometrics: Modern Approack. ISE International Student Edition. Wouth-Westerm, 2009, ISBN 9780324581621.
Xu, A. Marx, O. Mian, and Vreeken. Causal inference with heteroscedastic noise models. In Proceedings of the AAAI Workshop on Information Theoretic Causal Inference and Discovery (ITCI'22, Proceedings of the AAAI Workshop on Information Theoretic Causal Inference and Discovery, pages 9311-9323, 2022.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Sections and 6, well as the assumptions, specifically Assumption 1.
(c) Did you discuss any potential negative societal impacts of your work? [Yes] See Section 6.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set assumptions all theoretical results? [Yes] See Assump- tions 1-5.
(b) Did you include complete proofs of all theoretical results? [Yes] We provide references and intuitive justifications for the results in the main text and complete proofs in the Supplement.
3. If you experiments..

---

(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material as a URL)? [Yes] The code and instructions on how to execute it are provided in the supplemental material. The code includes a method for generating the synthetic data.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] All hyperparameters
listed in Section 5.
(c) Did you report error bars (e.g., with respect the random seed after running experi- ments multiple times)? [Yes] See Figures 2 and 3.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Section
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]] See Section 5.
(b) Did you mention the license of the assets? [Yes] See Section
(c) Did you include any new assets either the supplemental material or as URL? [Yes]
The code is included in the supplemental material.
(d) Did you discuss whether and how consent obtained from people whose data you're
using/curating? [N/A]
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
12


Reviews:
[{"review_id": "yDCHKcdqaXb", "reviewer": "Reviewer_hnmZ", "summary": "A conditional independence test is proposed based on partial correlation, but adapted to the case of heteroskedasticity. It is also studied how this test can be used in the PC algorithm. All results a rigorously proved, as well as illustrated using experiments.", "questions": "* On line 102, I was confused by the $[0, \\infty)$. This seems to be the domain of $t$, though as far as I can tell, this is always a discrete sampling index, not a continuous value.\n\n* On line 270-271, it is described that the case of \"parent and child being affected by the same kind of heteroskedasticity\" is being avoided.  I suppose this is simply because that could lead to much larger heteroskedasticity strengths than what is supposed to happen for a given value of the parameter $s$. Or is there a  more fundamental problem with this case?\n\nMinor:\n\n* Line 122-123: \"the\" missing before \"hat operator\"", "limitations": "Yes", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This contribution is of high theoretical quality. It is also of great practical relevance, as homoskedasticity assumptions are often violated in real datasets. There are currently some major limitations (which are clearly mentioned in the paper), mainly that the proposed partial correlation test only works if for each variable, there is at most one source of heteroskedasticity, and for theoretical consistency of the PC algorithm, that this source is not another variable. Improvement in these areas is listed as future work.\n\nThe paper is clearly written; I especially appreciate the illustrations of the problems in Figure 1. The application to causal discovery is original, though if I understand correctly, the test itself is quite similar to existing work on heteroskedasticity.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "wggSWbb4Wx", "reviewer": "Reviewer_9b5D", "summary": "In this work, authors consider heteroskedasticity in a structural causal model (SCM) and develop an adapted weighted least-squares (WLS) partial correlation variant as a conditional independence (CI) test which is able to deal with heteroskedastic noise with expert knowledge about the heteroskedastic relationships are given. Combined with the PC algorithm, the proposed CI test can also identify the causal relationships between variables under some assumptions regarding the heteroskedastic relationships. Both theoretical consistency analysis and simulation experiments demonstrate that the proposed CI test achieves better false positive rates and also improves upon the detection power as compared to the standard partial correlation CI test. ", "questions": "1. The authors did not mention the computational complexity of the proposed causal discovery algorithm (i.e., inject the proposed CI test into the PC algorithm), and it would be great if the authors could discuss the computational complexity of the proposed algorithm.\n\n2. Since only simulated experiments are conducted in the paper, the reviewer wonders whether the proposed algorithm works well on real-world datasets?", "limitations": "Authors have discussed the Limitations and societal impacts in the Conclusion Section.\n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. The authors propose an adapted WLS partial correlation test for conditional independence testing in the presence of heteroskedastic noise.\n\n2. The authors provide a detailed analysis of the WLS method's consistency with the proposed weight approximation method, and further inject the proposed CI test into the PC algorithm to learn the causal structure.\n\n3. Experiments on simulated datasets demonstrate that the proposed WLS partial correlation CI test achieves better performance than the OLS partial correlation CI test in the presence of heteroskedastic noise. Combined with the PC algorithm, the proposed WLS partial correlation CI test has a smaller false positive rate in identifying causal relationships.\n\nWeaknesses:\n\n1. The idea of using the inverse of error's variance is intuitive, however, it would be great if the authors could further discuss why WLS can deal with heteroskedastic noise and some possible better re-weighting strategies in detail.\n\n2. The authors formulate several assumptions in the paper, however, the heteroskedastic relationships require expert knowledge. The reviewer wonders if the authors could discuss further how to relax the heteroskedastic relationships assumption.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Uezy6JIzXmu", "reviewer": "Reviewer_vG5d", "summary": "This paper proposes a conditional independence test in a linear Gaussian setting but the variance of the error term may not be a constant but affected by a function of a subset of its parents or sampling index (e.g., individual-dependent variance or a time-series like setting). Hence, the authors designed heteroskedasticity-aware weighted partial correlation where the weights come from inverse of variances, which is estimated using a windowing approach. However, this test works with expert knowledge about heteroskedasticity (whether it is caused by sampling index or heteroskedastic parent). In a causal discovery setting, not knowing the parents of variables, it only assumes that the heteroskedasticy comes from sampling index.", "questions": "\n- What if we can use non-parametric conditional independence test? Will it be (certainly) even weaker than OLS?\n- Can we combine both heteroskedastic parent and sampling index? In such a case, can we just use “sampling index” approach? (or impossible to combine them?) \n- In estimating weights, no clear justification is given for sorting residuals with respect to the values of H. I understand that being closer in H values implies having similar variances. But I would like to see more formal statements on this. \n- I am not sure whether a similar idea is employed in a similar or different setting. Like, transforming data in a way to be more homoskedastic, and then applying ParCorr-OLS.\n", "limitations": "\n- Requires prior knowledge on the reasons for heteroskedasticity.\n", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nStrengths\n- well-motivated problem and a theoretically correct solution (under some assumptions.)\n- a simple solution to an important problem.\n\nWeaknesses\n\n- As the authors acknowledged, Assumption 1 (and 5 for PC algorithm) requires expert knowledge.\n\nIt seems hard for me to find weaknesses of the paper considering that the authors did a right thing to solve a problem and acknowledged limitations clearly that are not due to the authors but comes from the nature.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}]

Author rebuttal:

Decision: Accept