Testing of Reverse Causality Using Semi-Supervised Machine Learning

Abstract

Two potential obstacles stand between the observation of a statistical correlation and the design (and deployment) of an effective intervention, omitted variable bias and reverse causality. Whereas the former has received ample attention, comparably scant focus has been devoted to the latter in the methodological literature. Many existing methods for reverse causality testing commence by postulating a structural model that may suffer from widely recognized issues such as the difficulty of properly setting temporal lags, which are critical to model validity. In this article, we draw upon advances in machine learning, specifically the recently established link between causal direction and the effectiveness of semi-supervised learning algorithms, to develop a novel method for reverse causality testing that circumvents many of the assumptions required by traditional methods. Mathematical analysis and simulation studies were carried out to demonstrate the effectiveness of our method. We also performed tests over a real-world dataset to show how our method may be used to identify causal relationships in practice.

1 Introduction

A fundamental purpose of research in psychology—and many other disciplines in social sciences for that matter—is to understand causal relationships between variables. In particular, it is both theoretically and practically important to distinguish between the mere observation of associations (between variables) and cases where causality can be inferred. When controlled randomized experiments are impractical, how to properly do so has garnered considerable attention in multiple disciplines (e.g., Pearl, 1998) including psychology (Rogosa, 1980). The practical significance of this inquiry becomes apparent considering the pivotal role that causal inference plays in formulating effective intervention strategies. For example, research has consistently shown that individuals who occupy central positions in their social networks (i.e., centrality) tend to receive more favorable assessment of charisma from their connections (Balkundi et al., 2011). Yet, the mere identification of this correlation ¹ does not permit the conclusion that enhancing an individual’s social-network centrality would be an effective intervention for boosting their charisma. As a case in point, whereas some scholars view charisma as an outcome of social-network patterns (Pastor et al., 2002), others contend that charisma is what attracts followers and allows an individual to occupy a central position in the first place (Shils, 1965), suggesting that the aforementioned intervention might be less effective.

Two potential obstacles stand between (a) the observation of a robust correlation between two constructs X and Y, like the centrality-charisma correlation, and (b) the effectiveness of deploying an intervention on X to change Y, like an attempt to boost centrality in order to increase charisma (Bollen, 1989, p. 41). The first is omitted variable bias (Mauro, 1990), meaning that an unobserved confounder induces the co-variation of both X and Y. For example, individuals high on extraversion may excel on both social-network centrality and charisma. The other is reverse causality (Leszczensky & Wolbring, 2022), meaning that the causal influence flows in the reverse direction (i.e., $Y \to X$ ) from what is required by the intended intervention—e.g., charisma influences centrality but not the other way around. There is a substantial body of methodological research on assessing the magnitude of omitted variable bias in a relationship (e.g., Busenbark et al., 2022; Cinelli & Hazlett, 2020; Harring et al., 2017; Mauro, 1990). Yet comparatively little attention has been directed toward testing for reverse causality (Leszczensky & Wolbring, 2022). Admittedly, the issue may not be applicable in situations where a clear causal direction is obvious (e.g., cancer $\to $ age is implausible). Nonetheless, for psychological constructs, concerns on reverse causality are prevalent, especially when theoretical arguments could be made for both directions. The focus of this work is to develop a method for testing reverse causality in observational data, with an emphasis on panel data.

To empirically investigate the issue of reverse causality in observational data, some existing methods, like cross-lagged panel models (CLPM; Hamaker et al., 2015), commence by presuming a reciprocal relationship, before postulating a structural model based on the assumption and empirically estimating the model to assess the magnitude in both directions while considering auto-regressive effects. Some other methods infer the direction of causality by exploiting certain distributional and/or functional features assumed of the underlying data-generating process, such as non-normality in a linear model (e.g., Shimizu et al., 2006; von Eye & DeShon, 2012), statistical independence between input variables and additive noise (e.g., Hoyer et al., 2008; Rosenström et al., 2023), etc. These methods pose two unsolved issues. First, from a theoretical perspective, given that many psychological theories (and theories in other fields) stipulate unidirectional (instead of reciprocal) relationships, there exists substantial theoretical interest in testing the direction of causality rather than assuming it away with the adoption of a reciprocal model like CLPM. For example, the structural advantage theory of social networks (Brass, 1984; Burt, 1992) posits a unidirectional effect of one’s social network characteristics (e.g., centrality) on individual behaviors (e.g., charismatic leadership). Given the robust correlation reported in the literature for centrality and charisma (e.g., Balkundi et al., 2011), a well-powered empirical study that shows the absence of a reverse effect (i.e., charisma $\not \to $ centrality) would provide strong evidential value for supporting the structural advantage theory. Second, from a methodological perspective, formulating an appropriate model—or accurately specifying its distributional/functional features (e.g., whether it is linear)—can be particularly challenging when dealing with panel data (Hamaker, 2024; Lucas, 2023). For example, the validity of many existing panel models (or even the notion of Granger causality itself) is known to break down when the sampling frequency is improperly specified (Shojaie & Fox, 2022; Vaisey & Miles, 2017), when there are confounders that are unaccounted for (e.g., Hamaker et al., 2015), etc.

A promising avenue for addressing these issues of existing methods arises from causal learning (Peters et al., 2017), a branch of machine learning that injects causal inference into the design of learning algorithms. Central to this approach is a fundamental question: if we were to train a machine learning model that predicts Y from X (based on a limited number of training data points $\langle X, Y \rangle $ ), could knowledge about the probability distribution of X (i.e., $P(X)$ ) help us improve the predictive accuracy of the trained model? Schölkopf et al. (2012) show that the answer is positive if and only if the causal direction flows from Y to X. This suggests that, by demonstrating a machine learning model’s ability to capitalize on $P(X)$ for enhancing predictive accuracy toward Y, we would identify the existence of reverse causality. Even more importantly, this proposition would hold irrespective of the functional form of the X–Y relationship.

Whereas Schölkopf et al. (2012) establish the dependence between causal direction and the predictive accuracy of semi-supervised learning algorithms, they leverage this finding to understand why semi-supervised learning works over some datasets but not others, instead of developing a concrete method for testing reverse causality. In the current research, we address this gap by developing a novel method for testing reverse causality over panel data. The contribution of our work is two-fold. First, our method for reverse causality testing allows researchers to rigorously test their causal theories by addressing two issues that researchers face. On the one hand, our method allows researchers to directly test unidirectional relationships, instead of assuming these relationships are reciprocal in nature. On the other hand, our method does not require the specification of distributional or functional features such as the sampling frequency or the shape of the proposed relationship. By demonstrating the effectiveness of our method using simulation studies and a case study, we show its applicability for psychology researchers in theory building and testing. Second, our work pioneers the integration of advancements in causal learning into the methodological arsenal of psychology for causal inference. We present conceptual arguments and mathematical formalization that link reverse causality testing with the predictive accuracy of semi-supervised learning (Van Engelen & Hoos, 2020). In doing so, we enrich the understanding of how machine learning could contribute to the psychological methods literature (e.g., Sterner et al., 2023; Wilcox et al., 2023; Zimmer & Debelak, 2023) and open up future avenues of inquiry at the intersection of psychological research and computer science.

2 Literature review

In this section, we briefly review the existing literature on reverse causality testing and the machine learning method we propose to use (i.e., semi-supervised learning).

2.1 Reverse causality testing

Compared with the rich and growing literature on omitted variable bias—in psychology (Harring et al., 2017), sociology (Halaby, 2004) and economics (Wüthrich & Zhu, 2023)—researchers across disciplines have made relatively limited progress on the testing of reverse causality over panel data. As summarized by Leszczensky & Wolbring (2022), a common approach is to specify causal direction in a panel model by applying temporal lags on variables that represent the “cause” after partialing out auto-regressive effects. For example, a causal direction of $X \to Y$ would be reflected by setting a lagged (e.g., previous-wave) value of X and the contemporary value of Y as independent and dependent variables, respectively, in the panel model, suggesting that X has a causal, lagged, effect on Y. A variety of panel models follow this idea (Orth et al., 2021), such as lagged first-difference (LFD) models (Vaisey & Miles, 2017), CLPM (Hamaker et al., 2015), etc.

As these existing panel models rely on temporal lags to identify the causal direction, how to specify the amount of this temporal lag becomes a prominent question. Theoretically deriving the “correct” temporal lag is obviously even more challenging than discerning the causal direction, suggesting the need for methods to be robust to misspecified temporal lags. Unfortunately, Vaisey & Miles (2017) show that panel models such as LFDs can be highly sensitive to the misspecification of temporal lags. Leszczensky & Wolbring (2022, Figure 2) also demonstrate that, when a contemporaneous effect is mis-characterized as lagged in a CLPM, the model can produce highly biased estimates. More fundamentally, the very notion of Granger causality, which underpins all panel models, is known to break down with a mis-specified temporal lag (Shojaie & Fox, 2022).

