Assessing Mediational Processes Using Piecewise Linear Growth Curve Models with Individual Measurement Occasions

In this section, we first introduce a simple mediation model, which is represented in the diagram by Figure 1 (Baron and Kenny,, 1986). As shown in the figure, a simple mediation model contains three variables, a predictor $X$, a mediator $M$, and an outcome $Y$, with $X$ is proposed as influencing $Y$ through $M$. Note that in this setup, both $M$ and $Y$ are dependent variables. There are two possible pathways through which the predictor $X$ can impact $Y$: the path leads from $X$ to $Y$ without passing through $M$ (referred to as the direct effect of $X$ on $Y$), and the path from $X$ to $Y$ through $M$ (referred to as the indirect effect or mediated effect of $X$ on $Y$). The direct effect represents how $Y$ is influenced by $X$ directly, while the indirect effect represents how $Y$ is affected by $X$ through a causal sequence where $X$ affects $M$, which, in turn, impacts $Y$.

The primary research question when employing a mediation model is estimating and interpreting the direct and indirect effects. For example, in the simple mediation model, the estimation of these effects is through two linear regression models, one for each dependent variable, as shown below

		$\displaystyle M=\beta_{M0}+a\times X+\epsilon_{M},$
		$\displaystyle Y=\beta_{Y0}+c\times X+b\times M+\epsilon_{Y},$

where $\beta_{M0}$ and $\beta_{Y0}$ are regression constants, $\epsilon_{M}$ and $\epsilon_{Y}$ are errors of $M$ and $Y$, respectively. In addition, $a$, $b$, and $c$ are regression coefficients, among which $c$ quantifies the direct effect of $X$ on $Y$, while the product of the other two coefficients, $a\times b$, quantifies the indirect effect of $X$ on $Y$. According to Baron and Kenny, (1986), the most common situation in social science research is partial mediation, where the mediator only explains part of the relationship between the predictor and the outcome. Therefore, the total effect of $X$ on $Y$ is the sum of direct and indirect effects $c+a\times b$ (MacKinnon et al.,, 2000).

One challenge in the mediation analyses is how to obtain standard errors of mediated effects, based on which we can conduct a hypothesis test for statistically significant partial mediation and construct a $95\%$ confidence interval (CI) for the indirect effect. Generally, earlier studies have recommended employing the Delta method, for example, Sobel, (1982, 1986); MacKinnon et al., (2002) to approximate the standard error of the product of $a$ and $b$, or utilizing bootstrap techniques to calculate more accurate estimates (Hayes,, 2009; Shrout and Bolger,, 2002). Moreover, Cheung and Lau, (2008) has shown that the bias-corrected bootstrap CI performs best in the tests for the mediated effects compared to other methods. In the present study, we do not intend to develop a novel approach to obtaining standard errors of mediated effects or compare different approaches for constructing $95\%$ CIs. Instead, we utilize the Delta method to obtain standard errors of mediated effects as more accurate standard errors are out of the primary research interests of the current study. Note that the standard errors calculated using the Delta method can be automatically obtained from multiple software of the structural equation modeling (SEM) framework such as Mplus 8 (Muthén and Muthén,, 2017) and OpenMx (Boker et al.,, 2020).

Earlier studies, for example, Gollob and Reichardt, (1987) have described multiple fundamental problems of applying mediation models to cross-sectional data. On the one hand, the causal relationships often take time to unfold, implying that the magnitude of the causal effect usually depends on the elapsed time between measurement occasions of the predictor, mediator, and outcome. However, the utilization of cross-sectional data assumes that the causal effects are instantaneous and that the magnitude of a causal effect remains the same over time. On the other hand, using cross-sectional data leaves out measurements at previous times that are potential predictors of the constructed model. Accordingly, longitudinal data are favored for testing the mediation hypotheses. Earlier studies proposed applying the simple mediation model to analyze data where the predictor, mediator, and outcome are from multiple waves. Such applications focus on interindividual differences. However, a critical aspect of a longitudinal study is intraindividual change, which cannot be analyzed explicitly by the simple mediation model (Selig and Preacher,, 2009). Additionally, Cole and Maxwell, (2003); Maxwell and Cole, (2007) have demonstrated that the simple mediation model typically yields substantially biased estimates for longitudinal parameters.

In this section, we briefly review modeling frameworks to investigate mediation in longitudinal data, which allow for the examination of intraindividual change. Researchers have explored mediational processes in three frameworks: cross-lagged panel models (Cole and Maxwell,, 2003), latent growth curve models (Cheong et al.,, 2003), and latent change score models (MacKinnon,, 2008, Chapter 8). Cole and Maxwell, (2003) and MacKinnon, (2008, Chapter 8) have provided extensive overviews of the application of cross-lagged panel mediation models. At least three waves of all three constructs in a fully-specified model, the predictor, the mediator, and the outcome are needed to test how the predictor at the first wave (i.e., $t_{1}$) affects the mediator at the second wave (i.e., $t_{2}$), and in turn, influences the outcome at the third wave (i.e., $t_{3}$).

The second type of longitudinal mediation model stems from the latent growth curve modeling framework. In the scenario where there are a baseline predictor (i.e., a time-invariant predictor), a longitudinal mediator, and a longitudinal outcome, researchers employ linear growth curve models to depict the intraindividual change of the mediator and the outcome and examine how a baseline predictor affects the status or change of the mediator affects the change in the outcome. For example, Cheong et al., (2003) describes how an intervention program (the predictor) affects the change in the outcome through the change of the mediator using a parallel latent growth curve model. Soest and Hagtvet, (2011) showed that, other than the change in the longitudinal mediator process, its initial status can also be viewed as a mediator in a longitudinal mediation model with a note that the regression of the initial status is comparable to the cross-sectional mediation analysis. In addition, MacKinnon, (2008, Chapter 8) added that the predictor can also be a longitudinal variable and demonstrated how to extend the model to perform mediation analyses among three parallel longitudinal processes.

Earlier studies such as MacKinnon, (2008, Chapter 8) and Selig and Preacher, (2009) also extended the latent change score modeling framework (McArdle and Nesselroade,, 1994), where intraindividual change is captured by latent change scores, to examine longitudinal mediation. This modeling framework allows for the examination of mediated effects among the status of the predictor, mediator, and outcome. Its additional advantage is that it enables researchers to investigate how the change score of the predictor during the first time interval (i.e., $t_{2}-t_{1}$) influences the change score of the mediator during the second time interval (i.e., $t_{3}-t_{2}$), and in turn, affects the change score of the outcome during the third time interval (i.e., $t_{4}-t_{3}$). Selig and Preacher, (2009) have shown that the mediation effects can also be any combination of statuses and change scores of constructs. Earlier studies, for example, MacKinnon, (2008, Chapter 8) and Selig and Preacher, (2009), have provided a detailed discussion in terms of the theory of change for the three variables (i.e., the predictor, mediator, and outcome), the role of time, and the types of indirect effects of these three longitudinal mediation models.

