Estimating Latent Baseline-by-Treatment Interactions in Statistical Mediation Analysis

Abstract Statistical mediation analysis is used to uncover intermediate variables, known as mediators [M], that explain how a treatment [X] changes an outcome [Y]. Often, researchers examine whether baseline levels of M and Y moderate the effect of X on posttest M or Y. However, there is limited guidance on how to estimate baseline-by-treatment interaction (BTI) effects when M and Y are latent variables, which entails the estimation of latent interaction effects. In this paper, we discuss two general approaches for estimating latent BTI effects in mediation analysis: using structural models or scoring latent variables prior to estimating observed BTIs and correcting for unreliability. We present simulation results describing bias, power, type 1 error rates, and interval coverage of the latent BTIs and mediated effects estimated using these approaches. These methods are also illustrated with an applied example. R and Mplus syntax are provided to facilitate the implementation of these approaches.

In the social sciences, statistical mediation analysis is used to investigate the intermediate variables, known as mediators (M), by which a treatment (X) led to a change in an outcome (MacKinnon, 2008).Recent examples of mediation hypotheses include a program on adaptive regulation strategies (X) decreasing externalizing problems (Y) by improving emotion regulation (M) in adolescents (te Brinke et al., 2021), and an imagery rescripting treatment (X) decreasing nightmare distress (Y) by increasing self-efficacy (M) in individuals with nightmare disorder (Kunze et al., 2019).Often, researchers include baseline levels of M and Y in the mediation model (i.e., measures prior to the treatment or intervention) to account for confounding of the M to Y relation at posttest (Valente et al., 2020).Furthermore, with the inclusion of baseline measures, researchers could also test whether the effect of the program depends on the baseline levels of M and Y (Baron & Kenny, 1986;Morgan-Lopez & MacKinnon, 2006), which helps inform future implementation of prevention programs (MacKinnon, 2008).For example, in a randomized prevention program (X) to decrease the use of anabolic steroids (Y) in high school athletes by increasing their knowledge of the effects of steroids on health (M), the program was shown to be most effective for athletes who had high intentions of using steroids at baseline (Fritz et al., 2005).As such, baseline levels of M and Y could moderate the effect of the program on the mediator or the outcome, and this moderation can be assessed by including baseline-by-treatment interactions in the model (Morgan-Lopez & MacKinnon, 2006).
Previous studies suggest that mediation effects can be accurately estimated in single mediator models with one latent interaction using methods cited above (e.g., Cheung & Lau, 2017;Gonzalez & Valente, 2022).However, it is unclear if those same findings hold in two-wave mediation models with latent baseline-by-treatment interaction (BTI) effects.Estimating interaction effects in mediation models with two waves of data is more complex than in mediation models with a single wave of data because the inclusion of latent M and Y at pretest and posttest creates two latent interaction terms (instead of one), and four interaction effects (instead of one).Additionally, correlated item residuals often arise when the same latent variable (e.g., M or Y) is assessed with the same measure at pretest and posttest, which must be adjusted for to obtain accurate parameter estimates.In this paper, we describe and evaluate the performance of methods developed to estimate latent interaction effects in two-wave mediation models with BTI effects and provide the following contributions.First, we provide covariance expectations for a better understanding of the two-wave mediation model with BTIs and apply the methods to empirical data.Second, we use Monte Carlo simulations to assess the recovery of the mediation paths and interaction effects across methods.Third, we provide code to estimate this model with the discussed methods.Our goal is to provide researchers with resources to understand and estimate two-wave mediation models with latent BTIs and facilitate the implementation of these models.
The structure of the paper is the following.First, we introduce the statistical mediation model with BTIs.Then, we discuss five different ways to estimate latent BTI effects in the mediation model.Next, we present the results of a Monte Carlo simulation study on the bias, power, type 1 error rate, and interval coverage of the parameter estimates of the statistical mediation model with latent BTIs.Lastly, we demonstrate the methods with data from the JOBS II study.In the supplement, we include the covariance expectations for mediation models with BTIs, along with R and Mplus syntax for estimating those models.