Beyond panel models, other existing methods for reverse causality testing rely on certain distributional/functional features assumed of the underlying data-generating process. Some—like Direction Dependence Analysis (DDA; Li & Wiedermann, 2020; Pornprasertmanit & Little, 2012; von Eye & DeShon, 2012; Wiedermann & Li, 2018) and Linear Non-Gaussian Acyclic Models (LiNGAM; Shimizu et al., 2006)—exploit the non-normality of data distributions in a linear model to infer causal direction. In the case of DDA, for example, the causal direction may be inferred by comparing the degrees of departure ² from normality across different variables. For longitudinal data, extensions of these methods (e.g., Bauer et al., 2016; Geiger et al., 2015; Hyvärinen et al., 2010) leverage the non-normality of noise to identify causal direction in multivariate time series. Additionally, methods like Rosenström et al.’s (2023) directional analysis approach, additive noise models (Hoyer et al., 2008; Peters et al., 2014) and (more generally) post-nonlinear models (Zhang & Hyvärinen, 2009) infer causal direction by assuming statistical independence between the input variables and an additive noise component within linear or nonlinear models. Yet, like panel models’ reliance on properly specified temporal lags, these methods also hinge on accurately defining certain distributional or functional features of the relationship between variables—knowledge that may not be available a priori.

These issues of existing methods raise an important question: Can we test for reverse causality without constraining the functional form of the data-generating process? Doing so would not only circumvent the complexities of specifying the temporal lag, but also relax the linearity or additive noise assumptions that permeate existing methods. We seek to answer this question in the current work by integrating recent advances in machine learning and the broader causal inference literature.

2.2 Semi-supervised learning

The objective of semi-supervised learning is to approximate the function f that links independent variables X to a dependent variable Y such that $Y = f(X)$ . This is done by learning from two datasets containing i.i.d. samples. The first dataset, typically known as labeled set, provides $\langle x_i, y_i\rangle $ (i.e., paired X–Y values) for $n_1$ data points. The second dataset provides only the values of X, but not Y, for other $n_2$ data points, and is therefore known as the unlabeled set. In machine learning, it is generally assumed that $n_1$ is much smaller than $n_2$ , due mainly to the high cost of acquiring Y in practice. For example, a task for which semi-supervised learning has shown remarkable success is image classification (Xie et al., 2020). With this task, X is an image and Y is its category (e.g., landscape, portrait). It is virtually cost-less to collect millions of images from the web, but considerably more expensive to hire human workers to properly label the collected images. In this case, we may opt to manually label only a few observations of X, leading to the assumption that $n_1 \ll n_2$ .

The limited size of the labeled set means that, if we were to launch a canonical supervised learning algorithm (e.g., OLS or logistic regression), which can only learn from the labeled set, we would not be able to obtain an accurate prediction of Y. The uniqueness of semi-supervised learning lies in its ability to leverage the $n_2$ unlabeled, X-only, data points to significantly improve the predictive accuracy of the learned $\hat {f}$ . To this end, numerous methods have been proposed for semi-supervised learning (Van Engelen & Hoos, 2020). A famous example, which we will further elaborate later in the article, is self-training (e.g., Sohn et al., 2020), which can be readily integrated with many supervised learning algorithms such as logistic regression. With self-training, we start by running a supervised learning algorithm (e.g., logistic regression) over the $n_1$ labeled data points to generate an approximation of f, denoted by $\hat {f}_1$ , which predicts Y based on an input X. Then, we apply $\hat {f}_1$ over each of the $n_2$ unlabeled data points to predict its pseudo-label (i.e., an estimate of Y). Note that when a machine learning model is used for prediction, it may generate not only a point-estimate (e.g., binary prediction in logistic regression) but also a confidence level associated with the estimate (e.g., log-odds in logistic regression). Leveraging this, we select from the $n_2$ pseudo-labels those with confidence above a pre-determined threshold, add their X values paired with pseudo-labels (i.e.,) $\langle x_i, \hat {f}_1(x_i) \rangle $ to the labeled set, before running the supervised learning algorithm again to update our approximation of f. This process can continue iteratively until no more pseudo-labels can be added.

Whereas the efficacy of semi-supervised learning has long been established (Van Engelen & Hoos, 2020), there has been limited research probing why unlabeled data, which lack any information about Y, can bolster our understanding of the X–Y relationship. To this end, Buja et al. (2019) establish that the distribution of a variable X may convey rich information about the X–Y relationship when the relationship is not strictly linear. For the high-dimensional case where X comprises multiple variables, Niyogi (2013) attributes the success of semi-supervised learning to the concept of manifold regularization (Belkin et al., 2006). At its core, Niyogi (2013) posits that almost all machine learning algorithms predicate their predictions on a smoothness assumption: if two data points are similar in X, then they also tend to be analogous in Y. Unfortunately, defining “similar” in the context of multi-dimensional X is challenging—and requires structural insights into the multivariate distribution of X—because common similarity measures like Euclidean distance are known to become meaningless in high-dimensional spaces (Aggarwal et al., 2001). By providing a more refined depiction of the multi-dimensional distribution of X, unlabeled data allow us to more accurately estimate the similarity of two data points in X, thereby enhancing our predictive precision for Y (Niyogi, 2013).

While this line of research sheds light on the mechanics of semi-supervised learning, it does not address why semi-supervised learning achieves high predictive accuracy with certain datasets but performs poorly with others (Schölkopf et al., 2012). The missing piece—the key to understanding the conditions under which semi-supervised learning succeeds—lies in the causal effects driving the generation of the observed data, as we will elaborate in the next section.

3 Linking causality with semi-supervised learning

In this section, we draw upon insights from causal learning (Peters et al., 2017)—specifically Janzing & Schölkopf’s (2015) mathematical characterisation of Schölkopf et al.’s (2012) seminal finding on the working condition for semi-supervised learning—to elucidate how the predictive accuracy of semi-supervised learning over panel data may reveal the existence of reverse causality in the underlying data-generating process. Our explanation unfolds in three steps. First, we offer conceptual arguments through an illustrative example. Following this, we delve into mathematical formalization, analyzing a special case where the data-generating process features linear, cross-lagged, effects. Finally, we expand our discussions to justify the reasoning behind a test that can be applied to any arbitrary data-generating process. The specific computational algorithms for semi-supervised learning will be described in the next section.

3.1 Conceptual illustration

In what follows, we offer an intuitive explanation for the key insight of Schölkopf et al. (2012) through two observations: (1) $P(X)$ is useless for predicting Y if $X \to Y$ , and (2) if $Y \to X$ , $P(X)$ (when combined with a small labeled set) may lead to an accurate prediction model toward Y. The first is obvious: When $X \to Y$ fully characterizes the X–Y relationship, we can represent their relationship as a stochastic function f such as $Y = f(X)$ . Changing f clearly has no impact on $P(X)$ . Consequently, knowledge of $P(X)$ reveals no information about f, making it useless for the prediction of Y (Janzing & Schölkopf, 2010).

For the second observation, consider a simple example of $Y \to X$ depicted in Figure 1, where $Y \in \{0, 1\}$ is binary and $X = \alpha Y + \epsilon $ with $\epsilon \sim N(0,\sigma ^2)$ . In this case, altering the causal mechanism (i.e., $\alpha $ ) obviously modifies $P(X)$ . From the perspective of machine learning, this means that $P(X)$ now encapsulates certain cues regarding $\alpha $ , offering the potential for building an accurate prediction model toward Y. As can be seen from Figure 1, for this specific example, knowledge of $P(X)$ alone permits an exact inference ³ of $\alpha $ , which is all that is required to build a Bayes-optimal predictior of Y.

Figure 1 Illustrative example for $Y \to X$ .Note: Both panels depict the probability density function of $P(X)$ when $X = \alpha Y + \epsilon $ , where Y follows Bernoulli distribution with $p = 0.5$ and $\epsilon \sim N(0,1)$ . Note that, in either case, $P(X)$ follows a Gaussian mixture distribution with two equal-weight components, which are illustrated in red dashed lines. The mean difference between the two Gaussian components is always equal to $\alpha $ , suggesting that the functional relationship between X and Y can be precisely inferred from $P(X)$ .

More generally, Schölkopf et al. (2012) contend that $P(X)$ is “independent” of f if reverse causality does not exist (i.e., $Y \not \to X$ ). However, this independence may no longer hold when $Y \to X$ . This notion of “independence” can be formalized mathematically in at least three distinct ways: (1) through Kolmogorov complexity (Janzing & Schölkopf, 2010), (2) by examining the uncorrelatedness between p and the derivative of f (Janzing & Schölkopf, 2015), and (3) through the uncorrelatedness between p and the logarithm of the derivative of f (Daniusis et al., 2010). Whereas the mathematical formulation of Schölkopf et al.’s (2012) insight necessarily varies with this underlying notion of independence—e.g., see Janzing and Schölkopf (2015, Lemma 1) for one variation—the two observations discussed before hold true regardless of this specific mathematical formulation.

3.2 Mathematical analysis of cross-lagged effects

To mathematically illustrate the link between causality testing and semi-supervised learning, we use the CLPM (Hamaker et al., 2015) as an example of the underlying data-generating process. This analysis aims to emphasize two key observations. First, in the absence of reverse causality (i.e., $Y \not \to X$ ), no information about the coefficient representing the $X \to Y$ relationship can be inferred from a dataset containing only X. Second, when reverse causality exists (i.e., $Y \to X$ ), a standard least squares estimate for the coefficient representing $Y \to X$ can be obtained even without Y in the dataset.

In its simplest form, CLPM assumes a model of structural equations:

