Interaction Effects in Cross-Lagged Panel Models: SEM with Latent Interactions Applied to Work-Family Conflict,Job Satisfaction,and Gender

Time-series and panel data methods are used to study a variety of phenomena in organization science (Bliese et al., 2020). Of these, the cross-lagged panel model (CLPM) and its variants are a popular approach (e.g., Zyphur, Allison, et al., 2020; called panel vector autoregressions or panel VARs in econometrics; see Abrigo & Love, 2016; Love & Zicchino, 2006). These models estimate lagged relationships to allow inferences that are consistent with the fact that causality unfolds over time (Zyphur, Voelkle, et al., 2020), typically estimated as structural equation models (SEM) with fewer than 10 observations over time.

The CLPM can also be used to study conditions under which lagged effects vary by testing interactions. ¹ Such tests often use product terms such as xz for predictors x and z (Cortina et al., 2021; Preacher et al., 2006, 2016), which can be calculated as the products of lagged predictors. For example, Eby et al. (2015) used the product of lagged predictors to test the interaction between (past) organizational support and (past) mentoring to predict (current) organizational citizenship behaviors. Related examples include: 1) models that account for measurement error with latent two-way (e.g., Dormann & Zapf, 1999) and three-way interactions (e.g., Lian et al., 2014); 2) combining interactions and nonlinear terms (e.g., Lin et al., 2017); 3) tests for sub-group differences in lagged effects (e.g., Cieslak et al., 2007; Houkes et al., 2003; Zablah et al., 2016); and 4) lag-as-moderator models (LAM) wherein the amount of time between measurements is treated as a moderator of lagged effects (Selig et al., 2012).

However, recent work suggests caution is needed when testing interactions in panel data, because different types of interactions may be conflated. Specifically, panel data contain stable ‘between-level’ (B) variation (e.g., personality; Hamaker et al., 2015) and dynamic ‘within-level’ (W) variation due to time-varying factors (e.g., emotions; Zyphur, Allison, et al., 2020). The multilevel modeling literature notes the importance of separating these sources of variation when testing interactions to allow precisely matching theory to hypotheses while avoiding the conflation of B and W effects (Preacher et al., 2016). However, when testing interactions in CLPMs, researchers have overlooked the implications of W and B sources of variation. In turn, this may have led to interaction estimates being difficult to interpret or resulted in flawed conclusions pertaining to significance testing.

In this paper we address testing interactions in CLPMs by integrating the literatures on CLPMs, multilevel moderation, and latent interactions. We start with the multilevel nature of longitudinal data and how panel data allow testing three types of interactions: (1) purely within-level or WxW interactions among dynamic factors; (2) cross-level or BxW interactions among stable B and dynamic W factors; and (3) purely between-level or BxB interactions among stable factors (Preacher et al., 2016). By extension, we also describe how to incorporate nonlinear and higher-order interactions of various kinds, including but not limited to within-level or W² nonlinearity; between-level or B² nonlinearity; and combinations of these with cross-level interactions, which we call B²xW, BxW², or B²xW² cases.

We then show how CLPMs may conflate such effects, including when estimating latent interactions (e.g., Dormann & Zapf, 1999; Lian et al., 2014) or attempting to control for stable B factors (e.g., Carpenter & Fredrickson, 2001). To address this issue, we describe a computationally tractable SEM framework for testing interactions in CLPMs with a Bayes estimator and latent B and W components (Asparouhov & Muthén, 2020). We illustrate our method for estimating and interpreting CLPM interaction effects using a dataset with time-varying measures of job satisfaction and work-family conflict as well as a stable measure of gender. We also run simulations showing that our approach is less biased than alternatives. We provide code for implementing our approach in Mplus, as this software is widely known in organization research and elsewhere (see Cortina et al., 2021; Heck & Thomas, 2015; Preacher et al., 2010, 2016). We conclude by discussing how our general approach can be used in other types of longitudinal panel data models including the general cross-lagged panel model (GCLM; Zyphur, Allison, et al., 2020), and latent curve or latent growth models (LGM; Chan, 1998). In sum, by capitalizing on the multilevel nature of longitudinal data we offer guidance for testing different types of multilevel interaction effects in a familiar CLPM framework estimated as an SEM.

A Multilevel Understanding of Cross-Lagged Panel Models

Longitudinal data are inherently complex, involving multiple sources of variation. When emotions are measured over time, for example, the resulting data exhibit both stability due to time-invariant factors and instability due to occasion-specific factors (Beal, 2015; Gabriel et al., 2019), which implies a multilevel structure. From this perspective, stable factors are higher-level or between-level (e.g., between-person or merely B), whereas dynamic factors are lower-level or within-level (e.g., within-person or merely W) phenomena (Zyphur, Allison, et al., 2020). In brief, B components reflect stable factors such as time-invariant individual differences (e.g., personality), and W components reflect temporary deviations around stable levels, due to occasion-specific factors (e.g., random events). Thus, to represent W components necessitates repeated measures of a variable with W variation and based on these a stable B component in the same content domain can be treated as a latent variable reflected in the repeated measures. It is notable that for variables without W variation such as stable gender or race, a single stable measure of the B variable will suffice and drawing an analogy to multilevel discussions of ‘direct-consensus’ measures, it would also be possible to collect ‘direct-rating’ measures for B terms (e.g., participants’ reporting their ‘usual’ or ‘overall’ job satisfaction across a studied timespan). ²

We illustrate the multilevel nature of longitudinal data using two variables x and y, which might reflect work-family conflict and job satisfaction measured for a person i at an occasion t, where i = 1, 2, …, N and t = 1, 2, …, T. The multilevel logic of a stable B and dynamic W component can then be shown as follows (see Figure 1 for an SEM diagram):

x i t = x B i + x W i t

(1)

y i t = y B i + y W i t

(2)

OPEN IN VIEWER

Figure 1.

An SEM diagram distinguishing W and B components in a CLPM for x and y measured at three occasions.

In multilevel parlance, x i t and y i t would be Level-1 variables measured at the lowest level of analysis and thus containing latent B and W variation. For example, if x i t were work-family conflict measured over time, x B i would capture time-invariant stable factors impacting an individual’s average level of work-family conflict (i.e., a “trait-like” part; Hamaker et al., 2015, p. 102), whereas x W i t would capture time-varying dynamic factors associated with person i's experienced conflict at a specific time t. Indeed, we could use single-level SEM notation with a stable B term η i and a dynamic W term as ε i t (see Curran, 2003). However, in keeping with our multilevel conceptualization, we use B and W terms here and emphasize to the reader that these are actually latent.

In this context, there are two primary ways to decompose a variable, namely manifest group-mean centering (i.e., an ‘observed-means’ approach) versus latent group-mean centering (see Lüdtke et al., 2008; Preacher et al., 2010). In the observed-means approach, researchers estimate the B part for each sampled entity using an average of the observed variables over time, with the deviations from this average representing the W part. In the latent centering case, a model-based approach is used to estimate the B and W variance components without attempting to estimate the underlying B and W scores directly. We compare these approaches later with real data and in simulations, but substantial prior research suggests that latent centering is superior to observed centering because the latter fails to account for uncertainty in the B and W components (Preacher et al., 2016). This is particularly important in longitudinal panel data models with lagged effects, wherein observed-means centering is known to cause ‘Nickell’ bias in lagged effect estimates (see Nickell, 1981; Zyphur, Allison, et al., 2020).

To illustrate, we start by describing the dynamic W terms used for lagged effects in CLPMs, using notation from Zyphur, Allison, et al. (2020) where β coefficients have subscripts to indicate predictors and superscripts for outcomes:

x W i t = β W x x ⋅ x W i t − 1 + β W y x ⋅ y W i t − 1 + u W i t x

(3)

y W i t = β W y y ⋅ y W i t − 1 + β W x y ⋅ x W i t − 1 + u W i t y

(4)

In Eqs. 3 and 4, the autoregressive (AR) terms β W y y and β W x x indicate persistence in a process over time; cross-lagged (CL) terms β W y x and β W x y indicate unique predictions of one variable's future state by another variable's past state; and a time-specific ‘impulse’ u W i t captures occasion-specific events or other factors that cause unpredictable (often conceptualized as random) increases or decreases in x and y at any occasion. As is common in the time-series literature, we use the term impulse instead of error (or an estimated residual) to emphasize the comprehensive nature of this term, encompassing not just variance unaccounted for by prior occasions but also potentially occasion-specific environmental changes (see Zyphur, Allison, et al., 2020). Using our example, if x were work-family conflict and y were job satisfaction, the AR term β W x x would indicate persistence in work-family-conflict and the AR term β W y y would represent persistence in job satisfaction over time. Further, for CL effects, β W y x would indicate how past job satisfaction uniquely predicts future work-family-conflict and β W x y would indicate the same for work-family-conflict predicting job satisfaction. Lastly, the u W i t ‘impulse’ terms for x and y would represent occasion-specific random fluctuations in work-family conflict and job satisfaction that are not predicted by either job satisfaction or work-family conflict at the previous occasion. For example, the impulse term for job-satisfaction at a particular occasion may capture the effect of an individual unexpectedly receiving praise from their supervisor.

Equations 3 and 4 represent a familiar lagged-effects logic (for an associated logic of causality, although not the focus of our paper, see Granger, 1969; for a critical discussion see Zyphur, Allison, et al., 2020; Zyphur, Voelkle, et al., 2020). However, by substituting the lagged W terms from Eqs. 3 and 4 into the original multilevel decomposition in Eqs. 1 and 2, the model can be shown as a type of ‘fixed-effects’ CLPM (in econometrics terms) that implies stable B terms are held constant when examining W effects (see Figure 1):

x i t = x B i + β W x x ⋅ x W i t − 1 + β W y x ⋅ y W i t − 1 + u W i t x

(5)

y i t = y B i + β W y y ⋅ y W i t − 1 + β W x y ⋅ x W i t − 1 + u W i t y

(6)

In Eqs. 5 and 6, x B i and y B i represent the time-invariant latent means of x and y for person i. Notably, these are called ‘random intercepts’ in Hamaker et al. (2015) and ‘unit effects’ in Zyphur, Allison, et al. (2020), but we use the more common ‘fixed effects’ terminology from the econometrics literature to emphasize that these B terms are held constant when estimating W effects (see also Hamaker & Muthén, 2020). As we show later with an SEM, a fixed-effects specification that controls for stable B components is facilitated by allowing the latent B variables to freely covary (as in Hamaker et al., 2015; Zyphur, Allison, et al., 2020). ³

To continue, a more general structural specification allows regression among B terms when causal effects among them can be theoretically justified, such as when studying individual differences or organizational climate and culture. For example, consider a model for y B i that includes the B part of x as a predictor and an observed B predictor variable z B i (where z B i is a Level-2 variable in terms of its measurement at the stable B level):

y B i = β B x y ⋅ x B i + β B z y ⋅ z B i + ζ B i y

(7)

with β B x y and β B z y implying, respectively, B effects of x B i and z B i on y B i , and ζ B i is a B residual or disturbance term associated with y B i . Here, the stable B part of y i t (e.g., the latent mean of job satisfaction) is predicted by the stable B part of x i t (e.g., the latent mean of work-family conflict) and a stable B predictor z B i (e.g., the observed variable gender). Therefore Equation (7) represents a regression among model-estimated latent averages for each person, which might be called a type of means-as-outcomes multilevel model (Raudenbush & Bryk, 2002).

By substituting the B model for y from Equation (7) into the larger model for y in Equation (6), we can show a full B + W representation with both stable and dynamic parts as follows:

y i t = ( β B x y ⋅ x B i + β B z y ⋅ z B i + ζ B i y ) + ( β W y y ⋅ y W i t − 1 + β W x y ⋅ x W i t − 1 + u W i t y )

(8)

In Equation (8), the first three terms represent the B model and the latter three represent the W model. When using the example of job satisfaction ( y i t ) , work-family conflict ( x i t ) and gender ( z i ) , the first two predictors may be separated into dynamic W and stable B components, whereas gender is treated as time-invariant and therefore purely B. Such a decomposition into B and W parts is common in the multilevel literature (Preacher et al., 2010), and forms a basis for CLPMs when estimated in a multilevel framework, such as dynamic structural equation modeling (DSEM; Asparouhov et al., 2018; Zhou et al., 2021). A key point of this decomposition is that it implies fundamentally different kinds of effects: the W model contains lagged effects that are of primary interest to researchers who focus on AR dynamics and CL terms for causal inference, whereas the B model involves regression among stable factors. Our primary focus here is on the W model because the purpose of a CLPM is typically to estimate causal relationships by analyzing temporally ordered variables (Zyphur, Allison, et al., 2020). Yet, for the interested researcher, we formulate both B and W interaction models for completeness.

Latent Interactions in Cross-Lagged Models

