Bayesian Inference of Dynamic Mediation Models for Longitudinal Data

Abstract Mediation analysis is widely applied in various fields of science, such as psychology, epidemiology, and sociology. In practice, many psychological and behavioral phenomena are dynamic, and the corresponding mediation effects are expected to change over time. However, most existing mediation methods assume a static mediation effect over time, which overlooks the dynamic nature of mediation effect. To address this issue, we propose dynamic mediation models that can capture the dynamic nature of the mediation effect. Specifically, we model the path parameters of mediation models as auto-regressive (AR) processes of time that can vary over time. Additionally, we define the mediation effect under the potential outcome framework, and examine its identification and causal interpretation. Bayesian methods utilizing Gibbs sampling are adopted to estimate unknown parameters in the proposed dynamic mediation models. We further evaluate our proposed models and methods through extensive simulations and illustrate their application through a real data application.


Introduction
Mediation analysis is a family of statistical methods to investigate whether the effect of a predictor variable (X) on an outcome variable (Y) can be explained by a third variable (called a mediator M).It has wide applications in a variety of scientific disciplines, including but not limited to psychology (Bharadwaj & Bezborah, 2021;Vuorre & Bolger, 2018), prevention programs (Fang & Schinke, 2014), and epidemiology (Richiardi et al., 2013).As shown in Figure 1, a mediator (M) is a variable that transmits the causal effect of a predictor variable (X) to an outcome variable (Y).To test mediation effect, researchers typically employ the causal steps approach (Baron & Kenny, 1986;Judd & Kenny, 1981), which conducts a series of separate significance tests on the relations among X, M, and Y.In addition to the causal steps approach, several other conventional mediation methods (Freedman & Schatzkin, 1992;Sobel, 1982) have been proposed to estimate and test the mediation effect.A recent development in mediation analysis is the causal mediation method (Imai et al., 2010(Imai et al., , 2010;;), which provides definitions and estimators of direct and indirect effects based on the potential outcome framework (Pearl, 2005;Rubin, 1974).This method has been expanded to more complex issues, such as categorical data (Albert & Nelson, 2011) and multiple mediators (VanderWeele & Vansteelandt, 2014;Vansteelandt & Daniel, 2017).However, most of these methods focus on the mediation effect at a specific time, while the mediation analysis for longitudinal data remains relatively uncommon.
With the rapid advancement of data collection technology, longitudinal data collection is more and more common.In order to analyze such data, various longitudinal mediation models have been proposed, including the autoregressive (AR) model (Cole & Maxwell, 2003;Wu et al., 2018), the latent growth curve model (Cheong et al., 2003;Sullivan et al., 2021), and the latent difference score model (Simone & Lockhart, 2019).The auto-regressive model is based on structural equation modeling (SEM), and assumes that each variable in the model depends causally not only on prior variables at the same time but also on the same variable at previous time points.The latent growth curve model uses the latent growth factors of M and Y (e.g., intercepts and slopes) to represent the roles of variables M and Y, respectively, in mediation model.It investigates how X affects the growth of M, which subsequently affects the growth of Y.The latent difference score model uses latent variables to represent the difference between adjacent measurements of variables, and examines how the difference in X affects the difference of M, which in turn affects the difference in Y.
Despite the fast growing literature on longitudinal mediation analysis, most of the existing methods focus on how the variables (X, M, and Y) vary over time, but overlook how the mediation effect varies over time.However, in many psychological and behavioral activities, the mediation effect is expected to change over time.Our motivating example is a psychological research by Hardy et al. (2014), which investigated the mediation effect of daily spiritual experiences (M) on the relationship between daily religious activities (X) and daily moral emotions (Y) over 50 days.Because Hardy et al. (2014) found that there are no significant lagged effects among these variables, we adopted a concurrent association method to investigate the mediation effects at different time points.Similar methods were used in Blood andCheng (2011, 2012), where they estimated time-specific time effects at each time point, and applied the static mediation method to each time point separately.However, the estimates of mediation effects in Blood andCheng (2011, 2012) were inefficient with large variances due to a lack of consideration for the relationship across different time points.To address this issue, Huang and Yuan (2017) proposed a dynamic mediation analysis method that models and estimates the dynamic mediation effects using a non-parametric penalized spline approach.However, such smoothing method may not be effective for non-smooth dynamic trends, such as the empirical dynamic trend shown in Figure 2, based on the mediation effects of daily spiritual experiences (Hardy et al., 2014).
Inspired by the non-smooth dynamic trend illustrated in Figure 2, we propose a new dynamic mediation analysis method, which models the path parameters of the mediation process as AR processes over time.To be specific, we propose three types of models with common or individual path parameters to capture the dynamic trend of mediation effect.Our method offers researchers a different method to the penalized spline approach (Huang & Yuan, 2017), and is able to capture the non-smooth dynamic trend of mediation effect in longitudinal data.Moreover, compared to the penalized spline approach, the AR processes method is easier to comprehend and interpret in longitudinal data analysis.In order to define the mediation effects in our models, we utilize the potential outcome framework and discuss the definition, identification assumptions and causal interpretation of mediation effects, while there is no identification guarantee for the penalized spline approach.A Bayesian inference method utilizing Gibbs sampling is adopted to estimate model parameters, which is conceptually simple to make inferences for our complex models.Our simulation and real data results demonstrate the efficiency of the proposed method in capturing the dynamic nature of the mediation effect.
In the subsequent sections of this paper, we present our dynamic mediation models in Section 2. Section 3 defines the mediation effect and examines its identification and causal interpretation.Section 4 provides specific steps of Bayesian inference for our models.Sections 5 and 6 demonstrate the validity and advantages of our method through extensive simulations and a real data application.Section 7 concludes the paper with a discussion of extensions and future directions.

Dynamic Mediation Models
Let Y it , M it , and X it , respectively, denote the outcome, mediator, and predictor at time t ¼ 1, . . ., T for individuals i ¼ 1, . . ., N: We now propose three dynamic mediation models for longitudinal data.The first model assumes homogeneity in path parameters across individuals, resulting in shared non-individual path parameters.Therefore, common causal effects can be inferred, but individual causal effects cannot.The second and third models, on the other hand, allow for heterogeneity in path parameters among individuals.Thus, individual causal effects can be inferred, and the effect due to within-individual variation can be separated from the effect due to between-individuals variation.These models offer a valuable framework for examining mediation effects in longitudinal data.

Dynamic Mediation Model with Non-Individual Path Parameters
The first model is referred to as a dynamic mediation model with non-individual path parameters.This model allows the path parameters between the independent variable X and the mediator variable M, between the mediator variable M and the outcome variable Y, and between the independent variable X and the outcome variable Y, to change over time t.Below, we model the change trends of path parameters as first-order AR (i.e., AR(1)) processes: (1) where e 1it and e 2it are normal random errors: e 1it $ Nð0, r 2 1 Þ and e 2it $ Nð0, r 2 2 Þ, and a 0 , a 1 , b 0 , b 1 , c 0 , and c 1 denote the auto-regressive parameters, and v jt $ Nð0, d 2 j Þ for j ¼ 1, 2, 3 and t ¼ 1, . . ., T: The random errors fe 1it , e 2it , v 1t , v 2t , v 3t : i ¼ 1, . . ., N; t ¼ 1, . . ., Tg are independent.The corresponding causal diagram is presented in Figure 3.The AR(1) processes in model ( 1) can be easily extended to the p-order AR (i.e., AR(p)) processes.

Dynamic Mediation Model with Individual Path Parameters
The second model, called the dynamic mediation model with individual path parameters, allows the individual difference in path parameters: (2) where A it , B it , and C it are both individual-specific and timespecific path parameters, and e 1it , e 2it , e a, i , e b, i , and e c, i are independent random errors with e 1it $ Nð0, The time-specific path parameters a t , b t , and c t shared by all individuals are assumed to follow the AR(1) processes

Dynamic Mediation Model with Individual AR Processes
The third model is called the dynamic mediation model with individual AR processes, which assumes weaker associations among different individuals compared to the second model.In this model, the path parameters of different individuals follow different AR processes, while the associations among different individuals lies in the fact that the autoregressive parameters of these different AR processes are the same: where e 1it and e 2it are independent normal random errors: e 1it $ Nð0, r 2 1 Þ and e 2it $ Nð0, r 2 2 Þ: a 0 , a 1 , b 0 , b 1 , c 0 and c 1 are AR parameters, and v jit $ Nð0, d 2 j Þ are independent random errors for j ¼ 1, 2, 3, i ¼ 1, . . ., N, and t ¼ 1, . . ., T: Moreover, it is worth noting that the AR(1) processes in model (4) can also be extended to AR(p) processes.
This model comprises numerous effective parameters a it , b it , and c it , which are difficult to accurately estimate due to the weak associations among individuals.However, since the auto-regressive parameters a 0 , a 1 , b 0 , b 1 , c 0 , and c 1 are individual and time independent, our method is able to accurately estimate them.This is verified by the simulation results reported in Section 5.1.The AR process property indicates that these auto-regressive parameters determine the value of time-averaged causal mediation effect, as detailed in Theorem 3.1.Therefore, our method can make valid inference for the time-averaged causal mediation effect.
It is worth noting that although the path parameters are estimated using the concurrent data of the mediation model variables, the lag relationships are also modeled in our three dynamic mediation models through the auto-regressive association of the path parameters.Potentially, other configurations of the lag relationship models can be used, depending on the specific research questions to answer.In this paper, our focus is on modeling the auto-regressive structure of the mediation effects.

Definition, Identification, and Causal Interpretation of Mediation Effects
In this section, we first use the potential outcome framework (Pearl, 2005;Rubin, 1974) to define mediation effects for the three models we proposed, then discuss the associated identifiability problem, and provide a causal interpretation of the defined mediation effects.In the potential outcome framework, M it ðxÞ denote the potential value of a mediator at time t for individual i if the predictor X it is set to be x, and Y it ðx, mÞ denote the potential value of the outcome at time t for individual i if the predictor X it is set to be x and the mediator M it is set to be m.Moreover, we use X , M, and Y to denote the supports of the distributions of X it , M it , and Y it , respectively.It is worth noting that our models in Section 2 have not considered the existence of baseline covariates, but in this section we consider q baseline covariates Z i ¼ ðZ 1i , . . ., Z qi Þ T for individuals i ¼ 1, . . ., N to expand the application range of our identification theory.Specially, the number of baseline covariates q is zero in our models, and Z is used to represent the support of the distribution for Z i : Then, we can define the causal mediation effect (CME) for individual i at time t as the average causal mediation effect (ACME) at time t as and the time-averaged causal mediation effect (TCME) for individual i as Following the work by Imai et al. (2010); Imai et al. (2010), we propose the following assumptions to ensure the identification of mediation effect.
Assumption 1 The conditional probability density functions f ðX it ¼ xjZ i ¼ zÞ > 0 for all z 2 Z and x 2 X, and Assumption 2 X it ??Y it ðx, mÞ Z i ¼ z j for all z 2 Z, x 2 X and m 2 M, i.e., no unmeasured confounders between the predictor and the outcome; Assumption 3 X it ??M it ðxÞjZ i ¼ z for all z 2 Z and x 2 X, i.e., no unmeasured confounders between the predictor and the mediators; x 2 X and m 2 M, i.e., no unmeasured confounders between the mediators and the outcome; x 0 2 X and m 2 M, i.e., no exposure-induced confounders between the mediators and the outcome.
With these assumptions, we can derive the following identification, and causal interpretation for mediation effect: Theorem 3.1 (Identification, and causal interpretation for mediation effect) Under Assumptions 1-5, the average causal mediation effect ACME t ðx, x 0 Þ and the time-averaged causal mediation effect TCME i ðx, x 0 Þ are nonparametrically identifiable for any x 2 X.For the dynamic mediation model with non-individual path parameters in Equation ( 1), ACME t ðx, x 0 Þ ¼ a t b t ðx 0 À xÞ; for the dynamic mediation model with individual path parameters in Equations ( 2) and (3), ACME t ðx, x 0 Þ ¼ a t b t ðx 0 À xÞ; and for the dynamic mediation model with individual AR processes in Equation (4), TCME i ðx, x 0 Þ ¼ a 0 1Àa 1 b 0 1Àb 1 ðx 0 À xÞ: Theorem 3.1 ensures the identification of mediation effect and gives the causal interpretation for our three dynamic mediation models, respectively.The proof of Theorem 3.1 is provided in Online Supplementary Material -Appendix A.

Bayesian Inference
In this section, we develop a statistical inference method for model parameters and mediation effects.The inference for our complex dynamic mediation models poses challenges when adopting a frequentist perspective.In contrast, the Bayesian approach offers several notable advantages.First, the frequentist mediation method encounters substantial computational and estimation difficulties when applied to complex mediation models (Kenny et al., 2003).By contrast, Bayesian methods offer a promising alternative, circumventing such complexities.Second, the construction of confidence intervals in frequentist methods relies on the asymptotic normal distribution of the mediation effect under a large sample size, but this distribution is usually skewed even in the simple mediation model (MacKinnon et al., 2002).To solve this issue, some researchers utilized the bootstrap method to construct confidence intervals.However, the bootstrap method is computationally expensive for our complex multilevel models, and how to conduct bootstrap with our data structure is also tricky.From the Bayesian perspective, the estimation and construction of credible intervals are conceptually simple and straightforward based on the posterior samples.Note that credible intervals in Bayesian methods are analogous to confidence intervals in frequentist methods, although they differ on definitions.Third, Bayesian methods allow for the incorporation of prior information into the inference process, potentially enhancing efficiency (Yuan & MacKinnon, 2009).Therefore, we choose the Bayesian method, which has been used by Yuan and MacKinnon (2009) in simple mediation analysis, for our inference.
To apply the Bayesian method to our models, the first and foremost step is to specify the prior distributions of the unknown parameters.In this study, we focus on the use of uninformative priors.Taking the AR(1) dynamic mediation model with non-individual path parameters as an example, the following independent vague priors are assigned to the unknown model parameters fa where Nð0, 10 3 Þ denotes the normal distribution with mean 0 and variance 10 3 , UðÀ1, 1Þ denotes the uniform distribution on the interval ðÀ1, 1Þ, and IGð10 À3 , 10 À3 Þ denotes the inverse gamma distribution with shape parameter of 10 À3 and scale parameter of 10 À3 : If one has prior information from pilot studies or expert opinions, informative priors can be used to improve the efficiency of Bayesian estimates.
With the above priors, the joint posterior distribution of the model parameters is where represents the joint distribution of data ðY, M, XÞ given the parameters, corresponding to the proposed dynamic mediation models.Note that for the simplicity of presentation, we did not include the random-effects in the formulation.It is difficult to calculate the closed form of the joint posterior distribution function.Therefore, we adopt the Gibbs sampling method to obtain posterior samples of the model parameters, which can be implemented conveniently through the software OpenBUGS (Lunn et al., 2013).Based on the posterior samples, we then obtain the estimates (posterior means) and credible intervals of the model parameters and mediation effects.
To be specific, OpenBUGS is a powerful and flexible software that utilizes Markov chain Monte Carlo (MCMC) methods to conduct Bayesian analysis for complex statistical models.The R2OpenBUGS package allows users to run OpenBUGS from R easily.The first step in OpenBUGS is to specify the statistical model of interest, including the likelihood function and the prior distributions of parameters.The corresponding code is presented in Online Supplementary Material -Appendix B. With the specified model, OpenBUGS automatically generates multiple Markov chains from the posterior distribution.To assess the convergence of these Markov chains, we can first examine their trace plots.If the trace plots show that the parameter values have stabilized and do not show any significant trends or unexpected patterns, it suggests that convergence has been achieved.We can also use the Gelman-Rubin diagnostic statistics (Gelman & Rubin, 1992) to assess the convergence, with a value less than 1.1 indicating the chains have reached convergence.It should be noted that in simulations, it is hard to examine the trace plots for each repeated simulation, and hence, the Gelman-Rubin diagnostic statistics is often monitored as in the current study.In the subsequent simulations and real data analysis, we set the number of Markov chains as 3, the number of retained draws as 10,000, and the number of burn-in draws as 4000.Moreover, we found that in our dynamic mediation models, the thinning parameter had little impact on the inference results.Therefore, to save computation time, we set the thinning parameter as 1 (i.e., no thinning).
The method described above relies on the specification of order (p) of the models.We adopt the deviance information criterion (DIC) to choose an appropriate order p. DIC was first proposed by Spiegelhalter et al. (2002) and widely used in Bayesian model comparison and selection (e.g., Cain & Zhang, 2019).Specifically, DIC is a Bayesian generalization of the Akaike information criterion (AIC) and Bayesian information criterion (BIC), and OpenBUGS can calculate DIC directly for our models.Like AIC and BIC, models with smaller DIC values are preferred over models with larger DIC values.Moreover, how to deal with missing data is important in practice.In the Bayesian framework, missing data can be regarded as unknown parameters and sampled in Gibbs sampling, which is the default option in OpenBUGS.Therefore, our inference method can handle missing data.

Simulation Studies
In this section, we conducted three simulation studies to examine the performance of Bayesian inference for our dynamic mediation models.We first investigated the performance of the three dynamic mediation models in Section 5.1, then we compared the performance of our method with three existing methods in Section 5.2.Finally, in Section 5.3, we investigated the performance of our dynamic mediation models under missing data.All three simulation studies were conducted using R and OpenBUGS, and the code for simulation can be found in Online Supplementary Material -Appendix B. Simulations were repeated 1000 times for each simulation condition.

Performance of the Proposed Method
This simulation study was designed to evaluate the performance of the Bayesian estimation method for the three proposed models.First, we investigated the performance for the AR(1) dynamic mediation model with non-individual path parameters.To be specific, X it was generated from the standard normal distribution N ð0, 1Þ and M it and Y it were generated according to the model described in Equation (1).The following parameters were kept the same in all conditions: and r 2 2 were considered, which correspond to R 2 values ranging from 0.1 through 0.9 (See Online Supplementary Material -Appendix C and Table S.1).Moreover, various combinations of sample size N and number of measurements T were considered.
The evaluation metrics included estimation bias (average estimation bias), mean squared error (MSE) and coverage probability (CP) of 95% credible interval (CI) for parameters and d 2 3 : To be specific, the 95% CIs were constructed using the 2.5% and 97.5% quantiles of the posterior samples, and the CP of the 95% CI was approximated by the proportion of the CIs that cover the true parameters across 1000 simulations.The results for N ¼ 100 and T ¼ 100 are summarized in Table 1, and the results for other combinations of N and T are provided in Online Supplementary Material -Appendix D. The Bayesian inference algorithm was hard to converge when R 2 ¼ 0:1, so we only reported the results for R 2 ¼ 0:9, 0.7, 0.5, and 0.3.To evaluate the estimation performance of mediation effects, we calculated the estimation bias, MSE and CP of 95% CI for ACME t ðx, x þ 1Þ ¼ a t b t at time points t ¼ 10, 20, . . ., 100 (Table 2).
It can be concluded from Tables 1 and 2 that our method performs well for large R 2 (! 0:5), with small estimation biases and MSEs.The coverage probabilities are close to the nominal value 95%.As the value of R 2 decreases, the estimation biases and MSEs increase.This is expected since the R 2 value is a statistical measure of how well the outcomes are explained by our model, that is, a small value of R 2 corresponds to a small proportion of the variance explained by the predictors.
Second, we considered the AR(1) dynamic mediation model with individual path parameters in Equations ( 2) and (3).The following parameters were kept the same in all conditions: and d 2 3 ¼ 1: Various combinations of r 2 1 and r 2 2 were considered for a range of R 2 values (Table S.1 in Online Supplementary Material).Tables S.17-S.19 in Online Supplementary Material -Appendix D displays the estimation bias, MSE and CP of 95% credible interval for model parameters a Results for other combinations of N and T are also presented in Online Supplementary Material -Appendix D. Evidently, our method performs well with a large value of R 2 (! 0:5).To be specific, our method possesses small estimation biases and MSEs, and good coverage probabilities that are close to the nominal value.Moreover, similar to the above model, the estimation biases and MSEs decrease with the increase of R 2 .Third, we investigated the AR(1) dynamic mediation model with individual AR processes in Equation ( 4).The following parameters were kept the same in all conditions: and d 2 3 ¼ 1: Due to the complexity of this model, the Bayesian inference algorithm was hard to converge with small values of R 2 , so only r 2 1 and r 2 2 corresponding to R 2 ¼ 0:9 was considered.Presented in Tables S.26 and S S.26 in Online Supplementary Material that our method performs well for the estimation of model parameters, with small biases and MSEs.However, the estimates of CMEs have large MSEs (Table S.27 in Online Supplementary Material), indicating that our method is unable to accurately estimate the individual-specific CME it ðx, x þ 1Þ ¼ a it b it : This is not surprising since this model comprises a large number of effective parameters so that the variances of the estimates could be very large.On the other hand, the model parameters a 0 , a 1 , b 0 , b 1 , c 0 , and c 1 determine the time-averaged causal mediation effect TCME i ðx, x þ 1Þ, refer to Theorem 3.1.Therefore, under our model with weak associations among individuals, our method can make valid inference for the time-averaged causal mediation effect.
Based on the simulation results, we can make some recommendations for the choice of the proposed three models.When we have any prior information that the individual path parameters are not differential, we recommend the use of the first model with non-individual path parameters.With prior knowledge of slight individual difference, the second model with individual path parameters is suggested.If we have prior information on the existence of large individual difference, the third model with individual AR process is recommended.It is worth noting that under the third model, our method can only make valid inference for the time-averaged causal mediation effects In this study, we employ the DIC as a model selection tool to identify the most appropriate model.To evaluate the performance of DIC, we conducted a simulation study.The data were generated based on the AR(1) dynamic mediation model with non-individual path parameters with N ¼ 50 and T ¼ 50.Specifically, the parameter values, namely a 0 , a 1 , b 0 , b 1 , c 0 , c 1 , d 2 1 , d 2 2 , d 2 3 , were set consistently with those utilized in the above simulations.Additionally, we considered variances r 2 1 ¼ 0:3 and r 2 2 ¼ 1:8, resulting in an R 2 value of 0.9.After the data generation, we applied the three proposed models to the simulated dataset and selected the model with the lowest DIC value.The process was repeated 1000 times, and the frequencies of model selection were recorded in Table 3.The results reveal that the true model (model-1-AR(1)) obtained the highest number of selections (926 in 1000 simulations), thereby validating the efficacy of DIC as a model selection criterion.
Furthermore, the simulation results presented above reveal that our Bayesian inference method encountered challenges in achieving convergence under certain configurations with low R 2 values.For instance, in the AR(1) dynamic mediation model with non-individual path parameters, when R 2 ¼ 0:1, we conducted 1000 simulations and observed that only 20 simulations achieved convergence.Despite increasing the number of retained draws and burnin draws, the convergence performance remained unsatisfactory.Subsequently, focusing on the 20 convergent simulations, we assessed the performance metrics, as depicted in Tables 4 and 5. Notably, the results show that even when convergence was attained with R 2 ¼ 0:1, the inference performance remained poor, particularly for parameters c 0 , c 1 , and d 2 3 : In order to address the non-convergence issue, we investigated weak-informative priors as potential solutions.Specifically, the following weak-informative priors were considered: Through employing these weak-informative priors, convergence was achieved in 823 out of 1000 simulations, and inference performance also improved, as demonstrated in Tables 4 and 5. Therefore, when prior information regarding the model parameters is available, we recommend the utilization of weak-informative priors in practical applications.

Performance Comparison with Existing Methods
Simulations were conducted to compare the performance of our method and three existing methods.Since most existing methods assume that individuals are not differential, we focused on our AR(1) dynamic mediation model with nonindividual path parameters, and the corresponding method is abbreviated as DM.Three competing methods were considered, namely, the Bayesian static mediation method (SM, Yuan & MacKinnon, 2009), the relaxed Bayesian mediation method (RM, Yuan & MacKinnon, 2009), and the non-parametric Bayesian dynamic mediation analysis method (NDM, Huang & Yuan, 2017).The SM method is a single-level Bayesian mediation analysis method, which assumes that the mediation effect is static over time.The RM method assumes independent time-specific path parameters and uses the SM method to estimate the mediation effects at each time point separately.The NDM method assumes that the mediation effect is a smooth function of time, and adopts the penalized spline approach to estimate the smooth function.
We generated data under two different models.First, we generated data under the AR(1) dynamic mediation model with non-individual path parameters.X it was generated from the standard normal distribution N ð0, 1Þ, and all parameters were kept the same across conditions: 3, and r 2 2 ¼ 1:8, which correspond to R 2 ¼ 0:9: Under this setting, the true ACME t ðx, x þ 1Þ ¼ a t b t is time varying with large variations (Figure 4).Various sample sizes (N) and numbers of measurement occasions (T) were considered.As shown in Table 6, the DM method outperforms the other methods in estimating the mediation effect, which has minimal estimation biases and MSEs, and coverage probabilities close to the nominal value 95%.Note that the SM method assumes that the mediation effects are static but the true values are time varying to a large extent, so it is not surprising that its estimation performance is poor, with the largest MSEs and lowest coverage probabilities.Since the smooth method NDM cannot accurately capture the dynamic nature of AR processes, it is not surprising that the corresponding MSEs were large and the coverage probabilities were low.However, since it assumes dynamic mediation effects, the NDM method outperformed the SM method.The RM method assumes independent time-specific parameters, so it requires a large sample size N to accurately estimate the parameters.Moreover, the RM method ignores the correlation among time-specific parameters, so it is intuitive that no information from other time points can be used to estimate each time-specific parameter.The DM method resolves this issue by modeling the parameters as AR process over time, and it outperformed the RM method overall in terms of bias and MSE.
Next, we generated data under the simulation settings of Huang and Yuan (2017) to evaluate the robustness of our method.Huang and Yuan (2017) considered four different shapes of the mediation effect: one static mediation effect, and three types of smooth mediation effects.The results are presented in Tables S.28-S.31 in Online Supplementary Material -Appendix E. Under the static mediation effect, the SM method performed the best among all considered methods, which is reasonable since the SM method assumes the static mediation effect.The DM method outperformed the NDM and RM methods in terms of bias and MSE.Under the three types of smooth mediation effect, it is not surprising that the NDM method outperformed the other methods, while the DM method outperformed the SM and RM methods in terms of bias and MSE.Note that the DM method had inflated coverage probabilities in all situations, which is not surprising since the data were not generated STRUCTURAL EQUATION MODELING: A MULTIDISCIPLINARY JOURNAL from the model underlying the DM method.In conclusion, with data generated from the flat or smooth mediation effect settings, our DM method still performed well in terms of bias and MSE, albeit to inflated coverage probabilities.

Dealing with Missing Data
Missing data are common in longitudinal data analysis.Therefore, we conducted an additional study to evaluate the performance of our method in the presence of missing data.We focused on the AR(1) dynamic mediation model with non-individual path parameters.First, we used the data generated in Section 5.1 as complete data.Then, we randomly deleted 10% of X, M, Y to simulate missing data under the missing completely at random (MCAR) mechanism.After that, we analyzed the data using our Bayesian method.The results summarized in Tables 7 and 8 show that the missingness had a minor impact on the estimation results 3 under the AR(1) dynamic mediation model with non-individual path parameters with data missing rate ¼ 0.1.compared to no missing data situation.We also considered the other two proposed models.The corresponding results are similar and are not reported.

A Real Data Application
We applied our proposed method to the real data set analyzed in Hardy et al. (2014).With this data set, Hardy et al. (2014) used a static mediation model to examine relations among intra-individual variability in daily religious activities, daily spiritual experiences, and daily moral emotions.They found that the daily spiritual experiences mediate relations between daily religious activities and daily moral emotions.The original sample included 142 individuals between 18 and 69 years old, who completed daily online surveys for up to 50 days on their religious activities, spiritual experiences, and moral emotions.
We removed those participants who missed over 10 days of participation, resulting in 94 participants for the current study.
Based on the previous analysis, we chose daily religious activities as the independent variable (i.e., X), daily spiritual experiences as the mediator variable (i.e., M), and daily empathy as the outcome variable (i.e., Y).Then, we applied the three proposed dynamic mediation models to the real data.The DICs corresponding to various orders are present in Table 9.The orders corresponding to the smallest DICs (i.e.,37110,36490,and 35860) for the three models are 1, 1, and 3, respectively, and the corresponding models are denoted as model-1-AR(1), model-2-AR(1), and model-3-AR(3), respectively.Of these three selected models, model-3-AR(3) has the smallest DIC.However, it should be pointed out that model-3-AR(3) has much more effective parameters than the other two selected models (1050 vs. 72 and 258), potentially leading to over-parameterization.Furthermore, model-2-AR(2), model-2-AR(3), and model-3-AR(1-3) failed to achieve convergence when applied to the real data.Nevertheless, despite these considerations, the subsequent discussion encompasses the results of all three selected models.The analysis results by the SM, NDM, and RM methods are also included for comparison.
Figure 5 displays the estimated ACMEs for all of the six methods.Since the true ACME values are unknown for the real data set, there is no ground truth.As the RM method specifies independent time-specific parameters, the corresponding estimated ACME values are used as the baseline for comparison.Notably, the ACMEs in model-3 cannot be directly estimated due to the unique model structure.In this scenario, we employ the average values of the CMEs for the first 10 individuals as proxies for ACME estimates, although caution should be exercised regarding their reliability.As shown in Figure 5, the ACME estimates of our dynamic mediation models (model-1-AR(1) and model-2-AR(1)) share the same dynamic patterns as that of the RM method, while the SM method cannot capture the dynamic trend of ACME.Furthermore, the smooth change pattern of the NDM method does not match that of the RM method.
Figure 6 shows the 95% CI of the estimated ACMEs by model-1-AR(1).Evidently, the mediation effect is significant at some but not all measurement times, again illustrating the importance of evaluating the dynamics of the mediation effect.Focusing on the mediation effect at a specific time, our method is capable of detecting significant mediation effects.Among the three selected models, model-2-AR(1) and model-3-AR(3) allow the analysis of individual path parameters.The CMEs for the first 10 individuals estimated by model-2-AR(1) and model-3-AR(3) are displayed in Figure 7.Under model-2-AR(1), the estimated individual dynamic trends of the 10 individuals are similar (Figure 7a), while under model-3-AR(3), the 10 individuals have quite different estimated dynamic trends (Figure 7b).This is reasonable since compared to model-2-AR(1), model-3-AR(3) assumes weaker associations of path parameters among different individuals.In conclusion, our method reveals a dynamic nature of mediation effect in this real data.
Moreover, the estimators of R 2 values for M and Y can be obtained through the estimators of r 2 1 and r 2 2 , which are presented in Table 10.It is noteworthy that the estimated R 2 values are relatively modest, potentially attributable to the omission of potential covariates and confounders in our models.This omission results in the variance contribution of these variables being erroneously ascribed to the residual error.Consequently, future investigations should account for these factors, and explore their influence on the results in this real dataset.

Discussion
In this article, three dynamic mediation models are proposed to capture the dynamic nature of the mediation effect in longitudinal data, in contrast to most existing methods that assume a time-invariant mediation While the non-parametric Bayesian dynamic mediation analysis method (Huang & Yuan, 2017) models the path parameters as smooth functions of time, the proposed method models the path parameters as AR processes.Compared to the relaxed Bayesian mediation (RM) method that specifies independent time-specific parameters, our method utilizes the relationship across the time-specific parameters at different time points.Additionally, under our proposed models, the dynamic mediation effect can be easily defined under the potential outcome framework, which is identifiable and has a causal interpretation under certain conditions.A Bayesian inference method is developed to estimate unknown parameters in our proposed models.Compared to the frequentist methods for mediation inference, our Bayesian method enables straightforward construction of credible intervals for model parameters and mediation effects based on posterior distribution samples.Furthermore, our Bayesian method enables researchers to utilize prior information to improve inference efficiency.
Three simulation studies were conducted to examine the desired performance of the proposed method.The results demonstrated that compared to existing methods, the proposed method exhibits superior performance in estimating dynamic mediation effects.Moreover, our method can deal with missing completely random data reasonably well with a    slight efficiency loss.In the application to the real data, our method reveals the dynamic nature of mediation effect in this real data, illustrating the necessity of the dynamic mediation analysis.
The proposed method in this study is subject to certain limitations.First, DIC, which is utilized to choose an appropriate model, might not be effective (Spiegelhalter et al., 2014).The DIC values of models with different AR(p) processes are observed to be close to each other in our real data application, posing a challenge in selecting the most suitable model.Hence, it deserves further investigation to develop a new model selection criterion for our models.Second, our proposed method is currently designed for continuous independent variable, mediator, and outcome, and it is of interest to generalize our method to other types of variables (binary, survival, categorical, and so on) in the future.Third, our analysis of missing data is conducted from a Bayesian perspective, where missing data is treated as unknown parameters and estimated accordingly.This analysis strategy is valid only when the data is missing completely at random.Advanced strategies such as imputation or inverse probability weighting may be more appropriate for other types of missing data, warranting further investigation (Little & Rubin, 2019).Fourth, in this paper, we focused on the auto-regressive relationship for the mediation effect directly, but ignore the directly modeling of the lag effects of variables (due to the research question we focus).To address this limitation, it is suggested that future research could explore the generalization of our method to the auto-regressive mediation model (Cole & Maxwell, 2003;Wu et al., 2018), which considered the lag effects among variables through SEM.Finally, our models ignore the potential covariates and confounders that could influence the causal inference results.To mitigate this limitation, future studies should consider such factors and evaluate the sensitivity of the inference.

Funding
The work of SZ and HZ was supported by the National Natural Science Foundation of China [No. 72091212, 11771096].The work of ZZ was supported by a grant from the US Department of Education [R305D210023].However, the contents of paper do not necessarily represent the policy of the US Department of Education, and you should not assume endorsement by the Federal Government.

Figure 1 .
Figure 1.The causal diagram of a simple static mediation model.

Figure 2 .
Figure2.The time plot of the estimated mediation effects of the concurrent association method(Blood & Cheng, 2011, 2012).

Figure 3 .
Figure 3.The causal diagram of the AR(1) dynamic mediation model with non-individual path parameters.
.27 in Online Supplementary Material -Appendix D are the estimation bias, MSE and CP of 95% CI for model parameters a 0 , a 1 , b 0 , b 1 , c 0 , c 1 , d 2 1 , d 2 2 , d 2 3 , and CMEs of the first individual: CME 1t ðx, x þ 1Þ ¼ a 1t b 1t , with N ¼ 100 and T ¼ 100.Results for other combinations of N and T are presented in Online Supplementary Material -Appendix D. It can be concluded from Table

Figure 4 .
Figure 4.The time plot of the true ACME in one replication of simulation 5.2 with T ¼ 100.

Figure 6 .
Figure 6.The estimated ACMEs of our AR(1) dynamic mediation model with non-individual path parameters, and the shade shows the 95% credible interval.

Figure 7 .
Figure 7.The estimated CME it ðx, x þ 1Þ ¼ a it b it of different participants.(a) Dynamic mediation AR(1) model with individual differences; (b) dynamic mediation AR(3) model with individual AR processes.Due to the large number of participants, we only demonstrate the results of the first 10 participants, and the lines of different colors correspond to different participants.
2 M , R 2 value for M it ; R 2 Y , R 2 value for Y it ; model-1, dynamic mediation model with non-individual path parameters; model-2, dynamic mediation model with individual path parameters; model-3, dynamic mediation model with individual AR processes.

Table 2 .
Mean of absolute true value (MAT), bias, MSE (Â10 À2 ), and CP of the 95% CI (%) for the ACME t ðx, x þ 1Þ ¼ a t b t at various time points t 0 s under the AR(1) dynamic mediation model with non-individual path parameters.

Table 3 .
The selection frequencies of various dynamic mediation models with different AR(p) processes in 1000 simulations.Note.Model-1, dynamic mediation model with non-individual path parameters; model-2, dynamic mediation model with individual path parameters; model-3, dynamic mediation model with individual AR processes.

Table 5 .
Mean of absolute true value (MAT), bias, MSE (Â10 À2 ), and CP of the 95% CI (%) for the ACME t ðx, x þ 1Þ ¼ a t b t at various time points t 0 s under the AR(1) dynamic mediation model with non-individual path parameters (R 2 ¼ 0:1).

Table 6 .
Mean of absolute true value (MAT), bias, MSE, and CP of the 95% CI (%) of the ACME t ðx, x þ 1Þ ¼ a t b t at various time points t 0 s, under four different methods.Note.SM, Bayesian staticmediation analysis method; NDM, non-parametric Bayesian dynamic mediation analysis method; RM, relaxed Bayesian mediation analysis method; DM, AR(1) dynamic mediation method with non-individual path parameters.

Table 8 .
Mean of absolute true value (MAT), bias, MSE (Â10 À2 ), and CP of the 95% CI (%) for the ACME t ðx, x þ 1Þ ¼ a t b t at various time points t 0 s under the AR(1) dynamic mediation model with non-individual path parameters with data missing rate ¼ 0.1.

Table 10 .
The estimated R 2 values for different models.