(1) $$ \begin{align} (x_{it} - \kappa_i) &= \alpha (x_{i,t-1} - \kappa_i) + \beta (y_{i, t-1} - \omega_i) + u_{it},\end{align} $$

(2) $$ \begin{align} (y_{it} - \omega_i) &= \delta (y_{i,t-1} - \omega_i) + \gamma (x_{i, t-1} - \kappa_i) + v_{it}, \end{align} $$

where $x_{it}$ and $y_{it}$ are the values of X and Y at time t, respectively; $\kappa _i$ and $\omega _i$ represent the sample-specific random intercepts for X and Y, respectively; and $u_{it}$ and $v_{it}$ are i.i.d. random impulses. With this model, the cross-lagged parameter $\beta $ captures the (within-person) reverse causal effect $Y \to X$ .

As the sample-specific random intercepts $\kappa _i$ and $\omega _i$ are extraneous to our ensuing analysis, we assume that $\kappa _i = \omega _i = 0$ for all $i \in [1, n]$ (where n is the sample size). For the ease of understanding, we also assume that all coefficients $\alpha $ , $\beta $ , $\delta $ , and $\gamma $ remain unchanged between time $t-1$ and $t+1$ . These simplifying assumptions yield the following structural equations:

(3) $$ \begin{align} x_{it} &= \alpha x_{i,t-1} + \beta y_{i, t-1} + u_{it}, \end{align} $$

(4) $$ \begin{align} y_{it} &= \delta y_{i,t-1} + \gamma x_{i, t-1} + v_{it}. \end{align} $$

The crux of connecting semi-supervised learning with causality testing hinges on the following inquiry: Can any information about the regression coefficients in Equation 4 (i.e., $\delta $ and $\gamma $ for predicting Y) be inferred from observations of X at time $t-1$ , t, and $t+1$ (i.e., $x_{1,t-1}, \ldots , x_{n,t-1}$ ; $x_{1t}, \ldots , x_{nt}$ ; $x_{1,t+1}, \ldots , x_{n,t+1}$ ), but no observation of Y, assuming the other coefficients (i.e., $\alpha $ and $\beta $ ) as given?

When $\beta = 0$ , this is obviously impossible because Equation 3 would only contain an autoregressive term $\alpha x_{i,t-1}$ and the remainder $u_{it}$ , meaning that we would have no information about Y from any of the inputs X, $\alpha $ , and $\beta $ .

When $\beta \neq 0$ , however, knowledge of $\alpha $ and $\beta $ would allow us to make a probabilistic inference of $y_{i,t-1}$ from Equation 3 based on $x_{it}$ and $x_{i,t-1}$ and, similarly, $y_{it}$ based on $x_{i,t+1}$ and $x_{it}$ . Substituting these inferred values of $y_{i,t-1}$ and $y_{it}$ into Equation 4 would then let us make a probabilistic inference of the regression coefficients $\delta $ and $\gamma $ . To elucidate the reasoning behind this procedure, we rewrite Equation 3 as

(5) $$ \begin{align}\kern-3.5pt x_{i,t+1} &= \alpha x_{it} + \beta(\delta y_{i,t-1} + \gamma x_{i, t-1} + v_{it}) + u_{i,t+1}, \end{align} $$

(6) $$ \begin{align} &= \alpha x_{it} + \delta (\beta y_{i,t-1}) + \gamma \beta x_{i,t-1} + \beta v_{it} + u_{i,t+1}, \end{align} $$

(7) $$ \begin{align} &= \alpha x_{it} + \delta (x_{it} - \alpha x_{i,t-1} - u_{it}) + \gamma \beta x_{i,t-1} + \beta v_{it} + u_{i,t+1}, \end{align} $$

which can be further simplified to

(8) $$ \begin{align} (x_{i,t+1}-\alpha x_{i,t}) &= \delta (x_{it}-\alpha x_{i,t-1}) + \gamma \beta x_{i, t-1} + \epsilon_i, \end{align} $$

where $\epsilon _i = \beta v_{it} + u_{i,t+1} - \delta u_{it}$ is a remainder term that aggregates various random impulses and satisfies $E(\epsilon _i) = E(\epsilon _i x_{it}) = E(\epsilon _i x_{i,t-1}) = 0$ . Equation 8 can be used to obtain a standard least squares estimate for $\delta $ and $\gamma $ , i.e.,

(9) $$ \begin{align} \begin{bmatrix} \hat{\delta}\\ \hat{\gamma} \end{bmatrix} = ([\mathbf{x}_t - \alpha \mathbf{x}_{t-1}, \beta \mathbf{x}_{t-1}]^\top [\mathbf{x}_t - \alpha \mathbf{x}_{t-1}, \beta \mathbf{x}_{t-1}])^{-1} [\mathbf{x}_t - \alpha \mathbf{x}_{t-1}, \beta \mathbf{x}_{t-1}]^\top (\mathbf{x}_{t+1} - \alpha \mathbf{x}_t), \end{align} $$

where $\mathbf {x}_j$ represents the vector of X at time j. The variances of $\hat {\delta }$ and $\hat {\gamma }$ are

(10) $$ \begin{align} \mathrm{Var}(\hat{\delta}) &= \sigma^2 ((\mathbf{x}_t - \alpha \mathbf{x}_{t-1})^\top(\mathbf{x}_t - \alpha \mathbf{x}_{t-1}) - ((\mathbf{x}_t - \alpha \mathbf{x}_{t-1})^\top \mathbf{x}_{t-1})^2 (\mathbf{x}^\top_{t-1} \mathbf{x}_{t-1})^{-1})^{-1}, \end{align} $$

(11) $$ \begin{align} \mathrm{Var}(\hat{\gamma}) &= \sigma^2 (\beta^2 \mathbf{x}_{t-1}^\top \mathbf{x}_{t-1} - ((\mathbf{x}_t - \alpha \mathbf{x}_{t-1})^\top \mathbf{x}_{t-1})^2 ((\mathbf{x}_t - \alpha \mathbf{x}_{t-1})^\top(\mathbf{x}_t - \alpha \mathbf{x}_{t-1}))^{-1})^{-1}, \end{align} $$

where $\sigma ^2$ is the variance of $\epsilon _i$ . Note from Equation 11 that the variance of $\hat {\gamma }$ tends to infinity when $\beta $ approaches 0. This substantiates our prior assertion that estimating the coefficient for $X \to Y$ (i.e., $\gamma $ ) from X is only feasible when the coefficient for $Y \to X$ (i.e., $\beta $ ) is significant.

Given that $\delta $ and $\gamma $ can be inferred (at least asymptotically) using just three inputs— $\alpha $ , $\beta $ , and X—without any information on Y, a natural question is: which among the three inputs carries information about $\delta $ and $\gamma $ ? We can rule out $\alpha $ and $\beta $ because they are free parameters (independent of $\delta $ and $\gamma $ ) per CLPM. This leaves the only possibility to be that the distribution of X carries certain information about the $X \to Y$ relationship (i.e., $\delta $ and $\gamma $ ). Importantly, this information carriage is valid if and only if reverse causality $Y \to X$ is present (i.e., $\beta \neq 0$ ).

At this juncture, the connection between semi-supervised learning and reverse causality testing becomes evident. Recall that a semi-supervised learning algorithm aims to learn the $X \to Y$ relationship $f(\cdot )$ —i.e., $\gamma $ in the case of CLPM—using a small labeled set of $\langle X, Y \rangle $ pairs and a large unlabeled set consisting solely of X. Clearly, the only information revealed by the unlabeled set is the distribution of X. When reverse causality does not exist (i.e., $\beta = 0$ ), the distribution of X does not carry any information about $\gamma $ . In this case, the unlabeled set is useless for improving the prediction of Y, making the deployment of a semi-supervised learning algorithm futile. However, when reverse causality exists (i.e., $\beta \neq 0$ ), the distribution of X—and therefore the unlabeled set—does carry information about $\gamma $ (as shown in Equation 9), making semi-supervised learning potentially fruitful.

For a more specific example, consider the aforementioned self-training method for semi-supervised learning. When the X–Y relationship is unidirectional flowing from X to Y (i.e., $\beta = 0$ ), expanding the labeled set with pseudo-labels will not lead to a more accurate prediction of Y for the simple reason that even an infinitely large unlabeled set, which perfectly reveals $P(X)$ , still contains no information about the $X \to Y$ relationship. Put simply, we can test reverse causality by comparing the predictive accuracy pre and post self-training, with an increase in predictive accuracy suggesting the existence of reverse causality.

3.3 Generalized test

Whereas we developed the connection between semi-supervised learning and reverse causality testing through the structural model of CLPM, the connection indeed persists regardless of the underlying data-generating process. To understand why, consider an extension of Equations 3 and 4 to data-generating processes beyond CLPM. When the X–Y relationship is unidirectional flowing from X to Y, we have

(12) $$ \begin{align} x_{it} &= f_t(x_{i,t-1},u_{it}), \end{align} $$

(13) $$ \begin{align} y_{it} &= g_t(x_{i,t-1},y_{i,t-1},v_{it}), \end{align} $$

where $f_t$ and $g_t$ can be any arbitrary stochastic function.