The present study focuses on the second type of longitudinal mediation model, latent growth mediation models. The existing latent growth mediation models developed in Cheong et al., (2003) and MacKinnon, (2008, Chapter 8) assume that change patterns of each repeated variable take the linear functional form. Given that a process may show nonlinear change patterns with respect to time, Cheong et al., (2003) suggested employing a two-stage piecewise (also referred to as a bilinear spline, see Grimm et al., (2016, Chapter 12)) parallel growth model instead of the parallel linear growth model to evaluate nonlinear trajectories. Accordingly, in this present study, we developed two models to investigate mediational processes where the parallel bilinear spline growth model (PBLSGM) (Liu and Perera,, 2021) is employed to examine joint longitudinal processes. The first model can be utilized to evaluate a mediational process with a baseline predictor and subsequently two variables that are repeatedly measured over time, while the second model can be employed to investigate a mediational process with three longitudinal variables.

There are two statistical challenges to utilizing a linear-linear piecewise functional form. First, the change point or knot where the change of the rate-of-change occurs must be determined. Earlier studies have demonstrated that the knot can either be pre-specified by domain theories (Dumenci et al.,, 2019; Flora,, 2008; Riddle et al.,, 2015; Sterba,, 2014) or be estimated as an unknown parameter if the knowledge of the knot is unavailable or the estimation of the knot itself is of research interest. (Cudeck and du Toit,, 2003; Harring et al.,, 2006; Kwok et al.,, 2010; Kohli,, 2011; Kohli et al.,, 2013; Kohli and Harring,, 2013; Preacher and Hancock,, 2015; Liu,, 2019; Liu et al., 2019a, ; Liu et al., 2019b, ; Dominicus et al.,, 2008; McArdle and Wang,, 2008; Wang and McArdle,, 2008; Muniz Terrera et al.,, 2011; Kohli et al.,, 2015; Lock et al.,, 2018). In the two proposed models, we assume that the knot of each longitudinal process is unknown and to be estimated for two considerations. On the one hand, time is important in a longitudinal mediation model, and for each process, the knot is the transition time of the two stages. Accordingly, we want to avoid any unnecessary pre-specification of the knot. On the other hand, the knot can be viewed as a developmental milestone or an event in a longitudinal process. For example, when examining sleep behavior of infants over time, the knot could be the age of crawling since the sleep behavior is expected to be different pre- and post- crawling age (Grimm et al.,, 2016, Chapter 10). Similarly, when analyzing the recovery process from knee replacement surgery, the knot is viewed as the transition time from the recovery of surgical pain to the gradual recovery period (Dumenci et al.,, 2019). So ideally, the knot of the predictor process occurs first, followed by the knot of the mediator process, and then the knot of the outcome process. We want to obtain data-driven evidence for this temporal order from the proposed models.

Second, for a longitudinal process that takes a linear-linear piecewise functional form, we have three statuses, the initial status, the status at the knot (i.e., the change-point at which two linear segments join together), and the status at the end of the study. In addition, we have a slope for each of the two linear pieces, which captures the short-term change rate and the long-term change rate, respectively. Harring et al., (2006) has shown that the degrees of freedom of the linear-linear piecewise functional form is four; therefore, we cannot freely estimate all the statuses and slopes. Most existing studies with an unknown knot usually view the initial status, two slopes, and the knot as four free growth coefficients (Harring et al.,, 2006; Preacher and Hancock,, 2015; Liu et al., 2019a, ; Peralta et al.,, 2020; Liu and Perera,, 2021). In this project, we consider the two slopes, the knot and the measurement at the knot, as the free growth coefficients since the measurement at the knot is more relevant than the initial status to longitudinal mediation analysis given two reasons. First, the initial status of a longitudinal variable only provides information at baseline; therefore, the regression of the initial status of the longitudinal mediator on that of the predictor only estimates instantaneous effects, and so does the regression of the initial status of the longitudinal outcome (Soest and Hagtvet,, 2011). Second, as stated earlier, the knot is usually viewed as a milestone or an event in a longitudinal process; therefore, the status of the knot in one longitudinal variable may influence a long-term change in the other longitudinal process. Suppose we have two types of endpoints to evaluate the treatment effects of psychotherapy: one being a clinical indicator and the other a patient-reported outcome (PRO) related to illness burden and health-related quality of life. It is often expected to observe improvement in clinical indicators first and arrive at a plateau, leading to improved PROs. Therefore, the measurement at the plateau (i.e., the knot) of a clinical indicator is likely to affect the improvement in PROs later on.

Based on the above considerations, we view the knot as unknown for each longitudinal process in the two proposed models and consider the two slopes, the knot, and the status at the knot as the free growth coefficients. We then estimate the two slopes and the status of the knot at the individual level with an assumption that each process-specified knot is fixed across all individuals¹¹1Although multiple existing studies, for example, Preacher and Hancock, (2015); Liu et al., 2019a ; Peralta et al., (2020); Liu and Perera, (2021) have demonstrated that the variance of the knot can also be estimated, we assume that each process-specified knot is fixed across all individuals since the knot variance is out of the research interests of the present study. . Additionally, similar to Liu and Perera, (2021), we construct the proposed models in the framework of individual measurement occasions using the ‘definition variables’ approach (Mehta and West,, 2000; Mehta and Neale,, 2005; Sterba,, 2014) to avoid potential inadmissible estimation (Blozis and Cho,, 2008; Coulombe et al.,, 2015) and allow for different time structures across constructs (Liu and Perera,, 2021).

The remainder of this article is organized as follows. In the method section, we start from a bilinear spline growth model to estimate a fixed knot for a univariate longitudinal process. We then extend it to longitudinal mediation models and introduce the model specification of the proposed models. We define the mediational process as either the baseline covariate or the change in covariate influencing the change in the mediator, which, in turn, affects the change in the outcome. Next, we describe the model estimation and model evaluation that is realized by the Monto Carlo simulation. We then present simulation results and evaluate the proposed models in terms of non-convergence rate and the performance measures, which include the relative bias, the empirical standard error (SE), the relative root-mean-squared error (RSME), and the empirical coverage probability for a nominal $95\%$ confidence interval of each parameter. Next, in the application section, we demonstrate how to apply the proposed models to examine a data set of longitudinal records of reading, mathematics, and science abilities from ECLS-K: $2011$ (Tourangeau et al.,, 2013). Finally, discussions are framed regarding practical considerations, methodological considerations, and future directions.

In this section, we briefly describe a latent growth curve (LGC) model with a linear-linear functional form to estimate a fixed knot in the framework of individual measurement occasions. As shown in Figure 2, we specify a separate linear function for each of the two phases of the developmental process and express the measurement for the $i^{th}$ individual at the $j^{th}$ time in the framework of individual measurement occasions as

$$y_{ij}=\begin{cases}\eta^{[y]}_{0i}+\eta^{[y]}_{1i}t_{ij}+\epsilon^{[y]}_{ij}&t_{ij}\leq\gamma^{[y]}\\ \eta^{[y]}_{0i}+\eta^{[y]}_{1i}\gamma+\eta^{[y]}_{2i}(t_{ij}-\gamma^{[y]})+\epsilon^{[y]}_{ij}&t_{ij}>\gamma^{[y]}\\ \end{cases},$$

(1)

where $y_{ij}$ and $t_{ij}$ are the measurement and measurement occasion of the $i^{th}$ individual at wave $j$. In Equation (1), $\eta^{[y]}_{0i}$, $\eta^{[y]}_{1i}$ and $\eta^{[y]}_{2i}$ are the individual-level intercept, first slope and second slope, which are usually called ‘growth factors’ in the latent growth curve modeling framework; they three along with the fixed knot $\gamma^{[y]}$ together determine the functional form of the growth curve of $\bm{y}_{i}$.

