Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
56,446 characters · 18 sections · 53 citation commands
Testability of Reverse Causality Without Exogenous Variation
{\it Keywords:} Endogeneity; Reverse causality; Conditional Independence; Causal Dis-covery; Reproducing kernel Hilbert spaces.
Endogeneity is a central problem in econometric models which potentially invalidates estimates of causal effects. Reverse causality, where the dependent variable causes the independent variable, is one source of such endogeneity. The aim of this paper is to show that reverse causality is testable under mild assumptions and without relying on exogenous variation. We build on hoyer09anm who provide a link between nonlinear model structure and causality, namely, that a nonlinear relation between cause and effect leads to observable signals about causal direction using observational data. While their theoretical results are striking, their results assume homoskedastic errors and do not generalize to settings with additional control variables. We show that this assumption can be relaxed to allow for heteroskedasticty with respect to additional control variables, as is commonly the case in econometric applications. Our primary conditions for achieving testability are twofold: First, we necessitate a nonlinear relationship between the dependent variable and the regressor (with no restrictions on how controls are incorporated into the model). Second, we impose additively separable errors.
Following hoyer09anm, we also demonstrate how our testability result can be applied in empirical practice. Specifically, we show that identification of the causal direction is equivalent to a conditional independence test of covariates and error terms given control variables. We make use of conditional independence tests based on kernel mean embeddings, i.e., maps of probability distributions into reproducing kernel Hilbert spaces (RKHS) muandetetal16kernel. Intuitively, this corresponds to approximating conditional distributions with unconditional ones by weighting with an appropriate kernel, and evaluating their covariance in an RKHS. The method can detect nonlinear dependencies.
We consider two formal applications of our testability result. First, we explore testing for causal direction based on conditional independence. As already indicated in the related literature, achieving exact size control can be challenging within this framework, and we provide a detailed discussion of this issue below. Second, we conduct causal discovery, where we remain agnostic about the causal direction and compare test statistics for two rivaling models to gain insight into which one represents the true causal structure peters14.
In Monte Carlo simulations, we investigate the power of our approach to detect reverse causality. We see that the degree of nonlinearity increases the power to detect reverse causation. The procedure has surprisingly high accuracy in detecting the true causal direction even under moderate form of nonlinearities and can be powerful even in linear models under restrictions on the distribution of errors terms, which is a result described by shimizu06lingam. Furthermore, we provide an empirical illustration using data from the German Survey of Income and Expenditure. We show that our algorithm can infer from purely observational data that work experience is a causal driver for income, not vice versa. Substantively, this is not a surprising result; however, the fact that it can be inferred without exogenous variation is.
\paragraph{Related literature}
Our test rests on the idea that $ X $ causing $ Y $ implies an independence between the error of a regression of $ Y $ on $ X $, and $ X $. This idea goes back to engleetal83, who propose a definition of an exogenous relation in terms of conditional densities. In particular, they argue that, if a joint probability density of two random variables $ Y $ and $ X $ factorizes as $ f(Y,X) = f(Y|X)f(X) $ and the conditional density $ f(Y|X) $ is invariant to changes in the marginal density $ f(X) $, then $ X $ is called “super exogenous" (p. 278). Statistical tests for the notion of “super exogeneity" are proposed by faverohendry92,englehendry93,hendrysantos10. These tests rely on analyzing to what extent parameter values are sensitive to exogenous interventions on the purported cause. Thus, their results specifically rely on exogenous variation (e.g. in the form of instrumental variables) whereas the approach at hand does not require such variation.
The problem of identifying causal structure from non-experimental data is receiving considerable attention in the causal machine learning literature mooijetal16,eci17,scholkopf2021toward. In its bivariate form, the problem is concerned with deciding whether a variable $ X $ is causing $ Y $ or vice versa solely based on a non-experimental joint probability distribution of the two variables. Without making any assumptions regarding the true underlying data-generating process, no identification is possible.\footnote{Previous work shows that the causal direction cannot be identified without making further assumptions. peters12 proves that for every joint distribution of two variables, $ X $ and $ Y $, there is a model $Y = h(X,\varepsilon), \; \text{with} \; X\mathrel{\perp\mspace{-10mu}\perp}\varepsilon $ with $ h $ a measurable function and $ \varepsilon $ a real-valued noise variable. The roles of $ X $ and $ Y $ can be easily interchanged showing that the joint distribution itself does not identify the causal direction in this most general form.} shimizu06lingam show that non-Gaussianity of observed variables leads to identifiability. Subsequently, hoyer09anm show that nonlinearity of $ h $ can play a similar role as regards the identifiability of the causal direction as non-Gaussianity. If the true model is of a nonlinear form, one can infer the causal direction without making any assumptions about the distribution of the error.
hoyer09anm and mooijetal16 discuss inference of the causal direction between two random variables (cause and effect) from observational data. peters14 constitutes a theoretical extension of these methods to more than two variables. The paper at hand falls between these two strands as it accounts for more than two variables, yet its primary concern is the causal directionality between a subset of just two of them. The remaining variables $W$ serve as controls.
In this paper, we rely on RKHSs and kernel mean embeddings, which are not widely used in the econometrics literature, albeit with notable exceptions. carrasco00continuum_gmm and carrasco07 discuss the usefulness of RKHS theory given infinite number of moment conditions. singh19kernel_iv study the use of kernel methods in the context of instrumental variable (IV) methods. They use kernel mean embeddings of the conditional distribution of the covariates given the instrument to propose a nonlinear extension of linear IV implementations. zhang2020maximum propose kernelized moment restrictions to estimate nonlinear IV estimators. grunewalder12 analyze connections between kernel mean embeddings and vector-valued functions to analyze Markov decision processes. flaxman15 use kernel mean embeddings to analyze who cast their votes for Obama in the 2012 US presidential election.
This paper is also related to a strand of the literature, which make use of exogenous variations to detect endogeneity of regressors. The idea to make use of instrumental variables to detect endogeneity was originally proposed by hausman78specification. More recently, blundell07exog and breunig15 provide exogeneity tests using instrumental variables for nonparametric models with additively separable errors, feve18estimation and breunig18specification for models with nonseparable errors.\\
The remainder of the paper is organised as follows. Section (ref) establishes testability of reverse causality under nonlinear regression functions and introduces the RKHS test for independence. In Section (ref), we analyze finite sample power of our RKHS procedure in a Monte Carlo simulation study. Section (ref) provides an application of our method to empirical data. Appendix (ref) provides a proof of our main testability result. Appendix (ref) gives a review of the construction of RKHS. Appendix (ref) contains additional Monte Carlo simulation results.
We show how to test for reverse causality between two variables $ X $ and $ Y $ in the presence of additional covariates $ W $. First, we introduce the model, discuss how the model specification relates to the existing causal discovery literature, and derive testable implications. Second, we present the conditional independence test that is a central component of the test. Third, we present the implementation of the test.
Consider a model where observable continuous scalar variable $ X $ causes observable scalar variable $ Y $ in the presence of the vector of covariates $W$ (which we refer to in the following as the model):
where $\varepsilon$ are unobservable variables and $ \sigma(\cdot)$ some strictly positive function. In addition, we assume $E[\varepsilon]=0$ without loss of generality.
Note that the error $ U $ is additively separable. The additive separability of $ U $ precludes the dependence of marginal effects on unobservables except through a dependence via $W$. The model allows for heteroskedasticity of the error term $U$ with respect to the control variables $ W $. In particular, model equation (ref) implies $ U \mathrel{\perp\mspace{-10mu}\perp} X|W $, which is also known as conditional exogeneity whitechalak2010condtlexo. It corresponds to the unconfoundedness assumption in the treatment effects literature imbensrubin15 and is also closely related to the special regressor assumption lewbel14.
The main idea of this paper is to study the conditions under which this model is distinguishable form reversed analog without relying on exogenous information. The reverse model, where $ Y $ is causing $ X $, again in the presence of the vector of covariates $W$, is defined as
where $\widetilde{\varepsilon}$ are unobservable variables and $ \widetilde{\sigma}(\cdot)$ some strictly positive function.\\
We denote the probability density function of a random vector $V$ by $f_V$ and make the following assumptions.
The nonlinearity of regression functions (Assumption (ref)) can be used to make inference on the causal structure of a model is obtained first by hoyer09anm. We extend their work by allowing for additional control variables $W$ and also considering heteroskedasticity of the error term with respect to these covariates. Intuitively, nonlinearity of $h$ ensures that the error terms in the reverse model are not independent of the regressor, which provides power of the test. While linear models are used in many economic applications, they are typically seen as approximations of nonlinear relationships between dependent variable and regressors.
We are now in a position to formulate the main theorem of this paper.
The proof of this statement can be found in Appendix (ref). For the proof of the result, we build on hoyer09anm and extend their result to allow for control variables $W$ with additional form of heteroskedasticity. Intuitively, it is shown that causal and anticausal models can only exist simultaneously under very specific circumstances: if the joint distribution of $(Y,X,W)$ satisfies both a causal and a anticausal model, we can show that densities $\log f_{X|W}$ and $\log f_{\varepsilon}$ of the causal model have to satisfy the linear inhomogeneous differential equation (ref). The solutions of this differential equation restrict the log density of $ X $ given $W$ to lie in a (specific) three-dimensional space, although a priori the (generic) space of possible log marginal densities of $ X $ given $W$ is infinite-dimensional hoyer09anm. To achieve testability of reverse causality, we need to exclude those distributions that satisfy the differential equation, i.e., we need to exclude specific combinations of $\log f_{X|W}$, $\log f_{\varepsilon}$, and $h(\cdot)$.
The next result provides a more concrete formulation of Theorem (ref) and how it can be applied to model specification testing. This corollary follows immediately from Theorem (ref) and its proof is thus omitted.
Corollary (ref) allows identification of the causal direction from observational data by analyzing to what extent the independence of errors and covariates holds. The nonlinearity of $ h $ and the additive separability of the error term, $ U $, give the proposed test power.
Note that even when $h$ is linear ($Y = X + U$, simplifying by abstracting from $W$), it is only possible to rewrite this as a reverse model, $X = Y + \widetilde{U}$, with $Y$ independent of $\widetilde{U}$, if all variables are Gaussian. If at most one of the exogenous variables is non-Gaussian, residual and purported cause are dependent in the reverse model shimizu06lingam.
Algorithm (ref) shows detailed steps of the implementation of the test. In words, after some pre-processing (Step 1 and 2), we propose estimating a nonlinear model in both directions, i.e. with $ Y $ and $ X $ as dependent variables respectively (Step 3), calculating residuals for both models (Step 4) and testing for conditional independence using a kernel conditional independence test (KCI) introduced by zhang11conditional (Step 5). We discuss this test statistic, see eq. (ref), and its development in Section (ref); see the detailed discussion in Section (ref).
Alternatively, instead of imposing model (ref) as a maintained hypothesis, we can also remain agnostic about the true causal relationship and infer which of the models is the correct one. Formally, this requires assuming that either model (ref) or (ref) hold. This is formalized in the following corollary which follows immediately from Theorem (ref) and its proof is thus omitted.
Algorithm (ref) shows detailed steps of the implementation of the bivariate causal discovery which is motivated by (ref). Steps 1 to 4 are the same as in (ref). In step 5, we compute the test statistics corresponding to two conditional independence tests: one for model (ref) and one for (ref). The relative size of the resulting test statistics is informative about which model is the correct causal model (Step 6).
This section introduces the concept of Hilbert Space embeddings of probability distributions and their use for (un)conditional independence testing of random variables. Since this notion is not common in the econometrics literature and conditional independence testing forms a central part of the proposed algorithm, we discuss the procedure in detail. We first intuitively introduce important underlying concepts such as feature maps, reproducing kernel Hilbert spaces, etc. keeping technical details to a minimum before turning to how these constructs can help to formulate a conditional independence test. See Appendix (ref) for the formal statements.
To introduce the usefulness of a feature map, consider the following problem. Terms used loosely in this paragraph are precisely defined below. Imagine you want to distinguish between two groups of subjects each characterized by two dimensions, say $x=$ weight and $y=$ height, by using a linear classifier (i.e. a linear regression that serves as a boundary between the two classes). If the data looks like those in Figure (ref)(a), a linear classifier will perform poorly since there is no linear decision boundary that it could uncover. A solution to the problem lies in mapping the data from two-dimensional input space to a higher-dimensional feature space by introducing an additional feature $ z = x^2 + y^2 $ that complements existing features $ x $ and $ y $ (here the map is from a two-dimensional to a three-dimensional space; in practice the feature space will have many more dimensions). In this higher-dimensional space, there is a linear boundary that separates the two classes, see Figure (ref)(b). This example is adopted from dlp16.
Similarly to the linear classifier in Figure (ref)(a) that does not succeed in distinguishing between two classes that are separated by a nonlinear decision boundary in input space, the (linear) covariance between two random variables does not succeed in detecting nonlinear statistical dependencies. Mapping the data from input to feature space enables the exemplary classifier to linearly describe the decision boundary in feature space despite it being nonlinear in input space. Similarly, one can use the theory on reproducing kernel Hilbert spaces (RKHS) to construct a representation of marginal and conditional probability distributions in higher-dimensional feature space. The covariance operator between two random variables in that feature space is then informative about nonlinear dependencies in input space. In sum, any linear algorithm in high-dimensional feature space corresponds to a nonlinear algorithm in input space. Crucially, inner products between feature space representations can be estimated without knowing the exact feature representation itself (the so-called `kernel trick'). We now turn to a formal definition of a RKHS and kernel mean embedding of probability distributions.
In practice, instead of manually defining a set of appropriate features (such as $ z = x^2+y^2 $ in the previous example), flexible functions can be used to define the feature map. Formally, we define a feature map $ \Phi $ from input space $ \mathcal{X} $ to the space of functions $ \mathbb{R}^\mathcal{X} $:
where $ k $ is a positive-definite kernel\footnote{A positive-definite kernel is a kernel with an associated kernel matrix $ K $, which has entries $ K_{ij} := k(x_i,x_j) $, that is positive-definite.}, such as the Gaussian kernel, which is defined as
for arbitrary vectors $ v $ and $ v' $ and a bandwidth parameter $ \lambda>0$, where $ \left\lVert\cdot\right\rVert_{\ell_2} $ denotes the $\ell_2$ norm. Each data point can thus be richly represented by its similarity (defined by the kernel) to all other data points. It can be shown that the inner product of two such feature maps in an RKHS reduces to an evaluation of the kernel itself kernels01:
This result shows that the inner product of possibly infinite-dimensional feature representations, $ \langle \Phi(x),\Phi(x')\rangle $, can be evaluated through the kernel $ k $ without making the feature representation explicit (the so-called `kernel trick' in machine learning). Any algorithm or other data processing technique that relies on calculating inner products between data representations can be `kernelized,' i.e. transformed into a nonlinear algorithm by mapping the data into a higher-dimensional Reproducing Kernel Hilbert Space. The covariance, which can be defined as a dot product, falls into this category.
Instead of representing a specific data point by means of a feature vector, we subsequently intend to represent a whole probability distribution in terms of a higher-dimensional vector. One way to think about this procedure intuitively is to note that probability distributions can be characterized uniquely by an infinite sequence of their moments. Thus, the elements of the infinite-dimensional feature vector can be populated by moments of increasing order when embedding a probability distribution in the RKHS, which gives rise to a unique representation of the probability distribution.
The conditional independendece test we use relies on a characterization of conditional independence as a vanishing partial cross-covariance operator between two RKHSs. To get an intuition, consider an analogy to the characterization of conditional independence for jointly Gaussian variables in terms of vanishing partial correlation. First note that, for jointly Gaussian variables $ (Z_1,Z_2,Z_W) $, the conditional independence, $ Z_1 \mathrel{\perp\mspace{-10mu}\perp} Z_2 | Z_W $ can be characterized as the correlation between $ Z_1|Z_W $ and $ Z_2|Z_W $ being zero. Partial correlation is a linear concept defined by the orthogonality of linear maps of $ Z_1 $ and $ Z_2 $ on the space orthogonal to $ Z_W $. It can only characterize conditional independence for jointly Gaussian variables because of the linearity of the underlying maps. Intuitively, one can extend the results to apply to nonlinear dependence of arbitrarily distributed random variables if such maps can be described more flexibly. We have seen how maps of data into higher-dimensional RHKS enables the use of linear algorithms to study nonlinear relationships. This reasoning also underlies the following characterization of conditional independence for arbitrarily distributed random variables.
daudin80 establishes the equivalence
for properly chosen function spaces, based on function spaces $ \mathcal{F}_{\widetilde{X}} := \big\{{f} \in L^2_{\widetilde{X}} : \; \mathbb{E}[{f(\widetilde{X})}|W]= 0\big\} $ and $ \mathcal{F}_{\widetilde{U}} := \big\{{g}: \; g(\widetilde{U}) = \check{g}({U}) - \mathbb{E}[\check g(U)|W]\text{ where } \check{g} \in L^2_U \big\} $ where for any random variable $Z$ the Hilbert space $ L^2_Z =\{f:\, \mathbb{E}[f^2(Z)]<\infty\} $.
Throughout the remainder of this section, we consider continuous random variables $ X $, $ U $ and $ W $ with domains $ \mathcal{X} $, $ \mathcal{U} $ and $ \mathcal{W} $, and with positive definite kernels $ k_\mathcal{X} $, $ k_\mathcal{U} $, and $ k_\mathcal{W} $ defined on these domains. These give rise to RKHSs $ \mathcal{H}_\mathcal{X} $, $ \mathcal{H}_\mathcal{U} $ and $ \mathcal{H}_\mathcal{W} $ respectively. Further, we make use of the notation $ \widetilde{X} = (X,W) $ and $ \widetilde{U} = (U,W) $ and define $ k_{\widetilde{\mathcal{X}}} = k_\mathcal{X}\times k_\mathcal{W} $ and corresponding RKHS $ \mathcal{H}_{\widetilde{\mathcal{X}}} $. Denote the feature maps corresponding to RKHS $ \mathcal{H} $ with $ \phi_\mathcal{H} $. To reduce the complexity of Daudin's equivalence result, zhang11conditional show that restricting function spaces of $ {f} $ and $ {g} $ to lie in RKHSs $ \mathcal{H}_{\widetilde{\mathcal X}} $ and $ \mathcal{H}_{\mathcal U} $. This restrictions are sufficient to derive an estimable statistic in terms of the Hilbert Schmidt norm of a partial cross-covariance operator:
where $ \Sigma_{\widetilde{X}U} $ is the covariance operator $ \Sigma_{\widetilde{X}U}: \mathcal{H}_\mathcal{\widetilde{X}} \rightarrow \mathcal{H}_\mathcal{U} $ defined as
Following zhang11conditional, a vanishing Hilbert Schmidt norm of the partial cross-covariance operator characterizes conditional independence yielding the equivalence result:
where the Hilbert Schmidt norm of an operator $ A : \mathcal{H}_\mathcal{\widetilde{X}} \rightarrow \mathcal{H}_\mathcal{U}$ is defined as $ \left\lVertA\right\rVert_{HS} = \sqrt{\sum_{j,l \geq 1} \langle f_j,A e_l\rangle^2_{\mathcal{H}_{\mathcal U}}}$ where $ \{e_j\}_{j\geq 1} $ and $ \{f_j\}_{j\geq 1} $ are an orthonormal basis in $ \mathcal{H}_{\widetilde{\mathcal X}} $ and $ \mathcal{H}_{\mathcal U}$, respectively.
zhang11conditional build on the conditional independence characterization in (ref) and define the KCI test statistic:
where $ \widetilde{K}_{\widetilde{X} |W}$ and $\widetilde{K}_{U|W} $ are centralized kernel matrices defined as follows. The centralized kernel matrix for any variable $ Z $ is given by $ \widetilde{K}_Z = HK_Z H $ where $ K_Z $ is the uncentralized kernel matrix, i.e. a matrix whose $ (i,j) $ element is given by $ k(x_i,x_j) $ where $ k $ is the Gaussian kernel in eq. (ref), and $ H = \mathbf{I}_n-n^{-1} \mathbf{1}_n\mathbf{1}_n^\top $ where $ \mathbf{I}_n $ denotes the identity matrix of size $ n $ and $ \mathbf{1}_n $ a vector of ones of length $ n $. These centralized kernel matrices need to be adjusted to reflect the conditioning on $ W $. This is achieved using kernel ridge regression to derive a matrix $R_W$
where $ \lambda_R $ is a regularization parameter. Finally, the kernel matrix $ \widetilde{K}_{\widetilde{X}|W}$ can be expressed as
Similarly, the construction for $\widetilde{K}_{U|W}$ is analogous and yields
Following zhang11conditional, we normalize the data and choose the hyperparameters heuristically as follows. The bandwidth parameter $ \lambda $ is set at $ \lambda = 0.8 $ if sample size $ n \leq 200 $, $ \lambda = 0.3 $ if sample size $ n > 1200 $, and $ \lambda = 0.5 $ otherwise for the construction of $ \widetilde{K}_{\widetilde X} $ and $ \widetilde{K}_{U} $, and at half that value for $ \widetilde{K}_{W} $. The regularization parameter $ \lambda_R $ is set to $ 10^{-3}$. These parameters are deemed appropriate when the dimenstionality of $ W $ is small, say less than three. For higher-dimensional $ W $, zhang11conditional recommend choosing different $ \lambda $ and $ \lambda_R $ for $ \widetilde{K}_{\widetilde X} $ and $ \widetilde{K}_{U} $, which can be achieved using cross-validation.
zhang11conditional derive the asymptotic distribution and corresponding p-values of the KCI test statistic under $ H_0: $ conditional independence and show that it achieves pointwise asymptotic level strobl19. In sum, the idea of feature representations motivates the map of the distributions of $ X $ and $ U $ conditional on $ W $ into higher-dimensional spaces where linear correlations correspond to nonlinear dependencies in original space.
The theory implies an independence of errors and the covariate in the causal model and a dependence between the errors and the covariate in the anticausal model. The algorithm involves testing the independence between errors and covariate conditional on $ W $ in both causal and anticausal model. Ideally, the test would conclude with the following decisions: i) if independence can be rejected at a pre-specified significance level in one model but not in the other, one would conclude that the latter model represents the correct causal relation, ii) if independence is rejected in both models, one would conclude that the relation between $ X $ and $ Y $ is confounded, and iii) if independence cannot be rejected in either model, one would conclude that the test does not have sufficient power to decide on the causal direction. It is not possible to implement such a strategy in practice because the true errors are unobserved and the practitioner has to rely on estimated errors. Specifically, the practitioner does not have a sample of $ U := Y - h(X,W) $ in eq. (ref) at their disposal and, therefore, must rely on estimated errors $ \widehat{U} := Y - \widehat{h}(X,W) $, and the respective estimated errors of the model in eq. (ref), to investigate which model is correct. That these residuals are estimated and, in particular, that they depend on the estimated $ \widehat{h} $, poses a challenge that we discuss now.
mooijetal16,hoyer09anm propose randomly splitting the available data $\mathcal{D} = \{Y_i,X_i,W_i\}_{i=1}^{n} $ in training and test sets, denoted $\mathcal{D}^{tr} = \{Y_i,X_i,W_i\}_{i=1}^{n/2} $ and $\mathcal{D}^{te} = \{Y'_i,X'_i,W'_i\}_{i=(n/2)+1}^{n} $, respectively. $\mathcal{D}^{tr}$ is used to get an estimate $\widehat{h}$ of the true regression function $ h $. $\mathcal{D}^{te}$ is then used to get estimates $\widehat{\varepsilon}' := Y' - \widehat{h}(X',W')$ of the true errors $\varepsilon$. An error in the estimated $\widehat{h}$ induces a dependence of $\widehat{\varepsilon}'$ and $X'$ (conditional on $ W' $) even though $\varepsilon$ and $X$ are truly independent (conditional on $ W $). Consequently, conventional thresholds for the independence test tend to be too loose and would ideally incorporate the fact that $\widehat{h}$ is estimated. Specifically, for a conventional threshold of, say, $ \alpha^* = 0.05 $ the empirical rejection rate will be larger than $ \alpha^* $ in the causal model even though under $ H_0 $ we have that $ U \mathrel{\perp\mspace{-10mu}\perp} X | W $, which should lead to an empirical rejection rate roughly equal to $ \alpha^* $. To achieve an empirical size of $ \alpha^* $, one needs to use a threshold $ \alpha = \alpha^* \times \lambda_{\alpha} $ with $ 0 < \lambda_{\alpha} < 1 $. There are no theoretical results on how to choose $ \lambda_{\alpha} $ to account for the dependence of $\widehat{\varepsilon}'$ and $X'$.\footnote{Simulation studies, which are not replicated here, show that $ \lambda_{\alpha} $ depends on the type of distribution that the true error follows. Since there is no way for a practitioner to get a hold on that error distribution, it is impossible to propose rules of thumb, substantiated by simulation exercises, to indicate the level of $ \lambda_{\alpha} $ as a function of observable or estimable quantities.} However, mooijetal16 show that one can infer the correct directionality under additional assumption that the causal or the anticausal model exist. Therefore, the identifiability result in Theorem (ref), which states that either causal or anticausal model, but not both, can satisfy the independence of the error with the covariate, in combination with the existence assumption imposed in Corollary (ref), which states that either causal or anticausal model exist, allows us to infer the directionality with Algorithm (ref).
In particular, under the conditions of Corollary (ref), one can infer that the model with the lower KCI test statistic (i.e. a larger p-value of the conditional independence test) is the correct causal model, thereby circumventing the lack of theoretical guidance about an appropriate threshold. Making this assumption comes at a cost; namely, a procedure that relies on comparing two test statistics can never conclude that there is not enough information in the data to decide on the causal direction. In other words, such a procedure will never conclude that there is a lack of power to make a decision.
We investigate finite sample performance of our heteroskedasticity robust reverse causality test and its use in the causal discovery context with Monte Carlo experiments. The results are based on 500 Monte Carlo replications in each experiment and the sample size is varied with $ n \in \{250,500,1000\}$.
We simulate data for the following model:
and
where $ \phi $ denotes the standard normal probability density function, $\text{sgn}(\cdot)$ is the sign function, and $q,\rho, c_{\rho,q}$ are constants which vary in the experiments below. The dependent variable is generated by
where $ j \in \{1,2\} $ and the functions $\kappa_j$ are given by
Here, $ \tau $ controls the degree of nonlinearity between $ Y $ and $ X $. The linear case corresponds to $ \tau = 0 $. The parameter $\rho$ captures degree of the heteroskedasticity w.r.t. the control variable $W$. That is, $\mathrm{Var}(U|W=w)= (1+\phi(w))^\rho$ and hence, $ \rho = 0 $ corresponds to the homoskedastic case. We simulate data for $ \tau \in \{0, 0.25,0.5,0.75,1\} $ and $ \rho \in \{0,1\} $. We rely on the R package CondIndTests CondIndTests for the implementation of the HSIC independence test.
To explore the robustness of our results with respect to the distribution of the error $ {U} $, we run the simulation with errors drawn from sub- and super-Gaussian distributions. We choose the constant $c_{\rho,q}$ (via numerical approximation) such that the variance of $U$ is normalized to one under each choice of $\rho$ and $q$. We estimate both causal and reverse models with a Generalized Additive Model using smoothing splines.
We implement the test as described in Algorithm (ref) for the nominal level $\alpha=0.05$. Figure (ref) reports the empirical rejection probabilities of testing conditional independence of covariates and estimated residuals under the reverse model (ref). Overall, the empirical rejection rates are close to the nominal level for different choices of $\rho$, $q$, and $\tau$. This is remarkable as the testing problem is complex and builds on nonparametric procedures. As such, we cannot expect exact control over the significance level; see also the discussion in Section (ref). The test shows some degree of oversizing under super-Gaussian error distributions, indicating the complexity of the testing problem.
We illustrate the empirical power of the reverse causality test in Figure (ref). The power of the test, i.e., the probability of rejecting the independence of estimated residual and candidate cause in model (ref), increases sharply as the degree of nonlinearity $\tau$ increases. While our identifiability results rely on a nonlinear $h$, see Assumption (ref), it can be seen that non-Gaussianity can be a source of power similar to nonlinearity: with sufficiently many samples (rightmost column), the empirical rejection rates lie above 0.4 even in the linear case ($\tau=0$) as long as $q\neq1$. This is a well-known result in the causal discovery community, see e.g. shimizu06lingam.
In Figure (ref) we report empirical probabilities of correct classification of the causal direction. When the relationship between $ X $ and $ Y $ is linear, i.e. $ \tau = 0 $, and the error Gaussian, i.e. $ q=1 $, the algorithm performs at about chance level. This is consistent with the theory since the causal direction is not identifiable in the linear case.\footnote{Note that shimizu06lingam show that the direction is identifiable in the linear case with a non-Gaussian, homoskedastic error distribution, i.e. when $ q \neq 1 $ and $ \rho = 0 $. Though our simulation results suggest that the result also holds with heteroskedastic errors w.r.t. $ W $, we do not extended it formally.} As soon as the relation between cause and effect becomes nonlinear, i.e. $ \tau \neq 0 $, the accuracy of the reverse causal discovery algorithm increases. For instance, when $ n = 500 $ and the relation between cause and $ \tau = 1 $ the algorithm arrives at the correct conclusion in more than 95% of the Monte Carlo runs, regardless of the level of heteroskedasticity (parameterized by $ \rho $) and shape of the error distribution (parameterized by $ q $). We also see that our procedure has power to detect the correct causal directions in cases where the nonlinearity is less pronounced, i.e., $0<\tau<1$. Moreover, the performance of the algorithm is robust to changes in the specification of the functional form, which can be seen by comparing the $ \kappa_1 $ and $ \kappa_2 $ rows in Figure (ref). For a given $ \tau $, the test has more power for $\kappa_1$ because the nonlinear relation between $ X $ and $ Y $ is more pronounced than for $ \kappa_2 $. Overall we see that the accuracy of causal discovery increases with sample size, where this change is stronger when the regression function is given by $\kappa_2$. We show that the results remain robust to different error variances in Appendix (ref).Specifically, we consider a low noise regime where $c_{\rho,q}$ is chosen such that $\mathrm{Var}(U)=0.8$ and a high noise regime where $\mathrm{Var}(U)=1.2$.
In sum, two observations are worth stressing. First, the results show that the algorithm has surprisingly high power when cause and effect are related nonlinearly. Second, the performance of the algorithm does not suffer from heteroskedastic errors w.r.t. $ W $.
We use data from the 2013 Survey of Income and Expenditure (“Einkommens- und Verbrauchsstichprobe”, EVS), which is a voluntary survey of roughly 60,000 households in Germany, to test the proposed algorithm. We consider the following variables: income, expenditure, highest educational attainment of the main earner, highest professional training of the main earner, and age group of the main earner. We analyze the causal direction between income and work experience, which we proxy by age group.
Hump-shaped income profiles over the life-cycle are well-documented in labor economics heckman06earnings. It is interesting to test the algorithm for a cause-effect pair where the causal direction is a priori clear. Since work experience mechanically increases over the life-cycle, it can be credibly assumed not to be caused by income changes. Therefore, we analyze the directionality between income (Inc) and age where age can be interpreted as proxy for work experience (Exp). We posit the correct causal model to be
where experience Exp causes income Inc. Vice versa, the anticausal model in which income causes work experience is given as
where in each model $ W $ contains all remaining covariates as control (expenditure, highest educational attainment of the main earner, highest professional training of the main earner).
We aim to alleviate the problem that we are likely to omit many crucial confounding variables by splitting the data in $ n_q $ quantiles of the expenditure distribution. At least part of the omitted confounding factors can be assumed to be fixed within given quantiles as they collect individuals with roughly similar life-styles. This argument applies even more strongly as the expenditure distribution is split into a larger number of quantiles. On the other hand, the larger $ n_q $ the smaller the number of observations within each quantile and the lower the power of the test to identify the correct causal direction. Therefore, we show results for a set of $ n_q = \{4, \dots, 20\} $ quantiles.\footnote{The KCI test, which forms an important part of the algorithm, requires the inversion of $ n \times n $ matrices where $ n $ is the number of observations. Constraints on local computing power preclude running the test on the whole sample with roughly 60,000 observations or with $ n_q = \{1, 2, 3\} $.} For each number of quantiles $ n_q $, we run the causal discovery algorithm in each of these $ n_q $ quantiles and plot the share of quantiles in which the algorithm prefers either model (note that the $ x $-axis in Figure (ref) refers to the number of quantiles the income distribution is split in, not the quantiles as such). For example, the bar above $ n_q = 5 $ in Figure (ref) denotes that in 4 of the 5 quantiles, i.e. 80%, the algorithm concludes that experience causes income. Regardless of $ n_q $, the test always favours the model where work experience causes income in at least 50% of quantiles. The algorithm favours the causal model in more than 75% of the quantiles for most $ n_q $.
In sum, this application documents that our algorithm gives economically meaningful results in empirical applications.
Endogeneity is a common threat to causal identification in econometric models. Reverse causality is one source of such endogeneity. We build on work by hoyer09anm,mooijetal16 who have shown that the causal direction between two variables $ X $ and $ Y $ is identifiable in models with additively separable error terms and nonlinear function forms. We extend their results by allowing for additional control covariates $ W $ and heteroskedasticity w.r.t. them and, thus, provide a heteroskedasticity-robust method to test for reverse causality. In addition, we show how this test can be extended to a bivariate causal discovery algorithm by comparing the test statistics of residual and purported cause of two candidate models. We extend known results on causal identification and causal discovery to settings with heteroskedasticity with respect to additional control covariates.
An empirical application underscores the feasibility of the proposed algorithm. We analyze the causal link between income and work experience, as proxied by age, and show that our procedure provides evidence that the true causal direction is from work experience to income. Though this result is not substantively surprising because income cannot causally influence work experience, it is encouraging that our algorithm can distinguish between the causal directions without resorting to instruments or other sources of exogenous variation. This underscores the value of our proposed methodology to shed light on the causal structure of economic phenomena without resorting to exogenous variation.
\printbibliography[heading=none]