Two key insights emerge from the equations. First, the $X \to Y$ relationship is wholly captured by $g_t$ . Second, alterations in $g_t$ has no influence on the value (and distribution) of X. Combining these two insights, it is evident that the distribution of X does not carry any information about the $X \to Y$ relationship unless there exists a reciprocal $Y \to X$ relationship. This underscores that the nexus between semi-supervised learning and reverse causality testing remains intact, independent of assumptions about the data-generating process.

4 Reverse causality testing using semi-supervised learning

In this section, we develop our novel method for reverse causality testing. As established in the last section, a semi-supervised learning algorithm can effectively leverage the large unlabeled (i.e., X-only) data to improve predictive accuracy for Y only when reverse causality $Y \to X$ is at play. Building upon this insight, our method revolves around assessing the predictive accuracy of semi-supervised learning. In the passages that follow, we outline our methodological design in two steps. First, we describe the input and output of semi-supervised learning, explaining in detail how the output of semi-supervised learning is used for reverse causality testing. Then, we delve into the algorithmic design of semi-supervised learning.

4.1 Input and output of semi-supervised learning

The objective of our method is to detect the presence of reverse causality, i.e., $Y \to X$ , given a longitudinal panel dataset $D = \{\langle x_{i1}, y_{i1}, \ldots , x_{im}, y_{im} \rangle | i \in [1, n]\}$ , where n is the sample size and m is the number of waves. Our method requires at least three waves (i.e., $m \geq 3$ ) for identification. It imposes no assumptions about the lag, provided that the lag is not so extensive that $y_{i1}$ exerts no causal effect on $x_{im}$ , as this would render reverse causality unidentifiable within the given dataset. Our method also makes no assumption about the functional form of the relationships between X and Y. The sole assumption underlying our method is the working condition of causal learning, which states that $P(X)$ is independent ⁴ of f if $X \to Y$ fully describes the relationship between X and Y (Schölkopf et al., 2012). Also note that, if the purpose is to test the existence of $X \to Y$ instead, the only revision required is to swap X and Y in the input data.

4.1.1 Input

Recall from the literature review that a semi-supervised learning algorithm takes as input two datasets, a small labeled set $D_{\mathrm {lb}}$ , which consists of $n_1$ data points with paired X–Y values, and a large unlabeled set $D_{\mathrm {ul}}$ , which consists of $n_2$ data points ( $n_2 \gg n_1$ ) with X values only. To properly specify $D_{\mathrm {lb}}$ and $D_{\mathrm {ul}}$ based on the input data D, there are three issues to be addressed.

First, since semi-supervised learning algorithms generally require the prediction target Y to be a scalar (i.e., a single variable), we need to select (the Y value of) one wave as the prediction target for semi-supervised learning. Mathematically, any wave except the last one would work because, as discussed earlier, if $Y \to X$ exists with a lagged effect, then the distribution of X in the tth wave carries information about Y in all previous waves (i.e., 1 to $t-1$ ). In other words, semi-supervised learning could effectively boost the predictive accuracy toward Y for all but the last wave. Since we are interested in minimizing the number of waves required for identification, a proper choice is to select Y in the first wave, i.e., $\{y_{i1}\}$ , as the prediction target because otherwise data in the first wave would become useless for identification (as Y in a latter wave cannot have a causal effect on X in the first wave).

Second, we need to determine the variable composition of X. When $Y \to X$ exists, the distribution of X from the second wave onward is informative for predicting the value of Y from the first wave. We therefore include $\langle x_{i2}, \ldots , x_{im} \rangle $ as the predictor vector X.

Third, we also need to determine $n_1$ and $n_2$ , the sample sizes for $D_{\mathrm {lb}}$ and $D_{\mathrm {ul}}$ , respectively. Recall that our purpose is to test whether semi-supervised learning can effectively leverage the unlabeled set $D_{\mathrm {ul}}$ to enhance predictive accuracy. Clearly, assuming $Y \to X$ exists, the smaller $n_1$ and the larger $n_2$ is, the more likely we would be able to detect the effectiveness of semi-supervised learning. From this perspective, a natural choice is to make $n_1$ the minimum labeled-sample size required by the semi-supervised learning algorithm, and to include all other data in the unlabeled set (i.e., $n_2 = n - n_1$ ). For example, the aforementioned self-training algorithm generally requires $n_1 \geq m$ to avoid an initial degenerate solution. Hence, we generate the labeled set $D_{\mathrm {lb}}$ and the unlabeled set $D_{\mathrm {ul}}$ by first randomly permuting the order of all data points in D, before setting

(14) $$ \begin{align} D_{\mathrm{lb}} &= \{\langle x_{i2}, \ldots, x_{im}, y_{i1} \rangle | i \in [1, m]\}, \mbox{ and} \end{align} $$

(15) $$ \begin{align} D_{\mathrm{ul}} &= \{\langle x_{i2}, \ldots, x_{im} \rangle | i \in [m+1, n]\}. \end{align} $$

4.1.2 Output

Our method for reverse causality testing focuses on the accuracy of $f_{\mathrm {ssl}}$ , the output of semi-supervised learning. To determine whether the unlabeled set $D_{\mathrm {ul}}$ is useful for improving the predictive accuracy of $f_{\mathrm {ssl}}$ (which is only possible when reverse causality exists), we compare $f_{\mathrm {ssl}}$ against a baseline model $f_0$ , which is the initial model of semi-supervised learning generated from the small labeled set $D_{\mathrm {lb}}$ only before accessing the unlabeled set $D_{\mathrm {ul}}$ . A smaller error of $f_{\mathrm {ssl}}$ would serve as evidence against the null hypothesis of no reverse causality.

Recall from earlier discussions that $D_{\mathrm {lb}}$ and $D_{\mathrm {ub}}$ are drawn uniformly at random from the input data D. Since $D_{\mathrm {lb}}$ and $D_{\mathrm {ub}}$ are random samples, comparing the predictive errors of the machine learning models trained on them—i.e., $f_0$ and $f_{\mathrm {ssl}}$ , respectively—requires a statistical significance test. Specifically, we need to determine whether any observed reduction in predictive error of $f_{\mathrm {ssl}}$ relative to $f_0$ is statistically significant. To assess this, we repeatedly draw i.i.d. samples of $D_{\mathrm {lb}}$ from D, derive the corresponding $D_{\mathrm {ub}}$ , and compare the predictive error of $f_{\mathrm {ssl}}$ against $f_0$ for each sample. The statistical significance of the reduction in predictive error can then be evaluated based on these comparisons. There are three key issues worth discussing in this design: (1) how to measure the predictive error of a model, (2) how to assess the statistical significance of the reduction in predictive error, and (3) how to determine the number of resamples necessary for the test.

First, in terms of measuring the predictive error of a machine learning model, a prevalent practice in machine learning is to separate the testing dataset from the training data due to concerns of overfitting (Bishop & Nasrabadi, 2006). It is important to note that overfitting is not a concern in our case because the prediction target $y_{i1}$ for the vast majority of data (i.e., $i \in [m+1, n]$ ) is hidden from semi-supervised learning. As such, we can assess the predictive error of $f_{\mathrm {ssl}}$ (and the baseline $f_0$ ) directly over the input dataset D. Further, our method has no specific requirement on the metric for predictive error, as a reduction of any error metric indicates the existence of reverse causality. For example, when $y_{i1}$ is binary, we could use the total number of prediction errors as the metric. The cross-entropy metric (Bishop & Nasrabadi, 2006) can be used for categorical $y_{i1}$ . When $y_{i1}$ is continuous, the error metric could be the mean squared error, mean absolute error, etc.

Second, for comparing the predictive errors of two machine learning models, Demšar (2006) reviewed three types of statistical tests used in the literature: (1) a parametric test, like the paired t-test (Dietterich, 1998), (2) a nonparametric test that makes assumptions about the distribution of difference in predictive accuracy, like the Wilcoxon signed-rank test (Santafe et al., 2015), and (3) a nonparametric test that only requires predictive accuracy to be comparable (i.e., at least on an ordinal scale), like the binomial test, also known as the sign test (Salzberg, 1997).

Problems with parametric tests such as the paired t-test have been well documented in the literature (Dietterich, 1998). As summarized by Demšar (2006), these tests suffer from sensitivity to outliers and the issue of commensurability between different runs, making it possible for the total failure of semi-supervised learning on a single sample of $D_{\mathrm {lb}}$ to dominate the test result. Similarly, key assumptions of the Wilcoxon signed-rank test do not hold in our context. For instance, it assumes that the difference in predictive accuracy between the two models is symmetrically distributed around a central value, yet there is no evidence to support this symmetry in our setting. In fact, semi-supervised learning is likely to amplify the skewness in the labeled set, suggesting that the difference could follow a heavy-tailed rather than symmetric distribution, violating the assumptions of the Wilcoxon signed-rank test. This leaves the (exact) binomial test (i.e., sign test) as a suitable alternative, as it requires only that each pair of predictive accuracy scores be comparable (Salzberg, 1997) and is inherently robust to outliers (Demšar, 2006). In our case, we employ the binomial test (with $\pi _0 = 0.5$ ) to compare the number of runs (i.e., $D_{\mathrm {lb}}$ samples) where $f_{\mathrm {ssl}}$ outperforms $f_0$ against the number of runs where the reverse is true.

