2.1 Reverse causality testing

Compared with the rich and growing literature on omitted variable bias—in psychology (Harring et al., Reference Harring, McNeish and Hancock2017), sociology (Halaby, Reference Halaby2004) and economics (Wüthrich & Zhu, Reference Wüthrich and Zhu2023)—researchers across disciplines have made relatively limited progress on the testing of reverse causality over panel data. As summarized by Leszczensky & Wolbring (Reference Leszczensky and Wolbring2022), 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., Reference Orth, Clark, Donnellan and Robins2021), such as lagged first-difference (LFD) models (Vaisey & Miles, Reference Vaisey and Miles2017), CLPM (Hamaker et al., Reference Hamaker, Kuiper and Grasman2015), 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 (Reference Vaisey and Miles2017) show that panel models such as LFDs can be highly sensitive to the misspecification of temporal lags. Leszczensky & Wolbring (Reference Leszczensky and Wolbring2022, 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, Reference Shojaie and Fox2022).

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, Reference Li and Wiedermann2020; Pornprasertmanit & Little, Reference Pornprasertmanit and Little2012; von Eye & DeShon, Reference von Eye and DeShon2012; Wiedermann & Li, Reference Wiedermann and Li2018) and Linear Non-Gaussian Acyclic Models (LiNGAM; Shimizu et al., Reference Shimizu, Hoyer, Hyvärinen, Kerminen and Jordan2006)—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 departureFootnote ² from normality across different variables. For longitudinal data, extensions of these methods (e.g., Bauer et al., Reference Bauer, Schölkopf and Peters2016; Geiger et al., Reference Geiger, Zhang, Schoelkopf, Gong and Janzing2015; Hyvärinen et al., Reference Hyvärinen, Zhang, Shimizu and Hoyer2010) leverage the non-normality of noise to identify causal direction in multivariate time series. Additionally, methods like Rosenström et al.’s (Reference Rosenström, Czajkowski, Solbakken and Saarni2023) directional analysis approach, additive noise models (Hoyer et al., Reference Hoyer, Janzing, Mooij, Peters and Schölkopf2008; Peters et al., Reference Peters, Mooij, Janzing and Schölkopf2014) and (more generally) post-nonlinear models (Zhang & Hyvärinen, Reference Zhang and Hyvärinen2009) 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., Reference Xie, Luong, Hovy and Le2020). 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, Reference Van Engelen and Hoos2020). A famous example, which we will further elaborate later in the article, is self-training (e.g., Sohn et al., Reference Sohn, Berthelot, Carlini, Zhang, Zhang, Raffel, Cubuk, Kurakin and Li2020), 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, Reference Van Engelen and Hoos2020), 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. (Reference Buja, Brown, Berk, George, Pitkin, Traskin, Zhang and Zhao2019) 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 (Reference Niyogi2013) attributes the success of semi-supervised learning to the concept of manifold regularization (Belkin et al., Reference Belkin, Niyogi and Sindhwani2006). At its core, Niyogi (Reference Niyogi2013) 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., Reference Aggarwal, Hinneburg and Keim2001). 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, Reference Niyogi2013).

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., Reference Schölkopf, Janzing, Peters, Sgouritsa, Zhang and Mooij2012). 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.1 Conceptual illustration