Existing SEM software such as Mplus (Muthén and Muthén,, 2017) and the R package OpenMx (Boker et al.,, 2020) does not allow for a conditional statement of an unknown parameter such as $\gamma^{[y]}$ in Equation (1), so the piecewise functional form above cannot be specified directly. In order to unify pre- and post-knot expressions, we have to reparameterize growth factors, which can be realized through multiple approaches. Harring et al., (2006) proposed to reparameterize three growth factors in Equation (1) (i.e., $\eta^{[y]}_{0i}$, $\eta^{[y]}_{1i}$ and $\eta^{[y]}_{2i}$) to the average of the two intercepts, the average of the two slopes, and the half difference between the two slopes. Alternatively, Grimm et al., (2016, Chapter 11) suggested reexpressing the three growth factors to the measurement at the knot and two slopes. Additionally, Liu et al., 2019a reparameterized $\eta^{[y]}_{0i}$, $\eta^{[y]}_{1i}$ and $\eta^{[y]}_{2i}$ as the measurement at the knot, the average of the two slopes, and the half difference between the two slopes.

The reparameterized growth factors proposed in Harring et al., (2006); Liu et al., 2019a are no longer directly related to the underlying developmental process and therefore are less useful in the present study where we want to estimate regression coefficients between original growth factors. Accordingly, we follow the reparameterized approach in Grimm et al., (2016, Chapter 11), where all three reparameterized coefficients are still directly related to the growth patterns. We then write the repeated outcome as

	$\displaystyle y_{ij}$	$\displaystyle=(\eta^{[y]}_{0i}+\eta^{[y]}_{1i}\gamma^{[y]})+\eta^{[y]}_{1i}\min(0,t_{ij}-\gamma^{[y]})+\eta^{[y]}_{2i}\max(0,t_{ij}-\gamma^{[y]})+\epsilon^{[y]}_{ij}$
		$\displaystyle=\eta^{[y]}_{\gamma_{i}}+\eta^{[y]}_{1i}\min(0,t_{ij}-\gamma^{[y]})+\eta^{[y]}_{2i}\max(0,t_{ij}-\gamma^{[y]})+\epsilon^{[y]}_{ij},$

where $\eta^{[y]}_{\gamma_{i}}$ is the measurement at the knot of the $i^{th}$ individual.

In this section, we extend the univariate linear-linear latent growth curve model to analyze a mediational process with baseline predictor, longitudinal mediator, and longitudinal outcome. Suppose both the mediator and the outcome take the linear-linear functional form with an unknown fixed knot. In such situations, we can construct a model for a mediational process where the baseline predictor influences the change in the outcome directly and indirectly through its effect on the change in the mediator. The longitudinal mediation model can be specified as

$$\begin{pmatrix}\bm{m}_{i}\\ \bm{y}_{i}\end{pmatrix}=\begin{pmatrix}\bm{\Lambda}_{i}^{[m]}&\bm{0}\\ \bm{0}&\bm{\Lambda}_{i}^{[y]}\end{pmatrix}\times\begin{pmatrix}\bm{\eta}^{[m]}_{i}\\ \bm{\eta}^{[y]}_{i}\end{pmatrix}+\begin{pmatrix}\bm{\epsilon}^{[m]}_{i}\\ \bm{\epsilon}^{[y]}_{i}\end{pmatrix},$$

(2)

where $\bm{m}_{i}$ and $\bm{y}_{i}$ are a $J\times 1$ vector of the repeated measurements of the mediator and outcome of the $i^{th}$ individual, respectively (in which $J$ is the number of measurement occasions). With an assumption that the trajectories of mediator and outcome process take the linear-linear functional form with an unknown knot, $\bm{\eta}^{[u]}_{i}$ ($u=m,y$) is a $3\times 1$ vector of growth factors,

$$\bm{\eta}^{[u]}_{i}=\left(\begin{array}[]{rrr}\eta^{[u]}_{1i}&\eta^{[u]}_{\gamma_{i}}&\eta^{[u]}_{2i}\end{array}\right)^{T},$$

in which $\eta^{[u]}_{1i}$, $\eta^{[u]}_{\gamma_{i}}$, and $\eta^{[u]}_{2i}$ represent the slope of the first stage, the measurement at the knot, and the slope of the second stage, respectively. The corresponding factor loadings $\bm{\Lambda}^{[u]}_{i}$, a $J\times 3$ matrix, is expressed as

$\displaystyle\bm{\Lambda}^{[u]}_{i}=\left(\begin{array}[]{rrr}\min(0,t_{ij}-\gamma^{[u]})&1&\max(0,t_{ij}-\gamma^{[u]})\end{array}\right)$

$\displaystyle(j=1,\cdots,J).$

where the subscript $i$ indicates that we build the model in the framework of individual measurement occasions. Additionally, $\bm{\epsilon}^{[u]}_{i}$ is a $J\times 1$ vector of residuals of the $i^{th}$ individual.

To define a mediation model with a baseline predictor, a longitudinal mediator, and a longitudinal outcome, we regress the growth factors of the mediator process on the predictor for the $i^{th}$ individual

$$\bm{\eta}^{[m]}_{i}=\bm{\alpha}^{[m]}+\bm{B}^{[x\rightarrow{m}]}\times x_{i}+\bm{\zeta}^{[m]}_{i},$$

(3)

where $x_{i}$, which is either continuous or binary, is the baseline covariate of the $i^{th}$ individual, $\bm{\alpha}^{[m]}$ is a $3\times 1$ vector of growth factor intercepts of the mediator process (which is the mean vector of growth factors of the mediator if the covariate is centered), and $\bm{B}^{[x\rightarrow{m}]}$ is a $3\times 1$ vector of regression coefficients from the predictor to the first slope, to the measurement at the knot, and to the second slope of the mediator process, that is

$$\bm{B}^{[x\rightarrow{m}]}=\left(\begin{array}[]{rrr}\beta^{[x\rightarrow{m}]}_{1}&\beta^{[x\rightarrow{m}]}_{\gamma}&\beta^{[x\rightarrow{m}]}_{2}\end{array}\right)^{T}.$$

(4)

Additionally, we need to regress the growth factors of the outcome process on the predictor and the growth factors of the mediator process

$$\bm{\eta}^{[y]}_{i}=\bm{\alpha}^{[y]}+\bm{B}^{[x\rightarrow{y}]}\times x_{i}+\bm{B}^{[m\rightarrow{y}]}\times\bm{\eta}^{[m]}_{i}+\bm{\zeta}^{[y]}_{i},$$

(5)

where $\bm{\alpha}^{[y]}$ is a $3\times 1$ vector of growth factor intercepts of the outcome process, $\bm{B}^{[x\rightarrow{y}]}$ is a $3\times 1$ vector of regression coefficients from the predictor to the first slope, to the measurement at the knot, and to the second slope of the outcome process, that is

$$\bm{B}^{[x\rightarrow{y}]}=\left(\begin{array}[]{rrr}\beta^{[x\rightarrow{y}]}_{1}&\beta^{[x\rightarrow{y}]}_{\gamma}&\beta^{[x\rightarrow{y}]}_{2}\end{array}\right)^{T}.$$

(6)