Finally, two key factors influence the determination of the number of resamples required for the test. First, increasing the number of resamples enhances the statistical power of the test, which is especially critical given the well-documented limitation of the binomial test in terms of low statistical power (Demšar, 2006; Rainio et al., 2024; Salzberg, 1997). Second, a larger number of resamples incurs higher computational costs, as the semi-supervised learning algorithm must be executed for each sampled $D_{\mathrm {lb}}$ and $D_{\mathrm {ul}}$ . This computational burden can be particularly significant for resource-intensive algorithms, such as those based on deep learning (Goodfellow et al., 2016). To balance these considerations, we recommend following the default setting proposed by Demšar (2006) and performing 1,000 resamples. As demonstrated in the simulation studies, this setting offers a practical trade-off between statistical power and computational efficiency.

4.2 Design of semi-supervised learning algorithm

Numerous algorithms have been proposed for semi-supervised learning (Van Engelen & Hoos, 2020). We chose to implement the self-training paradigm discussed earlier. This choice is driven by two primary considerations. First, self-training is highly versatile, as it can be seamlessly integrated with a wide range of supervised learning algorithms as its base learner (Sohn et al., 2020). Second, the mathematical analysis of self-training is one of the few that have been thoroughly developed in the literature (Amini & Gallinari, 2002; Grandvalet & Bengio, 2004), providing a clear connection between the Bayes risk of its predictions and the foundational working condition of semi-supervised learning that has been discussed earlier in the article.

In the passages that follow, we describe the mathematical foundation and algorithmic design of self-training. Note that, our method represents a novel use of self-training for the purpose of reverse causality testing, and we did not make any change to its canonical algorithmic design. For the ease of discussion, we start with a simple setting where the prediction target $y_{i1}$ is binary (i.e., $y_{i1} \in \{0,1\}$ ) and the underlying learning algorithm is logistic regression. At the end of this section, we describe a generalization to continuous variables and any regression algorithm.

4.2.1 Mathematical foundation of self-training

Any algorithm aiming to learn a prediction model that estimates a binary $y_{i1}$ based on a predictor vector $\mathbf {x}_i = \langle x_{i2}, \ldots , x_{im} \rangle $ can be viewed as learning a function ⁵ $f(\mathbf {x}_i)$ that approximates $\Pr \{y_{i1} = 1 | \mathbf {x}_i\}$ . For example, logistic regression specifies $f(\mathbf {x}_i)$ as

(16) $$ \begin{align} f(\mathbf{x}_i) = \mathrm{sigmoid}(\beta_1 + \beta_2 \cdot x_{i2} + \cdots + \beta_m \cdot x_{im}), \end{align} $$

where are the coefficients to be estimated from training data.

Semi-supervised learning in general, and self-training in particular, starts by estimating from the small labeled set $D_{\mathrm {lb}}$ . Specifically, it does so by finding that maximizes the following log-likelihood function:

(17) $$ \begin{align} \ell(f) &= \log P(\mathbf{Y}_{\mathrm{lb}} | f, \mathbf{X}_{\mathrm{lb}}) \end{align} $$

(18) $$ \begin{align} &= \sum^{m}_{i=1} \left(y_{i1} \cdot \log(f(\mathbf{x}_i)) + (1-y_{i1}) \cdot \log(1-f(\mathbf{x}_i))\right), \end{align} $$

where $\mathbf {X}_{\mathrm {lb}}$ and $\mathbf {Y}_{\mathrm {lb}}$ represent the predictor (i.e., $x_{i2}, \ldots , x_{im}$ ) and prediction target (i.e., $y_{i1}$ ) portions of $D_{\mathrm {lb}}$ , respectively. In the Bayesian framework, these initial coefficient estimates form the maximum a posteriori (MAP) estimate under the uniform prior. As shown by Seeger (2000), this MAP estimate does not change if we merely add the unlabeled set $D_{\mathrm {ul}}$ into the observed data, because semi-supervised learning only works under the condition that the distribution of $\mathbf {x}_i$ reveals information about $y_{i1}$ . In other words, in order to proceed beyond this initial step and generate the MAP estimate for semi-supervised learning, we need to encode its working condition into the prior distribution used to calculate the MAP estimate.

In the Bayesian framework, a common method for deriving the prior distribution from a given constraint (e.g., the working condition for semi-supervised learning) is the principle of maximum entropy (Jaynes, 1957), which sets the prior distribution as the one that satisfies the given constraint while having the maximum information entropy (Cover & Thomas, 2006). Grandvalet & Bengio (2004) followed this principle to prove that, for a semi-supervised learning algorithm that uses both $D_{\mathrm {lb}}$ and $D_{\mathrm {ul}}$ , the MAP estimate for is the maximizer of the following criterion $C(f)$ ,

(19) $$ \begin{align} C(f) = \ell(f) - \lambda \cdot \sum^{n}_{i = m+1} \left(-f(\mathbf{x}_i) \cdot \log(f(\mathbf{x}_i)) - (1-f(\mathbf{x}_i)) \cdot \log(1-f(\mathbf{x}_i))\right), \end{align} $$

where $\lambda $ ( $\lambda> 0$ ) is the Lagrange multiplier that, roughly speaking, captures the amount of information about $y_{i1}$ that can be revealed by the distribution of $\mathbf {x}_i$ . Amini & Gallinari, 2002 analyzed self-training with logistic regression as the underlying learning algorithm, and proved that the resulting coefficient estimates maximizes $C(f)$ when $\lambda = 1$ . In other words, semi-supervised learning in general, and self-training in particular, can be viewed as approximating the MAP estimate of f under the working condition of semi-supervised learning. Specifically, it does so by finding f that maximizes $C(f)$ in Equation 19.

4.2.2 Algorithmic design of self-training

The algorithmic design of self-training can be readily derived from Equation 19. Note that the second term in the equation is proportional to the sum of

(20) $$ \begin{align} H(f(\mathbf{x}_i)) = -f(\mathbf{x}_i) \cdot \log(f(\mathbf{x}_i)) - (1-f(\mathbf{x}_i)) \cdot \log(1-f(\mathbf{x}_i)) \end{align} $$

for all data points in the unlabeled set $D_{\mathrm {ul}}$ . Recall that $f(\mathbf {x}_i)$ is the prediction from semi-supervised learning for $\Pr \{y_{i1}=1|\mathbf {x}_i\}$ . This makes $H(f(\mathbf {x}_i))$ in Equation 20 the entropy (Cover & Thomas, 2006) of the predicted distribution of $y_{i1}$ given $\mathbf {x}_i$ , meaning that it captures the amount of uncertainty in machine learning predictions. For example, we have $H(f(\mathbf {x}_i)) = 0$ , its minimum possible value, when semi-supervised learning is fully confident of its prediction (i.e., $f(\mathbf {x}_i) = 0$ or 1). In contrast, when $f(\mathbf {x}_i) = 0.5$ (i.e., maximum uncertainty), $H(f(\mathbf {x}_i))$ reaches its maximum value ⁶ of 1.

With this understanding, the goal of self-training—i.e., the maximization of $C(f)$ —pertaining to the unlabeled set is equivalent with minimizing prediction uncertainty (i.e., $H(f(\mathbf {x}_i))$ ) for the unlabeled data. Self-training does so in an iterative manner. As described earlier, the initial iteration estimates from the labeled set $D_{\mathrm {lb}}$ by applying the underlying learning algorithm, which in this case is the standard logistic regression. The estimated is then used to compute $f(\mathbf {x}_i)$ , and thereby $H(f(\mathbf {x}_i))$ , for all data points in the unlabeled set $D_{\mathrm {ul}}$ . The first iteration concludes by adding into the labeled set $D_{\mathrm {lb}}$ all unlabeled data points with prediction uncertainty below a pre-determined threshold h, using their predicted labels as if they were real. In other words, the labeled and unlabeled sets are updated as

(21)

(22) $$ \begin{align} D_{\mathrm{ul}} &:= \{\mathbf{x}_i | \mathbf{x}_i \in D_{\mathrm{ul}}, H(f(\mathbf{x}_i))> h\}, \end{align} $$

where is the predicted label for $\mathbf {x}_i$ —i.e., an indicator function that returns 1 if $f(\mathbf {x}_i) \geq 0.5$ and 0 otherwise. The updated $D_{\mathrm {lb}}$ and $D_{\mathrm {ul}}$ are then entered as input to the next iteration. The iterative process ends when no new X–Y pair is added to $D_{\mathrm {lb}}$ after an iteration. As can be seen from the description, self-training minimizes prediction uncertainty among unlabeled data using the idea of pseudo-labels, i.e., by promoting those with minimal prediction uncertainty to the labeled set, using their predicted labels as if they were real. These promoted data points, in turn, reduce prediction uncertainty for the remaining unlabeled data (Amini & Gallinari, 2002), pushing the coefficient estimates closer to their MAP values that maximize $C(f)$ .

4.2.3 Generalization to continuous variables

