Sample sizes in the 13 studies varied between 142 and 2,235 (M = 856.5). Data were collected in eight different countries (China, Finland, Italy, Norway, Portugal, Romania, The Netherlands, and USA) and the populations included, for example, school children, firefighters, and general practitioners. Percentage of male participants varied between 0 and 90% (M = 49.1%) and mean age at the first measurement varied between 12.7 and 48.6 years (M = 36.4 years).

The Maslach Burnout Inventory (MBI, Maslach et al., 1996) was the most commonly used instrument to measure burnout or exhaustion (in 7 of 13 studies). The predictors varied between the studies and were, consequently, measured with different instruments. However, the instruments appeared to have satisfactory homogeneity, with Cronbach’s alpha at the initial measurement varying between 0.66 and 0.92 (M = 0.81).

For each of the 23 effects in Table 1, we simulated data with the same sample size and six correlations between the predictor and burnout/exhaustion measured at two occasions. These correlations were reported in the reanalyzed studies and they are included in our descriptive dataset available at the Open Science Framework at https://osf.io/smy5n/. In each simulated dataset, the four variables (i.e., P₁, P₂, E₁, and E₂ in Figure 1) were drawn from a standard (M = 0, SD = 1) normal distribution. We did not include any missing values in the simulated data, meaning that procedures for handling missing data were not required. We used simulated data as the original data were not available to us. It is important to note that standardized regression effects are functions of correlations. The standardized effect of X₁ on Y₂ when adjusting for Y₁ is given by Equation 1. (Cohen et al., 2003) and the effect of X₁ on the Y₂-Y₁ difference score (all three variables standardized) is given by Equation 2 (Guilford, 1965). This means that if a simulated dataset has the same correlations between variables as an empirical dataset, regression effects estimated in the simulated dataset will be the same as if estimated in the empirical dataset. This is true even if the simulated dataset does not match the empirical dataset in some other regards, e.g., the distribution of the variables. Moreover, if the simulated dataset has the same sample size, the statistical significance of the regression effects will be the same as if estimated in the empirical dataset. (1) β X 1 , Y 2. Y 1 = r X 1 , Y 2 − r X 1 , Y 1 r Y 1 , Y 2 1 − r X 1 , Y 1 2 (2) β X 1 , Y 2 − Y 1 = r X 1 , Y 2 − r X 1 , Y 1 2 ( 1 − r Y 1 , Y 2 )

We fitted four complementary models to the simulated data: (1) A traditional CLPM where subsequent burnout/exhaustion was regressed on prior exhaustion/burnout and a prior score on the predictor and vice versa (Figure 1A); (2) A reversed CLPM where initial burnout/exhaustion was regressed on subsequent burnout/exhaustion and an initial score on the predictor (Figure 1B). This model was based on the logic that time-reversal should result in reversed signs of effects (Campbell and Kenny, 1999; Haufe et al., 2013). For example, if initial depression had an increasing effect on burnout, we should expect a negative effect of initial depression on concurrent burnout when adjusting for subsequent burnout. This negative effect would suggest that low initial depression had counteracted high initial burnout and allowed individuals to reach the same subsequent level of burnout as individuals who had a lower initial level of burnout but also a higher initial degree of depression; (3) A latent change score model (LCSM; McArdle, 2009; Ghisletta and McArdle, 2012; Kievit et al., 2018), where an initial score on the predictor predicted subsequent change in burnout/exhaustion and vice versa (Figure 1C). With a genuine increasing or decreasing effect of the predictor on burnout/exhaustion, this effect should be positive or negative, respectively; (4) A model of artifactualness, where initial and subsequent scores on the predictor and burnout/exhaustion were regressed on latent general predictor and burnout/exhaustion factors, respectively. These latent factors were, in turn, regressed on a second-order latent general negativity factor. Observed scores from the same occasion were allowed to correlate to account for presumed effects of various state factors, e.g., temporary mood (Figure 1D). A good fit of this model would indicate that data may have been generated by a model without any direct effects between the predictor and burnout/exhaustion, i.e., such effects may have been artifactual. In all analyses we used the lavaan default convergence criteria of 0.0001, meaning that model iterations stop when unscaled parameter values change less than 0.0001 (in absolute value)¹.

Simulations and analyses were conducted with R 4.4.0 statistical software (R Core Team, 2025) using the MASS (Venables and Ripley, 2002), lavaan (Rosseel, 2012) and osfr (Wolen et al., 2020) packages. Data and the analytic script are available at the Open Science Framework at https://osf.io/smy5n/.

The point that adjusted cross-lagged effects may be artifactual has been made before (Campbell and Kenny, 1999; Glymour et al., 2005; Eriksson and Häggström, 2014; e.g., Castro-Schilo and Grimm, 2018; Sorjonen et al., 2019; Lucas, 2023). However, the output of studies using the cross-lagged panel model, often including uncritical causal conclusions, does not seem to subside. Therefore, reiteration of this point is warranted.