Additionally, $\bm{B}^{[m\rightarrow{y}]}$ is a $3\times 3$ matrix of regression coefficients from the growth factors of the mediator process to those of the outcome process. The coefficient matrix can be expressed as

$$\bm{B}^{[m\rightarrow{y}]}=\left(\begin{array}[]{rrr}\beta^{[m\rightarrow{y}]}_{11}&0&0\\ \beta^{[m\rightarrow{y}]}_{1\gamma}&\beta^{[m\rightarrow{y}]}_{\gamma\gamma}&0\\ \beta^{[m\rightarrow{y}]}_{12}&\beta^{[m\rightarrow{y}]}_{\gamma 2}&\beta^{[m\rightarrow{y}]}_{22}\\ \end{array}\right).$$

(7)

In Equations (4), (6), (7), the superscript and subscript of $\beta$ together define the path of the corresponding coefficient. For example, $\beta^{[x\rightarrow{m}]}_{1}$ is the path coefficient from the baseline covariate to the first slope of the mediator process. Similarly, $\beta^{[m\rightarrow{y}]}_{1\gamma}$ is the coefficient from the first slope of the mediator process to the knot measurement of the outcome process. Additionally, $\bm{\zeta}^{[u]}_{i}$ is a $3\times 1$ vector of deviations of the $i^{th}$ individual from the mean values of the variable-specific growth factors. We list six possible paths of indirect effects of the predictor on the outcome process in the specified model in Table 1.

In this section, we extend the univariate linear-linear latent growth curve model to investigate a mediational process with a longitudinal predictor, mediator, and outcome. Suppose all three variables take a linear-linear functional form with an unknown fixed knot. In such situations, we can construct a longitudinal mediation model to explore how the change in covariate influences the change in the outcome directly and indirectly through its effect on the change in the mediator. The longitudinal mediation model can be specified as

$$\begin{pmatrix}\bm{x}_{i}\\ \bm{m}_{i}\\ \bm{y}_{i}\end{pmatrix}=\begin{pmatrix}\bm{\Lambda}_{i}^{[x]}&\bm{0}&\bm{0}\\ \bm{0}&\bm{\Lambda}_{i}^{[m]}&\bm{0}\\ \bm{0}&\bm{0}&\bm{\Lambda}_{i}^{[y]}\end{pmatrix}\times\begin{pmatrix}\bm{\eta}^{[x]}_{i}\\ \bm{\eta}^{[m]}_{i}\\ \bm{\eta}^{[y]}_{i}\end{pmatrix}+\begin{pmatrix}\bm{\epsilon}^{[x]}_{i}\\ \bm{\epsilon}^{[m]}_{i}\\ \bm{\epsilon}^{[y]}_{i}\end{pmatrix}.$$

(8)

which defines a longitudinal process for the predictor in addition to the mediator and outcome processes specified in Equation (2). In Equation (8), $\bm{x}_{i}$ is a $J\times 1$ vector of the repeated measures of the predictor for the $i^{th}$ individual. In addition, $\bm{\eta}^{[x]}_{i}$, $\bm{\Lambda}^{[x]}_{i}$, and $\bm{\epsilon}^{[x]}_{i}$ are the growth factors (a $3\times 1$ vector), the corresponding factor loadings (a $J\times 3$ matrix), and the residuals (a $J\times 1$ vector) of the covariate process, respectively. We then write the growth factors of the covariate process as deviations from the corresponding growth factor means

$$\bm{\eta}^{[x]}_{i}=\bm{\mu}^{[x]}_{\bm{\eta}}+\bm{\zeta}^{[x]}_{i},$$

(9)

where $\bm{\mu}^{[x]}_{\bm{\eta}}$ and $\bm{\zeta}^{[x]}_{i}$ are the mean vector of growth factors (a $3\times 1$ vector) and the deviations (a $3\times 1$ vector) of the covariate process. The growth factors of the mediator process are regressed on the those of the covariate, as shown below

$$\bm{\eta}^{[m]}_{i}=\bm{\alpha}^{[m]}+\bm{B}^{[x\rightarrow{m}]}\times\bm{\eta}^{[x]}_{i}+\bm{\zeta}^{[m]}_{i},$$

(10)

where $\bm{\alpha}^{[m]}$ is a $3\times 1$ vector of growth factor intercepts of the mediator, $\bm{B}^{[x\rightarrow{m}]}$ is a $3\times 3$ matrix of regression coefficients from the growth factors of the covariate process to those of the mediator process

$$\bm{B}^{[x\rightarrow{m}]}=\left(\begin{array}[]{rrr}\beta^{[x\rightarrow{m}]}_{11}&0&0\\ \beta^{[x\rightarrow{m}]}_{1\gamma}&\beta^{[x\rightarrow{m}]}_{\gamma\gamma}&0\\ \beta^{[x\rightarrow{m}]}_{12}&\beta^{[x\rightarrow{m}]}_{\gamma 2}&\beta^{[x\rightarrow{m}]}_{22}\\ \end{array}\right).$$

(11)

Similarly, the growth factors of the outcome process are regressed on those of the covariate and mediator

$$\bm{\eta}^{[y]}_{i}=\bm{\alpha}^{[y]}+\bm{B}^{[x\rightarrow{y}]}\times\bm{\eta}^{[x]}_{i}+\bm{B}^{[m\rightarrow{y}]}\times\bm{\eta}^{[m]}_{i}+\bm{\zeta}^{[y]}_{i},$$

(12)

where $\bm{\alpha}^{[y]}$ is a $3\times 1$ vector of growth factor intercepts of the outcome, $\bm{B}^{[x\rightarrow{y}]}$ ($\bm{B}^{[m\rightarrow{y}]}$) is a $3\times 3$ matrix of regression coefficients from the growth factors of the covariate (mediator) process to those of the outcome process. The coefficient matrix $\bm{B}^{[x\rightarrow{y}]}$ is defined as

$$\bm{B}^{[x\rightarrow{y}]}=\left(\begin{array}[]{rrr}\beta^{[x\rightarrow{y}]}_{11}&0&0\\ \beta^{[x\rightarrow{y}]}_{1\gamma}&\beta^{[x\rightarrow{y}]}_{\gamma\gamma}&0\\ \beta^{[x\rightarrow{y}]}_{12}&\beta^{[x\rightarrow{y}]}_{\gamma 2}&\beta^{[x\rightarrow{y}]}_{22}\\ \end{array}\right)$$

(13)

and $\bm{B}^{[m\rightarrow{y}]}$ has the same expression as such in Equation (7). Ten possible paths of indirect effects of the predictor process on the outcome process in the specified model are also listed in Table 1.

To simplify estimation, we make the following four assumptions. First, we assume that the covariate’s growth factors are normally distributed. Second, we assume that the mediator’s growth factors are normally distributed conditional on the baseline predictor (or growth factors of the predictor process). The third assumption is that the outcome’s growth factors are normally distributed conditional on the baseline covariate (or growth factors of the covariate process) and the growth factors of the mediator process. Accordingly, $\bm{\zeta}^{[u]}_{i}\sim\text{MVN}(\bm{0},\bm{\Psi}^{[u]}_{\bm{\eta}})$ ($u=x,m,y$), where $\bm{\Psi}^{[u]}_{\bm{\eta}}$ is a $3\times 3$ (unexplained) variance-covariance matrix of growth factors of the variable u. We also assume that the individual residuals $\bm{\epsilon}^{[u]}_{i}$ are identical and independent normal distributions over time, and the covariances between residuals are homogeneous over time. We define $\bm{I}$ as a $J\times J$ identity matrix, and for the first and second proposed models, the corresponding variance-covariance matrix of residuals can be expressed as