Based on previous work on multilevel moderation (Preacher et al., 2016), we can also begin drawing conclusions about interactions in CLPMs. An interaction exists when the effect of a predictor x on an outcome y varies across the levels of another variable z. To test this, x and z are usually multiplied and their product xz is used as a predictor of y. As noted above, lagged-effects models that test interactions use this or conceptually similar approaches with lagged predictors. The idea behind using product terms in CLPM is that “extending this method to autoregressive models is straightforward, and no modifications to the mathematics are required to accommodate lagged values. That is, rather than (or in addition to) controlling for information available in prior lags, use lag variables as moderators” (Hayes, 2015, p. 16). This belief about the generality of methods for testing interactions has motivated the use of related interaction testing methods in longitudinal data (e.g., Selig et al., 2012; Zablah et al., 2016). However, some modifications to these approaches are needed to ensure that different sources of multilevel variation and effects are not conflated.

First, product terms among observed variables will confound different types of multilevel interactions when observed variables have B and W parts—consider the implications of expanding the product term x i t y i t = ( x B i + x W i t ) ( y B i + y W i t ) = x B i y B i + x B i y W i t + y B i x W i t + x W i t y W i t ; where there are now purely B, purely W, and two cross-level interactions implied in what would be a single predictor x i t y i t (we return to this point in Equation. (24), below). Second, B and W components are not observed and thus are better conceptualized as latent variables (consider interactions among η i and ε i t terms in SEM; see Ozkok et al., 2019).

In sum, multilevel research shows that testing interactions with multilevel data requires (1) decomposing observed variables into B and W parts, and then (2) specifying latent interactions separately with these parts (Preacher et al., 2016; Zyphur et al., 2018). Complementing this, recent work on autoregressive confirmatory factor analysis (AR-CFA) shows how such effects can be tested with latent interactions (see Ozkok et al., 2019). To adapt this research for CLPMs, we start by describing three types of interactions with panel data: purely within-level (WxW) interactions among dynamic factors; cross-level (BxW) interactions among stable B and dynamic W factors; and purely between-level (BxB) interactions among stable factors. We then elaborate on other cases.

Purely Within-Level or WxW Interactions

Many studies that test interactions using CLPMs seek to examine how CL effects vary across the levels of another dynamic variable. For example, Zablah et al. (2016) investigate whether the effect of past customer satisfaction on future job satisfaction is moderated by past customer engagement. In such cases, lagged predictors and moderators are often described as time-varying (alluding to W effects) and discussions of the effects imply that time-varying moderating variables are of interest. This focus on the dynamic W parts of variables also exists in studies of time-varying moderators of AR terms, such as investigations of how time-varying stress moderates the persistence (i.e., AR term) of emotional states (see Koval & Kuppens, 2012). Thus, the literature has various examples that analyze W interactions between lagged predictors ( x W i t − 1 and z W i t − 1 ) as these interact to predict an outcome ( y W i t ) . However, Koval and Kuppens (2012) findings illustrate that—although it is often overlooked—the lagged predictors in a typical CLPM may interact to predict their own future values, such as the interaction among x W i t − 1 and y W i t − 1 predicting x W i t in Eq. 3 (along with predicting y W i t in Equation (4)). Therefore, for simplicity and novelty reasons, as well as for illustration, we use this bivariate case to present WxW interactions here. It is easy to extend this to trivariate CLPM models with another time-varying variable z for testing xz, yz, or xyz interactions, as needed.

With an interest in forming latent interactions among W factors we first present purely W moderation effects. We show the simplest WxW case that does not require any external moderators and instead forms latent interactions among the lagged W components of x and y from Eqs. 3 and 4, which are always available in CLPMs (see Figure 2a):

x W i t = β W x x ⋅ x W i t − 1 + β W y x ⋅ y W i t − 1 + β W x y x ⋅ x W i t − 1 ⋅ y W i t − 1 + u W i t x

(9)

y W i t = β W y y ⋅ y W i t − 1 + β W x y ⋅ x W i t − 1 + β W x y y ⋅ x W i t − 1 ⋅ y W i t − 1 + u W i t y

(10)

In Equations (9) and (10), the purely W interaction effects among the dynamic parts of x and y are represented by β W x y x and β W x y y (we show both predictors in subscripts to signify an interaction). Uniquely, these effects imply that the lagged relationship for each predictor-outcome pair is conditional on the lagged value of the other predictor, which has not been explored in most CLPMs. To explain, β W x y x and β W x y y imply that (i) an outcome variable's past moderates the CL effect of a predictor, such as in the case of past job satisfaction impacting the effect of work-family conflict on future job satisfaction; or (ii) that a predictor's past moderates the persistence (or ‘inertia’) of an outcome variable, such as previous work-family conflict moderating the AR effect of job satisfaction. These effects may be quite common, as in Koval and Kuppens (2012) who found that inducing stress caused a decrease in emotional inertia (i.e., a smaller AR effect). In terms of our example, consider that the effect of past work-family conflict on future job satisfaction may be weaker when past job satisfaction is high. This would imply that high job satisfaction provides a buffer, which mitigates the effects of work-family conflict on future satisfaction. We examine this case later, but again note that additional predictors can easily be included for more typical two-way WxW interactions.

OPEN IN VIEWER

Figure 2a.

CLPM for x and y at three occasions with latent WxW interactions.

OPEN IN VIEWER

Figure 2b.

CLPM for x and y at three occasions with latent BxW interactions.

OPEN IN VIEWER

Figure 2c.

CLPM for x and y at three occasions with latent BxB interactions.

OPEN IN VIEWER

Figure 2d.

CLPM for x and y at three occasions with latent WxW, BxW, and BxB interactions.

The point is that only the W part of each variable should form product terms to estimate a purely W interaction. This means that the stable B components are not only controlled, they also do not bias estimates of the interaction among x and y (a concern in the multilevel literature; see Preacher et al., 2016). Indeed, this would apply even when assessing WxW effects among variables that are relatively stable over time (e.g., personality), because these might still exhibit some W variation (e.g., across contexts; Fleeson & Gallagher, 2009; Tett & Burnett, 2003) and tests of any dynamic W parts should not be biased by their stable B parts. Finally, our exploration of WxW effects shows how estimating these has always been possible in CLPMs because, by design, these models always include at least two lagged predictors. This realization opens up new opportunities to use CLPMs to test interactions—consider that every published CLPM can be re-estimated to test WxW interactions.

Cross-Level or BxW Interactions

Next, consider that time-invariant B variables, such as the stable components of personality, gender, or the environment may moderate lagged relationships in the W model. This is a familiar type of cross-level interaction that in CLPMs involves a B variable moderating a lagged W effect. An example from the emotions literature is that stable trait neuroticism (B) is associated with stronger persistence of negative emotions over time (i.e., slower reversion to the mean), such that the W AR terms for negative emotions are larger at higher levels of the stable B variable neuroticism (Koval et al., 2016; Suls et al., 1998). Similarly, the stable B part of negative emotions measured over time may moderate its own W AR effect (for insight see Figure 3 in Hamaker et al., 2018). Cross-level interactions may be present elsewhere, such as in Zablah et al. (2016), if the lagged W effect of customer satisfaction ( x W i t − 1 ) on employee job satisfaction ( y W i t ) were moderated by the stable B part of customer engagement ( z B i ) , which these authors tested as a moderator.

OPEN IN VIEWER

Figure 3.

Graphic representation of the full model with estimates.

For this, the B and W parts of observed variables must be decomposed or else the product terms will conflate W and B variation when estimating interactions (Preacher et al., 2016). To illustrate a solution to this problem, below we study how stable B gender ( z B i ) moderates the W cross-lagged effect of work-family conflict ( x W i t − 1 ) on job satisfaction ( y W i t ) , an effect which was previously supported with models that likely conflated multilevel effects (e.g., Grandey et al., 2005). Furthermore, we show how the stable B parts of job satisfaction and work-family conflict can also moderate the W lagged effects among the same variables in a CLPM—again, this represents new opportunities as all published CLPMs can be re-analyzed to test for cross-level interactions among the modeled variables.

To formalize this familiar type of cross-level moderation we start with Eqs. 3 and 4, adding latent interactions among B and W terms (see Figure 2b). The equations are rather cumbersome, so for concision we show only the model for y. For clarity, we separate the purely W effects (including the WxW term) from the cross-level terms wherein a B variable moderates AR and CL effects in the W model. For the cross-level terms, Equation (11) describes a single B variable z moderating both W effects:

y W i t = β W y y ⋅ y W i t − 1 + β W x y ⋅ x W i t − 1 + u W i t y + z B i ( β B z W y y ⋅ y W i t − 1 + β B z W x y ⋅ x W i t − 1 )

(11)

In Equation (11), we use subscripts B and W on the cross-level interaction effect β B z W x y to indicate the stable B and dynamic W predictors involved. In this equation, the B variable z moderates: (1) the AR term, with the level of persistence in y differing across values of the B variable (i.e., β B z W y y ); and (2) a CL term, with the lagged effect of one W variable on the other differing across values of the B variable (i.e., β B z W x y ). Therefore, Equation (11) shows the logic of a stable B variable measured at Level-2, such as gender ( z B i ) , moderating the W effect of one dynamic variable, such as work-family conflict ( x W i t − 1 ) , on another dynamic variable, such as job satisfaction ( y W i t ) . However, because every CLPM involves at least two variables with repeated measurements, Equation. (11) can be expanded to include BxW interactions involving the latent B components (i.e., stable parts of x and y) in addition to an observed z B i with:

y W i t = β W y y ⋅ y W i t − 1 + β W x y ⋅ x W i t − 1 + u W i t y + y B i ( β B y W y y ⋅ y W i t − 1 + β B y W x y ⋅ x W i t − 1 ) + x B i ( β B x W y y ⋅ y W i t − 1 + β B x W x y ⋅ x W i t − 1 ) + z B i ( β B z W y y ⋅ y W i t − 1 + β B z W x y ⋅ x W i t − 1 )

(12)