Whereas the above description of self-training is based on a binary $y_{i1}$ and logistic regression being the underlying learning algorithm, the same iterative process can be adapted to support continuous $y_{i1}$ and any underlying learning algorithm. This adaption requires addressing two issues. One is the design of an appropriate uncertainty measure (i.e., $H(f(\mathbf {x}_i))$ in the binary case), and the other is the generation of pseudo-label (i.e., in the binary case). The reason why these two issues arise for continuous variables is because, unlike in the binary case where $f(\mathbf {x}_i)$ , as an estimate of $\Pr \{y_{i1} = 1 | \mathbf {x}_i\}$ , inherently captures uncertainty about the predicted $y_{i1}$ , a point-estimate for a continuous $y_{i1}$ contains no such uncertainty information. Therefore, instead of deriving both the uncertainty measure and the pseudo-label from $f(\mathbf {x}_i)$ itself, we may have to resort to other information in order to do so in the continuous case.

A well-known method for addressing both issues in semi-supervised learning is co-training (Blum & Mitchell, 1998). With this method, instead of generating a single prediction of $f(\mathbf {x}_i) \approx y_{i1}$ , we generate two predictions, $f_1$ and $f_2$ , based on two (slightly) different subsets of variables in $\mathbf {x}_i$ . Then, the difference between $f_1$ and $f_2$ (i.e., $|f_1 - f_2|$ ) is a natural measure of uncertainty, while the mean of the two (i.e., $(f_1 + f_2)/2$ ) can be used as pseudo-label for expanding the labeled set.

More specifically, recall that our method includes in $\mathbf {x}_i$ a total of $m - 1$ variables $\langle x_{i2}, \ldots , x_{im} \rangle $ , where $m \geq 3$ because we require a minimum of three waves for identification. A natural choice is to associate $f_1$ with the first $m - 2$ variables $\mathbf {x}^{\mathrm {F}}_i = \langle x_{i2}, \ldots , x_{i,m-1} \rangle $ , and $f_2$ with the last $m - 2$ variables $\mathbf {x}^{\mathrm {L}}_i = \langle x_{i3}, \ldots , x_{im} \rangle $ . This way, we allow a divergence of two predictions (to allow the uncertainty estimate) while minimizing the number of predictors withheld from either. With this design, the updates of labeled and unlabeled set in each iteration become

(23) $$ \begin{align} D_{\mathrm{lb}} &:= D_{\mathrm{lb}} \cup \{\langle \mathbf{x}_i, (f_1(\mathbf{x}^{\mathrm{F}}_i) + f_2(\mathbf{x}^{\mathrm{L}}_i))/2 \rangle | \mathbf{x}_i \in D_{\mathrm{ul}}, |f_1(\mathbf{x}^{\mathrm{F}}_i) - f_2(\mathbf{x}^{\mathrm{L}}_i)| \leq h\}, \end{align} $$

(24) $$ \begin{align} D_{\mathrm{ul}} &:= \{\mathbf{x}_i | \mathbf{x}_i \in D_{\mathrm{ul}}, |f_1(\mathbf{x}^{\mathrm{F}}_i) - f_2(\mathbf{x}^{\mathrm{L}}_i)|> h\}, \end{align} $$

whereas everything else in the iterative process follows directly from the binary case. Clearly, this design for continuous $y_{i1}$ is compatible with any underlying learning algorithm for generating $f_1$ and $f_2$ , e.g., linear regression (Stine, 1985), support vector machines (SVMs; De Brabanter et al., 2010), neural networks (Heskes, 1996), etc. The pseudocode of our algorithm for continuous variables is available in Table 1.

Table 1 Pseudocode for reverse causality testing with continuous variables

4.3 Transparency and openness

The complete code implementation of our algorithm is publicly available at https://github.com/calearn/revc (Python) and https://github.com/calearn/revc_m (MATLAB).

5 Simulation studies

We conducted two simulation studies. The main study evaluates the statistical power of our method in identifying reverse causality, and a followup study assesses the Type I error rate of our method in the absence of reverse causality, because existing methods such as the random-intercept CLPM (RI-CLPM; Hamaker et al., 2015) are known to generate Type I errors.

5.1 Data-generating process

In the main simulation study, we aimed to delineate the primary factors influencing the statistical power of our method. To achieve this, we employed a straightforward structural model-Equations 1 and 2—as the data-generating process, establishing a clear ground truth for the causal direction. To ensure a comprehensive analysis, we created multiple levels for three key parameters in the data-generating process: the total number of observations N and the cross-lagged parameters $\beta $ (i.e., $Y \to X$ ) and $\gamma $ (i.e., $X \to Y$ ). For the sample size N, we created four levels: 100, 250, 500, and 1,000. For the cross-lagged parameters $\beta $ and $\gamma $ , we created six levels for each: 0, 0.1, 0.2, 0.3, 0.4, and 0.5. In total, the design for the main simulation study consists of $4$ (N) $\times\ $ $6$ ( $\beta $ ) $\times\ $ $6$ ( $\gamma $ ) = $144$ unique conditions.

We adopted a standard parameter setup (e.g., Hamaker et al., 2015) to configure the other fixed parameters for the data-generating process in the main simulation study. Specifically, we followed Hamaker et al. (2015) to set the autoregressive parameters for both variables to $\alpha = \delta = 0.5$ . All random impulses $u_{it}$ and $v_{it}$ were generated from a Gaussian distribution $\mathcal {N}(0, 0.5^2)$ , while the sample-specific random intercepts were fixed at $\kappa _i = \omega _i = 0$ for all $i \in [1,\ N]$ . To generate the initial (Wave 1) values of $x_{i1}$ and $y_{i1}$ , we deliberately used different distributions ⁷ to emphasize that our method does not rely on specific distributional assumptions for the input data. For each $i \in [1,\ N]$ , $x_{i1}$ was sampled uniformly at random from the interval $[-3,\ 3]$ , while $y_{i1}$ was drawn from a Gaussian distribution $\mathcal {N}(0,\ 1)$ .

In the followup simulation study, we followed Lucas (2023), which demonstrates the spurious cross-lagged effects generated by CLPM, in adopting the widely used Stable Trait Autoregressive Trait and State (STARTS) model (Kenny & Zautra, 1995) as the underlying data-generating process. We also followed Lucas (2023) in the parameter setup, specifically by setting (for both X and Y) the stability parameter as 0.5, the random intercept variance as 1, and variance of autoregressive component as 1. We also added a measurement error of variance 0.3 to both X and Y. Since the focus of the followup study is on Type I errors, the cross-lagged parameter for reverse causality ( $Y \to X$ ) was always set to zero, while the parameter for $X \to Y$ (represented as $\gamma $ for consistency with the main study) was varied from 0.1 to 0.5. Like in the main study, we created four levels, 100, 250, 500, and 1,000, for the sample size N. In total, the design for the followup simulation study consists of $4$ (N) $\times\ $ $5$ ( $\gamma $ ) = $20$ unique conditions.

5.2 Algorithmic implementations

Recall from earlier discussions that our method applies to both discrete and continuous Y, and can be used with any supervised learning algorithm as its base learner. To demonstrate the versatility of our method, we implemented two variants of it. The first was designed for continuous Y and implemented in MATLAB R2023b (with statistics and machine learning toolbox). We selected a simple design, i.e., OLS regression, as the underlying learning algorithm. We set the number of resamples for binomial test to be 1,000, leading to a computational overhead of about 20 seconds per simulation run on a laptop computer with Apple M2 CPU and 8GB RAM. For the number of waves taken as input, we set $m = 3$ (i.e., the first three waves), the minimum value that satisfies the identification requirement of our method while providing a conservative estimate of its statistical power. For the uncertainty threshold h (i.e., the maximum difference in prediction between the two models for an unlabeled data point to be added to the labeled set, see Equation 24), after testing a wide range of threshold values, we found that the output of our method is insensitive to the threshold setting as long as it is neither too large—so as to allow all unlabeled data to enter the labeled set at once—nor too small to permit any unlabeled data point into the labeled set. Due to this finding, we set the threshold to a constant of $h = 0.1$ , which is about 10% of the standard deviation of the label (i.e., Y).

Since the OLS regression algorithm relies on a linear model and our data-generating process is also linear, concerns may arise regarding a potential unfair advantage due to their inherent consistency. To address these concerns, we implemented a variant of our method using a nonlinear learning model. Specifically, for the underlying base learner in self-training, we employed the SVM (Cortes & Vapnik, 1995) algorithm with a polynomial kernel of degree 3 (i.e., a cubic kernel SVM).

Given the significantly higher computational overhead of SVM, we implemented this variant in Python using the scikit-learn (Pedregosa et al., 2011) machine learning library. We adhered to the default settings of scikit-learn’s LIBSVM implementation (Chang & Lin, 2011) for all SVM parameters, including the kernel configuration (i.e., in function sklearn.svm.SVC). For the uncertainty threshold in the self-training algorithm, we used the built-in options for the criterion parameter in sklearn.semi_supervised.SelfTrainingClassifier. Apart from the use of SVM as the base learner, all other design elements were identical to the OLS implementation, except for an additional input data normalization process required for the SVM implementation. Specifically, Y was converted to binary values (using median dichotomization), and X was normalized to the range $[-1,\ 1]$ (through min-max scaling) to account for SVM’s sensitivity to input feature scaling (Chang & Lin, 2011; Tax & Duin, 2004).

To accommodate the computational demands of the nonlinear SVM algorithm and the large number of simulation runs required for the main study, this implementation was executed on Amazon Web Services (AWS) Batch using Elastic Container Service (ECS) clusters configured with 256 virtual CPUs, provisioned and scaled automatically using AWS Fargate. With the same number of resamples (1,000) as the OLS implementation, the SVM variant completed each simulation run in approximately 30 seconds on the AWS cluster.