$$\begin{pmatrix}\bm{\epsilon}^{[m]}_{i}\\ \bm{\epsilon}^{[y]}_{i}\end{pmatrix}\sim\text{MVN}\bigg{(}\bm{0},\begin{pmatrix}\theta^{[m]}_{\epsilon}\bm{I}&\theta^{[my]}_{\epsilon}\bm{I}\\ &\theta^{[y]}_{\epsilon}\bm{I}\end{pmatrix}\bigg{)},$$

and

$$\begin{pmatrix}\bm{\epsilon}^{[x]}_{i}\\ \bm{\epsilon}^{[m]}_{i}\\ \bm{\epsilon}^{[y]}_{i}\end{pmatrix}\sim\text{MVN}\bigg{(}\bm{0},\begin{pmatrix}\theta^{[x]}_{\epsilon}\bm{I}&\theta^{[xm]}_{\epsilon}\bm{I}&\theta^{[xy]}_{\epsilon}\bm{I}\\ &\theta^{[m]}_{\epsilon}\bm{I}&\theta^{[my]}_{\epsilon}\bm{I}\\ &&\theta^{[y]}_{\epsilon}\bm{I}\end{pmatrix}\bigg{)},$$

respectively.

The parameters in the first proposed model include the mean and variance of the covariate, the intercepts ($\bm{\alpha}^{[u]}$), and unexplained variance-covariance ($\bm{\Psi}^{[u]}_{\bm{\eta}}$) of the growth factors of the mediator (outcome) process, the coefficients from the covariate to the mediator (outcome) growth factors and those from the mediator growth factors to the outcome growth factors, the mediator (outcome) residual variance, and the residual covariance. We define $\bm{\Theta}_{1}$ as

$$\bm{\Theta}_{1}=\{\mu_{x},\Phi_{x},\bm{\alpha}^{[u]},\bm{B}^{[x\rightarrow{m}]},\bm{B}^{[x\rightarrow{y}]},\bm{B}^{[m\rightarrow{y}]},\bm{\Psi}^{[u]}_{\bm{\eta}},\theta^{[u]}_{\epsilon},\theta^{[my]}_{\epsilon}\}\ (u=m,y)$$

to list the parameters specified in the first proposed longitudinal mediation model, where $\bm{B}^{[x\rightarrow{m}]}$ and $\bm{B}^{[x\rightarrow{y}]}$ are those defined in Equations (4) and (6), respectively.

The parameters in the second proposed model include the mean vector ($\bm{\mu}^{[x]}_{\bm{\eta}}$) and variance-covariance matrix ($\bm{\Psi}^{[x]}_{\bm{\eta}}$) of the growth factors of the predictor process, the intercepts ($\bm{\alpha}^{[u]}$) and unexplained variance-covariance ($\bm{\Psi}^{[u]}_{\bm{\eta}}$) of the growth factors of the mediator (outcome) process, the coefficients between growth factors of three longitudinal processes, the residual variances and covariances. We define $\bm{\Theta}_{2}$ as

$$\bm{\Theta}_{2}=\{\bm{\mu}^{[x]}_{\bm{\eta}},\bm{\Psi}^{[x]}_{\bm{\eta}},\bm{\alpha}^{[u]},\bm{B}^{[x\rightarrow{m}]},\bm{B}^{[x\rightarrow{y}]},\bm{B}^{[m\rightarrow{y}]},\bm{\Psi}^{[u]}_{\bm{\eta}},\theta^{[x]}_{\epsilon},\theta^{[u]}_{\epsilon},\theta^{[xm]}_{\epsilon},\theta^{[xy]}_{\epsilon},\theta^{[my]}_{\epsilon}\}\ (u=m,y)$$

and list the parameters specified in the second longitudinal mediation model, where $\bm{B}^{[x\rightarrow{m}]}$ and $\bm{B}^{[x\rightarrow{y}]}$ are those defined in Equations (11) and (13), respectively.

We estimate $\bm{\Theta}_{1}$ and $\bm{\Theta}_{2}$ using full information maximum likelihood (FIML) to account for the individual measurement occasions and the potential heterogeneity of individual contributions to the likelihood function. In this work, we constructed the proposed longitudinal mediation models using the R package OpenMx with CSOLNP optimizer (Pritikin et al.,, 2015; Neale et al.,, 2016; Boker et al.,, 2020; Hunter,, 2018).