Here, Equation (12) shows how the stable B parts of all variables can, in theory, be used to specify and test cross-level moderation (e.g., higher stable levels of work-family conflict x B i might sensitize people to the same variable's dynamic changes over time, thus increasing the CL effect β W x y of time-varying work-family conflict x W i t − 1 on job satisfaction y W i t via the cross-level interaction β B x W x y ). Of course, the WxW interaction effect in Equation (10) could also be moderated by a stable B variable, such as x B i ⋅ x W i t − 1 ⋅ y W i t − 1 with a coefficient β B x W x y y , but this type of three-way BxWxW interaction is beyond the scope of our current paper.

Here, we focus on interpreting the effects in Equation (12) for readers who prefer either a traditional interaction/moderation depiction of coefficients or a multilevel random-slope depiction of cross-level interaction effects (as in Raudenbush & Bryk, 2002). We begin with a traditional depiction wherein ‘main effects’ and interactions are grouped as a compound coefficient on each primary lagged W predictor in the model, showing interactions separately for AR and CL terms as follows:

y W i t = ( β W y y + β B y W y y ⋅ y B i + β B x W y y ⋅ x B i + β B z W y y ⋅ z B i ) y W i t − 1 + ( β W x y + β B y W x y ⋅ y B i + β B x W x y ⋅ x B i + β B z W x y ⋅ z B i ) x W i t − 1 + u W i t y

(13)

The model in Equation (13) is identical to the one in Equation (12), but the implied modification of the AR and CL terms is made clearer with a familiar moderated regression style of presentation, wherein moderation by B terms is shown as modifying the coefficients within parentheses on lagged predictors.

Alternatively, the same BxW terms can be shown in a multilevel random-slope style as follows, with arbitrary but conceptually useful equations for the AR and CL coefficients:

y W i t = A R i ⋅ y W i t − 1 + C L i ⋅ x W i t − 1 + u W i t y

(14a)

A R i = β W y y + β B y W y y ⋅ y B i + β B x W y y ⋅ x B i + β B z W y y ⋅ z B i

(14b)

C L i = β W x y + β B y W x y ⋅ y B i + β B x W x y ⋅ x B i + β B z W x y ⋅ z B i

(14c)

In Equations (14a) to (14c), each slope is treated as varying (as in Raudenbush & Bryk, 2002; Yuan et al., 2014), but its variation is treated as fixed in the predictors rather than being given a random residual term. In other words, the ‘random’ part of a random slope is not present. This reduces model complexity and facilitates estimation particularly with small-T cases, like the CLPM wherein such a random variance would be estimated with low levels of precision. Instead, only the relevant fixed effects are estimated so that the AR and CL terms in the W model are treated as varying across the B variables, but not ‘randomly’ varying as in a random-slope model—note that we avoid the ‘Level-1’ and ‘Level-2’ nomenclature here to stay focused on the W and B effects (as in Preacher et al., 2010, 2016; Zhang et al., 2009).

This formulation clarifies the typical ‘slopes-as-outcome’ interpretation of multilevel interactions wherein the AR and CL coefficients β W y y and β W x y reflect the average effects across all N units of analysis when the interacting B predictors are mean-centered (i.e., y B i , x B i , z B i ). However, when the interacting B predictors are not mean-centered, the AR and CL coefficients are treated like random-slope intercepts, such that an average β W y y is equal to β W y y − β B y W y y ⋅ M y B i − β B x W y y ⋅ M x B i − β B z W y y ⋅ M z B i ; which shows why mean-centered predictors are often useful. Conveniently, in our SEM approach all latent B predictors are mean-centered by default because the mean structure will be accounted for by observed-variable intercepts a t . When using observed-variable B predictors such as gender, these can either be manually mean-centered or the implied averages can be derived using the equation for M as noted (assuming that the overall averages of AR and CL terms are of interest).

As our approach thus far makes clear, not only is it possible to conceptualize W interactions in CLPMs, but it is also possible to formulate cross-level BxW interactions because panel data imply B components that may serve as moderators of W lagged effects. Seen this way, the bulk of the literature on multilevel moderation appears relevant to CLPMs including for interpreting effects. For example, when B moderators are mean-centered the ‘main effects’ β W y y and β W x y are the average AR and CL effects across all N panels as in Equations (3) and (4). Yet, these may vary across the levels of stable B variables—here y B i , x B i , and z B i —which can be easily estimated by latent interactions among B and lagged W variables.

Purely Between-Level or BxB Interactions

Many researchers propose effects among variables that are conceptualized as having strong stable components, such as personality and job performance (see Barrick & Mount, 1991; Tett et al., 1991). In these cases, it is possible to test interactions among B variables. For example, Zablah et al. (2016) could have focused their theorizing on B effects, such that the stable B part of customer satisfaction predicts the stable B part of job satisfaction, and this effect is moderated by the stable B part of customer engagement. This model would describe a purely B form of moderation, which we refer to as a BxB interaction. This type of interaction can involve observed or latent variables at the B level. For example, gender could be treated as an observed B moderator ( z B i ) of the relationship between the stable B part of work-family conflict ( x B i ) and the stable B part of job satisfaction ( y B i ) .

With this kind of approach, the B model for y in Equation (7) could be shown as involving a B moderator z B i to represent a BxB interaction as follows (see Figure 2c):

y B i = β B x y ⋅ x B i + β B z y ⋅ z B i + β B x z y ⋅ x B i ⋅ z B i + ζ B i y

(15)

In Equation (15), we use the term β B x z x to indicate a purely B interaction effect among the stable components of x and z. In our example, the stable B part of work-family conflict x B i may interact with gender z B i to predict the stable B part of job satisfaction y B i . Notably, even if z B i is observed (such as gender, self-reported without any measurement error assumed), a latent interaction needs to be specified if one component (in this case x B i ) is latent.

Nonlinearity and Higher-Order Interactions

With these three types of interactions, more complex and interesting pictures of dynamic systems are available for theory testing (see Figure 2d for all interactions in one model). Yet, our approach also allows incorporating nonlinearity and higher-order interactions. Consider Lian et al.'s (2014) three-way interaction wherein abusive supervision, self-control capacity, and intention to quit at t = 1 were modeled as interacting to predict organizational deviance at t = 2. Their model focused on W relationships but did not account for stable B components. With our approach, building on Eqs. 9 and 10, it is possible to model pure WxWxW interactions of three variables by forming latent W interactions in a CLPM. A WxWxB interaction could also be modeled by using the B part of a variable and forming a three-way interaction with two W lagged predictors.

As another example, consider Lin et al. (2017) who studied a type of nonlinear moderation using squared predictors multiplied by a moderator. Specifically, perceived underemployment at t = 1 and task crafting at t = 2 had an inverted U-shaped relationship (modelled by including a squared term for underemployment), which was moderated by organizational identification at t = 1. Using our approach (assuming a sufficient number of measurement occasions) it would be possible to model the purely W nonlinear effect of underemployment on task crafting as a latent WxW interaction of underemployment with itself (what we call a W 2 model term), which could then be moderated by organizational identification's W or B part using a higher-order latent interaction—all while controlling for the stable B parts of the variables modeled.

Although the possibilities for testing nonlinear and higher-order interactions are essentially unlimited, we point out five potential cases that should be of primary interest and rather intuitive for most researchers. Notably, our examples here contain a variety of different types of interactions and nonlinear effects, not all of which may be of interest in any given application, but which we show for completeness. The first is a purely W quadratic effects or W 2 case that may include AR and CL terms using squared W lagged predictors as follows:

x W i t = β W x x ⋅ x W i t − 1 + β W x 2 x ⋅ x W i t − 1 2 + β W y x ⋅ y W i t − 1 + β W y 2 x ⋅ y W i t − 1 2 + u W i t x

(16)

y W i t = β W y y ⋅ y W i t − 1 + β W y 2 y ⋅ y W i t − 1 2 + β W x y ⋅ x W i t − 1 + β W x 2 y ⋅ x W i t − 1 2 + u W i t y

(17)

In Equations. (16) and (17), the coefficients β W x 2 x and β W y 2 y imply nonlinear AR effects, and coefficients β W y 2 x and β W x 2 y imply nonlinear CL effects. In the context of AR effects, this implies nonlinear persistence in a variable over time. For example, intense negative emotions may show higher inertia (i.e., greater persistence in the form of larger AR effects) because they are more difficult to down-regulate than mild negative emotions. Alternatively, consider the CL case where a W 2 term implies nonlinear lagged effects of a predictor. As we show below, these can be investigated by forming latent interactions among lagged W terms with themselves.

Next is the case of a quadratic B effect, which we refer to as B 2 , that can be shown by elaborating on the B model for y from Equation (7) as follows:

y B i = β B x y ⋅ x B i + β B x 2 y ⋅ x B i 2 + β B z y ⋅ z B i + ζ B i y

(18)

In Equation (18), β B x 2 y represents a quadratic effect of the stable B part of x. As an example, from Lin et al. (2017), if they had longitudinal panel data this could be used to test a nonlinear effect of the stable B part of perceived underemployment ( x B i ) on the stable part of task crafting ( y B i ). This relationship could also be moderated by the stable B part of organizational identification ( z B i ).

Now consider three unique cross-level moderation cases that combine the logic of the two previous W 2 and B 2 cases. We will refer to these as B 2 xW, Bx W 2 , and B 2 x W 2 cases. In the first B 2 xW case, a quadratic B term moderates a linear W relationship. We start with the purely B 2 (quadratic) relationship in Eq. 18, adding a purely W linear relationship between y and x from Eq. 12 to form the final B 2 xW relationship as follows:

y W i t = β W y y ⋅ y W i t − 1 + β W x y ⋅ x W i t − 1 + u W i t y + y B i ( β B y W y y ⋅ y W i t − 1 + β B y W x y ⋅ x W i t − 1 ) + x B i ( β B x W y y ⋅ y W i t − 1 + β B x W x y ⋅ x W i t − 1 ) + x B i 2 ( β B x 2 W y y ⋅ y W i t − 1 + β B x 2 W x y ⋅ x W i t − 1 ) + z B i ( β B z W y y ⋅ y W i t − 1 + β B z W x y ⋅ x W i t − 1 )

(19)

In Equation (19), the nonlinear cross-level B 2 xW moderation is represented by β B x 2 W y y and β B x 2 W x y . For example, if the stable part of work-family conflict ( x B i ) were posited to have a curvilinear moderating effect on the W effect of work-family conflict ( x W i t − 1 ) on job satisfaction ( y W i t ), this can be modeled in our approach using the term β B x 2 W x y in Equation (19).

In the second case of nonlinear Bx W 2 moderation, consider a nonlinear W 2 relationship that is moderated by the stable B parts of modeled variables. We start by writing purely W 2 relationship from Equation (17) and add cross-level interaction effects as follows:

y W i t = β W y y ⋅ y W i t − 1 + β W y 2 y ⋅ y W i t − 1 2 + β W x y ⋅ x W i t − 1 + β W x 2 y ⋅ x W i t − 1 2 + u W i t y + y B i ( β B y W y y ⋅ y W i t − 1 + β B y W y 2 y ⋅ y W i t − 1 2 + β B y W x y ⋅ x W i t − 1 + β B y W x 2 y ⋅ x W i t − 1 2 ) + x B i ( β B x W y y ⋅ y W i t − 1 + β B x W y 2 y ⋅ y W i t − 1 2 + β B x W x y ⋅ x W i t − 1 + β B x W x 2 y ⋅ x W i t − 1 2 ) + z B i ( β B z W y y ⋅ y W i t − 1 + β B z W y 2 y ⋅ y W i t − 1 2 + β B z W x y ⋅ x W i t − 1 + β B z W x 2 y ⋅ x W i t − 1 2 )

(20)

In Equation (20), the cross-level B moderation effects of the nonlinear W 2 terms are β B y W y 2 y , β B y W x 2 y , β B x W y 2 y , β B x W x 2 y , β B z W y 2 y , and β B z W x 2 y . These terms imply that nonlinear W effects may vary across the levels of the B parts of the variables x, y, and z. For example, if a curvilinear effect of work-family conflict ( x W i t − 1 2 ) on job satisfaction ( y W i t ) were theorized to be moderated by gender ( z B i ), this could be tested as β B z W x 2 y , which might show that family-work conflict has a stronger effect on job satisfaction at higher levels of family-work conflict, but only for women.

In the third case of nonlinear cross-level moderation, B 2 x W 2 , we extend Equation (20) to the case of a purely W quadratic relationship between x and y moderated with a quadratic effect of the B part of a predictor variable z:

y W i t = β W y y ⋅ y W i t − 1 + β W y 2 y ⋅ y W i t − 1 2 + β W x y ⋅ x W i t − 1 + β W x 2 y ⋅ x W i t − 1 2 + u W i t y + y B i ( β B y W y y ⋅ y W i t − 1 + β B y W y 2 y ⋅ y W i t − 1 2 + β B y W x y ⋅ x W i t − 1 + β B y W x 2 y ⋅ x W i t − 1 2 ) + x B i ( β B x W y y ⋅ y W i t − 1 + β B x W y 2 y ⋅ y W i t − 1 2 + β B x W x y ⋅ x W i t − 1 + β B x W x 2 y ⋅ x W i t − 1 2 ) + z B i ( β B z W y y ⋅ y W i t − 1 + β B z W y 2 y ⋅ y W i t − 1 2 + β B z W x y ⋅ x W i t − 1 + β B z W x 2 y ⋅ x W i t − 1 2 ) + z B i 2 ( β B z 2 W y y ⋅ y W i t − 1 + β B z 2 W y 2 y ⋅ y W i t − 1 2 + β B z 2 W x y ⋅ x W i t − 1 + β B z 2 W x 2 y ⋅ x W i t − 1 2 )

(21)

In Equation (21), the nonlinear cross-level moderation effects of interest are β B z 2 W y 2 y and β B z 2 W x 2 y . This model allows for a quadratic cross-level moderating effect of z B i 2  on the quadratic W relationship between work-family conflict x W i t − 1 2 and job satisfaction y W i t .

In all cases above, our approach decomposes the dynamic W and stable B parts of observed variables before forming product terms. This yields unbiased estimates of same-level and cross-level effects in CLPMs, while helping researchers focus on theory and hypothesis tests that are sensitive to interactions at different levels of analysis. As a beneficial byproduct, any heteroskedasticity that would otherwise exist in residuals can be accounted for by the interaction effects in the model (Raudenbush & Bryk, 2002). Before continuing, however, we would emphasize that simply because so many interactions and nonlinearities can be modeled does not imply that they should be. Indeed, a benefit of ‘unconflating’ B and W effects is that researchers need only model effects that are theoretically relevant.

Confounding in Tests of Interaction Effects

Given our presentation thus far, it is now straightforward to show how confounding occurs when variables are not disaggregated into B and W parts. For example, when the product of x and y is included in a typical CLPM, the net result can be shown as follows for the y variable (again, temporarily ignoring intercepts or an ‘occasion effect’ α t ):

y i t = β y y y i t − 1 + β x y ⋅ x i t − 1 + β x y y ⋅ x i t − 1 ⋅ y i t − 1 + e i t y

(22)

where the product of x and y is lagged and the coefficient β x y y is meant to indicate the interaction effect or degree of moderation for either the AR or CL term (as is well known, either x or y can be understood as the moderator of the other's effect in such a case).

However, the multilevel decomposition in Equations (1) and (2) suggests that this formulation has various potential shortcomings, which we help illustrate by using a compound residual term e i t y rather than the separate B and W parts ζ B i y and u W i t y , respectively. To further show the problems here, we can substitute the multilevel B and W counterparts in Eq. 22 as follows:

y i t = β y y ( y B i + y W i t − 1 ) + β x y ( x B i + x W i t − 1 ) + β x y y ( x B i + x W i t − 1 ) ( y B i + y W i t − 1 ) + e i t y

(23)

In Equation (23), the first two terms β y y ( y B i + y W i t − 1 ) and β x y ( x B i + x W i t − 1 ) exemplify the well-known multilevel problem of conflating potential B and W effects (for general insight see Preacher et al., 2010; Zhang et al., 2009; for the CLPM case specifically see Hamaker et al., 2015; Zyphur, Allison, et al., 2020; Zyphur, Voelkle, et al., 2020).

By not decomposing y and x into their level-specific components, interaction terms that conflate B and W parts can be shown by expanding the β x y y product term as follows:

y i t = β y y ( y B i + y W i t − 1 ) + β x y ( x B i + x W i t − 1 ) + β x y y ⋅ x B i ⋅ y B i + β x y y ⋅ x B i ⋅ y W i t − 1 + β x y y ⋅ y B i ⋅ x W i t − 1 + β x y y ⋅ x W i t − 1 ⋅ y W i t − 1 + e i t y

(24)

As we already showed above, Eq. 24 makes it apparent that the single coefficient β x y y actually represents four separate interactions: a purely between-level or BxB term ( x B i ⋅ y B i ), two cross-level or BxW terms ( x B i ⋅ y W i t − 1 ) and ( y B i ⋅ x W i t − 1 ), and a purely within-level or WxW term ( x W i t − 1 ⋅ y W i t − 1 ). By extension, as implied by Eq. 24, adding nonlinearity and/or three-way interactions would exacerbate the interpretational difficulties and degrees of confounding. In what follows, we illustrate our proposed latent interaction CLPM approach by applying it to a longitudinal panel dataset. Following this, we report the results of simulations that lend support for our proposed approach with latent interactions by comparison with observed-means centering and an uncentered approach (i.e., raw variables).

Method

Participants

Data were taken from the Swiss Household Panel (SHP) survey, which is a nationally representative random sample measured annually, beginning in 1999 with around 5,000 Swiss households and adding two additional samples in 2004 and 2013 (see Voorpostel et al., 2018). Individual household members are surveyed via telephone on a broad range of topics including health, employment, and leisure activities. The SHP has been used extensively in fields such as public health (e.g., Knecht et al., 2011), sociology (e.g., Oesch & Lipps, 2013), and economics (e.g., Frei & Sousa-Poza, 2012).

In the present study, we used data collected in the years 2004 through 2007. We chose 2004 as the onset of our study because this year marked the first that all study variables were measured annually, and the first year that included a second SHP sample, which increased our sample size. Although the SHP's sample included around 10,000 working people in 1999, due to attrition this number decreased to 6,000 in 2004. The addition of a second sample in 2004 increased the number of working participants to over 13,000. Furthermore, the waves 2004 to 2007 consistently demonstrated response rates of more than 70%.

We included all working individuals that responded to at least one of the four waves surveyed between 2004 and 2007. This brought our sample size to 7,748 individuals. The participants in our sample (51.6% women, 48.4% men) were on average 39.28 (SD = 14.50) years old and had 12.98 (SD = 3.20) years of education.

Measures

Descriptive statistics are reported in Table 1. We used measures of job satisfaction and work-family conflict. Past research shows that behavior-based types of work-family conflict (i.e., behavioral incompatibilities of different roles) have the strongest relationship to job satisfaction, particularly composite-based satisfaction (averaged across facets) rather than global measures (Bruck et al., 2002).

OPEN IN VIEWER

Table 1.

Descriptive Statistics and Correlations of Study Variables.

Variable		Mean	σ2	1	2	3	4	5	6	7	8	9	10	11	12
1	J9	7.825	1.533
2	J8	7.826	1.503	.634
3	J7	7.908	1.451	.568	.630
4	J6	7.966	1.507	.528	.559	.622
5	C9	3.802	7.620	−.288	−.236	−.201	−.189
6	C8	3.827	7.476	−.191	−.231	−.183	−.177	.546
7	C7	3.721	7.766	−.197	−.198	−.245	−.202	.489	.539
8	C6	3.652	8.218	−.183	−.183	−.176	−.219	.426	.462	.476
9	H9	13.731	6.918	−.059	−.065	−.078	−.043	.326	.294	.282	.245
10	H8	13.735	7.025	−.061	−.062	−.066	−.048	.281	.323	.317	.260	.834
11	H7	13.769	6.955	−.056	−.058	−.068	−.039	.267	.298	.347	.280	.781	.864
12	H6	13.804	6.723	−.071	−.077	−.088	−.062	.260	.280	.311	.316	.723	.808	.852
13	G	0.506	0.250	.020	.016	.035	.042	−.111	−.111	−.139	−.101	−.349	−.364	−.364	−.368

Note: N = 7,748 from the Swiss Household Panel (SHP) with variables named in accordance with the SHP waves; J9-J6 = job satisfaction measured in 2007–2004; C9-C6 = work-family conflict measured in 2007–2004; H9-H6 = hours worked measured in 2007–2004; G = gender coded 0 for men and 1 for women; σ² = estimated population variance.

Job Satisfaction

Participants indicated the extent to which they were satisfied with their job by responding to six items created for the SHP (0 = not at all satisfied to 10 = completely satisfied). The lead-in question was “Can you indicate your degree of satisfaction for each of the following points?” The six items were “your job in general”, “the level of interest in your tasks”, “the amount of your work”, “the income you get from your job”, “your work conditions”, and “the atmosphere between you and your colleagues.” For simplicity, we computed observed-variable mean scores across these six job satisfaction items. However, alternative higher-order latent variable models would also be possible given the general flexibility of SEMs (see Mulder & Hamaker, 2021). The alpha coefficients for the four waves were 0.78, 0.78, 0.79, and 0.79, respectively.

Work-Family Conflict

Participants rated the extent to which their work negatively impacted their family life by a single-item measure: “How strongly does your work interfere with your private activities and family obligations more than you would want this to be?” (0 = not at all to 10 = very strongly).

Hours Worked

Participants reported the number of hours worked per week, which we included as an observed variable in our model with the same W and B decomposition to control for hours worked per week when interpreting all effects. Before analysis we multiplied this variable by 4 to improve convergence when setting observed variances to small values as noted in the next section.

Gender

In every wave, participants indicated their gender using a dummy variable (0 = man, 1 = woman). Responses across all participants were invariant across all years, allowing this to serve as a purely B observed variable in our CLPM as described next.

Specification and Estimation

We first offer a baseline model decomposing W and B parts of observed variables (as shown in Equations. (1) and (2)). This model includes lagged effects for the W parts (as in Equations (5)and (6)); a B model with stable B job satisfaction predicted by stable B work-family conflict, gender, and hours worked (as in Eq. 7); stable B work-family conflict predicted by gender and hours worked; and gender predicting B hours worked. Details on this SEM specification can be found in Hamaker et al. (2015), including the addition of a time-varying intercept α t for each observed variable by default (typical in SEM). To elaborate on this, we estimated four other models that included interaction effects: within-level (WxW; as in Equations (9) and (10)); cross-level (BxW; as in Equations (12) to (14)); between-level (BxB; as in Equation (15)); and all simultaneously (WxW, BxB, and BxW). For comparison, we also estimated a traditional CLPM without the B and W decomposition so that product terms are formed with observed-variable lags, which we term the ‘conflated model’ below. All Mplus files and a summary of the results can be found in our Online Appendices.

To estimate our models, we considered using latent moderated structural equations or LMS (e.g., Preacher et al., 2016), but this is computationally infeasible in high dimensions (Zyphur et al., 2018). This problem is addressed with the Bayes estimator in Mplus, which provides optimal estimates and faster compute times using parallel processing (Asparouhov & Muthén, 2020). To do this requires first checking compute times for 400 iterations with different numbers of physical (rather than logical or ‘virtual’) processing cores, which we did for up to 48 cores using Amazon Web Services (AWS) ‘c5.24xlarge’ EC2 instances (with 48 physical and 96 virtual cores). We found that, in a Windows environment, 4 or 8 cores were optimal for models with fewer interactions, but with many interactions 16 cores was best. However, in the same AWS platform, we found that in a Linux environment estimation scaled much better with higher core counts. For example, estimating the WxW model below with 60,000 total iterations (as we note below, thinning every other iteration and requesting 30,000 iterations), we observed the following estimation times with our comparatively large model: 2 cores = 25.92 h; 4 cores = 13.33 h; 8 cores = 7 h; 16 cores = 3.58 h; 32 cores = 1.83 h; 48 cores = 1.36 h. As such, we recommend using a Linux environment rather than Windows, particularly given the Linux instance is roughly half price on AWS (∼$4.50 hly vs. Windows at ∼$9.00 hly).

For all models we checked convergence using posterior scale reduction (PSR) values < 1.05, at least doubling iterations after this to ensure stable convergence (Asparouhov & Muthén, 2010; Zyphur & Oswald, 2015). For this we were conservative and estimated more than will typically be required for other researchers using similar data, for example thinning every other iteration with 30,000 iterations set, for a total of 60,000 in the WxW case. To aid in convergence, it is notable that decomposing the W and B parts of observed variables requires accounting for all variance in an observed variable, which implies zero residual variances. This causes slow convergence with Bayes estimators, so we set residual variances to roughly 1–3% of the total observed variation to improve convergence (Asparouhov & Muthén, 2020), while allowing the decomposed W and B variances to be freely estimated.

Results for Real-Data Example

All model results are shown in Table 2 in a side-by-side format to allow easy comparisons of parameter estimates across models. We present results for each model separately, starting with the baseline model to help draw comparisons with those that follow, including interaction effects and finally the conflated model that does not decompose B and W parts with observed-variable products for interactions. Standardized effects can be computed using formulas presented in Asparouhov and Muthén (2020), but for brevity we present raw coefficients, as is often done in multilevel literature testing interactions. For brevity, Figure 3 contains results only for the full model with all effects except for the paths involving the control variable (hours worked) and with only two occasions of measurement shown to aid readability. Also, all reported p-values are two-tailed and calculated using Bayesian posteriors SDs (rather than their frequentist counterpart SEs).

OPEN IN VIEWER

Table 2.

Summary of Model Estimates.

Variables(1)	Baseline Model	WxW Model	BxB Model	BxW Model	Full Model	Conflated Model Variables	Conflated Model
	Estimates	Estimates	Estimates	Estimates	Estimates		Estimates
W Effects						B + W Effects
AR effect CW	0.086***	0.087***	0.088***	0.206***	0.214***	AR effect for C	0.465***
AR effect JW	0.203***	0.197***	0.204***	0.395***	0.395***	AR effect for J	0.639***
AR effect Hw	0.493***	0.507***	0.498***	0.557***	0.558***	AR effect for H	0.816***
CL effect CW→JW	0.006	0.005	0.006	0.071***	0.066***	CL effect C→J	0.006
CL effect JW→CW	−0.079	−0.064	−0.075	0.229***	0.207***	CL effect J→C	−0.202***
CL effect HW→CW	0.140***	0.156***	0.146***	0.236***	0.239***	CL effect H→C	0.136***
CL effect HW→JW	−0.005	−0.002	−0.005	0.027**	0.023**	CL effect H→J	−0.007
CL effect CW→HW	0.020	0.025*	0.021*	0.063***	0.065***	CL effect C→H	0.020***
CL effect JW→HW	0.046	0.053	0.046	0.121***	0.115***	CL effect J→H	0.005
WxW Interactions						Interactions
CW*JW→CW		0.044**			0.076***	C*J→C	−0.001
CW*JW→JW		0.030**			0.066***	C*J→J	−0.003
B Effects						Effects of Gender
CB→JB	−0.196***	−0.200***	−0.176***	−0.392***	−0.391***
GB→JB	0.000	−0.002	−0.011	−0.174***	−0.182***	GB→J	−0.010
HB→JB	0.027**	0.026**	0.029**	−0.003	0.002
GB→CB	0.077	0.072	0.074	−0.080	−0.095	GB→C	0.034*
HB→CB	0.357***	0.354***	0.356***	0.269***	0.260***
GB→HB	−2.000***	−1.999***	−2.106***	−1.987***	−1.987***	GB→H	−0.407***
BxB Interaction
CB*GB→JB			−0.045		−0.047***
BxW Interactions
CB*CW→CW				−0.179***	−0.218***
JB*CW→CW				−0.451***	−0.534***
GB*CW→CW				−0.132***	−0.130***
CB*JW→CW				0.121**	−0.021
JB*JW→CW				0.098	−0.175
GB*JW→CW				0.072	0.013
CB*CW→JW				−0.077***	−0.096***
JB*CW→JW				−0.239***	−0.283***
GB*CW→JW				−0.041**	−0.045***
CB*JW→JW				0.365***	0.373***
JB*JW→JW				0.837***	0.837***
GB*JW→JW				0.110**	0.096**

Note: All estimates are unstandardized.

⁽¹⁾

W = Within; B = Between; J = Job satisfaction; C = Work-family conflict; G = Gender; H = Hours worked; Arrows represent direction.

*p < 0.05, **p < 0.01, ***p < 0.001. Source: Swiss Household Panel (SHP).

Baseline Model. For this model, job satisfaction (J), work-family conflict (C), and hours worked (H) are decomposed into stable B and dynamic W parts with latent variables, whereas gender (G) is an observed B variable (to benefit the reader these variable labels are also in our online Mplus files). We use H as a main-effect control (i.e., no interactions) but decompose it into W and B parts to properly account for separate W and B associations with J and C. Notably, typical fit statistics are not available in the presence of latent interactions, but baseline models without interactions have fit equivalent to interaction models (Asparouhov & Muthén, 2020). The baseline model showed adequate fit, but with typical problems of absolute misfit using χ² metrics coupled with very good fit for other indices as follows: χ² difference for observed versus model-generated data shows a 95% CI ranging from 21.296 to 104.514 (excluding zero, indicating poor fit) with the number of parameters pD = 53.855; posterior predictive p-value or PPP = .001 (smaller than .05, indicating poor fit); but Bayes-equivalent RMSEA = .006 (smaller than .08, indicating good fit); and both CFI and TLI = .999 (larger than .95, indicating good fit); and DIC = 218442.103.

The following CLPM was implemented in Mplus (see the online files Baseline CLPM.inp and .out) to compute baseline W and B parts and slope coefficients as follows:

J i t = α t J + β B ( C ) J ⋅ C B i + β B ( G ) J ⋅ G B i + β B ( H ) J ⋅ H B i + β W ( J ) J ⋅ J W i t − 1 + β W ( C ) J ⋅ C W i t − 1 + β W ( H ) J ⋅ H W i t − 1 + ζ B i J + u W i t J

(25)

C i t = α t C + β B ( G ) C ⋅ G B i + β B ( H ) C ⋅ H B i + β W ( C ) C ⋅ C W i t − 1 + β W ( J ) C ⋅ J W i t − 1 + β W ( H ) C ⋅ H W i t − 1 + ζ B i C + u W i t C

(26)

H i t = α t H + β B ( G ) H ⋅ G B i + β W ( H ) H ⋅ H W i t − 1 + β W ( J ) H ⋅ J W i t − 1 + β W ( C ) H ⋅ C W i t − 1 + ζ B i H + u W i t H

(27)

The AR effects for work-family conflict, job satisfaction, and hours worked were significant with β W ( J ) J  = .203 (p < .001), β W ( C ) C  = .086 (p < .001), β W ( H ) H  = .493 (p < .001); but most CL effects were not: β W ( C ) J  = .006 (p = .354), β W ( J ) C  = −.079 (p = .076), β W ( H ) J  = −.005 (p = .73), β W ( C ) H  = .02 (p = .092), β W ( J ) H  = .046 (p = .19), except for one significant CL effect of hours worked on work-family conflict β W ( H ) C  = .14, (p < .001). The AR results imply that roughly 20% of a dynamic W perturbation to job satisfaction and 49% for hours worked persists from one occasion to the next, but only about 8% for work-family conflict. For CL effects, typical significance tests support only a direct effect of hours worked on work-family conflict.

Estimates for B terms supported a negative effect of average work-family conflict on average job satisfaction as β B ( C ) J  = −.196 (p < .001), but not for gender β B ( G ) J < 0.001 (p = .099). Gender was negatively associated with hours worked with β B ( G ) H  = −2.0 (p < .001) but not average work-family conflict β B ( G ) C  = .077 (p = .230). This indicates that women (Gender = 1) seem to work fewer hours but do not experience more work-family conflict. Average hours worked predicted average work-family conflict β B ( H ) C  = .357 (p < .001) and average job satisfaction β B ( H ) J  = .027 (p = .01). Taken together these B coefficients suggest that as a person works more, they experience more work-family conflict but also more job satisfaction, on average. Gender did not predict average job satisfaction or work-family conflict. If the B coefficients were causal, this would suggest that greater work-family conflict causes reduced job satisfaction, women work fewer hours, and hours worked increases work-family conflict, but working more also ironically increases job satisfaction.

In sum, the baseline model effects are important because researchers who use CLPMs are interested in knowing what dynamic effects are present in the W parts of their models. The B model is also relevant for those who are comfortable with their model specification and any inferences that follow from it. However, given the possibility of estimating different types of latent interactions, this typical CLPM may omit important moderation effects.

Within-Level Interaction (WxW). Our first interaction model includes W latent interactions of work-family conflict and job satisfaction as in Eqs. 9–10. This model allows studying whether AR and CL terms among the W dynamic variables depend on levels of the lagged variables. As we would expect, AR and CL estimates were very similar to the baseline model—because the W latent variables are mean-centered by default, the baseline model will reflect the average AR and CL terms in a sample as will any WxW interaction model. Specifically, AR terms were β W ( J ) J  = .197 (p < .001), β W ( C ) C  = .087 (p < .001), and β W ( H ) H  = .507 (p < .001), whereas CL terms were β W ( C ) J  = .005 (p = .464), β W ( J ) C  = −.064 (p = .12), β W ( H ) C  = .156, (p < .001), β W ( H ) J  = −.002 (p = .894), β W ( C ) H  = .025 (p = .042), and β W ( J ) H  = .053 (p = .08). Compared to the baseline model, including the WxW interaction terms reduced some uncertainty by making the CL effect of work-family conflict on hours worked significant, suggesting that a dynamic increase to work-family conflict may cause a slight increase in hours worked. Furthermore, the B estimates were also very similar to those obtained under the baseline model—which is expected because the B and W model parts are orthogonal and therefore adding WxW interactions should not markedly affect B estimates—with β B ( C ) J  = −.20 (p < .001), β B ( G ) J  = −.002 (p = .954), β B ( H ) J  = .026 (p = .008), β B ( G ) C  = .072 (p = .246), β B ( H ) C  = .354 (p < .001), and β B ( G ) H  = −1.999 (p < .001).

Moving on to the WxW interactions, we found statistically significant effects. First, when predicting work-family conflict, we found a positive interaction effect β W ( J C ) C  = .044 (p = .016) such that greater levels of lagged work-family conflict or job satisfaction increased the effect of the other on future work-family conflict. Figure 4 offers an example of how to visualize such interactions. This means that as work-family conflict increases, the CL effect of job satisfaction on conflict ( β W ( J ) C ) increases, and as job satisfaction increases, the AR effect of work-family conflict ( β W ( C ) C ) increases. Similarly, when predicting job satisfaction, we found a positive interaction effect β W ( J C ) J  = .03 (p = .02). This is plotted in Figure 5, which shows that greater levels of lagged job satisfaction increase the CL effect of conflict ( β W ( C ) J ), and greater levels of lagged conflict increase the AR effect of job satisfaction ( β W ( J ) J ). Although either variable could be theoretically seen as a moderator, our point is that these interaction effects are of potential interest to researchers who can always estimate them in CLPMs. Our results here show that the dynamic W parts of work-family conflict and job satisfaction appear to have implications for the other's dynamic effects.

OPEN IN VIEWER

Figure 4.

Plot of The WxW model's interaction effect of work-family conflict (C) and Job satisfaction (J). Note: Plotted using only AR and CL effects for simplicity, showing how the W part of C and the W part J can interact to influence the future state of the W part of C (from the model estimates: C _{Wt +1} = .087 * C _Wt − .064 * J _Wt  + .044 C _Wt * J _Wt ).

OPEN IN VIEWER

Figure 5.

Plot of The WxW model's interaction effect of work-family conflict (C) and Job satisfaction (J). Note: Plotted using only AR and CL effects for simplicity, showing how the W part of C and the W part J can interact to influence the future state of W part of J (from the model estimates: J _{Wt +1} = .197 * J _Wt  + .005 * C _Wt  + .03 C _Wt * J _Wt ).

Between-Level Interactions (BxB). Here we include B interactions as in Eq. 15. As Table 2 shows, AR and CL effects were again essentially unchanged compared to the baseline because the B and W models are orthogonal to each other, but the CL term for work-family conflict on hours worked became significant, as in the WxW model. As Table 2 shows, for the AR terms β W ( J ) J  = .204 (p < .001), β W ( C ) C  = .088 (p < .001), and β W ( H ) H  = .498 (p < .001), whereas for the CL terms β W ( C ) J  = .006 (p = .334), β W ( J ) C  = −.075 (p = .082), β W ( H ) C  = .146 (p < .001), β W ( H ) J  = −.005 (p = .698), β W ( C ) H  = .021 (p = .044), and β W ( J ) H  = .046 (p = .166).

Again, as in the WxW model, here the B effects were very similar to the baseline, with average hours worked associated with increased job satisfaction, β B ( H ) J  = .029 (p = .002), and average work-family conflict, β B ( H ) C  = .356 (p < .001). Being a woman was associated with fewer hours worked, β B ( G ) H  = −2.106 (p < .001). The key difference in this model is the B interaction among average work-family conflict and gender predicting job satisfaction, which was negative and significant, β B ( C G ) J  = −.045 (p = .008). By including this interaction with the gender variable, the significant effect of average work-family conflict on average job satisfaction for men (G = 0) was β B ( C ) J  = −.176 (p < .001). In turn, we can compute this coefficient for women by adding the interaction effect to net −.045 + (−.176) = −.221. This is sensible given that the estimated coefficient in the baseline model, which averages across men and women, is approximately the average of these two coefficients at −.196. The implication is that, on average over time, women appear to have a stronger association between work-family conflict and job satisfaction (see e.g., Grandey et al., 2005).

Finally, the estimated effect of gender on average job satisfaction was similar yet slightly smaller, β B ( G ) J  = −.011 (p = .694). However, in this model, this represents the effect of gender on average job satisfaction when average work-family conflict = 0, which is the grand-mean of the latent B work-family variable. However, in our model all latent variables have means of zero. Thus, this effect has a similar interpretation as in the baseline model. Of course, although not statistically significant, the interaction effect could imply that average work-family conflict changes the effect of gender on job satisfaction, making job satisfaction even lower for women by −.045 as when work-family conflict increases by 1 unit. Furthermore, as we noted previously, we are reticent to give causal meanings to the B coefficients, but we offer this as an example.

Cross-level Interactions (BxW). In this model, the latent B and W interaction effects among work-family conflict, gender, hours worked, and job satisfaction are tested based on Eqs. 12–14. As noted, in these models the interpretation of the AR and CL coefficients is equivalent to a ‘slopes-as-outcomes’ model in the multilevel context (e.g., Raudenbush & Bryk, 2002), with W estimates such as an AR term β W ( J ) J being equal to the estimated average W effect across all individuals when all B interacting variables are mean-centered. The stable B parts of work-family conflict ( C B i ) and job satisfaction ( J B i ) are mean-centered by default because means are accounted for by observed-variable intercepts α t . However, gender is coded 0 = man and 1 = woman, so all W estimates of AR and CL terms are for men, and cross-level interactions with gender indicate the difference in AR and CL terms for women (compared to estimates for men). By accounting for such cross-level interactions, we expect more efficient estimates and therefore smaller p-values for W terms.

Starting with AR terms, we estimated β W ( J ) J  = .395 (p < .001), β W ( C ) C  = .206 (p < .001), and β W ( H ) H  = .557 (p < .001), which is similar to our baseline estimates for hours worked but roughly double for job satisfaction and work-family conflict, indicating greater dynamic persistence (i.e., slower regression to the mean). For women, the cross-level interactions with gender for job satisfaction was β B ( G ) W ( J ) J  = .11 (p = .002), implying slower regression to the mean in job satisfaction for women as their AR term is .11 + .395 = .505. But, for work-family conflict the estimate was β B ( G ) W ( C ) C  = −.132 (p < .001), which implies an AR term for women of −.132 + .206 = .075, meaning only 7.5% of a temporary increase to work-family conflict persists from one year to the next. In sum, women tend to regress to their long-run average job satisfaction more slowly than men, but women tend to regress to their long-run average for work-family conflict faster than men after experiencing a temporary change (i.e., temporary changes to work-family conflict last longer for men).

For average (stable B) job satisfaction and work-family conflict the story is simpler, as greater levels of both appeared to increase persistence of dynamic W perturbations for job satisfaction, respectively β B ( J ) W ( J ) J  = .837 (p = < .001) and β B ( C ) W ( J ) J  = .365 (p < .001). This indicates that higher levels of average B job satisfaction and work-family conflict appear to be associated with slower regression to the mean for W job satisfaction. However, the story is reversed for W work-family conflict, with average (stable B) job satisfaction and work-family conflict associated with reduced persistence of W dynamic work-family conflict, β B ( J ) W ( C ) C  = −.451 (p < .001) and β B ( C ) W ( C ) C  = −.179 (p < .001), implying faster regression to the mean.

In terms of CL effects, accounting for cross-level interactions leads all of these to become significant compared to the baseline model, with all terms interpreted as the average CL effect across the sample for men (when gender = 0) as follows: β W ( C ) J  = .071 (p < .001), β W ( J ) C  = .229 (p < .001), β W ( H ) C  = .236 (p < .001), β W ( H ) J  = .027 (p = .004), β W ( C ) H  = .063 (p < .001), β W ( J ) H  = .121 (p < .001). Similar to the AR effects, most estimates were larger compared to the baseline model, which can be put into some context by considering cross-level interactions.

In terms of gender: women might have a larger W effect of job satisfaction on work-family conflict but this interaction is not significant, β B ( G ) W ( J ) C  = .072 (p = .25); and women show a significantly smaller W effect of work-family conflict on job satisfaction, β B ( G ) W ( C ) J  = −.041 (p < .001), meaning that (relative to men) women's job satisfaction appears to be less affected by a temporary W increase in work-family conflict. In terms of average (stable B) work-family conflict: a higher average level of work-family conflict seems to amplify the effect of a temporary increase of job satisfaction on work-family conflict, β B ( C ) W ( J ) C  = .121 (p = .01); and an increase in average work-family conflict seems to diminish the impact of work-family conflict's effect on job satisfaction levels over time, β B ( C ) W ( C ) J  = −.077 (p < .001), suggesting greater average work-family conflict may make people less sensitive to dynamic changes in work-family conflict in terms of its effect on future job satisfaction. In terms of average (stable B) job satisfaction: higher average job satisfaction was related with a weaker W effect of work-family conflict on job satisfaction, β B ( J ) W ( C ) J  = −.239 (p < .001), which implies that higher average job satisfaction may buffer people against the negative effects of dynamic changes to work-family conflict. In contrast, average job satisfaction was not associated with the W effect of job satisfaction on work-family conflict, β B ( J ) W ( J ) C  = .098 (p = .348).

In terms of the B model, cross-level interactions appear to have reduced uncertainty in these estimates, leading to smaller p-values. However, cross-level interactions also markedly change the point estimates, which was initially vexing given that cross-level interactions should impact only W rather than B parameters (Preacher et al., 2016). Specifically, the B estimates were: β B ( C ) J  = −.392 (p < .001), which is twice as large as in the baseline model; β B ( G ) J  = −.174 (p < .001), which is negative and significant rather than approximately zero in the baseline model; β B ( H ) J  = −.003 (p = .726), which is no longer positive and significant as in the baseline model; β B ( G ) C  = −.08 (p = .224), which is reversed in sign but non-significant as in the baseline model; β B ( H ) C  = .269 (p < .001), which is similar to the baseline model estimate; and β B ( G ) H  = −1.987 (p < .001), which is similar to the baseline estimate.

Importantly, the reason for these changes appears to be the use of the B parts of the job satisfaction and work-family conflict variables for BxW interactions. In brief, the B parts of these or any other variables may correlate with the B variation in the W effects implied by cross-level interactions—in the multilevel context this would be typical covariation among a random intercept (stable B part) and a random slope for a W coefficient (implying a potential cross-level interaction). By estimating BxW interactions using the B parts of the observed variables, this controls for the B part of job satisfaction and work-family conflict associated with BxW interactions. Thus, when introducing BxW interactions we eliminate any B variance in the observed variables that is correlated with the ‘random slope’ implied by the BxW interactions, which in turn adjusts the B parameters. Therefore, researchers may want to interpret B parameters separately from models that include such BxW interactions.

Full Model (WxW, BxB, and BxW). To provide the reader with a sense for how estimating all CLPM interaction effects simultaneously differs from models that estimate them separately, we also estimated a ‘full’ model that includes all types of interactions. To save space we refer the reader to Table 2 to compare each previously discussed model's results with the full model, with a few notable similarities and differences. First, results in the full model were very similar to the BxW model, with substantial differences in B coefficients compared to baseline, WxW, and BxB models for the reasons we have just described. Again, researchers interested in B parameters may want to avoid interpreting these B terms in models that also include latent cross-level interactions (i.e., BxW terms).

Second, the full model contains two WxW interactions whereas the BxW does not, and therefore it is relevant to compare the WxW terms against those in the full model. Consistent with accounting for significant variance with the addition of many interaction terms in the full model, the WxW interactions in the full model have smaller p-values (the estimates are also somewhat larger, but this may be merely due to variance caused by the estimator).

Conflated Model. Finally, we present a traditional CLPM with work-family conflict, gender, hours worked, and job satisfaction with interactions specified without decomposing the variables into B and W parts. We could have used a maximum likelihood estimator, but we used a Bayes estimator with comparable fit statistics (parameter estimates were almost identical across the two types of estimators). Model fit statistics worsened compared to the baseline model: χ² difference for observed versus model-generated data shows a 95% CI from 2069.879 to 2180.863 (not encompassing zero, indicating poor fit) with the number of parameters pD = 76.362; posterior predictive p-value or PPP < .001 (smaller than .05, indicating poor fit); but Bayes-equivalent RMSEA = .06 (smaller than .08, indicating good fit); CFI = .973 and TLI = .957 (larger than .95, indicating good fit); and DIC = 298804.233 (larger than the baseline model's DIC of 218442.103, indicating worse fit for this model).

AR effects were all significant and substantially larger than in the baseline model: β J J  = .639 (p < .001), much larger than the baseline model's .086; β C C  = .465 (p < .001), which can be compared to .203; and β H H  = .816 (p < .001), which again is much larger than the estimate obtained in the baseline model (.493). The issue here is that all stable B (co)variance in this typical CLPM must be accounted for by AR and CL terms in the model, which will often substantially overestimate AR terms in a positive direction (see Hamaker et al., 2015; Zyphur, Voelkle, et al., 2020).

The CL terms were also very different, including two that became significant in the conflated model. First, job satisfaction predicting work-family conflict β J C  = −.202 (p < .001) was over twice the baseline model's estimate (of −.079). Second, work-family conflict predicting hours worked β C H  = .02 (p < .001) was similar in magnitude but with a much smaller p-value. The other CL terms, which were similar both in magnitude and statistical significance, were: β C J  = .006 (p = .792); β H C  = .136 (p < .001); β H J  = −.007 (p = .054), β J H  = .005 (p = .67).

In terms of gender, a stable B variable, we used this to predict all occasions past the first (i.e., t > 1). Although this specification is uncommon (e.g., Shin & Konrad, 2017), it is the correct way to control for a B variable by using it to predict all dependent variables in the model. Results show β G J  = −.01 (p = .608), which is almost identical to the baseline, but β G C  = .034 (p = .012) is about half of the baseline model estimate and statistically significant. Last, β G H  = −.407 (p < .001) is less than one quarter of the baseline estimate and significant.

Finally, the product term estimates were not significant, with β C J J  = −.003 (p = .34) and β C J C  = −.001 (p = .92). These estimates are not very illuminating given the fact that the product terms conflate WxW, BxB, and cross-level BxW interaction effects, as we have noted. Estimating interactions with gender would result in the same problem except we would then be conflating only BxW and any BxB interaction effects (because gender was a purely B variable). Indeed, the fact that we are unable to estimate any such interactions with the latent B parts of job satisfaction and work-family conflict further clarifies the limitations of traditional CLPMs—apart from the difficulty controlling for the stable B parts as we and others have pointed out (see Hamaker et al., 2015; Zyphur, Allison, et al., 2020; Zyphur, Voelkle, et al., 2020). For further details, we point the reader to our Online Appendices which include the Mplus code and output for all models.

Simulations

To provide the reader with some preliminary insights into how CLPMs with latent interactions perform, we also ran simulations for the two cases that we believe will be most common in future research using our approach: purely within (WxW) models and (BxW) models. In the WxW case, this allowed us to compare our approach to observed-means centering and an uncentered approach with raw variables (this uncentered approach is equivalent to the ‘conflated model’ described previously and in Table 2 as a traditional CLPM; for all simulation results see Table 3). Alternatively, in the BxW case, we compared our approach to only observed-means centering because it would seem incoherent to propose a cross-level interaction but use uncentered variables to test this (see results in Tables 4a and 4b). As we will now describe, our proposed CLPM approach outperformed the alternatives, and for some parameters drastically so.

OPEN IN VIEWER

Table 3.

Results of Simulation Study 1 Examining WxW Interaction Effects.

	Proposed CLPM Approach					Observed-Means Approach					Uncentered Predictors Approach
	Estimates	SE	Bias in %	95% Coverage	Proportion of Significant Coefficients	Estimates	SE	Bias in %	95% Coverage	Proportion of Significant Coefficients	Estimates	SE	Bias in %	95% Coverage	Proportion of Significant Coefficients
WxW Interaction Effect = −0.3
W effects
AR effect X	0.502	0.061	−0.4	.957	.998	−0.373	0.032	174.6	.000	1.000	0.747	0.024	−49.4	.000	1.000
AR effect Y	0.500	0.061	0.0	.951	.998	−0.372	0.033	174.4	.000	1.000	0.748	0.024	−49.6	.000	1.000
CL effect YW→XW	0.299	0.052	0.3	.949	.992	0.044	0.033	85.3	.000	.314	0.175	0.024	41.7	.000	1.000
CL effect XW→YW	0.298	0.051	0.1	.947	.992	0.048	.033	84.0	.000	.330	0.175	0.024	41.7	.000	1.000
WxW interaction
XW*YW →XW	−0.300	0.046	0.0	.943	.998	−0.036	0.050	88.0	.022	.294	−0.079	0.015	73.7	.000	0.990
XW*YW →YW	−0.297	0.046	1.0	.943	.998	−0.038	0.051	87.3	.018	.276	−0.079	0.016	73.7	.000	0.994
B covariance
XB - YB	0.315	0.101	−5.0	.951	.927	0.707	0.111	−135.7	.024	1.000	—	—	—	—	—
Number of converged replications				494					500					500
WxW Interaction Effect = 0
W effects
AR effect X	0.489	0.090	2.2	.941	.988	−0.373	0.034	174.6	.000	1.000	0.725	0.026	−45.0	.000	1.000
AR effect Y	0.488	0.090	2.4	.946	.995	−0.373	0.036	174.6	.000	1.000	0.725	0.027	−45.0	.000	1.000
CL effect YW→XW	0.276	0.071	8.0	.939	.955	0.039	0.035	87.0	.000	.216	0.172	0.026	42.7	.000	1.000
CL effect XW→YW	0.277	0.071	7.7	.932	.962	0.043	0.036	85.7	.000	.252	0.173	0.026	42.3	.000	1.000
WxW interaction
XW*YW →XW	0.001	0.061	—	.927	.073	0.000	0.053	—	.932	.068	0.000	0.018	—	.928	.072
XW*YW →YW	−0.001	0.061	—	.951	.049	−0.002	0.055	—	.920	.080	0.000	0.018	—	.950	.050
B covariance
XB - YB	0.384	0.167	−28.0	.915	.594	0.620	0.103	−106.7	.100	1.000	—	—	—	—	—
Number of converged replications				426					500					500
WxW Interaction Effect = 0.3
W effects
AR effect X	0.505	0.060	−1.0	.946	.998	−0.372	0.032	174.4	.000	1.000	0.749	0.024	−49.8	.000	1.000
AR effect Y	0.505	0.060	−1.0	.926	.998	−0.372	0.033	174.4	.000	1.000	0.748	0.024	−49.6	.000	1.000
CL effect YW→XW	0.300	0.051	0.0	.948	.992	0.044	0.033	85.3	.000	.278	0.176	0.024	41.3	.000	1.000
CL effect XW→YW	0.299	0.051	0.3	.956	.992	0.047	0.033	84.3	.000	.316	0.176	0.024	41.3	.000	1.000
WxW interaction
XW*YW →XW	0.292	0.046	2.7	.946	.998	0.032	0.049	89.3	.014	.290	0.079	0.015	73.7	.000	.994
XW*YW →YW	0.294	0.046	2.0	.960	1.000	0.032	0.050	89.3	.016	.258	0.079	0.016	73.7	.000	.992
B covariance
XB - YB	0.318	0.101	−6.0	.950	.914	0.710	0.111	−136.7	.024	1.000	—	—	—	—	—
Number of converged replications				499					500					500

Note. Parameter estimates are based on 500 replications. The AR effects were fixed at 0.5. The CL effects were fixed at 0.3. The B covariance was fixed at 0.3. Each replication has N = 300. Bias was calculated using the formula (Population parameter – Parameter estimate) / (Population parameter). 95% Coverage indicates the proportion of replications whose 95% confidence intervals include the population parameter. Proportion of Significant Coefficients indicates the proportion of replications for which the null hypothesis is rejected at p = .05. This value represents the power when the population parameter is nonzero and the Type-I error rate when the population parameter is zero.

OPEN IN VIEWER

Table 4a.

Results of Simulation Study 2 Examining BxW Interaction Effects for N = 300 (with T = 3/T = 4 Results Shown).

	Proposed CLPM Approach					Observed-Means Approach
	Estimates	SE	Bias in %	95% Coverage	Proportion of Significant Coefficients	Estimates	SE	Bias in %	95% Coverage	Proportion of Significant Coefficients
BxW Interaction Effect = –0.3
W effects
AR effect X	0.499/0.498	0.078/0.054	0.2/0.4	.943/.946	1.000/1.000	−0.371/−0.136	0.035/0.033	174.2/127.2	.000/.000	1.000/.998
AR effect Y	0.505/0.502	0.073/0.049	−1.0/−0.4	.959/.958	1.000/1.000	−0.349/−0.097	0.033/0.031	169.8/119.4	.000/.000	1.000.908
CL effect YW→XW	0.291/0.293	0.059/0.041	3.0/2.3	.936/.960	.992/1.000	0.039/0.069	0.036/0.033	87.0/77.0	.000/.000	.206/.508
CL effect XW→YW	0.295/0.296	0.058/0.042	1.7/1.3	.943/.966	.992/1.000	0.054/0.084	0.033/0.029	82.0/72.0	.000/.000	.366/.810
BxW interaction
XW*YB →YW	−0.302/−0.303	0.073/0.047	−0.7/−1.0	.941/.958	1.000/1.000	−0.082/−0.124	0.025/0.023	−72.7/−58.7	.000/.000	0.910/.996
B covariance
XB - YB	0.306/0.304	0.127/0.107	−2.0/−1.3	.941/.944	.697/.859	0.621/0.660	0.107/0.102	−107.0/−120.0	.100/.032	1.000/1.000
Number of converged replications				488/498					500/500
BxW Interaction Effect = 0
W effects
AR effect X	0.492/0.496	0.090/0.060	1.6/0.8	.924/.945	.997/1.000	−0.371/−0.136	0.035/0.033	174.2/127.2	.000/.000	1.000/1.000
AR effect Y	0.510/0.506	0.089/0.060	−2.0/−1.2	.918/.916	.997/1.000	−0.373/−0.134	0.035/0.032	174.6/126.8	.000/.000	1.000/.998
CL effect YW→XW	0.277/0.291	0.071/0.046	7.7/3.0	.921/.933	.952/1.000	0.041/0.072	0.036/0.033	86.3/76.0	.000/.000	.226/.574
CL effect XW→YW	0.285/0.296	0.701/0.046	5.0/1.3	.946/.954	.983/1.000	0.043/0.072	0.035/0.031	85.7/76.0	.000/.000	.212/.628
BxW interaction
XW*YB →YW	−0.003/−0.005	0.096/0.054	—/—	.952/.933	.048/.067	0.000/−0.002	0.026/0.025	—/—	.930/.910	.070/.090
B covariance
XB - YB	0.319/0.282	0.181/0.151	−6.3/6.0	.898/.887	.414/.484	0.620/0.653	0.106/0.102	−106.7/−117.7	.096/.028	1.000/1.000
Number of converged replications				353/417					500/500
BxW Interaction Effect = 0.3
W effects
AR effect X	0.498/0.501	0.079/0.054	0.4/−0.2	.944/.932	.998/1.000	−0.372/−0.136	0.035/0.033	174.4/127.2	.000/.000	1.000/.998
AR effect Y	0.505/0.501	0.073/0.049	−1.0/−0.2	.944/.918	.998/1.000	−0.348/−0.099	0.033/0.031	169.9/119.8	.000/.000	1.000/.916
CL effect YW→XW	0.292/0.294	0.059/0.041	2.7/2.0	.963/.930	.994/1.000	0.040/0.069	0.036/0.033	86.7/77.0	.000/.000	.234/.528
CL effect XW→YW	0.295/0.296	0.059/0.042	1.7/1.3	.944/.954	.996/1.000	0.053/0.084	0.033/0.029	82.3/72.0	.000/.000	.358/.768
BxW interaction
XW*YB →YW	0.303/0.300	0.072/0.047	−1.0/0.0	.963/.930	1.000/1.000	0.081/0.120	0.025/0.023	73.0/60.0	.000/.000	.862/.994
B covariance
XB - YB	0.305/0.302	0.128/0.108	−1.7/−0.7	.944/.964	.649/.836	0.622/0.662	0.107/.103	−107.3/−120.7	.088/.024	1.000/1.000
Number of converged replications				481/499					500/500

Note. Parameter estimates are based on 500 replications. The first value in a cell indicates the findings for T = 3. The second value in a cell indicates the findings for T = 4. The AR effects were fixed at 0.5. The CL effects were fixed at 0.3. The B covariance was fixed at 0.3. Each replication has N = 300. Bias was calculated using the formula (Population parameter – Parameter estimate) / (Population parameter). 95% Coverage indicates the proportion of replications whose 95% confidence intervals include the population parameter. Proportion of Significant Coefficient indicates the proportion of replications for which the null hypothesis is rejected at p = .05. This value represents the power when the population parameter is nonzero and the Type-I error rate when the population parameter is zero.

OPEN IN VIEWER

Table 4b.

Results of Simulation Study 2 Examining BxW Interaction Effects for N = 2,000 (with T = 3/T = 4 Results Shown).

	Proposed CLPM Approach					Observed-Means Approach
	Estimates	SE	Bias in %	95% Coverage	Proportion of Significant Coefficient	Estimates	SE	Bias in %	95% Coverage	Proportion of Significant Coefficient
BxW Interaction Effect = –0.3
W effects
AR effect X	0.501/0.499	0.030/0.021	−0.2/0.2	.927/.950	1.000/1.000	−0.371/−0.134	0.013/0.013	174.2/126.8	.000/.000	1.000/1.000
AR effect Y	0.500/0.499	0.028/0.019	0.0/0.2	.941/.940	1.000/1.000	−0.348/−0.098	0.013/0.012	169.6/119.6	.000/.000	1.000/1.000
CL effect YW→XW	0.300/0.300	0.023/0.016	0.0/0.0	.949/.954	1.000/1.000	0.040/0.069	0.013/0.012	86.7/77.0	.000/.000	.832/.998
CL effect XW→YW	0.299/0.299	0.022/0.016	0.3/0.3	.943/.940	1.000/1.000	0.052/0.081	0.013/0.012	82.7/73.0	.000/.000	.982/1.000
BxW interaction
XW*YB →YW	−0.301/−0.300	0.024/0.017	−0.3/0.0	.947/.958	1.000/1.000	−0.083/−0.122	0.009/0.009	72.3/59.3	.000/.000	1.000/1.000
B covariance
XB - YB	0.302/0.300	0.048/0.040	−0.7/0.0	.927/.928	1.000/1.000	0.617/0.657	0.041/0.038	−105.7/−119.0	.000/.000	1.000/1.000
Number of converged replications				491/499					500/500
BxW Interaction Effect = 0
W effects
AR effect X	0.501/0.499	0.035/0.024	−0.2/0.2	.933/.952	1.000/1.000	−0.371/−0.134	0.013/0.013	174.2/126.8	.000/.000	1.000/1.000
AR effect Y	0.501/0.499	0.035/0.024	−0.2/0.2	.914/.952	1.000/1.000	−0.373/−0.135	0.013/0.012	174.6/127.0	.000/.000	1.000/1.000
CL effect YW→XW	0.296/0.300	0.028/0.018	1.3/0.0	.944/.958	1.000/1.000	0.042/0.071	0.013/0.012	86.0/76.3	.000/.000	.842/.998
CL effect XW→YW	0.296/0.299	0.028/0.018	1.3/0.3	.925/.932	1.000/1.000	0.040/0.069	0.013/0.012	86.7/77.0	.000/.000	.826/1.000
BxW interaction
XW*YB →YW	−0.002/0.000	0.022/0.016	—/—	.944/.952	.056/.048	−0.001/0.000	0.010/0.009	—/—	.914/.926	.086/.074
B covariance
XB - YB	0.303/0.297	0.077/0.060	−1.0/1.0	.880/.888	.889/.974	0.614/0.684	0.041/0.038	−104.7/−128.0	.000/.000	1.000/1.000
Number of converged replications				359/455					500/500
BxW Interaction Effect = 0.3
W effects
AR effect X	0.500/0.500	0.030/0.021	0.0/0.0	.924/.948	1.000/1.000	−0.371/−0.134	0.013/0.013	174.2/126.8	.000/.000	1.000/1.000
AR effect Y	0.501/0.498	0.029/0.019	−0.2/0.4	.922/.944	1.000/1.000	−0.348/−0.099	0.013/0.011	169.6/119.8	.000/.000	1.000/1.000
CL effect YW→XW	0.301/0.301	0.022/0.016	−0.3/−0.3	.939/.950	1.000/1.000	0.041/0.069	0.013/0.012	86.3/77.0	.000/.000	.838/.998
CL effect XW→YW	0.300/0.299	0.022/0.016	0.0/0.3	.947/.946	1.000/1.000	0.052/0.081	0.013/0.012	82.7/73.0	.000/.000	.976/1.000
BxW interaction
XW*YB →YW	0.301/0.300	0.025/0.017	−0.3/0.0	.930/.948	1.000/1.000	0.082/0.122	0.009/0.009	72.7/59.3	.000/.000	1.000/1.000
B covariance
XB - YB	0.297/0.301	0.048/0.040	1.0/−0.3	.930/.930	1.000/1.000	0.615/0.656	0.041/0.038	−105.0/−118.7	.000/.000	1.000/1.000
Number of converged replications				489/499					500/500

Note. Parameter estimates are based on 500 replications. The first value in a cell indicates the findings for T = 3. The second value in a cell indicates the findings for T = 4. The AR effects were fixed at 0.5. The CL effects were fixed at 0.3. The B covariance was fixed at 0.3. Each replication has N = 300. Bias was calculated using the formula (Population parameter – Parameter estimate) / (Population parameter). 95% Coverage indicates the proportion of replications whose 95% confidence intervals include the population parameter. Proportion of Significant Coefficient indicates the proportion of replications for which the null hypothesis is rejected at p = .05. This value represents the power when the population parameter is nonzero and the Type-I error rate when the population parameter is zero.

Each simulation condition was replicated 500 times. The true model always had a B variance of 1.0 and a model-implied W variance of 1.0 (adjusting W residuals appropriately to account for W effects), giving an ICC(1) = .50 as is often found in longitudinal data. These unit variances imply estimated AR, CL, and interaction effects can be interpreted as being standardized. In all simulation conditions, AR terms were fixed to .5 to reflect a reasonable degree of persistence and CL terms were fixed at .3 to reflect modest cross-lagged effects. We also specified a covariance of .3 among B components to reflect a moderate correlation among B factors. To improve convergence with the Bayes estimator, we allowed for a small observed-variable residual of .03, so the total observed variance for any variable was 2.03. To account for this, we fixed observed-variable residuals to .03 when estimating our latent CLPM. To ensure that this small variance did not impact the performance of the alternative models, for both the observed-means centering and uncentered approaches we specified a latent variable for each observed variable with variance fixed at .03, but uncorrelated with all other variables, thus adding no estimated parameters and not impacting model df or fit. For both the WxW and the BxW simulations, we varied the interaction effects in three conditions: fixed to −0.3, 0, and 0.3 reflecting moderately negative, zero, and moderately positive interaction effects, respectively. We chose these values to be in line with effect sizes found in prior work (e.g., Ernst Kossek & Ozeki, 1998, reported a meta-analyzed relationship of around -.3 between work-family conflict and job satisfaction).

In Simulation Study 1, we generated three waves of data (T = 3), as would be fairly common in applied research, and we set the sample size to N = 300 with a WxW interaction among the W parts of two observed variables x and y. This Simulation Study 1 thus has a 1 × 1 × 3 design: 1 (Sample size: N = 300) by 1 (Occasion: T = 3) by 3 (Interaction effect sizes: −.3, 0, .3). Results shown in Table 3 support the latent means approach, showing good statistical coverage and only comparatively small levels of bias for all parameters across all conditions. Alternatively, the observed-means approach underestimates AR terms, the uncentered approach overestimates AR terms, both of these methods underestimate CL terms to different degrees, and only in the no-interaction condition do these approaches offer good estimates of the interaction effect—in the −.3 and .3 interaction cases, the estimates are markedly shrunk towards zero.

Simulation Study 2 examines BxW interactions with a single interaction among the B part of y and the W part of x to predict the future W part of y (shown as X_W*Y_B →Y_W in Tables 4a and 4b), contrasting our proposed latent CLPM approach with an observed-means centering approach. In this simulation, we set T = 3 and N = 300 (shown in Table 4a), but because cross-level interactions are known to be sensitive to the number of lower-level observations per higher-level unit (see Raudenbush & Bryk, 2002), we further varied the number of waves to T = 4. Also, for comparison, we included conditions with the sample size of N = 2,000 (see Table 4b). ⁴ Simulation Study 2 thus has a 2 × 2 × 3 design: 2 (Sample size: N = 300 and N = 2,000) by 2 (Occasion: T = 3 and T = 4) by 3 (Interaction effect sizes: −.3, 0, .3). Results show that the latent CLPM approach outperforms the observed-means centering method in all cases, with a similar pattern for the observed-means approach: underestimated AR terms; CL terms shrunk towards zero; and cross-level interactions markedly shrunk towards zero. Notably, the larger sample size of N = 2,000 in Table 4b does not appear to improve estimates very much, but the larger T = 4 sample size does improve estimates for the observed-means case—this is to be expected given the longitudinal panel data modeling literature on the topic (see Nickell, 1981; Zyphur, Allison, et al., 2020; Zyphur, Voelkle, et al., 2020).

Before continuing, a relevant finding in both studies was that fewer latent CLPMs converged when interaction effects were zero, particularly with smaller samples (see Tables 3, 4a, and 4b). This may point to inherent instabilities with very small interaction effects, but in our view, this more likely suggests that a larger number of iterations may be required in the presence of very small interaction effects and small samples—we used a maximum of 20,000 iterations in our simulation runs. Requesting more iterations or thinning iterations should therefore alleviate these issues—something that is more feasible when estimating only one data set instead of 500, as was the case in our simulations. Conveniently, estimation becomes faster as sample sizes decrease, and therefore increasing the number of iterations should not present computational issues in practice. In sum, our results suggest taking a latent mean centering approach over others for estimating interactions in CLPMs. In both simulation studies, we find that our proposed approach outperforms both an observed-means as well as uncentered method with raw variables, providing less biased estimates with greater coverage to improve statistical inference. The full simulation files and results can be found in our Online Appendices.

OPEN IN VIEWER

Table 5.

Overview of Software Packages for Estimating Latent Interactions in SEMs.

Software package name	Latent interaction testing procedures in alternative software and their capabilities	Key reference(s)
Mplus	Both latent moderated structural equations (LMS) and Bayesian estimation allowed. The LMS approach runs into difficulties with more than only a few dimensions because it requires numerical integration, which is computationally heavy. Bayesian estimation (as used in our approach here) is significantly faster, less biased, and has much better convergence rates particularly in high dimensions.	Asparouhov & Muthén (2020); Muthén and Muthén (1998–2017)
blavaan and JAGS (R)	The blavaan package may soon be an appealing approach to latent interactions in SEM given their plans for this development, but currently it does not support this feature (Merkle & Rosseel, 2018). In addition, researchers have noted that JAGS seems to have difficulty in handling such complex models (Depaoli et al., 2016).	Depaoli et al. (2016); Merkle & Rosseel (2018)
gllamm and gsem (Stata)	The gllamm package in Stata allows the implementation of our approach and Mplus code logic. However, researchers have noted the long computation times and instability of gllamm when estimating complex models. More precisely, according to Grilli and Rampichini (2006), computation time for interaction effects is “approximately proportional to the product of the number of quadrature points for all latent variables used. For example, if there are two random effects at level 2 (a random intercept and slope) and 8 quadrature points are used for each random effect, the time will be approximately proportional to 64” (p. 4–5). Like Mplus, Stata is a commercial software and thus may not be available to all researchers. The recent gsem package in Stata, which appears to have partially supplanted gllamm, does not appear to allow for same-level latent interactions, which our approach requires.	Grilli and Rampichini (2006); Rabe-Hesketh et al. (2005); Rabe-Hesketh and Skrondal (2012)
lavaan (R)	Cortina et al. (2021) provide R code for implementing a latent moderated structural equations (LMS) method with fully and partially latent approaches. However, they note that latent procedures, specifically moderation, involve “complex or impractical” methods that make researchers abandon latent interaction testing (p. 16). Cheung and colleagues (2021) suggest that lavaan cannot be used for LMS. Similarly, Cortina et al. (2021) recommend the use of Mplus for LMS.	Cheung et al. (2021); Cortina et al. (2021)
nlsem (R)	This package is useful for nonlinear latent interactions. However, the estimation of complex models can be very slow particularly in high dimensions (Umbach et al., 2017).	Umbach et al. (2017)
OpenMx (R)	This package can be useful in simple multilevel SEMs, however, more complex models are currently not supported. Specifically, Neale et al. (2016) state that “the OpenMx development team is working toward a much more general solution that would accommodate cross-classified models as well as large and complex data. We are still working on syntax for general purpose multilevel models that is both comprehensive and simple to understand” (p. 547).	Neale et al. (2016)
xxM (R)	This package is designed for multilevel SEM analysis. However, Pritikin et al. (2017) report challenges associated with estimating complex models. Discussions in the OpenMx forum indicate that there are no recent developments regarding the implementation of latent interactions but interested readers can watch these packages as developers update them with new estimation algorithms and coding methods for latent interactions (AdminNeale, 2020a, 2020b; Pritikin, 2020).	AdminNeale (2020a, 2020b); Pritikin (2020); Pritikin et al. (2017)

Discussion

This paper describes how to estimate interactions among stable (between-level or B) and/or time-varying (within-level or W) parts of latent variables in CLPMs. We propose that prior research examining interaction effects in CLPMs overlooks the potential importance of separating different types of (latent) interaction effects. As a consequence, interpretations of interaction effects in past CLPMs and conclusions with regards to hypothesis testing may have been flawed. At the same time, some researchers may have concluded that interaction effects did not exist when, in fact, they did if B and W components had been properly separated. To address this issue and assist researchers interested in testing interaction effects in CLPMs, we offer a latent interactions approach that can be used to detect previously unconsidered forms of moderation involving B and/or W parts of predictors. As noted above, although not all of the many types of interactions and nonlinear effects that our framework allows may be of interest, the convenient feature of our proposed method is that any one or more of these effects can be precisely specified and estimated to test theory, while excluding others that are not theoretically relevant.

Conveniently, our approach can be easily applied to other methods that have a similar process of forming latent W terms, including the general cross-lagged panel model or GCLM wherein W terms are referred to as ‘impulses’ u (see Zyphur, Allison, et al., 2020); the ‘xtdpdml’ panel data model (see Allison et al., 2017; Moral-Benito et al., 2019; Williams et al., 2018); as well as other methods (e.g., Bollen & Brand, 2010). In general, whenever a lagged-effects structure exists that separates stable B and dynamic W components, we recommend formally delineating WxW, BxB, and BxW interactions—although in other models the B terms may need to be rescaled (for discussion see Zyphur, Allison, et al., 2020).

Furthermore, it is notable that WxW and BxW interactions can be applied to indirect effects using AR and CL terms (for general discussion see Preacher et al., 2016; Zyphur et al., 2018). As Zyphur, Allison, et al. (2020) note, the logic of an ‘impulse response’ can be used to model the long-run effects of a 1-unit change in a predictor variable as effects ‘flow’ through a system over time along AR and CL paths. The implied impulse responses at low and high levels of a moderating variable involved in an interaction can be plotted to gain insight into potential long-run effects at different values of the W or B components of variables involved in an interaction. The approach we take here can also be extended to additional cases, including latent-indicator models to account for measurement error or multiple-group models (see Mulder & Hamaker, 2021), which can be easily implemented in SEM (see Little, 2013), as well as indirect-effects in the B model.

In terms of our applied example, prior research has come a long way in attempts to tease apart causal associations and potential interactions/moderation effects. However, much of this work has not accounted for potential stable B versus dynamic W components (e.g., Spector et al., 2007). For example, Grandey et al. (2005) found that an annual change in job satisfaction was predicted by work-family conflict for women only. Alternatively, using our data and our B and W decomposition with latent interactions we found an effect of past work-family conflict on future job satisfaction, which was actually stronger for men than for women, after accounting for hours worked. Using the framework we advance here, future studies can better interrogate the causal relationships among these variables including potential higher-order interactions and forms of nonlinearity.

Limitations and Future Research

We note four key limitations of our proposed method. First, although our method is computationally tractable, it can still be burdensome to estimate many interaction effects in large datasets. To speed up estimation, we used an AWS EC2 instance with up to 48 physical cores (96 virtual cores), which could still take a few hours to converge—along with relevant checks for convergence as recommended by Asparouhov and Muthén (2020). Run on a Linux instance rather than Windows allowed using all 48 cores to speed up processing. However, discovering this fact by running Windows and Linux instances on AWS for quite some time, including thousands of simulation runs, cost our research team roughly $3,000 USD—an expense future researchers will not need to incur if using cheaper and much faster Linux instances on AWS. Furthermore, our sample size was comparatively large with almost 8,000 people measured at 4 occasions for 3 separate variables, resulting in 12 latent W variables and 3 latent B variables. Hopefully our efforts will have removed much of the required trial and error, and costs, associated with latent interaction CLPMs.

Furthermore, future research may be able to draw from advances in the field of Bayesian statistics to reduce convergence time (e.g., Hamiltonian Monte Carlo; see Gelman et al., 2013) and utilize upcoming CPUs with even greater core counts and higher performance. Also, in practice, researchers can exclude any latent interactions that are not theoretically relevant, thus reducing computation times. In fact, for Simulation Study 1, it took about 50 h to estimate a total of 1,500 models involving two WxW interaction effects (T = 3, N = 300) with only 2 time-series variables x and y using our proposed CLPM approach on a standard office computer. Assuming this sample size will be relatively standard, this is a modest 2 min of estimation time per model.

The second limitation is that our approach is currently feasible only in Mplus using a wide (single-level) SEM framework and is thus limited to approximately 10 occasions of measurement (T = 10). Mplus is currently the only available software package, to our knowledge, that can deal effectively with many latent interactions in SEMs. For the interested reader, we present commonly-used software packages and their limitations with regards to our proposed approach in Table 5. We hope that future work will provide researchers with alternatives to proprietary software. Regarding the number of measurement occasions, we suppose that for T > 10, there may be too many observed variables to ensure reasonable computation times. In such cases, researchers may opt for a multilevel framework such as DSEM (see Asparouhov et al., 2018), which currently allows cross-level BxW interactions but not yet latent WxW or BxB interactions (see Asparouhov et al., 2018; Asparouhov & Muthén, 2019; Hamaker & Muthén, 2020).

Relatedly, a third limitation is that our proposed approach makes it burdensome to estimate ‘random slopes’ for AR and CL terms, as in multilevel models (including DSEM). Although BxW interactions typically involve a ‘slopes-as-outcomes’ model, where a W effect is allowed to have a B variance (see Raudenbush & Bryk, 2002), the model input and logic in the single-level case (such as ours) is rather complicated. Although it is technically possible to specify a latent variable to capture B variation in AR and CL terms, results of our second simulation study suggest that estimates of this B variation will probably not be trustworthy given the small T involved (contrast T = 3 and T = 4 results in Tables 4a and 4b). Therefore, we have not pursued this approach here and we do not recommend it. Instead, we recommend our approach wherein BxW terms are treated as fixed, which the reader can selectively interpret as a slopes-as-outcomes model that excludes a B residual for the slope.

The fourth limitation is the fact that models that include BxW interactions with the stable B part of a longitudinal variable may bias estimates of purely B regressions, and likewise B regressions may bias the estimation of BxW interactions. This may occur with few occasions of measurement because of the inherent uncertainty in B estimates that are involved in these model components with small T. Therefore, if researchers are interested in testing hypotheses associated with purely B effects, then these should likely be evaluated in models that exclude BxW interactions among the B variables involved. However, we would point out that this does not seem like a very serious problem given that the entire point of lagged-effects models, such as CLPMs, is to avoid having to make causal inferences based on purely B ‘effects’. Instead, we recommend focusing on AR and CL terms in the W part of a CLPM and estimating WxW and BxW interactions as needed to test theory (while allowing for covariance among B parts to hold them constant, i.e., a ‘fixed-effects’ method).

In terms of future research, we propose that a key direction for innovation is the development of standardization methods for latent interaction effects in CLPMs. Typically, obtaining standardized effect estimates when testing interactions is done by standardizing variables before forming product terms. However, as we have decomposed observed variables into latent B and W parts for analysis, standardizing observed variables before analysis is not possible. Instead, it is theoretically possible to standardize the latent B and W variables in our model by imposing a unit variance and zero mean on them. To this end, Asparouhov and Muthén (2020) suggest a standardization approach that could overcome the need to impose such a model specification. When doing so, previous work has pointed out that it may be important to account for B variation in W variances used for standardization (see Schuurman et al., 2016). Future research should examine how this process can be facilitated in the context of latent interaction effects like those we have presented here.

Supplemental Material

sj-rar-1-orm-10.1177_10944281211043733 - Supplemental material for Interaction Effects in Cross-Lagged Panel Models: SEM with Latent Interactions Applied to Work-Family Conflict, Job Satisfaction, and Gender

Supplemental material, sj-rar-1-orm-10.1177_10944281211043733 for Interaction Effects in Cross-Lagged Panel Models: SEM with Latent Interactions Applied to Work-Family Conflict, Job Satisfaction, and Gender by Ozlem Ozkok, Manuel J. Vaulont and Michael J. Zyphur, Zhen Zhang, Kristopher J. Preacher, Peter Koval, Yixia Zheng in Organizational Research Methods