5.3 Simulation results

For the main study, Table 2 presents the statistical power achieved by our method under a significance level of $\alpha = .05$ for all simulation settings with $\beta> 0$ . For settings where $\beta = 0$ (highlighted in gray), the table reports the Type I error rates. Each cell in the table is based on 1,000 independent runs of our method. Note that results for the nonlinear SVM-based implementation are included only for $N = 1,000$ , as the SVM-based implementation required more than 500 samples to consistently converge to reliable predictions.

Table 2 Type I error rate and statistical power of our method in the main simulation study

We can draw three key observations from the table. First, regarding Type I error rates (highlighted rows with $\beta = 0$ ), our method consistently remains below 0.075 across all simulation conditions, aligning with Bradley’s (1978) liberal robustness criterion. This demonstrates the robustness of our method against finding spurious reverse causal effects.

Second, note that the top six rows of Table 2, along with the highlighted rows, represent simulation conditions where X and Y exhibit a unidirectional relationship. The statistical power achieved by our method in these scenarios highlights its effectiveness in identifying the direction of a unidirectional effect—a critical use-case for our method when competing theories stipulate different causal directions between X and Y. For example, Rows 3–6 show that, when $\beta \geq 0.2$ , the OLS implementation of our method achieves statistical power exceeding 0.97 when $N \geq 500$ , while the SVM implementation achieves the same when $N = 1,000$ .

Third, the statistical power of our method is influenced by three key factors. One factor is N, the input sample size. Larger sample sizes significantly improve statistical power. For example, when $\beta = 0.3$ and $\gamma = 0$ , the OLS implementation achieves a statistical power of only 0.44 with $N = 100$ but exceeds 0.99 when $N \geq 500$ . Another factor is $\beta $ , the cross-lagged parameter representing the strength of the reverse causal effect. When $\gamma = 0$ (i.e., unidirectional effect), larger $\beta $ values correspond to higher statistical power. For instance, with $N = 250$ and $\gamma = 0$ , the OLS implementation achieves a power of 0.19 at $\beta = 0.1$ , increasing to 0.77 at $\beta = 0.3$ , and exceeding 0.99 for $\beta \geq 0.4$ . The SVM implementation exhibits a similar trend. The final factor is $\gamma $ , the strength of the causal effect $X \to Y$ . The impact of $\gamma $ on statistical power varies depending on other parameters. For example, when $\beta = 0.5$ and $N = 1,000$ , a stronger forward effect of $\gamma = 0.5$ can obscure reverse causality, reducing power from 1.00 at $\gamma = 0$ to 0.86 at $\gamma = 0.5$ . Conversely, when $\gamma $ and $\beta $ are smaller yet closer in magnitude (e.g., $\gamma = 0.2$ , $\beta = 0.2$ , $N = 250$ ), $\gamma $ could potentially enhance power, increasing it from 0.30 at $\gamma = 0$ to 0.67 at $\gamma = 0.2$ . This increase can be attributed to the technical design of the two-step self-training algorithm in our method. Specifically, in situations where the causal effect $X \to Y$ is nonexistent (i.e., $\gamma = 0$ ), the initial step of self-training, which learns from the labeled set only, is destined to yield highly inaccurate predictions. In such scenarios, particularly with smaller sample sizes, it becomes challenging for the subsequent step (of learning from the unlabeled set) to significantly enhance predictive accuracy even in the presence of reverse causality. Therefore, a larger $\gamma $ value, which improves the accuracy of the initial step, also tends to boost the statistical power.

Interestingly, the negative impact of $\gamma $ on the statistical power of our method appears more pronounced in the SVM implementation compared to the OLS implementation. For example, with $N = 1,000$ , the OLS implementation maintains a statistical power of at least 0.70 for $\beta \geq 0.2$ , regardless of $\gamma $ . In contrast, the SVM implementation’s power drops from 1.00 at $\gamma = 0$ to 0.43 at $\gamma = 0.5$ for $\beta = 0.5$ . We attribute this discrepancy to SVM’s sensitivity to feature scaling (Tax & Duin, 2004) and vulnerability to outliers (Debruyne, 2009), both of which are exacerbated at larger $\gamma $ .

Table 3 summarizes the results of the follow-up simulation study comparing the Type I error rates of RI-CLPM and the OLS implementation of our method. ⁸ Following Mulder & Hamaker (2021), RI-CLPM was implemented using default settings in Mplus. ⁹ The results reveal that RI-CLPM frequently identifies spurious reverse causal relationships, with Type I error rates reaching 0.88 when $N = 1,000$ and $\gamma = 0.5$ . In general, RI-CLPM’s Type I error rate increases with larger N and higher $\gamma $ . In contrast, our method maintains a Type I error rate below 0.05 across all conditions, again aligning with Bradley’s (1978) liberal robustness criterion and demonstrating an advantage of our method over CLPM and its variants, which are known to find spurious cross-lagged effects (Lucas, 2023).

Table 3 Type I error rate of RI-CLPM and our method under STARTS model

6 A case study

We applied our method over a real-world panel dataset to demonstrate its value for reverse causality testing in practice. Specifically, we examined the relationship between work-family conflict (WFC)—i.e., “a form of inter-role conflict in which the role pressures from the work and family domains are mutually incompatible in some respect” (Greenhaus & Beutell, 1985, p. 77)—and job satisfaction (JAS)—i.e., the “overall evaluative judgement one has about one’s job” (Judge et al., 2017, p. 357). We chose this focal relationship for two main reasons. First, like many constructs in industrial-organizational psychology, it is often practically difficult, if not ethically dubious, to manipulate WFC or JAS in a randomized controlled trial. Second, whereas a robust correlation between WFC and JAS has been widely recognized (Allen et al., 2020; Amstad et al., 2011), there are ongoing debates about the causal direction between the two constructs. Some posit a unidirectional effect of WFC negatively affecting JAS (e.g., Allen et al., 2020; Amstad et al., 2011; Kossek & Ozeki, 1998). Others contend that the causal influence flows in the opposite direction, as higher JAS leads to greater work-life balance and, correspondingly, lower WFC (Landolfi et al., 2022). Yet others suggest the existence of reciprocal effects, with WFC and JAS affecting each other over time (e.g., Demerouti et al., 2004). Given these varied viewpoints, we sought to apply our method to empirically test the existence of causal effect in either direction.

6.1 Data

We drew from the Swiss Household Panel (SHP Group, 2023; Voorpostel et al., 2023) to examine the WFC–JAS relationship. The SHP began in 1999 with a nationally representative sample of 5,074 Swiss Households, introducing supplementary samples in 2004 (addition of 2,537 households), 2013 (addition of 3,989 households), and 2020 (addition of 4,380 households; Voorpostel et al., 2023). All household members aged 14 and above are surveyed annually via telephone and written surveys. We mirrored the timeframe used in prior work on CLPM and analyzed data collected in Waves 6 through 9 (i.e., the annual surveys between 2004 and 2007; Ozkok et al., 2022). We excluded participants who did not work or who did not complete a single survey in our timeframe, resulting in N = 7,748. Participants (51.55% women, 48.45% men) were on average 39.28 (SD = 14.50) years old. On average, participants had 12.98 (SD = 3.20) years of education.

In the case study, WFC was measured by a single item on a 11-point Likert scale (0 = not at all, 10 = very strongly). The item read “How strongly does your work interfere with your private activities and family obligations more than you would want this to be?” JAS was measured by a six-item measure created for the SHP on a 11-point Likert scale (0 = not at all satisfied, 10 = completely satisfied). The lead-in to the items was “Can you indicate your degree of satisfaction for each of the following points?” Sample items include “your job in general” and “the amount of your work.” The internal consistency for the four waves was .78, .78, .79, and .79, respectively.

6.2 Results

In terms of algorithmic implementations, we used the same implementations (OLS and SVM) as the simulation study, and set the number of resamples to 10,000 and 100, respectively. Note that we intentionally set the number of resamples for the OLS implementation to be much higher than required in order to support a detailed analysis described later. With either implementation, we first tested the presence of WFC $\to $ JAS by setting the self-training algorithm to predict WFC from JAS, before testing the presence of JAS $\to $ WFC by setting the self-training algorithm to predict JAS from WFC. For the OLS implementation, self-training was effective for predicting WFC from JAS in 5,743 out of 10,000 resamples (one-tailed binomial exact test with $\pi _0 = 0.5$ : $p < 10^{-6}$ ), and was effective for predicting JAS from WFC in 5,006 out of 10,000 resamples (one-tailed binomial exact test with $\pi _0 = 0.5$ : $p = 0.4562$ ). For the nonlinear SVM implementation, the results were consistent with the OLS implementation. Specifically, for the prediction of WFC from JAS, self-training was effective in 67 out of 100 iterations, leading to a p-value of $p < 0.0005$ . On the flip side, when predicting JAS from WFC, self-training was effective only in 37 out of 100 runs, ¹⁰ leading to a p-value of $p = 0.9967$ . In sum, the results suggest that work-family conflict influences JAS, but not vice versa.