In addition to the parameters that can be estimated directly, in practice, we are also interested in estimating the conditional mean vector and variance-covariance matrix of the growth factors of the mediator (outcome) process and indirect effects as well as the total effects listed in Table 1. When fitting the model using the R package OpenMx, we can specify the additional parameters in the function mxAlgebra() (Boker et al.,, 2020). We need to provide the algebraic expressions of the conditional mean vector and variance-covariance matrix of growth factors of the mediator (outcome) process and those in Table 1; then OpenMx is capable of computing the point estimates along with their standard errors of these additional parameters (that is realized by the Delta Method). We provide the expressions of the conditional mean vector and variance-covariance matrix of growth factors of the mediator (outcome) process for the two models in Appendix Appendix A and Appendix Appendix B, respectively. Additionally, we provide OpenMx code in the online appendix (https://github.com/Veronica0206/Extension_projects) and a demonstration to show how to build the proposed models and estimate these additional parameters.

We provide the detailed simulation designs for the two longitudinal mediation models in Table S1 and Table S2. The parameters of the most interest in the proposed model are the coefficients $\bm{B}^{[x\rightarrow{m}]}$, $\bm{B}^{[x\rightarrow{y}]}$ and $\bm{B}^{[m\rightarrow{y}]}$, based on which we obtain direct, indirect, and total effects that the baseline predictor (or predictor process) has on the outcome process. The conditions hypothesized to influence the estimation of these coefficients and other model parameters include sample size, the number of repeated measurements, the knot locations, the ratio of indirect effects to direct effects, shapes of trajectories, and measurement precision. Accordingly, we fixed the factors, for example, the mean vectors and variance-covariance matrices of the growth factors of longitudinal processes, that presumably do not affect the model performance meaningfully.

The primary objective of the simulation study is to evaluate how the proportion of mediated effects affects the proposed models. Accordingly, for the first mediation model, we fixed $\bm{B}^{[x\rightarrow{y}]}$ so that the predictor account for $13\%$ variability of the growth factors of the outcome process (i.e., the medium level according to Cohen, (1988, Chapter 9)). Additionally, we considered three levels of $\bm{B}^{[x\rightarrow{m}]}$ so that the predictor explains zero, medium (i.e., $13\%$), and substantial (i.e., $26\%$) variability of the growth factors of the mediator process. We then manipulated $\bm{B}^{[m\rightarrow{y}]}$ to further adjust the proportion of mediated effects (see the setup of coefficients for each condition in Table S1 in the Online Supplementary Document). For the elements in $\bm{B}^{[m\rightarrow{y}]}$, we considered the diagonal coefficients and off-diagonal coefficients separately since the diagonal elements are related to the immediate effects²²2The immediate effects could be, for example, the coefficient from the first slope of the mediator process to the first slope of the outcome process or that from the knot measurement of the mediator process to the knot measurement of the outcome process., while the off-diagonal elements are related to the delayed effects³³3The delayed effects could be, for example, the coefficients from the first slope of the mediator process to the knot measurement of the outcome process. in a mediational process. For all three levels of $\bm{B}^{[x\rightarrow{m}]}$, we considered three levels of indirectly immediate effects and two levels of indirectly delayed effects. In addition, for $\bm{B}^{[x\rightarrow{m}]}$ where the predictor explains $26\%$ variability of the growth factors of the mediator process, we considered two additional conditions to examine how varying path coefficients with fixed mediated effects affect model performance. Among all conditions, the maximum mediated effects from $x$ to the outcome growth factors are provided in the footnote of Table S1. Note that in the condition with the maximum mediated effects, the sum of the indirect effects achieved $48\%$ of the total effect, which is comparable to other methodological work such as Cheong et al., (2003).

Similarly, in the second longitudinal mediation model, we fixed $\bm{B}^{[x\rightarrow{y}]}$ so that the growth factors of the predictor process account for $13\%$ variability of the corresponding immediate growth factor and $3.25\%$ variability of each delayed growth factor of the outcome process. We considered two levels of $\bm{B}^{[x\rightarrow{m}]}$ so that the growth factors of the predictor process account for $13\%$ ($26\%$) variability of the corresponding immediate growth factor and $3.25\%$ ($6.50\%$) variability of each delayed growth factor of the mediator process (see the setup of coefficients for each condition in Table S1 in the Online Supplementary Document). For each level of $\bm{B}^{[x\rightarrow{m}]}$, we assessed two levels of non-zero coefficients for the elements in $\bm{B}^{[m\rightarrow{y}]}$. Among all conditions, the maximum mediated effects from the growth factors of the predictor process to those of the outcome processes are provided in Table S2. Note that in the condition with the maximum mediated effects, the sum of the indirect effects achieved $51\%$ of the total effect, which is comparable to other methodological work such as Cheong et al., (2003).

Another important factor in the simulation design to evaluate longitudinal models is the number of repeated measures since the proposed model is employed to analyze a longitudinal data set. Intuitively, the model should perform better with more repeated measurements. We also realized that the transition time to the second stage of the covariate (mediator) process, ideally, should occur no later than that of the mediator (outcome) process in longitudinal mediation models. We then decided to test two different levels of the number of measurements: six and ten. We selected six as the minimum number of repeated measurements to ensure the proposed model was fully identified⁴⁴4Bollen and Curran, (2005, Chapter 4) has shown that the latent growth model with linear-linear functional form can be identified with at least five waves with a knot at the midway of the study duration, though no studies provided information on identification for the model with an unknown knot.. For the conditions with six repeated measures, we set the transition time to the second stage of the mediator and outcome process at halfway through the study duration ($\mu^{[u]}_{\gamma}=2.5$, $u=x,m,y$). We considered the other condition, ten measurement occasions, for two reasons. On the one hand, we wanted to examine whether an increasing number of repeated measures would improve model performance. More importantly, this condition allowed us to evaluate the model performance under two scenarios: (1) the knot of each longitudinal process is in the middle of the study duration (i.e., $\mu^{[u]}_{\gamma}=4.5$, $u=x,m,y$); (2) the transition time of the mediator (covariate) process occurs earlier than that of the outcome (mediator) process (i.e., $\mu^{[m]}_{\gamma}=3.5$ and $\mu^{[y]}_{\gamma}=5.5$ for the first model, while $\mu^{[x]}_{\gamma}=3.5$, $\mu^{[m]}_{\gamma}=3.5$ and $\mu^{[y]}_{\gamma}=5.5$ for the second model). In addition, to account for individual measurement occasions, we allowed a time window with width $(-0.25,+0.25)$ around each wave, which is considered a ‘medium’ deviation as in Coulombe et al., (2015). Moreover, we examined three common change patterns, two levels of residual variances ($\theta^{[u]}_{\epsilon}=1$ or $2$) and set the residual correlation as $0.3$. We also examined the models at two levels of sample size, $n=200$ and $n=500$.

We evaluate the proposed longitudinal mediation models through the performance measures, including the relative bias, the empirical standard error (SE), the relative root-mean-square error (RMSE), and the empirical coverage probability for a nominal $95\%$ confidence interval of each parameter of the proposed models. Table 2 lists the definitions and estimates of these four performance measures. Specifically, the relative bias and empirical SE quantify whether the model, on average, targets the population values and whether estimates are precise, respectively. The relative RMSE integrates the bias and the precision metric into one measure. The coverage probability tells how well the interval estimate covers the corresponding population value.

We decided to replicate the simulation study $S=1,000$ by an empirical approach following Morris et al., (2019). Specifically, we run a pilot simulation study, the standard error of bias of all parameters except the (unexplained) variance of knot measurement of each longitudinal process (i.e., $\psi_{\gamma\gamma}^{[u]}$, $u=x,m,y$) was less than $0.15$. Therefore, we need at least $900$ replications to keep the Monte Carlo standard error of the bias below $0.005$⁵⁵5$\text{Monte Carlo SE(Bias)}=\sqrt{Var(\hat{\theta})/S}$ (Morris et al.,, 2019).. Out of more conservative consideration, we decided to proceed with $S=1,000$.

We carried out the following general steps for the simulation study of the proposed mediation models:

1. 1.

   Generate the baseline predictor (or growth factors of the predictor process), growth factors of the mediator process, and those of the outcome process simultaneously with the prespecified mean vectors and variance-covariance matrices (details are provided in Appendix Appendix C) using the R package MASS (Venables and Ripley,, 2002),

2. 2.

   Generate a scaled and equally-spaced time structure with $J$ waves $t_{j}$ and obtain individual measurement occasions: $t_{ij}\sim U(t_{j}-0.25,t_{j}+0.25)$ by allowing a time-window with width $(-0.25,0.25)$ around each wave,

3. 3.

   Calculate factor loadings for each individual of the (predictor or) mediator or outcome process from individual measurement occasions and the process-specific knot location,

4. 4.

   Calculate the values of the repeated measurements for each longitudinal process from corresponding growth factors, factor loadings, as well as residual variances and covariance,

5. 5.

   Implement the proposed models on the generated data set, estimate the parameters, and construct corresponding $95\%$ Wald CIs,

6. 6.

   Repeat the steps as mentioned above until having $1,000$ convergent solutions.

We first examined the convergence⁶⁶6The convergence was defined as arriving at OpenMx status code $0$ that indicates a successful optimization until up to $10$ runs with different sets of initial values (Neale et al.,, 2016). rate under each condition for each model before evaluating how the proposed longitudinal mediation models performed. The proposed model converged satisfactorily. For the first model, where we have baseline covariate and longitudinal mediator and outcome, $365$ conditions out of $396$ conditions reported a $100\%$ convergence rate. Only $31$ condition(s) reported up to five non-convergence replications. For the second model, in which we have longitudinal covariate, mediator, and outcome, 166 conditions out of a total 180 conditions reported a $100\%$ convergence rate. Only $14$ condition(s) reported up to five non-convergence replications.

We present the performance metrics, including the relative bias, empirical standard error (SE), relative root-mean-square error (RMSE), and empirical coverage probability (CP) for a nominal $95\%$ confidence interval of each parameter for the proposed models. For each parameter of interest, given the size of the conditions and parameters, we first obtained its four performance measures over $1,000$ replications under each condition in the simulation design. We then summarized these metrics across conditions as the corresponding median and range⁷⁷7For relative biases/RMSEs, we removed the conditions with coefficients with zero population values when summarizing the related summary statistics.. We provide the summary of the four performance metrics for the two models in Table S3-Table S5. These tables show that both models are capable of estimating parameters unbiasedly and precisely and exhibiting target confidence interval coverage in general.

Both models yielded unbiased point estimates along with small empirical standard errors generally. Specifically, for both models, the magnitude of relative biases of the growth factor means was under $0.004$, suggesting that we are able to have unbiased mean trajectories from both models. Additionally, for the first model, the magnitude of relative biases of the coefficients from the covariate to the mediator (outcome) growth factors was under $0.03$. Yet some coefficients from mediator growth factors to outcome growth factors, for example, $\beta^{[m\rightarrow{y}]}_{12}$ and $\beta^{[m\rightarrow{y}]}_{22}$ exhibit some bias (i.e., the relative biases were greater than $10\%$). These biased coefficients further led to biased estimates for the mediated effects. For the second model, the magnitude of the coefficients from covariate growth factors to mediator (outcome) growth factors except $\beta^{[x\rightarrow{m}]}_{12}$ ($\beta^{[x\rightarrow{y}]}_{12}$) was under $0.08$. Some coefficients from the mediator growth factors to the outcome growth factors exhibit some bias.

For the first model, we provide Figures 3 and 4 to further examine the relative bias pattern for the coefficients and indirect effects with some bias greater than $10\%$, respectively. From Figure 3, we noticed that the estimates of $\beta^{[m\rightarrow{y}]}_{12}$ and $\beta^{[m\rightarrow{y}]}_{22}$ were satisfactory in general (i.e., less than $10\%$ relative bias). Most factors we considered in the simulation design, such as the trajectory shapes and knot locations, did not affect the relative bias meaningfully. These biased estimates were generated under the conditions with a smaller sample size ($n=200$), a shorter study duration (i.e., six measurements), and a larger residual variance (i.e., $\theta^{[u]}=2$). There are two additional findings from this figure. First, the first model tended to overestimate the path coefficient from the first slope of the mediator to the second slope of the outcome, yet to underestimate the path coefficient from the second slope of the mediator to the second slope of the outcome. Second, the magnitude of relative biases was relatively small under the conditions where there was a relatively greater relationship between the mediator and the outcome. From Figure 4, we also noticed that the mediated effect involved $\beta^{[m\rightarrow{y}]}_{12}$ was likely to be overestimated, while that involved $\beta^{[m\rightarrow{y}]}_{22}$ was likely to be underestimated. In addition, we noticed that varying path coefficients with fixed mediated effects might slightly affect the performance of the first model.

For the second model, we plot Figures 5 and 6 to further examine the relative bias pattern for the coefficients and indirect effects with some bias greater than $10\%$, respectively. From Figure 5, we noticed that the estimates of these coefficients were satisfactory in general (i.e., less than $10\%$). Some factors we considered in the simulation design, such as the trajectory shapes, the knot locations, and the magnitude of indirect effects, did not affect the relative bias meaningfully. These biased estimates were generated under the conditions with a smaller sample size ($n=200$), a shorter study duration (i.e., six measurements), and a larger residual variance (i.e., $\theta^{[u]}=2$). Similar to Model 1, Model 2 was likely to overestimate path coefficients from a first slope to a second slope, but to underestimate the coefficients between two second slopes. From Figure 6, we noticed that the mediated effects involved path coefficients from a first slope to a second slope were tended to be overestimated. In contrast, those involved path coefficients between second slopes were tended to be underestimated.

Moreover, estimates obtained from the proposed models were precise: the magnitude of empirical standard errors of the slope- or knot-related parameters was under $0.20$, although this value of the parameters related to the knot measurement could achieve $0.57$. The relatively large empirical SEs of the parameters related to the knot measurement were due to the large scale of their population values (the population value of intercept means was around $100$).

The proposed models are capable of estimating parameters accurately. For both models, the magnitude of relative RMSEs of growth factor means was under $0.09$, and for knot was under $0.05$. The relative RMSE magnitude of the path coefficients, indirect effects, and direct effects was relatively large. In addition, both models performed well regarding empirical coverage since the coverage probabilities of all parameters were around $0.95$. We noticed that the coverage probabilities of indirect effects could achieve $100\%$, which is greater than the nominal coverage probability ($95\%$). As stated earlier, we constructed a $95\%$ Wald confidence interval with the SE obtained by the Delta Method and generated by mxAlgebra() automatically for each indirect effect, which should be very similar to the Sobel SE (Sobel,, 1982, 1986; Baron and Kenny,, 1986). Accordingly, the issue of overestimated SE is in line with an earlier study (Cheung,, 2009) and within our expectations to see.

Ideally, as stated earlier, the knot of the covariate process, indicating a milestone or an event in the process, occurs first, followed by the knot of the mediator process and then that of the outcome process. Therefore, we can specify the paths from the covariate knot to the other two knots and the path from the mediator knot to the outcome knot if we observe this temporal order. Otherwise, these paths between the knots may not be reasonable since we usually do not expect a later event predict an earlier one. Therefore, in this section, we first constructed a univariate latent growth curve model with the linear-linear functional form to estimate the fixed knot for each longitudinal process and then decide on possible paths for the proposed models. The estimated knot of the development of reading, mathematics, and science ability was around $93$, $100$, and $100$ months, respectively, suggesting that we can construct longitudinal models with the paths between these knots.

In this section, we built the first longitudinal mediation model to assess how the baseline attentional focus affects the development of reading ability and then the development of mathematics ability. As shown in Figure 7(a), the estimates obtained from the first model lead to a model trajectory that sufficiently captures the smooth line of observed trajectories for each ability. Table 3 presents the estimates of parameters of interest for the first model. Post-knot development of both abilities slowed down substantially. The transition to the slower rate occurred earlier in reading ability ($93$ months) than in mathematics ability ($99$ months). Baseline attentional focus positively affects the pre-knot development and the knot measurement of reading ability, and the knot measurement of mathematics ability. For example, the knot measurement of reading ability and mathematics ability increased $5.897$ and $2.277$ with one standardized unit increase in baseline attentional focus. Moreover, a child who developed more rapidly in reading ability tended to develop more rapidly in mathematics; a child who had higher scores at the knot of reading development tended to have higher mathematics scores at its knot. Additionally, pre-knot development of reading affects the knot measurement of mathematics development positively (estimate was $2.871$ with the p-value of $0.0453$).

In terms of mediated effects, the baseline attentional focus is associated with the pre-knot development, knot measurement, and post-knot development of the mathematics process indirectly through the corresponding counterpart of the reading process. In addition, we noticed that the effects on the post-knot development of both abilities of the baseline attentional focus were trivial.

In this section, we built the second longitudinal mediation model to evaluate how the development of reading ability affects the development of mathematics ability and then the development of science ability. From Figure 7(b), the estimates obtained from the second model also lead to a model-implied trajectory that sufficiently captures the smooth line of observed values for each disciplinary subject. We list the estimates of parameters of interest for the second model in Table 4, where we can see post-knot development of science ability slowed down slightly. The change to the slower rate occurred earliest in reading ability ($93$ months), followed by mathematics ability ($99$ months), and then science ability ($99$ months). In general, a child who developed more rapidly in reading ability tended to develop more rapidly in mathematics (science) at an early stage, and a child who had higher scores at the knot of reading development tended to have higher mathematics (science) scores at the corresponding knot. We also noticed that pre-knot development of mathematics impacted the knot measurement of science IRT scores negatively, which seems counter-intuitive. We will provide possible explanations in the Discussion section.

In terms of mediated effects, the pre-knot development, knot measurement, and post-knot development of the reading ability positively affected the corresponding values of science ability through the counterpart of mathematics ability. Additionally, pre-knot development of reading ability influenced the knot measurement of mathematics scores and the knot measurement of science scores. The knot measurement of reading ability also impacted the post-knot development of mathematics skills and the post-knot development of science skills. In terms of total effects, all three growth factors of reading development affected the corresponding immediate and delayed growth factors of science development. For example, the total effects of the pre-knot development of reading ability on the pre-knot development, knot measurement, and post-knot development of science ability were $0.339$, $7.347$, and $-0.330$, respectively.

In this section, we provide a set of recommendations for empirical researchers who are interested in employing the proposed models. First, we recommend building a univariate growth curve model before constructing the proposed models, as we did in the Application section. These univariate growth curve models allow us to assess whether the proposed functional form can capture the underlying change patterns of each longitudinal process. More importantly, we obtain the estimated knot of each longitudinal process, which allows for examining whether the temporal order of the covariate knot, mediator knot, and outcome knot is satisfied. If we have the temporal order as in the Application section, it is reasonable to add the paths between these knots. Otherwise, we may want to remove the paths to avoid a logical fallacy. For example, if the outcome knot occurred earlier than the mediator knot, we may not want to have the path from the mediator knot to the outcome knot in longitudinal mediation models.

Second, we constructed full models in the simulation study to examine model performance in estimating each possible parameter within research interests. However, we suggest only using the models with all possible paths for exploratory purposes. We recommend constructing a parsimonious model for a more confirmatory analysis due to two considerations. First, it helps avoid the potential collinearity issue in a confirmatory model, which we encountered in the Application section. In the second example, we found that the direct effects on the knot measurement of science ability of the pre-knot development of reading and mathematics ability were $9.455$ and $-8.667$, respectively. Although pre-knot development of mathematics ability might negatively affect the knot measurement of science ability, a more reasonable explanation is that the pre-knot development of reading and mathematics abilities were highly correlated. Therefore, the estimate of either individual effect was invalid, although the bundle of the pre-knot development of the two abilities can still predict the knot measurement of science ability. Accordingly, the paths in a more confirmatory model, where the individual effect is of interest, need to be carefully determined by specific research questions. More importantly, a parsimonious model with fewer paths may help model identification by decreasing the likelihood of the appearance of a not positive-definite sample covariance matrix, a common error in all models in the SEM framework. One underlying mechanism is related to the proposition regarding the positive definite property of the symmetric block matrix, as demonstrated in Appendix Appendix C).