In what follows, we offer an intuitive explanation for the key insight of Schölkopf et al. (Reference Schölkopf, Janzing, Peters, Sgouritsa, Zhang and Mooij2012) 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, Reference Janzing and Schölkopf2010).

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 inferenceFootnote ³ 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. (Reference Schölkopf, Janzing, Peters, Sgouritsa, Zhang and Mooij2012) 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, Reference Janzing and Schölkopf2010), (2) by examining the uncorrelatedness between p and the derivative of f (Janzing & Schölkopf, Reference Janzing and Schölkopf2015), and (3) through the uncorrelatedness between p and the logarithm of the derivative of f (Daniusis et al., Reference Daniusis, Janzing, Mooij, Zscheischler, Steudel, Zhang and Schölkopf2010). Whereas the mathematical formulation of Schölkopf et al.’s (Reference Schölkopf, Janzing, Peters, Sgouritsa, Zhang and Mooij2012) insight necessarily varies with this underlying notion of independence—e.g., see Janzing and Schölkopf (Reference Janzing and Schölkopf2015, 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., Reference Hamaker, Kuiper and Grasman2015) 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.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 independentFootnote ⁴ of f if $X \to Y$ fully describes the relationship between X and Y (Schölkopf et al., Reference Schölkopf, Janzing, Peters, Sgouritsa, Zhang and Mooij2012). 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.2 Design of semi-supervised learning algorithm

Numerous algorithms have been proposed for semi-supervised learning (Van Engelen & Hoos, Reference Van Engelen and Hoos2020). 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., Reference Sohn, Berthelot, Carlini, Zhang, Zhang, Raffel, Cubuk, Kurakin and Li2020). Second, the mathematical analysis of self-training is one of the few that have been thoroughly developed in the literature (Amini & Gallinari, Reference Amini and Gallinari2002; Grandvalet & Bengio, Reference Grandvalet and Bengio2004), 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 functionFootnote ⁵ $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 (Reference Seeger2000), 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, Reference Jaynes1957), which sets the prior distribution as the one that satisfies the given constraint while having the maximum information entropy (Cover & Thomas, Reference Cover and Thomas2006). Grandvalet & Bengio (Reference Grandvalet and Bengio2004) 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, Reference Amini and Gallinari2002 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, Reference Cover and Thomas2006) 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 valueFootnote ⁶ 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, Reference Amini and Gallinari2002), 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, Reference Blum and Mitchell1998). 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, Reference Stine1985), support vector machines (SVMs; De Brabanter et al., Reference De Brabanter, De Brabanter, Suykens and De Moor2010), neural networks (Heskes, Reference Heskes1996), 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.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., Reference Hamaker, Kuiper and Grasman2015) to configure the other fixed parameters for the data-generating process in the main simulation study. Specifically, we followed Hamaker et al. (Reference Hamaker, Kuiper and Grasman2015) 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 distributionsFootnote ⁷ 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 (Reference Lucas2023), 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, Reference Kenny and Zautra1995) as the underlying data-generating process. We also followed Lucas (Reference Lucas2023) 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, Reference Cortes and Vapnik1995) 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., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot and Duchesnay2011) machine learning library. We adhered to the default settings of scikit-learn’s LIBSVM implementation (Chang & Lin, Reference Chang and Lin2011) 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, Reference Chang and Lin2011; Tax & Duin, Reference Tax and Duin2004).

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 (Reference Bradley1978) 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, Reference Tax and Duin2004) and vulnerability to outliers (Debruyne, Reference Debruyne2009), 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.Footnote ⁸ Following Mulder & Hamaker (Reference Mulder and Hamaker2021), RI-CLPM was implemented using default settings in Mplus.Footnote ⁹ 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 (Reference Bradley1978) 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, Reference Lucas2023).

Table 3

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

6.1 Data

We drew from the Swiss Household Panel (SHP Group, 2023; Voorpostel et al., Reference Voorpostel, Tillmann, Lebert, Kuhn, Lipps, Ryser, Antal, Dasoki and Wernli2023) 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., Reference Voorpostel, Tillmann, Lebert, Kuhn, Lipps, Ryser, Antal, Dasoki and Wernli2023). 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., Reference Ozkok, Vaulont, Zyphur, Zhang, Preacher, Koval and Zheng2022). 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.

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., Reference Peters, Janzing and Schölkopf2017), in particular Schölkopf et al.’s (Reference von Eye and DeShon2012) 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., Reference Schölkopf, Janzing, Peters, Sgouritsa, Zhang and Mooij2012) or improve (Kügelgen et al., Reference Kügelgen, Mey and Loog2019) 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