We recommend fitting complementary models to data, e.g., a time-reversed model where an initial score on the outcome is regressed on a subsequent score on the outcome in addition to an initial score on the predictor. However, the time-reversed model is just as susceptible to bias and artifactual findings as the traditional cross-lagged panel model. It is also possible, despite our arguments above, for the original and the time-reversed effects to have the same sign even in the presence of true causal effects (Sorjonen et al., 2024a). Moreover, positive and negative effects of a predictor on a latent change score of an outcome may be due to influence by unmodeled state factors rather than indicating true increasing and decreasing effects, respectively. With this in mind, we would not recommend claiming a decreasing effect of depression on burnout despite a negative effect being revealed by a latent change score model (Figure 2C). Consequently, with our recommendation to analyze complementary models we do not claim that they deliver infallible evidence of causality. However, it is our conviction that considering several fallible pieces of information make for better conclusions than considering just one fallible piece of information. We do not believe that it would be preferable to consider results only from the traditional cross-lagged panel model, which is susceptible to bias and artifactual findings, than from several complementary models. As an analogy, we do not believe that it would be tenable to argue that prosecutors should interview only one witness per case, instead of several, because human perception and memory is error-prone.

We used simulated rather than empirical data as the empirical data were not available to us. Somebody might consider this as a major limitation. However, we analyzed data with four complementary structural equation models (SEM, Figure 1) and SEM uses sample sizes and covariances or correlations (for standardized parameters) as input and estimates parameter values that minimize the difference between empirical covariances/correlations and covariances/correlations predicted by the defined model. Consequently, two datasets, empirical or simulated, with the same sample sizes and covariances/correlations between variables would yield very similar results. Therefore, the points made in the present study of questionable prospective effects on burnout/exhaustion cannot be explained away by our use of simulated data.

As mentioned above, SEM and path models, including the CLPM (Figure 1A), the reversed CLPM (Figure 1B), the LCSM (Figure 1C), and our model of artifactualness (Figure 1D), try to minimize differences between predicted and empirical covariances/correlations and fit indices reflect how good models are at minimizing this difference. However, correlations do not prove causality (Reichenbach, 1971). This means that statistically significant parameter values in SEM and path models, including directed regression effects, do not prove causality if the model has been fitted on observational (i.e., non-experimental) data, not even if the model has a good fit. For example, a negative effect of initial depression on subsequent latent change in burnout in a LCSM (effects 11.1–11.3 in Table 2) does not prove that depression has a true causal decreasing effect on burnout. Moreover, several alternative models may fit the same data well. This means that a well-fitting model may be the true data generating model but it does not have to be as other alternatives are possible. For example, a good fit of our model of artifactualness (Figure 1D) suggests that data may have been generated without any direct effects between predictors and burnout/exhaustion, but it does not prove that data must have been generated without such effects.

There were some discrepancies between our estimated cross-lagged effects on subsequent burnout/exhaustion when adjusting for initial burnout/exhaustion (column “A” in Table 2) and effects reported in the original studies (column “b” in Table 1). These discrepancies were presumably mainly due to adjustment for additional variables in the original studies rather than due to differences between the empirical and our simulated data. For example, the cross-lagged effect of initial work-to-family facilitation on subsequent exhaustion reported by Innstrand et al. (2008) (b = −0.24, p < 0.05, row 4.2 in Table 1) was estimated while adjusting for initial work-to-family conflict, family-to-work conflict, and family-to-work facilitation in addition to initial exhaustion. Our result (b = −0.01, p > 0.05, row 4.2 in Table 2) suggested that if Innstrand et al. (2008) had fitted a simpler model to their data, only adjusting for initial exhaustion, they would not have found a statistically significant cross-lagged effect of initial work-to-family facilitation on subsequent exhaustion. However, this does not mean that the effect presented by Innstrand et al. (2008) therefore, should be assumed to prove a true causal decreasing effect of work-to-family facilitation on exhaustion. “The effect was truly causal” does not follow deductively from “the effect was affected by adjustment for covariates.”

Moreover, the reversed effects in column B in Table 2 tend to be stronger than the “non-reversed” effects in column A. Adjustments for covariates tend to move regression effects closer to zero (although exceptions are possible) and we see no reason to assume that this effect should be stronger on the reversed compared with the non-reversed effects. Hence, there is no reason to assume that adjustments for covariates would “correct” the discrepant findings by moving the reversed effects to the other side of zero (i.e., changing the sign of the effect from positive to negative or vice versa) while, at the same time, allowing the non-reversed effects to remain on “the right” side of zero. The same reasoning can be applied on the discrepant effects on latent change scores in column C in Table 2. Consequently, the points made in the present study of questionable prospective effects on burnout/exhaustion cannot be explained away by us not adjusting for covariates.

Our selection of studies to reanalyze was not systematic. Instead, the set of selected studies could be characterized as a convenience sample. Consequently, it is possible, although hardly likely, that most studies using cross-lagged panel models to estimate prospective effects on burnout/exhaustion would, differently from the studies reanalyzed here, withstand closer scrutiny. However, even if that would be the case, the main methodological point of the present study, that effects in cross-lagged panel models do not prove causality, would still be valid.