Third, the models proposed and the analyses performs in the two examples in the Application section were pedagogical. We analyzed how the baseline attentional focus affects mathematics development through reading development by constructing the first model. Although we can maintain the temporal order, we cannot claim causality. The alternative pathway that is attentional focus affects the mathematics ability and then influences the reading ability is also plausible. Similarly, an alternative pathway for the second model could be that the development of mathematics development affects science development via the mediation of reading development. Moreover, in the proposed models, only the relationship between those between-construct growth factors is captured by unidirectional paths; that is, the development during the first phase of one ability is not considered a predictor of the development during the second phase of the same ability. Therefore, Equations A.2, A.4, B.2 and B.4 are helpful to obtain the covariance (correlation) between such within-construct growth factors if it is of research interest. In addition, similar to all mediation models, other confounders, such as learning approach or socioeconomic status, could also explain the relationship.

Last, researchers should be aware of the multiple statistical tests conducted when fitting a longitudinal mediation model and should consider adjusting hypothesis tests to control the Type I error rate or false discovery rate. Furthermore, decisions about the importance of insights should not only be based on p-values but also consider effect sizes, prior evidence, and alternative explanations (Wasserstein et al.,, 2019). For instance, in the first example, the effects size of the baseline attentional focus on the post-knot development of reading ability was $-0.047$. Although it is statistically significant at the $0.05$ level, this finding may not be meaningful in practice due to the trivial effect size. Additionally, one should carefully interpret the effect of the pre-knot slope of the covariate to the post-knot slope of the mediator (outcome), especially under conditions with fewer measurement occasions and less precise measurement where the estimates could exhibit some bias greater than $10\%$ based on our simulation study.