To further inspect why semi-supervised learning succeed in predicting WFC but not JAS, we examined how its effectiveness depends on the small labeled set, which is the only information source about Y that semi-supervised learning receives. Recall from the design of self-training that, if the distribution of X reveals signals about the X–Y relationship, we would expect self-training to be more effective when the labeled set is a more representative sample of the full dataset. This is because, when the labeled set is a severely biased sample, say featuring only a single value of Y, then self-training would stand no chance in becoming effective, as it would not even know what the other values of Y might be.

To inspect how the effectiveness of semi-supervised learning varies with the bias of labeled sample, we leveraged the 10,000 resample runs executed for the OLS implementation. Specifically, we sorted all runs according to the distributional distance (as measured by Kolmogorov–Smirnov test statistic; Daniel, 1990) between the labeled set and the full dataset, before stratifying them into 10 equi-sized bins (i.e., each containing 1,000 runs) according to the sorted order and calculating the average effectiveness of semi-supervised learning for each bin as measured by the fraction of runs that improve predictive accuracy minus the fraction of those that reduce accuracy. As discussed before, if self-training were ineffective, we would expect the average effectiveness measure to hover around zero for all ten bins. In contrast, if self-training were effective, we would expect the average effectiveness to be high when the distributional distance is small, and gradually decrease when the distributional distance becomes larger.

As can be seen from Figure 2, when self-training is used to predict WFC (i.e., to test the presence of WFC $\to $ JAS), there is a clear, negative, correlation between its effectiveness and the distributional distance. In contrast, such correlation disappears in the prediction of JSA (i.e., testing of JAS $\to $ WFC), as the effectiveness of self-training hovers around zero regardless of the distributional distance. This clearly shows that, whereas the distribution of JSA reveals valuable signals for predicting WFC—meaning that WFC may have a causal effect on JAS—the distribution of WFC carries no information for predicting JAS, suggesting that the causal direction is unlikely to flow from JAS to WFC.

Figure 2 Relationship between effectiveness of self-training and representativeness of labeled set.Note: WFC, work-family conflict; JSA, job satisfaction. Source: Swiss Household Panel (SHP).

7 General discussion

In this section, we first discuss the research implications of our new method. We then review the limitations of our method and the future research needed to address them.

7.1 Research implications

For testing reverse causality, the main research implication of our semi-supervised learning based method is its ability to identify causal directions without presuming distributional or functional features of the data-generating process. Our method achieves this by leveraging recent advances in machine learning that directly link the causal direction with the effectiveness of semi-supervised learning. A unique feature of our method is that it is an umbrella algorithm independent of the underlying learning algorithm, which can be any algorithm for supervised learning. As such, researchers can freely choose a learning algorithm that fits the type and scale of their data before using our method for causal identification.

More broadly, the main research implication of our work is its transfer of insights of causal learning (Peters et al., 2017), in particular Schölkopf et al.’s (2012) seminal finding that links causal direction with the effectiveness of semi-supervised learning, into the methodological arsenal of causal inference in psychology. Whereas most existing work in causal learning—perhaps owing to the disciplinary focus of machine learning—leverages this link to explain (Schölkopf et al., 2012) or improve (Kügelgen et al., 2019) the performance of machine learning algorithms, our contribution is to demonstrate that the exact same link can be used to develop a concrete method for identifying reverse causality, addressing long-standing challenges in the analysis of longitudinal data in psychology and other disciplines.

7.2 Limitations and future directions

7.2.1 Capability of semi-supervised learning

The statistical power of our method depends on the capability of the semi-supervised learning algorithm in approximating the functional relationship between X and Y. On the one hand, this allows our method to generalize beyond linearity to allow nonlinear X–Y relationships. On the other hand, it also brings about the limitation that, if the semi-supervised learning algorithm being used cannot effectively approximate a nonlinear relationship between X and Y, our method may generate Type II errors.

In the main simulation study in our current work, we employed a linear model (i.e., CLPM) as the underlying data-generating process. This choice allowed us to concentrate on the overarching design of our method rather than on fine-tuning the underlying learning algorithm, given that even a simple algorithm (like OLS) would likely suffice for a linear relationship. Nonetheless, although the theoretical foundations of our method extend readily to more complex, nonlinear relationships between X and Y, we did not evaluate the performance of our method in the presence of such relationships. As a natural next step, future research could provide a comprehensive empirical assessment of the statistical power of our method across a wider range of data-generating processes, encompassing both linear and nonlinear scenarios beyond those considered in this study. Additionally, exploring various underlying learning algorithms and other semi-supervised learning designs will be valuable for capturing the inherent complexity of nonlinear relationships between X and Y.

In the longer term, recent advancements in machine learning offer two key insights. First, there are semi-supervised learning algorithms that offer universal approximation (Goodfellow et al., 2016), meaning that they are theoretically capable of approximating any arbitrary function in a Euclidean space. Examples include the use of our method (i.e., self-training) with deep neural networks (Hornik et al., 1989) or variational Gaussian processes (Tran et al., 2016) as the underlying learning algorithm. Second, unfortunately, one would have to restrict the type of relationship between X and Y in order to translate such theoretical feasibility into any practical guarantee. To understand why, consider a stylized example where Y is the encrypted value (i.e., ciphertext) of X based on a secret key. Theoretically, it is feasible for a machine learning algorithm to eventually learn how to predict Y from X. In practice, doing so constitutes a brute-force ciphertext-only attack against the encryption algorithm, which is commonly believed to be practically infeasible (Goldreich, 2004).

7.2.2 Leveraging other advances in machine learning

A central premise of our method is the link between causal direction and the predictive accuracy of machine learning algorithms. When the causal direction flows only from X to Y (i.e., $X \to Y$ ), the distribution of X would reveal no information about the X–Y relationship, rendering semi-supervised learning ineffective.

Following the same logic, causal direction could also affect the effectiveness of machine learning algorithms besides semi-supervised learning. For example, a common generalizability issue facing supervised learning, covariate shift (Sugiyama et al., 2007), arises when a machine learning model degrades in predictive accuracy because of the distributional differences in predictor variables (i.e., X) between the training dataset and the dataset actually in need of prediction. As discussed before, it would be impossible for covariate shift to arise if the causal direction flows only from X to Y because, in this case, a change of X’s distribution should have no bearing on the X–Y relationship. In machine learning, Kügelgen et al. (2019) leveraged this property to address the degradation of predictive accuracy caused by covariate shift. Similarly, future research could examine the use of covariate shift or, more broadly, the generalizability of supervised learning models to help identify the causal relationship between X and Y.

7.2.3 Integration of our method with existing methods

A key feature of our method is robustness to model misspecification, as it imposes no restrictions on the functional form of the data-generating process (e.g., linearity or the addition of independent noise). This flexibility minimizes the risk of spurious discoveries, such as the incorrect inference of reverse causality, caused by inappropriate model specifications or dependency structures. However, a notable limitation of our method lies in its relatively weak statistical power. For instance, our simulations demonstrate that it may struggle to detect reverse causality when the sample size is small. As such, our method may be better suited to larger datasets (e.g., with $N> 500$ ).

These characteristics make our method a valuable complement to existing approaches, such as additive noise models and DDA, which can be sensitive to model misspecifications (e.g., Schultheiss & Bühlmann, 2024; Thoemmes, 2015). For example, while existing methods may not be ideally suited for exploratory investigations into causal directions between variables (Wiedermann & von Eye, 2015), our method can serve this purpose effectively. Researchers could first employ our approach as an exploratory tool to identify promising relationships. Subsequently, they may leverage existing methods to conduct confirmatory analyses of those relationships supported by robust theoretical formulations but not identifiable through our approach (e.g., in cases where neither $X \to Y$ nor $Y \to X$ is detected).

Our method can also be seamlessly integrated with panel models such as CLPM and its variants. To illustrate, applying our method to infer the causal direction between X and Y could result in one of four possible outcomes: (1) the detection of $Y \to X$ but not $X \to Y$ , (2) the detection of $X \to Y$ but not $Y \to X$ , (3) the detection of both $X \to Y$ and $Y \to X$ , suggesting reciprocal relationships or the potential existence of an unobserved confounder, and (4) the detection of neither $X \to Y$ nor $Y \to X$ . With the first two outcomes, our model informs the model specification for CLPM, as researchers could choose to remove the cross-lagged effect inconsistent with the outcome of our method. An important future direction for research is to study whether causal learning could be used not only for identifying the causal direction (as in our work) but also for estimating the temporal lag of an effect, which would further improve the model specification for CLPM. The third outcome could prompt researchers to consider panel models that allow for certain types of latent confounders—e.g., the RI-CLPM (Hamaker et al., 2015), which models time-invariant confounders, or Latent Curve Models with Structured Residuals (LCM-CR; Curran et al., 2014), which models time-varying confounders with patterns such as autoregressive structures. Finally, as discussed earlier, researchers facing the fourth outcome could consider the use of existing confirmatory methods, e.g., additive noise models and DDA, to identify the causal direction.

8 Conclusion

In conclusion, our work proposes a novel method that integrates machine learning and reverse causality testing. Through mathematic analysis, simulation studies, and a case illustration, we demonstrates the effectiveness of this method. We hope this approach inspires future research to more strongly embrace the advancements in computer science and machine learning to enrich the methodological toolkit of psychological sciences.