Two-Wave Mediation Model
In this paper, we focus on the simplest mediation model for experimental studies with pretest and posttest data (Valente & MacKinnon, 2017).Although there are other ways to capture change in mediation models (e.g., Howe, 2019;Montoya, 2019;Valente & MacKinnon, 2017), the two-wave mediation model allows researchers to control for baseline levels of M and Y, which in turn mitigates confounding of the posttest M to Y relation (Valente et al., 2020).Twowave mediation models (some with other covariates) have been used to evaluate interventions for drug use (e.g., MacKinnon, 2001), internalizing symptoms (e.g., Jensen et al., 2014;Perrino et al., 2014), depressive symptoms (e.g., Chaudoir et al., 2021), antisocial behavior (e.g., Smith et al., 2014), and substance use and school discipline (e.g., Gonzales et al., 2012).
A mediation model with a randomized treatment, continuous pretest and posttest measures, and BTIs can be expressed as: In Equation ( 1) and ( 2), X represents a binary randomized treatment/control group indicator; M 1 and M 2 represent the mediator at pretest and posttest, respectively; and Y 1 and Y 2 represent the outcome at pretest and posttest, respectively.Note that in randomized treatments, X does not correlate with pretest measures M 1 and Y 1 in expectation.The a-path is the relation between X and M 2 , s m is the temporal stability of M, c y is the cross-lag path between M 2 and Y 1 , b is the relation between M 2 and Y 2 , c 0 is the relation between X and Y 2 , s y is the temporal stability of Y, and c m is the cross-lag path between Y 2 and M 1 : Furthermore, the four coefficients h 1 , h 2 , h 3 and h 4 represent the BTIs, where the relation between X and M 2 or Y 2 is moderated by the baseline scores.In the absence of BTIs, the mediated effect is defined by the product of ab.We can construct confidence intervals and test for statistical significance of ab using the distribution of the product of two random variables (Tofighi & MacKinnon, 2011), bootstrapping (MacKinnon et al., 2004), or deriving a Bayesian credible interval for ab (Yuan & MacKinnon, 2009).If h 1 or h 2 are significant, then there is a baseline-by-treatment interaction on the mediator.As such, the a-path varies as a function of the pretest measures, indicating moderated mediation (Hayes, 2015;MacKinnon, 2008).If h 3 or h 4 are significant, then there is a baseline-by-treatment interaction on the outcome, so the c'-path varies as a function of the pretest measures.
To ascribe a causal interpretation to the mediated effect, some assumptions must hold.For example, we assume that we have specified the appropriate functional form and temporal precedence between the variables.We also assume that the scores on all variables have high reliability and can be ascribed a valid interpretation (Gonzalez & MacKinnon, 2021).Furthermore, there are four no-confounding assumptions needed to identify the mediated effect (Imai et al., 2010;Pearl, 2001;Valeri & VanderWeele, 2013 In randomized treatment scenarios, assumptions 1 and 3 are typically satisfied, but assumptions 2 and 4 are not satisfied because individuals self-select their value on M 2 : Controlling for pretest levels of M and Y makes meeting assumptions 2 and 4 more tenable because M 1 and Y 1 form a joint confounder of the M 2 to Y 2 relation (Mayer et al., 2014;Valente et al., 2020) via the model stabilities and the cross-lags, which we control for in the model.Sensitivity analyses can be used to determine how likely it is that assumptions 2 and 4 are met (Cox et al., 2013;Fritz et al., 2016;Liu & Wang, 2021).Lastly, as mentioned above, the model presented in Equations ( 1) and (2) has linear relations, continuous M and Y, and no XM 2 interaction, so traditional mediated effects estimated from regression models and the causal mediation estimates from the potential outcomes framework are equivalent (MacKinnon et al., 2020;Valeri & VanderWeele, 2013).

Statistical Mediation Model with Latent Variables
As mentioned above, there are many situations in which the mediator or the outcome are measured with error (Gonzalez & MacKinnon, 2021).One way to adjust for measurement error is to specify latent variables for M and Y at each time point.The relation between items and latent variables can be represented using the linear factor model: where M Ti and Y Ti represent the latent scores for individual i at time¼T, which are assessed by a set of j items in vectors m Tij and k Tij : In this case, s mTj and s y Tj represent the item intercepts, k mTj and k y Tj are the factor loadings, and e mTij and e y Tij are the unique scores for individual i on item j at time T. For our purposes, we assume longitudinal scalar measurement invariance between pretest and posttest measures in the mediation model (e.g., Georgeson et al., 2021), and estimate correlated residuals between like-indicators of the same construct across time.An example of the model is presented in Figure 1.Because the model now includes latent variables, the BTIs are necessarily latent interactions.
Given that M 1 , M 2 , Y 1 , and Y 2 are not observed, we consider two general approaches to estimate the latent BTI effect.The first approach is to compute observed scores for M 1 , M 2 , Y 1 , and Y 2 , compute observed BTIs, and then adjust for unreliability and correlated errors.The second approach is to use structural models to estimate the latent BTI effects.We discuss five methods across these two general approaches next.

Scoring Latent Variables with Reliability Adjustment
As mentioned above, when measurement error in M and Y is ignored, interaction effects and mediation paths are attenuated (Gonzalez & Valente, 2022) (Cox & Kelcey, 2021;Croon, 2002;Hayes & Usami, 2020) to the two-wave mediation model and extend the methods when noted.The first step is to obtain scores for M and Y.One option to score the latent variable is to use unit-weighted summed scores (Method 1, referred to as summed scores with reliability adjustment; MacKinnon, 2008).For example, summed score M 1s is calculated by adding the item responses of M 1 : Another option is to use factor scores (Method 2, referred to as factor scores with reliability adjustment; Bartlett, 1937;Lawley & Maxwell, 1971;Skrondal & Laake, 2001).Factor scores can be thought of as weighted summed scores where the weights are determined by the parameter estimates from the linear factor model.In this paper, we use Bartlett scores for Method 2. The weights (i.e., scoring matrix) are estimated by using the expression, where K is a j x p matrix of factor loadings, where j are the indicators and p are the latent variables, and H e is the j x j covariance matrix of the unique scores.To obtain the B w weights to score M 1 and M 2 , K and H e are obtained from a correlated two-factor model with correlated residuals, specifying longitudinal scalar invariance.Once obtained, the item responses are multiplied by the weights and added to estimate factor scores M 1f and M 2f : A similar procedure is used to estimate Y 1f and Y 2f : The second step is to estimate the full covariance matrix of observed scores and interactions and correct specific submatrices of the full covariance matrix.For Method 2, we propose to apply the versions of Croon's (2002) correction to models with interaction effects (Cox & Kelcey, 2021) and correlated residuals (Hayes & Usami, 2020) simultaneously.For example, to correct the 2 Â 2 submatrix R M f , comprised of the variances of M 1f and M 2f and their covariance, we need the scoring weights for M, B wM , and the submatrix of the variances and covariances of the residuals for the j indicators of M, H eM : As such, one can use: 1 Similar steps are used to correct R Y f : To correct the submatrices in the full covariances matrix for Method 1, we adapt Equation ( 7) to summed scores by replacing B wM with an p x j summing matrix S M : S M is populated by 1s on the entries indicating which items load on each variable and 0s elsewhere.For example, to correct submatrix R M s , comprised of the variances of M 1s and M 2s and their covariance, one can use: 2 and R Y sc are obtained, then the variances of the observed interaction terms can be corrected.For interaction terms, their variance can be implied by the variances and covariances of their components (Bohrnstedt & Marwell, 1978).For example, when X and M 1s are centered, the implied variance of XM 1s is, For our specific model, Equation ( 9) simplifies because X is a perfectly reliable, randomized variable, so it does not covary with the pretest measures of M 1s (i.e., Cov X, M 1s ð Þ¼ 0), and its variance does not have to be corrected for measurement error.To correct Var XM 1s ð Þ, we use the corrected variance of M 1s (which is an entry in R M sc ) in Equation ( 9) instead of the observed variance of M 1s : Similar steps would follow to correct Once the variances and covariances of the scores and the interaction terms are corrected, the final step is to fit the mediation model with interactions (Equations ( 1) and ( 2)) to the full corrected covariance matrix using any structural equation modeling program.Note that by fitting the model to the full covariance matrix with maximum likelihood, we assume multivariate normality.

Structural Models for Latent BTIs
Another option to estimate latent BTI effects is to use structural models.Below we apply three approaches to the twowave mediation model that have been shown to recover parameters well in single mediator models with a latent interaction (Gonzalez & Valente, 2022).

Unconstrained Product Indicator (UPI) Approach
Product indicator methods estimate latent interaction effects by specifying a latent interaction variable (e.g., a latent variable for XM 1 ), whose indicators are products of the indicators of the variables involved (Kenny & Judd, 1984;Marsh et al., 2012).There are many ways to specify latent interactions with product indicators (e.g., Ayt€ urk et al., 2021;Jaccard & Wan, 1995;J€ oreskog & Yang, 1996;Kenny & Judd, 1984;Marsh et al., 2004;Ping, 1996).In this paper, we use the unconstrained product indicator (UPI) approach (Method 3; Marsh et al., 2004).For our specific model, indicators of the XM 1 latent variable were the product of each M 1 indicator and X, the mean of latent XM 1 was fixed to Cov X, M 1 ð Þ¼0 (because X is a randomized treatment), and the variance of latent XM 1 was fixed to Var X ð ÞVar M 1 ð Þ per Equation ( 9).The latent variable for XY 1 can be similarly specified.We chose the UPI approach because there are no constraints placed on the estimation of the factor loadings and residual variances of the indicators of the latent interaction terms, which differs from other product indicator methods.Note that the UPI approach relies on multivariate normality.

Latent Moderated Structural Equations (LMS)
Latent Moderated Structural Equations (LMS) is an indirect approach to estimate latent interaction effects (Method 4; Klein & Moosbrugger, 2000) and is available to Mplus users (Muth en & Muth en, 1998/2017).LMS approximates the likelihood of the model with latent interactions using a mixture of normal distributions.A mixture of distributions is used because the conditional distribution of the outcome given the predictors is nonnormal due to the nonlinear relations present in the model (i.e., the linear-by-linear interactions; Kelava et al., 2011).We can represent this model by separating Equations ( 1) and (2) into linear and nonlinear effects: The model is estimated using the EM algorithm (Dempster et al., 1977, Ng & Chan, 2020).Details of the procedure are highly technical, so we refer interested readers to Kelava et al. (2011) and Klein and Moosbrugger (2000).Note that LMS also relies on multivariate normality.

1
The general term in the right hand side of the equation is but Bartlett scores have the property that B w K ¼ I, which simplifies the expression.Regression scores do not have this property.

2
In Equations ( 7) and ( 8) we are summing the residual variances to estimate the error variance of the observed score.By subtracting the error variances from the observed score variance, what remains is the reliable variance.

Bayesian Mediation
Bayesian approaches have been shown to yield accurate parameter estimates from mediation models (Method 5; Enders et al., 2013;Mio cevi c, 2019;Yuan & MacKinnon, 2009) and handle latent interactions (Asparouhov & Muth en, 2021;Mio cevi c et al., 2018).In general, we can sample from posterior distributions of the model parameters in a mediation model with latent variables (with no interactions) using Markov-chain Monte Carlo (MCMC) estimation (Mio cevi c, 2019).An intermediate step in the estimation of the statistical mediation model with latent variables using a Gibbs sampler is the sampling of a latent variable score for each latent variable (Mio cevi c, 2019).As such, we can incorporate latent BTIs in a Bayesian mediation model by multiplying in each iteration the sampled value of the latent M 1 or Y 1 score by X and regress the posttest mediator and outcome on those values.

Present Study
Examining baseline-by-treatment interactions in statistical mediation analysis can provide insight into how and for whom an intervention works (MacKinnon, 2008).However, it is unknown whether the methods that have been shown to recover parameter estimates well in a single mediator model with one latent interaction (Cheung & Lau, 2017;Gonzalez & Valente, 2022) will also recover parameters well in a more complex model, such as the two-wave mediator model with two latent interactions and four BTI effects.Continuing, we are not aware of any resource that generally describes how to incorporate latent BTI effects in mediation analysis and that provides syntax to guide researchers on the estimation of these theoretically relevant models to intervention research.Below, we conduct a simulation study to examine the bias, power, type 1 error rate, and interval coverage of the estimates in the two-wave mediation model and the latent BTIs, illustrate the methods with an applied intervention example and provide syntax in the supplement to fit these models.

Simulation Factors
True covariance matrices for the two-wave mediation model with baseline-by-treatment interactions were derived and are shown in the supplement.Datasets were generated in the R statistical environment by varying the factors found in Table 1.
Factors that did not vary in the simulation were the effect size of the c'-path (small), the stabilities of M and Y (at r ¼ .7), the cross-lag effects (zero effect), and the correlation of M 1 and Y 1 at baseline (r ¼ .575).Overall, there were 96 conditions, with 500 replications per condition.
All path effect sizes were obtained from Valente and MacKinnon (2017) and approximately match Cohen's f 2 benchmarks.Residual covariance values were specified to be similar to Cohen's effect sizes of small and medium correlation, and the sample sizes studied are achievable in social science applications and have been found to recover parameter estimates well in a single mediator model with one latent interaction (Gonzalez & Valente, 2022).Items for M 1 , Y 1 , M 2 , and Y 2 were simulated to be unidimensional and adhere to strict invariance across time.All item intercepts were zero, factor loadings for M 1 , Y 1 , M 2 , and Y 2 were drawn from a uniform random distribution with limits between .5 and .7,and the residual variances were the complement of the factor loadings so that the total variance of each item at baseline was equal to 1.Note that all the indicators for the conditions presented here are multivariate normally distributed, so they meet the distributional assumptions of the methods examined.In the supplement, we show that the stability paths are biased when discrete, nonnormally distributed indicators are treated as continuous.Therefore, we do not recommend using these methods in these scenarios.

General Procedure
Each dataset was analyzed with each of the five methods to estimate the latent BTI effect: summed scores with reliability adjustment, factor scores (Bartlett) with reliability adjustment, unconstrained product indicator (UPI) approach, latent moderated structural equations (LMS), and Bayesian mediation analysis.Note that all longitudinal latent variable models fit to the data were specified to be scalar invariant and correlated residuals were prespecified.For the summed scores and factor score methods, we used the parameter estimates from the longitudinal factor model for M 1 and M 2 and the longitudinal factor model for Y 1 and Y 2 to correct the covariance matrices as discussed above.Then, the statistical mediation model with observed BTIs was fit to both corrected covariance matrices using lavaan in R (Rosseel, 2012).Furthermore, the LMS approach (using the XWITH argument for latent interactions, random effect estimation, and numerical integration) and UPI approach (with constraints explained above) were estimated in Mplus and processed in R using the MplusAutomation package (Hallquist & Wiley, 2018).Finally, the Bayesian mediation model was estimated in JAGS (Plummer, 2003) in R using the R2jags package (Su & Yajima, 2020). 3Similar to Gonzalez and Valente (2022), we specified diffuse prior [4 for M and 8 for Y] [6 for M and 6 for Y] Residual covariance (2) 0.10, 0.30 Note: S is small effect size; M is medium effect size.

3
We imposed longitudinal scalar invariance constraints by using the same parameter in different parts of the model.For example, a single posterior distribution represented the relation between M 1 and the item m 11 and M 2 and the item m 21 .Also, we modeled the correlated residuals like the specific factors in bifactor models: the factor had two indicators, which were the same item across time, and these factors were orthogonal to all other variables.distributions for model parameters to mimic situations with limited prior beliefs about the parameters: factor loadings had normal distributions N(1, 20) truncated at 0 (i.e., the loadings were positive), item intercepts had normal distributions N(0, 20), residual variances had inverse gamma distributions IG(1,1), latent variable means for exogenous variables had a normal distribution N(0,1), and regression coefficients and equation intercepts had normal distributions N(0, 1000).Pilot analyses suggested that the chains for each parameter achieved a potential scale reduction (PSR) factor <1.1 with at most 5,000 iterations.As such, we used three chains which each sampled 5,000 iterations (1,000 of which were burn-in iterations).

Monte Carlo Outcomes
We examined three bias outcomes.Raw bias was computed in each replication by the difference between the estimated parameter and the data-generating value. 4We also report relative bias (i.e., raw bias divided by the true value of the estimate) for conditions in which the true value of the parameter is nonzero, and standardized bias (i.e., raw bias divided by the standard deviation of the parameter across replications per condition) for conditions in which the true value of the parameter is zero.Bias outcomes were analyzed using a mixed ANOVA as a function of the simulation factors as between-group factors, the latent variable method as a within-group factor, and all possible interactions.Conditions with relative bias >.10 (Kaplan, 1988) or standardized bias >.40 (Collins et al., 2001), and predictors with partial-g 2 >.005 were further discussed.We also estimated 95% nonparametric percentile bootstrap confidence intervals (with 500 bootstrap draws) for model parameters and the mediated effect across all methods (MacKinnon et al., 2004) except Bayesian mediation, in which equal-tail credible intervals based on the posterior distributions were used (Yuan & MacKinnon, 2009).Power was defined as the proportion of replicates in which the interval did not contain zero for parameters with nonzero true values, and type 1 error was the proportion of replicates in which the interval did not contain zero for parameters with zero true values.Note that power was estimated for the a-path, b-path, c-path, stability estimates, and the BTI effects, while type 1 error was only estimated for the BTI effects and cross-lag effects.Finally, interval coverage was defined as the proportion of replicates where the derived interval contained the true value of the parameter.Power, type 1 error rates, and interval coverage were analyzed with regression trees (Gonzalez et al., 2018), where splits had to yield an increase of at least R 2 ¼ .005 on the outcome to be interpreted.

Parameter Bias and Interval Coverage
Table 2 shows the parameter bias for the mediated effect and the BTI coefficients.In the ANOVAs predicting parameter bias (i.e., raw, relative, and standardized bias) from all simulation factors and methods used to estimate the latent interaction effect, there were three parameters with factors that had effect sizes with g 2 > .005.There were main effect for the method used on the raw bias of the a-path (g 2 ¼ 0.007), c-path (g 2 ¼ 0.007), cross-lag effect of Y (g 2 ¼ 0.011), and on the raw bias (g 2 ¼ 0.041) and relative bias (g 2 ¼ 0.041) of the stability of Y. Also, there were significant interactions between the method and sample size and the method and the number of items on the raw bias (g 2 ¼ 0.014 and 0.030, respectively) and relative bias (g 2 ¼ 0.014 and 0.030, respectively) of the stability of M, and a significant interaction between the method and sample size for the standardized bias on the cross-lag effect of Y (g 2 ¼ 0.007).However, across methods, all the relative bias estimates for nonzero parameters were below .05,and all the standardized bias estimates for zero parameters were below .40.Furthermore, Table 2 shows the interval coverage across methods.Regression trees predicting coverage from the simulation factors and the BTI estimation method did not yield any splits, so there were no important differences in interval coverage across simulation factors.All interval coverage estimates are within the range of .925and .975.

Statistical Power and Type 1 Error Rate
Regression trees for the analysis of statistical power and type 1 error by the simulation factor are shown in the supplement.For statistical power, the effect size of the parameter and sample size were significant predictors, as expected.Furthermore, the h1 coefficient (e.g., the relation between the XM 1 interaction and M 2 ) had higher power in conditions with six items than with four items at N ¼ 500.Recall that all BTI effects had a small effect size.Power was .310and .430for conditions with four and six items in the mediator, respectively.We presume that difference in power might be because more items yield more stable latent variables.Also, the c-path had higher power in conditions when the a-path had a small effect size (power ¼ .750)than in conditions with a medium effect size (power ¼ .650).Although Table 3 might suggest that there are some power differences across methods, the method to estimate the latent interaction effect was not an important predictor in the trees.On the other hand, results suggest that there were no significant predictors for type 1 error rate for the BTI estimates and crosslag effects.Table 4 shows the type 1 error rates across methods, and all the estimates were within the bounds of .025and .075.

Summary
Our results suggest that, as long as the indicators are multivariate normally distributed, summed scores with reliability adjustment, factor scores with reliability adjustment, UPI, LMS, and Bayesian mediation provide unbiased parameter estimates in the mediation model, the 95% confidence/credible interval coverages perform appropriately, and the type 1 error rates were close to .05.Furthermore, methods did not differ significantly on parameter bias, power, type 1 error rates, and interval coverage.As such, similar to the findings in the single mediator model with one latent interaction, the 4 Rescaled true values for the corrected summed score model were derived to serve as reference values for raw bias.Summed scores are on a different metric than the latent variable models, similar to Georgeson et al. (2021).
studied methods (which include our proposed extension to the scoring approaches to adjust for measurement error and the effect of correlated residuals simultaneously) recover the parameter estimates well.

Illustration
Data were obtained from The Job Search Intervention Study (JOBS II), which is a randomized field intervention study that, among other aims, examined if a job training intervention for unemployed workers improved symptoms of depression (outcome) by enhancing job search self-efficacy (mediator), controlling for baseline covariates (Vinokur et al., 1995;Vinokur & Schul, 1997).The original sample size was N ¼ 1,801, but to simplify the analyses, we used only participants who had complete item-level data for the mediator and outcome at both timepoints.Therefore, the sample size used for this illustration was N ¼ 1126, with N ¼ 347 in the control group and N ¼ 779 in the treatment group.Depressive symptoms were measured by the eleven five-category items on the Hopkins Symptoms Checklist (coefficient a ¼ .807at pretest and a ¼ .884at posttest), which asked about the frequency of several symptoms, such as loneliness and hopelessness, experienced by the respondent in the past two weeks.Higher scores represented more depressive symptoms.The job search selfefficacy measure had six five-category items (coefficient a ¼ .868at pretest and a ¼ .886at posttest), which asked how confident respondents were about several aspects of the job process, such as completing a job application and a resume.Higher scores represented higher job search self-efficacy.Preliminary analyses suggest that the mediator and outcome measures are largely unidimensional and meet the longitudinal scalar invariance assumptions (see supplement).Here, we test a statistical mediation model with latent baseline-by-treatment interactions to investigate if the effect of the training on job search self-efficacy and depressive symptoms depends on the baseline levels of those constructs.The methods examined were summed scores with reliability adjustment, factor scores with reliability adjustment, UPI, LMS, and Bayesian mediation, and used either 95% percentile bootstrap confidence intervals or 95% equal-tail credible intervals for statistical inference.The specification of the models, including the prior distributions of parameters in Bayesian mediation, closely follow the specification described for the simulation.

Results
Table 5 shows the parameter estimates of the mediation model from Equations ( 1) and ( 2).The estimates were similar across models, so we will only interpret the estimates for the factor scores with reliability adjustment.

Discussion
Estimating baseline-by-treatment interactions can reveal how and for whom an intervention works, which is valuable information for the dissemination and implementation of interventions.In this paper, we discussed five different methods (e.g., summed scores with reliability adjustment, factor scores with reliability adjustment, unconstrained product indicator approach, latent moderated structural equations, and Bayesian mediation) for estimating the statistical mediation model with baseline-by-treatment interactions when the mediator and outcome are latent variables.We also evaluated the estimation properties of the approaches and illustrated the methods using an applied example and presented code in the supplement to facilitate testing for latent baseline-by treatment interactions.Our simulation results suggest that the five approaches examined led to unbiased effects, similar power, interval coverage of approximately .95, and type 1 error rates close to .05.As such, each of these methods estimated the mediation effects and the latent BTI effects accurately.However, across methods, the estimates of power to detect small effect sizes were well below .80 even for a condition with N ¼ 500, which might be a limiting factor in applied research.Larger sample sizes might be needed to detect small effects using this model.
Given that all methods provided accurate results, additional simulations may reveal situations in which they provide different answers.Determining which latent interaction approach to use with any given model remains an open question (Kelava & Brandt, 2022).As such, researchers estimating latent BTI effects in mediation models could follow current recommendations from the literature.In a recent review of work in latent interactions, Kelava and Brandt (2022) recommend using the LMS approach in situations in which there are normally-distributed latent variables and indicators.Kelava and Brandt (2022) mention that LMS is more efficient than the UPI approach (there is a loss of degrees of freedom when the estimation of the latent interaction term is unconstrained) and the Bayesian mediation approach (where efficiency depends on the specified prior), and that more research is needed to understand the performance of the scoring approaches.Other considerations include amount of background, analytical flexibility, computational time, and programming needed for their implementation.For example, the specification of the LMS procedure is relatively simple and took about 1-3 minutes to estimate a single model, and the procedure is available in Mplus.If access to Mplus is a limitation, then researchers could use any of the other methods presented-Bayesian mediation can be used in any program that does MCMC estimation, while the others can be estimated in any structural equation modeling program that takes a covariance matrix as an input.However, estimating the Bayesian mediation model would require researchers to be familiar with Bayesian inference and specifications of prior distributions, which might pose a challenge, along with long computation times.For the model in the illustration, the estimation and parallelization per MCMC chain took around 50 minutes in JAGS using a 13in MacBook Pro 2020 with 32GB of RAM.On the other hand, the coding to estimate the model using the UPI approach or summed scores and factor scores with reliability adjustment could be cumbersome, but the estimation time is almost immediate.As such, we advise researchers to review the code we provide in the appendix to gain familiarity with each of these approaches.

Extensions to the Model
In applied settings, researchers often use extensions of the two-wave model examined here to test specific theories.Depending on the model and situation, these extensions might or might not be easily accommodated by the procedures examined.These extensions require more simulation work to determine their performance, but we list how the methods examined could handle these extensions.For example, it may be of interest to include other covariates in Equations ( 1) and (2) to increase statistical power (e.g., Gonzales et al., 2012).In cases in which the covariates are observed and perfectly reliable, the extensions are simpler.If using summed scores and factor scores with reliability adjustment, the covariates would be included in the covariance matrix that is analyzed and additional adjustments to the covariance matrix would not be necessary.For the UPI, LMS, and Bayesian methods, the observed covariate could be specified and included in the corresponding equations.
In cases in which the covariates are latent or have measurement error, adding a covariate adds another layer of complexity.The methods with summed scores and factor scores with reliability adjustment would need to score the covariate, correct its variance for measurement error, and include it in the analyzed covariance matrix.For the UPI, LMS, and Bayesian mediation, a latent covariate could be specified and included in the corresponding equations.Furthermore, the methods examined could be extended to situations in which M and Y are measured at baseline and at two additional timepoints (i.e., three-wave experimental designs).A third wave adds a layer of complexity to the two-wave model because M3 and Y3 are two additional latent variables to estimate.Perhaps in this design, the mediated effect of most interest is the X !M2 !Y3 relation, where each path can be moderated by baseline levels of M and Y.The methods with summed scores and factor scores with reliability adjustments would correct the 3 Â 3 variance-covariance submatrices of M and Y for measurement error and correlated residuals.If interested, researchers can take the steps to correct the XM1 and XY1 effects in the former Y2 equation and correct the XM2 and XY2 effects in the Y3 equation.For the UPI, LMS, and Bayesian methods, latent M3 and Y3 could be specified and included in the models.
Finally, the methods examined could be extended to situations in which X is a continuous variable (e.g., an intervention fidelity indicator).Although the procedures would stay the same when X is perfectly reliable, the procedures would differ when X has measurement error or is represented by a latent variable.The methods with summed scores and factor scores with reliability adjustments would need to correct the variance of X for measurement error, and consequently the variances of XM1 and XY1, in the analyzed covariance matrix.For the UPI, LMS, and Bayesian methods, a (single-indicator or multiple-indicator) latent variable could be specified for X and latent interactions for XM1 and XY1 and then the researcher would proceed as described.Note that we assume multivariate normality for all of these effects.Beyond the modeling complexity, it is important to note that the nonrandomized nature of X opens the door to unmeasured confounders in the X !M2 and X !Y2 relations.

Limitations and Future Directions
An important limitation of this study is that the results we presented rely on the correct specification of the models.A future direction is to investigate how longitudinal noninvariance (Kelava & Brandt, 2022), unmodeled multidimensionality (Gonzalez & MacKinnon, 2018), distribution of the latent variables (Cham et al., 2012), or structural misspecifications (Kelava & Brandt, 2022) affect the estimation of the statistical mediation model with latent BTIs.Specifically, it would be important to study whether ignoring these aspects leads to spurious observed baseline-by-treatment interactions.Additionally, we could have studied how sensitive our results are to different specifications of the models.For example, one could study how sensitive our results are to prior distributions used in the Bayesian mediation model.In our study, diffuse priors were able to recover the parameters well, but this might not always be the case.Furthermore, although previous research suggests that varying stability does not affect the accuracy of the mediated effect estimate in two-wave mediation models without interactions (Valente & MacKinnon, 2017), one could examine if these findings also hold in two-wave mediation models with latent BTIs.Also, it would be interesting for future studies to increase the number of replications per simulation condition to obtain more precise estimates of our Monte Carlo outcomes.Finally, there are other methods to handle latent BTIs that we did not study, such as two-stage least squares methods (Bollen & Paxton, 1998) or others discussed by Kelava and Brandt (2022).Overall, we encourage researchers to estimate statistical mediation models with latent baseline-by-treatment interactions using any of the methods discussed in this paper for accurate estimation of the effects and to help with the implementation and dissemination of interventions.

Figure 1 .
Figure 1.Statistical mediation model with latent baseline-by-treatment interactions and correlated residuals.

Table 2 .
Bias estimates and interval coverage of the parameters from the mediation model with latent BTIs.

Table 3 .
Statistical power of the mediation model parameters across methods.Power estimates for the stability of M and Y are not shown because the power was almost 1 across conditions.(s) is for small effect size and (m) is for medium effect size.

Table 4 .
Type 1 error for parameters from the mediation model across methods.

Table 5 .
Parameter estimates with 95% bootstrap confidence intervals [lower limit, upper limits] from a statistical mediation model with latent BTIs in the JOBS II data.