There are several directions to consider for future studies. First, the present study allows for employing the bilinear spline functional form to approximate nonlinear change patterns. This functional form is versatile and valuable for exploring longitudinal processes where different change rates correspond to different stages. More importantly, the estimated knot, which indicates a milestone or an event of a longitudinal process, is especially useful in a mediation model because it helps avoid a logical fallacy and ensures consistency between theory and statistical methods. Although the linear-linear piecewise functional form is sufficient in most cases, such as a developmental process with earlier and later stages in psychological phenomena or a recovery process with short-term and long-term recovery periods in biomedical fields, a functional form with more linear segments may also prove useful in practice. The proposed models can also be extended accordingly.

Second, we constructed the Wald confidence interval (CI) with the SE obtained by the Delta Method for the indirect effects and the total effects in the simulation studies and real-world data analyses. In the simulation studies, we noticed that the coverage probabilities could achieve $100\%$, greater than the nominal coverage probability ($95\%$). Another possible future direction is to construct other types of confidence intervals, such as the bootstrap CIs. Additionally, it is worth conducting simulation studies to evaluate these alternative methods for constructing CIs.

Third, we may not want to include all possible paths to construct a parsimonious model and avoid the potential collinearity issue. Although we recommend building the reduced model driven by specific research questions, it is still worth developing hypothesis testing for comparing the full model and reduced models, especially in an exploratory stage where even empirical researchers only have vague assumptions about causal relationships. In addition, the development of longitudinal mediation models that allow for penalization and regularization, either through $L_{1}$ norm (i.e., least absolute shrinkage and selection operator, LASSO (Santosa and Symes,, 1986)), or $L_{2}$ norm (i.e., ridge adjustment (Hoerl and Kennard, 1970b, ; Hoerl and Kennard, 1970a, )), or both (i.e., elastic net regularization (Zou and Hastie,, 2005)), to automatically select the most important paths and/or shrink the estimates for highly correlated paths is an alternative solution to the issue of collinearity and overparameterization. For a specific research question, after the most important paths and then the most important mediated effects are identified, another possible future direction is to perform power analysis with these identified paths and more levels of sample size and provide some guidelines for future studies.

To summarize, in this article, we propose two longitudinal mediation models to explore multiple nonlinear longitudinal processes. We can evaluate how the baseline covariate or the covariate process affects the outcome process through the mediator process with the proposed models. As discussed above, the proposed methods can be further extended in practice and further examined in methodology.