Obtaining Interpretable Parameters From Reparameterized Longitudinal Models: Transformation Matrices Between Growth Factors in Two Parameter Spaces

Abstract

This study proposes transformation functions and matrices between coefficients in the original and reparameterized parameter spaces for an existing linear-linear piecewise model to derive the interpretable coefficients directly related to the underlying change pattern. Additionally, the study extends the existing model to allow individual measurement occasions and investigates predictors for individual differences in change patterns. We present the proposed methods with simulation studies and a real-world data analysis. Our simulation study demonstrates that the method can generally provide an unbiased and accurate point estimate and appropriate confidence interval coverage for each parameter. The empirical analysis shows that the model can estimate the growth factor coefficients and path coefficients directly related to the underlying developmental process, thereby providing meaningful interpretation.

Keywords

linear spline growth models unknown knots individual measurement occasions time-invariant covariates simulation studies

1. Introduction

Longitudinal studies of change are popular in various disciplines to evaluate individual growth over time. If a process under investigation is followed for a long enough time duration, it is likely to exhibit some degree of nonlinear change in which the curve has a nonconstant relationship to recorded time. One possible model to examine nonlinear individual change patterns is a piecewise linear latent growth model (Chou et al., 2004; Cudeck & Harring, 2010; Harring et al., 2006; Kohli, 2011; Kohli & Harring, 2013; Kohli et al., 2013; Sterba, 2014), also known as the linear spline growth model (Grimm et al., 2016, Chapter 11). By using at least two linear pieces, the linear spline (or piecewise linear) functional form can approximate more complex underlying change patterns. Previous studies have demonstrated broad application potential for linear spline trajectory models. For example, researchers have employed the model to assess postsurgical rehabilitation (Dumenci et al., 2019; Riddle et al., 2015), the learning process of a specific task (Cudeck & Klebe, 2002), the onset of alcohol abuse (Li et al., 2001), and intellectual development (Marcoulides, 2018).

One major challenge of analyzing linear piecewise change patterns is to decide the location of the change-point or “knot” that indicates the occurrence of the change in the developmental process. Although empirical researchers can specify the knot location a priori driven by domain knowledge, for example, Dumenci et al. (2019), Flora (2008), Riddle et al. (2015), and Sterba (2014), the knot is unknown in most cases. Marcoulides (2018) proposed to fit a pool of candidate models with knots at different locations and select the “best” one using the Bayesian information criterion (BIC). By building a multiphase mixed-effects model, Cudeck and Klebe (2002) demonstrated that the change-point can also be a random coefficient that allows individuals to transition from one phase to another at different time points. Earlier studies have demonstrated how to estimate knots using frequentist (Cudeck & du Toit, 2003; Harring et al., 2006; Kohli, 2011; Kohli & Harring, 2013; Kohli et al., 2013; Kwok et al., 2010; Preacher & Hancock, 2015) and Bayesian (Dominicus et al., 2008; Kohli, Hughes, et al., 2015; Lock et al., 2018; McArdle & Wang, 2008; Muniz Terrera et al., 2011; Wang & McArdle, 2008) mixed-effects models and growth models.

One intuitive longitudinal model with a linear spline functional form is a bilinear spline growth model (or linear-linear piecewise model) with an unknown knot. It can approximate complex nonlinear change patterns to answer specific research questions about when one developmental stage ends and the other one starts, as well as the development rate of each phase (Kohli, Hughes, et al., 2015; Kohli, Sullivan, et al., 2015; Sullivan et al., 2017). The piecewise linear-linear growth curve model can be fit in the mixed-effects modeling framework and the latent growth curve (LGC) modeling framework to estimate a knot. In the mixed-effects modeling framework, Kohli et al. (2016) estimated the knot with both fixed and random effects along with coefficients of the initial status and the rate of change and extended the model to the finite mixture modeling framework by implementing an R routine fitPMM (Zopluoglu et al., 2014) constructed on the marginal maximum likelihood estimation (MLE) procedure described in Du Toit and Cudeck (2009). Note that the marginal MLE procedure was proposed because the linear-linear piecewise growth model with a random knot is a nonlinear random coefficient model (NRCM), of which the joint likelihood function does not have a closed form and therefore cannot be calculated explicitly (Du Toit & Cudeck, 2009; Kohli et al., 2016).

In the LGC framework, Harring et al. (2006) developed a linear-linear piecewise model with an unknown fixed knot with a unified functional form of two linear pieces through reparameterization, which is realized by reexpressing growth factors as linear combinations of their original forms. Kohli (2011), Kohli et al. (2013), and Grimm et al. (2016, Chapter 11) derived the growth coefficients that are directly related to the underlying change patterns from the estimated coefficients. By implementing this model, Kohli et al. (2013) examined the development in each stage and estimated a fixed knot for procedural learning task research.

Preacher and Hancock (2015) extended the model with an unknown fixed knot to estimate the knot’s variance simultaneously, in which the knot location is allowed to vary and then is viewed as an additional growth factor besides the intercept and two slopes in the LGC modeling framework. More importantly, they considered incorporating obesity status to predict varying knot locations. They demonstrated that the linear-linear functional form with an unknown random knot is useful by reexamining a data set containing longitudinal records of plasma phosphate. To specify the NRCM in the structural equation modeling (SEM) framework, Preacher and Hancock (2015) linearized the model by reparameterizing growth factors through the Taylor series expansion (Browne & du Toit, 1991). One drawback of this linearized model lies in that it is fit in a reparameterized framework. Some reparameterized coefficients were no longer directly related to the specific research questions answered by employing such functional form and therefore lacked meaningful, substantive interpretation.

Some existing methodological studies, for example, Bauer (2003) and Curran (2003), have shown that the LGC modeling framework is mathematically equivalent to the mixed-effects model in most cases. However, researchers implemented different approximation methods on the NRCMs in the two frameworks: the marginal MLE procedure and the Taylor series expansion for mixed-effects and LGC modeling framework, respectively, as stated above. Specifically, the marginal MLE is utilized during the optimization process. The Taylor series expansion is employed to rewrite the repeated outcome as a linear combination of growth factors (i.e., prior to the optimization process), so that the normal theory of MLE is still applicable, which allows for manageable extensions, for example, adding covariates to account for variability in growth factors.¹ Based on this consideration, we propose inverse-transformation matrices for the model proposed by Preacher and Hancock (2015) to derive the growth factor coefficients in the original frame so that all coefficients are related to the specific research interest answered by utilizing the linear-linear piecewise functional form and interpretable, although these coefficients can be directly estimated in the mixed-effects modeling framework. In addition, we extend the model proposed by Preacher and Hancock (2015) by allowing time-invariant covariates (TICs) to predict individual differences in the initial status and rate-of-change of each stage, so that these covariates may account for the heterogeneity in trajectories instead of that only in knot locations. More importantly, we conducted simulation studies to assess how well the Taylor series approximation, which is necessary for reparameterization, works under various conditions.

Last, we propose the model in the framework of individual measurement occasions (Cook & Ware, 1983; Finkel et al., 2003; Mehta & West, 2000). This can occur when time is measured precisely (e.g., exact date and time of measurement) or self-initiated or randomly assigned responses. For instance, in an adolescent smoking study, participants were asked to complete a questionnaire on pocket computers immediately after smoking (Hedeker et al., 2006). Earlier studies, for example, Preacher and Hancock (2015) and Sterba (2014), have demonstrated that one possible way to fit growth models with individual measurement occasions is the definition variable approach, in which the “definition variables” are defined as observed variables that are employed to adjust model parameters to individual-specific values (Mehta & Neale, 2005; Mehta & West, 2000). In our case, these individual-specific values are individual measurement occasions.

The remainder of this article is organized as follows. In the Method section, we first present the model specification of the bilinear spline growth model with TICs to estimate a knot with its variability (the full model). Next, we introduce how to reparameterize growth factors and the corresponding path coefficients from the TICs to the growth factors to make them estimable in the LGC modeling framework and transform them back after estimation for interpretation purposes. Following that, we propose one parsimonious model to estimate the knot without variability (the reduced model). We then describe the model estimation as well as the Monte Carlo simulation design for model evaluation. In the Result section, we evaluate the model performance with regard to the nonconvergence rate, the proportion of improper solutions, as well as the performance measures, including the relative bias, the empirical standard error (SE), the relative root mean squared error (RMSE), and the empirical coverage probabilities for a nominal 95% confidence interval (CI) of each parameter of interest. In the Application section, by applying the proposed model to a longitudinal data set of mathematics achievement scores from the Early Childhood Longitudinal Study, Kindergarten Class of 2010–2011 (ECLS-K: 2011), we give a set of feasible recommendations for real-world practice. Finally, a broad discussion is presented regarding practical considerations, methodological considerations, and future directions.

2. Method

2.1. Model Specification

This section briefly describes a linear-linear growth curve model with an unknown random knot and TICs used to analyze nonlinear change patterns, estimate the knot with its variability, and predict individual differences in trajectories. As shown in Figure 1, the linear-linear piecewise latent growth model, which is an extension of the linear growth model, specifies a separate linear function for each of the two stages of the developmental process for each individual. We can express the measure for the ith individual at the $j th$ time point in the framework of individual measurement occasions as

y_{i j} = \{\begin{array}{l} η_{0 i} + η_{1 i} t_{i j} + ε_{i j} & t_{i j} \leq γ_{i} \\ η_{0 i} + η_{1 i} γ_{i} + η_{2 i} (t_{i j} - γ_{i}) + ε_{i j} & t_{i j} > γ_{i} \end{array},

where $y_{i j}$ and $t_{i j}$ are the measurement and measurement time of the $i th$ individual at time j, respectively. In Equation 1, $η_{0 i}$ , $η_{1 i}$ , $η_{2 i}$ , and $γ_{i}$ are individual-level intercept, first slope, second slope, and knot, which are usually called “growth factors” and all together determine the functional form of the growth curve of $y_{i}$ .

Figure 1.

Within-individual change over time with bilinear spline functional form.

Note that we cannot fit the linear-linear piecewise change pattern specified in Equation 1 in the LGC modeling framework for two reasons. First, the outcome variable $y_{i j}$ does not have the same model expression pre- and postknot. Second, the model in Equation 1 specifies a nonlinear relationship between the $y_{i j}$ and the growth factor (or random coefficient in the mixed-effects modeling framework) $γ_{i}$ and then cannot be estimated in the SEM framework directly (Grimm et al., 2016, Chapter 12). As shown in Tishler and Zang (1981), Seber and Wild (2003), Grimm et al. (2016, Chapter 11), and Liu (2019), the expression of repeated outcome in Equation 1 before and after the knot can be unified by reparameterization (see Online Appendix A.1 for the detailed derivation of reparameterization). Following Browne and du Toit (1991), Grimm et al. (2016, Chapter 12), and Liu (2019), we then write the repeated outcome as a linear combination of all four growth factors by the Taylor series expansion (see Online Appendix A.2 for the detailed derivation of Taylor series expansion). Accordingly, the model specified in Equation 1 can be written as a standard LGC model with reparameterized growth factors

y_{i} \approx Λ_{i}^{'} η_{i}^{'} + ε_{i},

where y _i is a $J \times 1$ vector of the repeated outcomes for the ith person (in which J is the number of measurements), $η_{i}^{'}$ is a $4 \times 1$ vector of individual-level reparameterized growth factors, and $Λ_{i}^{'}$ , which is a function of measurement occasions and the unknown knot, is a $J \times 4$ matrix of factor loadings.² Note that the subscript i in $Λ_{i}^{'}$ indicates that we construct the proposed model in the framework of individual measurement occasions. Additionally, $ε_{i}$ is a $J \times 1$ vector of residuals for the $i th$ person. For the $i th$ individual, the reparameterized growth factors (the measurement at the knot, the mean of two slopes, the half-difference of two slopes, and the knot) and their corresponding factor loadings are

η_{i}^{'} = {(\begin{array}{r} η_{0 i}^{'} & η_{1 i}^{'} & η_{2 i}^{'} & δ_{i} \end{array})}^{T} = {(\begin{array}{r} η_{0 i} + γ_{i} η_{1 i} & \frac{η_{1 i} + η_{2 i}}{2} & \frac{η_{2 i} - η_{1 i}}{2} & γ_{i} - μ_{γ} \end{array})}^{T},

and

\begin{array}{l} Λ_{i}^{'} = (\begin{array}{r} 1 & t_{i j} - μ_{γ} & | t_{i j} - μ_{γ} | & - μ_{η_{2}}^{'} - \frac{μ_{η_{2}}^{'} (t_{i j} - μ_{γ})}{| t_{i j} - μ_{γ} |} \end{array}) & (j = 1, \dots, J), \end{array}

respectively, where $μ_{γ}$ is the knot mean and $δ_{i}$ is the deviation from the knot mean of the ith individual. We provide the detailed deviation of the reparameterized growth factors and corresponding factor loadings in Online Appendix A.1 and A.2. The reparameterized growth factors of the ith individual can be further regressed on TICs,

η_{i}^{'} = α^{'} + B^{'} X_{i} + ζ_{i}^{'},

where $α^{'}$ is a $4 \times 1$ vector of the reparameterized growth factor intercepts (which is the mean vector of reparameterized growth factors if the TICs are centered), $B^{'}$ is a $4 \times c$ matrix of regression coefficients (in which c is the number of the TICs) from the TICs to the reparameterized growth factors, and X _i is a $c \times 1$ vector of covariates of the ith individual. Note that the covariates can be either continuous or binary. Additionally, $ζ_{i}^{'}$ is a $4 \times 1$ vector of deviations of the ith subject from the reparameterized growth factor means. To simplify estimation, we assume that the reparameterized growth factors are normally distributed conditional on the individual-level covariates X _i . Accordingly, $ζ_{i}^{'} \sim MVN (0, Ψ_{η}^{'})$ , where $Ψ_{η}^{'}$ is a $4 \times 4$ unexplained variance-covariance matrix of reparameterized growth factors. We also assume that individual residuals, $ε_{i}$ in Equation (2), are independent and identically normally distributed over time, that is, $ε_{i} \sim MVN (0, θ_{ε} I)$ , where I is a $J \times J$ identity matrix. Note that in Equation 5, the intercepts (i.e., $α^{'}$ ) and coefficients of reparameterized growth factors (i.e., $B^{'}$ ) as well as the unexplained growth factor variances (i.e., $Ψ_{η}^{'}$ ) are no longer directly related to the underlying developmental process and therefore lack meaningful, substantive interpretation. In the next section, we then demonstrate how to derive the coefficients in the original parameter space, so that all coefficients are interpretable.

2.2. Transformation and Inverse-Transformation Between Two Parameter Spaces

When fitting a model, especially a complex model, in the SEM framework, selecting a proper set of starting values can improve the likelihood of convergence and accelerate the computational process. Generally, descriptive statistics and visualization are tools for researchers to decide suitable starting values for the parameters of interest. However, it may not be straightforward for the reparameterized coefficients. Accordingly, the transformation from parameters in the original frame to those in the reparameterized setting help decide appropriate starting values. More importantly, some reparameterized coefficients are not straightforward to relate to developmental theory and then need to be transformed back to be interpretable after estimation. So the inverse transformation that helps derive the coefficients in the original frame is also needed.

By using the multivariate delta method (Lehmann & Casella, 1998, Chapter 1), we can transform the point estimates (that are the MLEs as we use full information maximum likelihood [FIML]) of coefficients in the two parameter-spaces. Suppose $f : R^{4} \to R^{4}$ is a function, which takes a point $η_{i} \in R^{4}$ as input and produces the vector $f (η_{i}) \in R^{4}$ (i.e., $η_{i}^{'} \in R^{4}$ ) as output. By the multivariate delta method,

η_{i}^{'} = f (η_{i}) \sim N (f (μ_{η}), \nabla_{f} (μ_{η}) Ψ_{η} \nabla_{f}^{T} (μ_{η})),

where $μ_{η}$ and $Ψ_{η}$ are the mean vector and variance-covariance matrix of original growth factors, respectively. As Equation 3, f is defined as

f (η_{i}) = {(\begin{array}{r} η_{0 i} + γ_{i} η_{1 i} & \frac{η_{1 i} + η_{2 i}}{2} & \frac{η_{2 i} - η_{1 i}}{2} & γ_{i} - μ_{γ} \end{array})}^{T} .

Similarly, suppose $h : R^{4} \to R^{4}$ is a function, which takes a point $η_{i}^{'} \in R^{4}$ as input and produces the vector $h (η_{i}^{'}) \in R^{4}$ (i.e., $η_{i} \in R^{4}$ ) as output. By the multivariate delta method,

η_{i} = h (η_{i}^{'}) \sim N (h (μ_{η}^{'}), \nabla_{h} (μ_{η}^{'}) Ψ_{η}^{'} \nabla_{h}^{T} (μ_{η}^{'})),

where $μ_{η}^{'}$ and $Ψ_{η}^{'}$ are the mean vector and variance-covariance matrix of reparameterized growth factors, respectively. With a straightforward algebraic calculation from Equation 3, we can express $η_{i}$ as a linear combination of $η_{i}^{'}$ and define h as

h (η_{i}^{'}) = {(\begin{array}{r} η_{0 i}^{'} - γ_{i} η_{1 i}^{'} + γ_{i} η_{2 i}^{'} & η_{1 i}^{'} - η_{2 i}^{'} & η_{1 i}^{'} + η_{2 i}^{'} & δ_{i} + μ_{γ} \end{array})}^{T} .

Based on Equations 6 and 7, we can make the transformation between the growth factor means of two parameter spaces by $μ_{η}^{'} \approx f (μ_{η})$ and $μ_{η} \approx h (μ_{η}^{'})$ , respectively. We can also define the transformation matrix $\nabla_{f} (μ_{η})$ and $\nabla_{h} (μ_{η}^{'})$ between the variance-covariance matrix of two parameter spaces as

\begin{array}{l} Ψ_{η}^{'} & \approx \nabla_{f} (μ_{η}) Ψ_{η} \nabla_{f}^{T} (μ_{η}) \\ = (\begin{array}{r} 1 & μ_{γ} & 0 & μ_{η_{1}} \\ 0 & 0.5 & 0.5 & 0 \\ 0 & - 0.5 & 0.5 & 0 \\ 0 & 0 & 0 & 1 \end{array}) Ψ_{η} {(\begin{array}{r} 1 & μ_{γ} & 0 & μ_{η_{1}} \\ 0 & 0.5 & 0.5 & 0 \\ 0 & - 0.5 & 0.5 & 0 \\ 0 & 0 & 0 & 1 \end{array})}^{T} \end{array},

and

\begin{array}{l} Ψ_{η} & \approx \nabla_{h} (μ_{η}^{'}) Ψ_{η}^{'} \nabla_{h}^{T} (μ_{η}^{'}) \\ = (\begin{array}{r} 1 & - μ_{γ} & μ_{γ} & 0 \\ 0 & 1 & - 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array}) Ψ_{η}^{'} {(\begin{array}{r} 1 & - μ_{γ} & μ_{γ} & 0 \\ 0 & 1 & - 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array})}^{T}, \end{array}

respectively.

When regressing growth factors on the individual-level covariates as we did in Equation 5, we also need to rereparameterize the growth factor intercepts and path coefficients. We center the covariates to simplify the calculation. On the one hand, the vector $α^{'}$ is equivalent to the mean vector of the reparameterized growth factor $μ_{η}^{'}$ when all covariates are centered. Accordingly, we can make the transformation between α and $α^{'}$ by $α^{'} = μ_{η}^{'} \approx f (μ_{η})$ and $α = μ_{η} \approx h (μ_{η}^{'})$ . On the other hand, through centering the covariates, we simplify the transformation and inverse-transformation matrices between the original path coefficients B and the reparameterized path coefficients $B^{'}$ as $\nabla_{f} (μ_{η})$ and $\nabla_{h} (μ_{η}^{'})$ , respectively. The detailed derivation is shown in Online Appendix A.3.

We obtain the SEs around the points estimates for the parameters in the original frame (i.e., the derived parameters) as follows. We first have the point estimates and the corresponding SEs of the coefficients in the reparameterized frame (i.e., the free parameters) directly. Note that these SEs are based on Fisher Information since we employ FIML. The inverse-transformation function and matrix described earlier define the relationships between the free parameters and the derived parameters. There are two common approaches to get the SEs around the point estimates for the derived parameters: an analytical approach and the delta method. Note that these two methods are equivalent if all the defined relationships are linear. We present how to calculate the SEs for the derived parameters in Online Appendix A.4. Additionally, SEM computer programs, such as Mplus and the R package OpenMx, are capable of providing the point estimates and corresponding SEs for the derived parameters automatically when we define the relationships.

When fitting the model using software such as the R package OpenMx, which allows for matrix calculation on those estimates from the model by the function mxAlgebra() (Boker et al., 2020), we only need to provide $h (μ_{η}^{'})$ and $\nabla_{h} (μ_{η}^{'})$ , and the package will derive the coefficients in the original setting and generate their point estimates with corresponding SEs automatically. Other SEM software such as Mplus can also calculate new parameters to be derived from those estimated automatically by specifying their relationship in the MODEL CONSTRAINT command. However, we need to derive the expression for each cell of the mean vector and the variance-covariance matrix of the original growth factors manually (instead of having software to do this) because Mplus does not allow for matrix algebra. All these expressions are in Online Appendix A.5.

2.3. Model Estimation

For the ith individual, the expected mean vector and the variance-covariance structure of the repeated measurements of the model given in Equations 2 and 5 can be expressed as

μ_{i} = Λ_{i}^{'} (α^{'} + B^{'} μ_{X}),

and

Σ_{i} = Λ_{i}^{'} Ψ_{η}^{'} Λ_{i}^{' T} + Λ_{i}^{'} B^{'} Φ B^{' T} Λ_{i}^{' T} + θ_{ε} I,

where $μ_{X}$ and $Φ$ are the mean vector ( $c \times 1$ ) and the variance-covariance matrix ( $c \times c$ ) of the TICs, respectively.

The parameters in the model given in Equations 2 and 5 include the mean vector and unexplained variance-covariance matrix of reparameterized growth factors, the reexpressed path coefficients as well as the means and the variance-covariance structure of covariates. Then, we can obtain the parameters in the original setting, as shown in Section 2.2. $Θ_{1}$ and $Θ_{1}^{'}$ are shown below

\begin{array}{l} Θ_{1} & = {α, Ψ_{η}, B, μ_{X}, Φ, θ_{ε}} \\ = {μ_{η_{0}}, μ_{η_{1}}, μ_{η_{2}}, μ_{γ}, ψ_{00}, ψ_{01}, ψ_{02}, ψ_{0 γ}, ψ_{11}, ψ_{12}, ψ_{1γ} {,ψ}_{22}, ψ_{2 γ}, ψ_{γ γ}, B, μ_{X}, Φ, θ_{ε}} \end{array},

and

\begin{array}{l} Θ_{1}^{'} & = {α^{'}, Ψ_{η}^{'}, B^{'}, μ_{X}, Φ, θ_{ε}} \\ = {μ_{η_{0}}^{′} {,μ}_{η_{1}}^{′} {,μ}_{η_{2}}^{′} {,μ}_{γ} {,ψ}_{00}^{′} {,ψ}_{01}^{′} {,ψ}_{02}^{′} {,ψ}_{0γ}^{′} {,ψ}_{11}^{′} {,ψ}_{12}^{′} {,ψ}_{1γ}^{′} {,ψ}_{22}^{′} {,ψ}_{2γ}^{′} {,ψ}_{γγ}, B^{'}, μ_{X}, Φ, θ_{ε}} \end{array},

list the parameters in the original and reparameterized frames, respectively.

We then use FIML to estimate $Θ_{1}^{'}$ due to the potential heterogeneity of individual contributions to the likelihood. The log-likelihood function of each individual and that of the overall sample can be expressed as

log l i k_{i} (Θ_{1}^{'} | y_{i}) = C - \frac{1}{2} ln | Σ_{i} | - \frac{1}{2} {(y_{i} - μ_{i})}^{T} Σ_{i}^{- 1} (y_{i} - μ_{i}),

and

log l i k (Θ_{1}^{'}) = \sum_{i = 1}^{n} log l i k_{i} (Θ_{1}^{'} | y_{i}),

respectively, where C is a constant, n is the number of individuals, and $μ_{i}$ and $Σ_{i}$ are the mean vector and the variance-covariance matrix of y _i , which have been defined in Equations 8 and 9, respectively. We construct the full model using the R package OpenMx with the optimizer CSOLNP (Boker et al., 2020; Hunter, 2018; Neale et al., 2016; Pritikin et al., 2015), with which we are able to fit the model and obtain interpretable coefficients efficiently, as shown in Section 2.2.

2.4. Reduced Model

With an assumption that the knot of each individual is roughly similar, we can fix the between-individual differences in the knot to 0 and construct a reduced model to estimate a knot without considering variability. We can specify the reduced model as

y_{i} = Λ_{i}^{'} η_{i}^{'} + ε_{i},

η_{i}^{'} = α^{'} + B^{'} X_{i} + ζ_{i}^{'},

where $η_{i}^{'}$ is a $3 \times 1$ vector of reparameterized growth factor (the measurement of the knot, the mean of two slopes, and the half-difference of two slopes), and $Λ_{i}$ , a function of time points and a fixed change-point, is a $J \times 3$ matrix of factor loadings. We can further express the reparameterized growth factors and corresponding factor loadings as

η_{i}^{'} = {(\begin{array}{r} η_{0 i}^{'} & η_{1 i}^{'} & η_{2 i}^{'} \end{array})}^{T} = {(\begin{array}{r} η_{0 i} + {γη}_{1 i} & \frac{η_{1 i} + η_{2 i}}{2} & \frac{η_{2 i} - η_{1 i}}{2} \end{array})}^{T},

and

\begin{array}{l} Λ_{i}^{'} = (\begin{array}{r} 1 & t_{i j} - γ & | t_{i j} - γ | \end{array}) & (j = 1, \dots, J) . \end{array}

The growth factor intercepts $α^{'}$ , path coefficients $B^{'}$ , and the deviation $ζ_{i}^{'}$ of the ith individual from the reparameterized growth factor means also reduce to be a $3 \times 1$ vector, $3 \times c$ matrix, and $3 \times 1$ vector, respectively. The (inverse-)transformation functions and matrices are also reduced accordingly. Specifically, we only need the first three entries of the functions f and h and then the first three rows and the first three columns in the matrices $\nabla_{f} (μ_{η})$ and $\nabla_{h} (μ_{η}^{'})$ , since only three growth factors need to be reparameterized.

For the ith individual, the expected mean vector and the variance-covariance matrix of the repeated outcomes of this model are

μ_{i} = Λ_{i}^{'} (α^{'} + B^{'} μ_{X}),

and

Σ_{i} = Λ_{i}^{'} Ψ_{η}^{'} Λ_{i}^{' T} + Λ_{i}^{'} B^{'} Φ B^{' T} Λ_{i}^{' T} + θ_{ε} I,

respectively. For this reduced model, $Θ_{2}$ and $Θ_{2}^{'}$ are defined as

\begin{array}{l} Θ_{2} & = {α, γ, Ψ_{η}, B, μ_{X}, Φ, θ_{ε}} \\ = {μ_{η_{0}} {, μ}_{η_{1}} {, μ}_{η_{2}} {, γ, ψ}_{00} {, ψ}_{01} {, ψ}_{02} {, ψ}_{11} {, ψ}_{12} {, ψ}_{22}, B, μ_{X}, Φ, θ_{ε}} \end{array},

and

\begin{array}{l} Θ_{2}^{'} & = {α^{'}, γ, Ψ_{η}^{'}, B^{'}, μ_{X}, Φ, θ_{ε}} \\ = {μ_{η_{0}}^{′} {, μ}_{η_{1}}^{′} {, μ}_{η_{2}}^{′} {, γ, ψ}_{00}^{′} {, ψ}_{01}^{′} {, ψ}_{02}^{′} {, ψ}_{11}^{′} {, ψ}_{12}^{′} {, ψ}_{22}^{′}, B^{'}, μ_{X}, Φ, θ_{ε}}, \end{array}

and they list the parameters in the original and reparameterized settings, respectively. By replacing $Θ_{1}^{'}$ in Equations 10 and 11 with $Θ_{2}^{'}$ and updating $μ_{i}$ and $Σ_{i}$ as such defined in Equations 12 and 13, we have the likelihood function of each individual and that of the overall sample. We build the reduced model using the R package OpenMx with the optimizer CSOLNP and employ FIML to estimate the parameters. We provide OpenMx code for the full model and its reduced version, as well as a demonstration in the first author’s GitHub website (https://github.com/Veronica0206/Dissertation_projects). For the researchers who are interested in using Mplus, we also provide Mplus 8 syntax for both models.

3. Model Evaluation

We evaluated the proposed method using a Monte Carlo simulation study with two goals. The first goal is to assess how the approximation based on the Taylor series expansion affects the performance measures, including the relative bias, the empirical SE, the relative RMSE, and the empirical coverage probability (CP) of a nominal 95% CI of each parameter. We provide the definitions and estimates of these performance metrics in Table 1. The second goal is to examine how the reduced model performed sufficiently well compared to the full model.

Table 1.

Performance Metric: Definitions and Estimates

Criteria	Definition	Estimate
Relative bias	$E_{\hat{θ}} (\hat{θ} - θ) / θ$	$\sum_{s = 1}^{S} ({\hat{θ}}_{s} - θ) / S θ$
Empirical standard error	$\sqrt{Var (θ^)}$	$\sqrt{\sum_{s = 1}^{S} {({\hat{θ}}_{s} - \bar{θ})}^{2} / (S - 1)}$
Relative root mean squared error	$\sqrt{E_{\hat{θ}} {(\hat{θ} - θ)}^{2}} / θ$	$\sqrt{\sum_{s = 1}^{S} {({\hat{θ}}_{s} - θ)}^{2} / S} / θ$
Coverage probability	$P r ({\hat{θ}}_{low} \leq θ \leq {\hat{θ}}_{upper})$	$\sum_{s = 1}^{S} I ({\hat{θ}}_{low, s} \leq θ \leq {\hat{θ}}_{upper, s}) / S$

Note. $θ$ = the population value of the parameter of interest

$\hat{θ}$ = the estimate of $θ$ ; S = the number of replications and set as $1, 000$ in our simulation study

$s = 1, \dots, S$ = indexes the replications of the simulation

${\hat{θ}}_{s}$ = the estimate of $θ$ from the $s th$ replication

$\bar{θ}$ = the mean of ${\hat{θ}}_{s}$ ‘s across replications

$I ()$ = an indicator function.

In the simulation study, we followed Morris et al. (2019) and decided the number of replications S = 1,000 by an empirical approach. We conducted a pilot simulation run and found that SEs of all coefficients except the unexplained intercept variance (i.e., $ψ_{00}$ ) were less than 0.15. As the most important performance metric in the simulation study was the (relative) bias, and since we wanted to keep its Monte Carlo SE ³ lower than 0.005, we needed at least 900 replications. We then decided to proceed with S = 1,000 for more conservative considerations.

3.1. Design of Simulation Study

We list all simulation conditions in Table 2. As mentioned earlier, the parameters of the most interest in the full model are the knot, its variance, and the path coefficients to the knot. The conditions hypothesized to influence the estimation of these knot parameters, along with other model parameters, include the number of individuals, the number of repeated measurements, the knot locations, the knot variances, shapes of trajectories, measurement precision, and the proportion of growth factor variances explained by TICs. Accordingly, we did not investigate some conditions that presumably would not affect the model performance meaningfully. For instance, the mean vector and variance-covariance matrix of the first three growth factors (i.e., the intercept and two slopes) usually change with the measurement scale and time scale in practice, and there is little theoretical reason to expect different time scales to have an impact on estimation (Biesanz et al., 2004). Accordingly, in the simulation design, we fixed the mean of intercept and first slope and only changed that of the second slope. We also fixed the variance-covariance matrix of the first three growth factors and kept the index of dispersion ( $σ^{2} /μ$ ) of each at one-tenth scale, guided by previous simulation studies for growth models (Bauer & Curran, 2003; Kohli, 2011; Kohli, Sullivan, et al., 2015). Further, the growth factors were set to be positively correlated to a moderate degree ( $ρ = .3$ ).

Table 2.

Simulation Design for the Models

Variables	Conditions
Fixed conditions
Mean of the intercept	$μ_{η_{0}} = 100$
Variance of the intercept	$ψ_{00} = 25$
Mean of the first slope	$μ_{η_{1}} = - 5$
Variance of slopes	$ψ_{11} = ψ_{22} = 1$
Correlations of growth factors	$ρ = .3$
Manipulated conditions
Partially crossed design
Time (t_j )	Six scaled and equally spaced $(j = 0 \dots, J - 1, J = 6)$
Time (t_j )	10 Scaled and equally spaced $(j = 0, \dots, J - 1, J = 10)$
(Mean of) the knot	$μ_{γ}$ at $t = 2.5$ for $J = 6$
(Mean of) the knot	$μ_{γ}$ at $t = 3.5$ or 4.5 or 5.5 for $J = 10$
Individual $t_{i j}$	$t_{i j} \sim U (t_{j} - Δ, t_{j} + Δ) (j = 0, \dots, J - 1; Δ = 0.25)$
Sample size	$n = 200$
Sample size	$n = 500$
Full factorial design
Standardized difference d_z	$μ_{η_{1}} - μ_{η_{2}} = \pm 1.6$ ( $1.0 -$ unit standardized difference)
	$μ_{η_{1}} - μ_{η_{2}} = \pm 2.4$ ( $1.5 -$ unit standardized difference)
	$μ_{η_{1}} - μ_{η_{2}} = \pm 3.2$ ( $2.0 -$ unit standardized difference)
Variance of knots	$ψ_{γγ} = 0$ ( $s d (γ) / (t_{j} - t_{j - 1}) = 0$ )
	$ψ_{γγ} = 0.09$ ( $s d (γ) / (t_{j} - t_{j - 1}) = 0.3$ )
	$ψ_{γγ} = 0.36$ ( $s d (γ) / (t_{j} - t_{j - 1}) = 0.6$ )
Coefficients ( $B$ )	TICs explain 13% variability of growth factors
Coefficients ( $B$ )	TICs explain 26% variability of growth factors
Residual variance	$θ_{ε} = 1$
Residual variance	$θ_{ε} = 2$

Note. TICs = time-invariant covariates.

For a model to analyze longitudinal data, the factor that we are most interested in is the number of repeated measures. Generally, the model should perform better if the study duration is longer, which we wanted to investigate by the simulation study. Another important factor for the linear-linear trajectory is the knot location. Intuitively, the model should perform the best if the knot is in the middle of study duration; it is our interest to test this conjecture as well. To this end, we chose two different levels of the number of measurements: six and 10. We selected six as the minimum number of repeated outcomes to make the full model fully identified.⁴ For the conditions with six repeated measures, we set the knot at halfway of study duration ( $μ_{γ} = 2.5$ ). The other level 10 was considered for two reasons. First, we wanted to evaluate whether a more extended study duration would improve model performance. Second, 10 measurements allowed us to place knots at different locations, say the halfway knot ( $μ_{γ} = 4.5$ ) as well as a left-shifted knot ( $μ_{γ} = 3.5$ ) and right-shifted knot ( $μ_{γ} = 5.5$ ) from the middle point of study duration. Around each wave, we allowed a time-window $(- 0.25, + 0.25)$ , which is a “medium” deviation as in Coulombe et al. (2015), for individual measurement occasions.

In order to fit the full model in the SEM framework and obtain interpretable coefficients, the approximation was introduced by the Taylor series expansion. As shown in Online Appendix A.2, one necessary condition to apply this approximation is that the knot variance approaches 0. We then investigated three levels of between-person differences in the knot to evaluate how the introduced approximation affects the model performance and compare the full model versus the reduced model. We set the knot standard deviation (SD) as 0, 0.3, and 0.6 as the zero, medium, and large difference in the knot location. Note that the magnitude of the knot SD is relative to the scale of measurement occasions.

Additionally, we examined how the shape of trajectories quantified by the standardized difference between two slopes affects the full model. We considered both positive and negative differences ( $η_{1 i} < η_{2 i}$ and $η_{1 i} > η_{2 i}$ , respectively) to assess model performance in the two cases as shown in Online Appendix, Figure A1. For each case, we considered three levels of standardized difference between two slopes to evaluate the models. We also considered two levels of effects of exogenous variables on endogenous growth factors,⁵ two levels of sample size, and two levels of measurement precision as shown in Table 2.

3.2. Data Generation and Simulation Step

For each condition that is listed in Table 2, the general steps of the simulation study for the models are outlined as follows:

Generate data for growth factors and TICs simultaneously using the R package MASS (Venables & Ripley, 2002; details are provided in Online Appendix A.6).

Generate the time structure with J scaled and equally spaced waves t_j and obtain individual measurement occasions: $t_{i j} \sim U (t_{j} - Δ, t_{j} + Δ)$ by allowing disturbances around each wave.

Calculate factor loadings, which are functions of individual measurement occasions and the knot, for each individual.

Calculate values of the repeated measurements based on growth factors, factor loadings, and residual variances.

Implement the full model and the reduced model on the simulated data, estimate the parameters, and construct corresponding 95% Wald’s CIs.

Repeat the above steps until after obtaining $1, 000$ convergent solutions.

The simulation study was conducted on a Linux Beowulf Cluster with R Version 3.4 and OpenMx Version 2.12.2.

4. Results

4.1. Model Convergence and Proper Solution

In this section, we first examined the convergence rate and component fit measures of each condition. The convergence⁶ is defined as to achieve OpenMx status code 0, which suggests a successful optimization, until up to 10 attempts with different sets of starting values (Neale et al., 2016). Convergence rate of each condition was investigated. Based on our simulation studies, the full model and its reduced version converged satisfactorily (the convergence rate of the full model achieved at least 95% while that of the reduced model was 100%). Of a total of 576 conditions, 490 conditions reported 100% convergence rate and 67 conditions reported convergence rate of 99% to 100%. The worst scenario in terms of the nonconvergence rate across all conditions was $46 / 1, 046$ , indicating that we needed to repeat the steps described in Section 3.2 $1, 046$ times to have $1, 000$ replications with a convergent solution. It occurred under the condition with the smaller sample size (i.e., $n = 200$ ), the shorter study duration (i.e., $J = 6$ ), less precise measurements (i.e., $θ_{ε} = 2$ ), the largest knot SD (i.e., $s d (γ) = 0.6$ ), and the smallest positive standardized difference between two slopes (i.e., $d_{z} = 1$ ).

We also conducted diagnostics to investigate “improper solutions” (referred to the estimates that are impossible in the population) following Bollen and Curran (2005, Chapter 2), such as estimates of growth factor variances less than 0 and/or estimates of correlations between growth factors beyond $[- 1, 1]$ . Table 3 demonstrates the number of improper solutions yielded by the full model under conditions with 10 repeated measures, the knot midway through the total time of observation, and 13% explained variance of the latent growth factors. The improper solutions include negative knot variances and its out-of-range (i.e., out of $[- 1, 1]$ ) correlations with any other growth factors. From Table 3, when the proposed model was correctly specified (i.e., under the conditions with a knot SD 0.3 or 0.6), the model is capable of producing proper solutions, although the number of improper solutions was relatively large if the model was overspecified under the conditions where the knot SD was 0.

Table 3.

Number of Improper Solutions Among 1,000 Convergent Replications Under Conditions With 10 Repeated Measures and Midway Knot, 13% Explained Variance

			$θ_{ε} = 1$		$θ_{ε} = 2$
			$n = 200$	$n = 500$	$n = 200$	$n = 500$
Small d_z	Positive d_z	$s d (γ) = 0$	$555 / / 48$ ^a	$516 / / 22$	$550 / / 40$	$507 / / 29$
		$s d (γ) = 0.3$	$117 / / 35$	$7 / / 11$	$278 / / 48$	$168 / / 21$
		$s d (γ) = 0.6$	$4 / / 3$	$0 / / 0$	$58 / / 25$	$7 / / 6$
	Negative d_z	$s d (γ) = 0$	$535 / / 37$	$542 / / 27$	$568 / / 44$	$531 / / 27$
		$s d (γ) = 0.3$	$106 / / 49$	$24 / / 14$	$279 / / 59$	$134 / / 39$
		$s d (γ) = 0.6$	$0 / / 8$	$0 / / 0$	$67 / / 45$	$8 / / 17$
Medium d_z	Positive d_z	$s d (γ) = 0$	$563 / / 49$	$553 / / 28$	$583 / / 39$	$573 / / 22$
		$s d (γ) = 0.3$	$17 / / 14$	$1 / / 0$	$141 / / 27$	$40 / / 11$
		$s d (γ) = 0.6$	$0 / / 1$	$0 / / 0$	$9 / / 2$	$0 / / 0$
	Negative d_z	$s d (γ) = 0$	$568 / / 42$	$542 / / 26$	$576 / / 50$	$513 / / 32$
		$s d (γ) = 0.3$	$19 / / 17$	$1 / / 0$	$130 / / 64$	$22 / / 32$
		$s d (γ) = 0.6$	$0 / / 0$	$0 / / 0$	$4 / / 23$	$0 / / 0$
Large d_z	Positive d_z	$s d (γ) = 0$	$555 / / 44$	$555 / / 36$	$562 / / 45$	$557 / / 19$
		$s d (γ) = 0.3$	$1 / / 0$	$0 / / 0$	$27 / / 19$	$0 / / 0$
		$s d (γ) = 0.6$	$0 / / 0$	$0 / / 0$	$0 / / 0$	$0 / / 0$
	Negative d_z	$s d (γ) = 0$	$548 / / 35$	$514 / / 30$	$555 / / 36$	$523 / / 33$
		$s d (γ) = 0.3$	$1 / / 2$	$0 / / 0$	$37 / / 39$	$2 / / 6$
		$s d (γ) = 0.6$	$0 / / 0$	$0 / / 0$	$0 / / 2$	$0 / / 0$

^a $555 / / 48$ indicates that among 1,000 convergent replications, we have 555 and 48 improper solutions due to negative knot variances and out-of-range correlations of the knot with other growth factors, respectively.

For the scenarios including knots with variability, we observed that the proper solution rate has positive associations with the magnitude of standardized difference between two slopes, the measurement precision, the sample size, and the knot SD. The sign of d_z , the shifted knot locations, or the proportion of explained variance of growth factors only affected the number of improper solutions slightly, although reduced follow-up time inflated this number. When such improper solutions occurred, we replaced the full model with its reduced version for the model evaluation.

4.2. Performance Measures

In this section, we present simulation results in terms of performance measures, including the relative bias, empirical SE, relative RMSE, and empirical CP for a nominal 95% CI for each parameter of interest. In general, the estimated coefficients from the full model were unbiased and accurate, with target coverage probabilities. We first provide the summary statistics (specifically, median and range) for each performance metric of each parameter across the conditions with 10 repeated measurements and the midway knot. Then, we discuss how the simulation conditions affect these performance metrics.

Tables 4 and 5 present the median (range) of the relative bias and empirical SE, respectively, for each parameter of interest across all conditions with 10 repeated measurements and the halfway knot for the full model and its reduced version. We first calculated each parameter’s relative bias/empirical SE across $1, 000$ replications under each condition and then summarized each parameter’s relative biases/empirical SEs across all conditions above using the median (range) of relative bias/empirical SE. As shown in Tables 4 and 5, both models generated unbiased point estimates with small empirical SEs to a large extent, and the full model performed better than its reduced version as the ranges of relative biases of the parameters from the full model were narrower than those from the reduced model. Specifically, for the full model, the magnitude of relative biases of the growth factor means was under 0.03, that of path coefficients to the intercept and two slopes was under .12, and that of unexplained variances of the intercept and two slopes was under .07. From Table 4, the full model may produce biased estimates for the unexplained knot variances and the path coefficients to the knot: The median values of relative biases of $β_{1γ}$ , $β_{2γ}$ , and $ψ_{γγ}$ were $- 0.2712$ , $- 0.2709$ , and 0.0403, respectively.

Table 4.

Median (Range) of Relative Bias of Each Parameter Across Conditions With 10 Repeated Measures and Midway Knot

Parameters	Reduced Model	Full Model
Parameters	Median (Range)	Median (Range)
Mean vector
$μ_{η_{0}}$	0.0000 ( $- 0.0005$ , 0.0006)	0.0000 ( $- 0.0006$ , 0.0006)
$μ_{η_{1}}$	0.0000 ( $- 0.0098$ , 0.0100)	0.0000 ( $- 0.0099$ , 0.0100)
$μ_{η_{2}}$	0.0000 ( $- 0.0063$ , 0.0280)	0.0000 ( $- 0.0062$ , 0.0273)
$μ_{γ}$	0.0000 ( $- 0.0014$ , 0.0009)	0.0000 ( $- 0.0017$ , 0.0006)
Path coefficient of x ₁
$β_{10}$	0.0030 ( $- 0.0823$ , 0.0773)	0.0027 ( $- 0.0345$ , 0.0367)
$β_{11}$	0.0023 ( $- 0.3069$ , 0.3055)	0.0010 ( $- 0.1113$ , 0.1052)
$β_{12}$	0.0016 ( $- 0.2984$ , 0.2970)	0.0002 ( $- 0.0940$ , 0.0946)
$β_{1γ}$	— ^a	$- 0.2712$ (NA,^b NA)
Path coefficient of x ₂
$β_{20}$	0.0006 ( $- 0.0697$ , 0.0690)	0.0007 ( $- 0.0216$ , 0.0250)
$β_{21}$	0.0006 ( $- 0.2901$ , 0.2950)	0.0000 ( $- 0.0882$ , 0.0933)
$β_{22}$	$- 0.0017$ ( $- 0.2888$ , 0.2901)	$- 0.0019$ ( $- 0.0901$ , 0.0878)
$β_{2γ}$	—	$- 0.2709$ (NA, NA)
Unexplained variance
$ψ_{00}$	$- 0.0107$ ( $- 0.0409$ , 0.0212)	$- 0.0107$ ( $- 0.0235$ , 0.0018)
$ψ_{11}$	$- 0.0070$ ( $- 0.0753$ , 0.1692)	$- 0.0138$ ( $- 0.0657$ , 0.0182)
$ψ_{22}$	$- 0.0080$ ( $- 0.0703$ , 0.1638)	$- 0.0130$ ( $- 0.0647$ , 0.0132)
$ψ_{γγ}$	—	0.0403 ( $- 0.4151$ , NA)

^a— indicates that the relative biases are not available from the reduced model. ^b NA indicates that the bounds of relative bias is not available. The model performance under the conditions with zero population value of difference in knot is of interest where the relative bias of those knot coefficients would go infinity.

Table 5.

Median (Range) of Empirical Standard Error of Each Parameter Across Conditions With 10 Repeated Measures and Midway Knot

Parameters	Reduced Model	Full Model
Parameters	Median (Range)	Median (Range)
Mean vector
$μ_{η_{0}}$	0.2576 (0.1903, 0.3550)	0.2578 (0.1903, 0.3547)
$μ_{η_{1}}$	0.0548 (0.0387, 0.0755)	0.0548 (0.0387, 0.0755)
$μ_{η_{2}}$	0.0552 (0.0400, 0.0755)	0.0552 (0.0400, 0.0762)
$μ_{γ}$	0.0418 (0.0173, 0.0854)	0.0412 (0.0173, 0.0922)
Path coefficient of x ₁
$β_{10}$	0.2761 (0.1957, 0.3734)	0.2768 (0.1967, 0.3735)
$β_{11}$	0.0579 (0.0412, 0.0800)	0.0583 (0.0412, 0.0806)
$β_{12}$	0.0574 (0.0412, 0.0819)	0.0583 (0.0412, 0.0819)
$β_{1γ}$	— ^a	0.0387 (0.0100, 0.0990)
Path coefficient of x ₂
$β_{20}$	0.2725 (0.1980, 0.3744)	0.2720 (0.1980, 0.3751)
$β_{21}$	0.0574 (0.0412, 0.0806)	0.0583 (0.0412, 0.0806)
$β_{22}$	0.0578 (0.0412, 0.0800)	0.0582 (0.0412, 0.0812)
$β_{2γ}$	—	0.0400 (0.0100, 0.1054)
Unexplained variance
$ψ_{00}$	1.6912 (1.1734, 2.3690)	1.6902 (1.1741, 2.3501)
$ψ_{11}$	0.0748 (0.0490, 0.1122)	0.0747 (0.0490, 0.1095)
$ψ_{22}$	0.0741 (0.0490, 0.1140)	0.0740 (0.0510, 0.1091)
$ψ_{γγ}$	—	0.0387 (0.0000, 0.1732)

^a— indicates that the empirical standard errors are not available from the reduced model.

We plotted the relative bias under each condition with 10 repeated measures and the midway knot for $ψ_{γγ}$ , $β_{1γ}$ , and $β_{2γ}$ in Figure 2 to further investigate the relative bias pattern for the unexplained knot variance and the path coefficients to the knot. We observed how the conditions we set up in the simulation design affected these estimates from the figure. First, under the conditions where the knot SD was 0.3 and/or the standardized difference between two slopes was medium or large (i.e., $d_{z} = 1.5$ or $d_{z} = 2.0$ ), relative biases were smaller. Second, the unexplained knot variance, $ψ_{γγ}$ , was underestimated under the conditions with 0.6 knot SD, whereas it could be overestimated under the conditions with 0.3 knot SD, and the path coefficients to the knot (i.e., $β_{1γ}$ and $β_{2γ}$ ) were underestimated under all investigated conditions. Third, other conditions, such as the sample size or the proportion of growth factor variance explained by the covariates, only influenced the relative biases slightly.

Figure 2.

Relative biases of path coefficients and residuals of knot under conditions with 10 repeated measures and midway knot.

From Table 5, estimates from the full model and its reduced version were precise: The magnitude of empirical SEs of slope or knot parameters were around 0.10, although those values of intercept parameters were relatively large: the empirical SE of the intercept mean and that of the unexplained intercept variance achieved 0.35 and 2.35, respectively.

Table 6 lists the median (range) of relative RMSE of each parameter for both models under the conditions with 10 repeated measures and the midway knot, which combines bias and precision to examine the point estimates holistically. From the table, both models estimated parameters accurately. The magnitude of relative RMSEs of growth factor means was under 0.05, and the median RMSE value of path coefficients to the intercept or slopes was around 0.20.

Table 6.

Median (Range) of Relative Root Mean Squared Error (RMSE) of Each Parameter Across Conditions With 10 Repeated Measures and Midway Knot

Parameters	Reduced Model	Full Model
Parameters	Median (Range)	Median (Range)
Mean vector
$μ_{η_{0}}$	0.0026 (0.0019, 0.0035)	0.0026 (0.0019, 0.0035)
$μ_{η_{1}}$	$- 0.0128$ ( $- 0.0178$ , $- 0.0078$ )	$- 0.0128$ ( $- 0.0178$ , $- 0.0078$ )
$μ_{η_{2}}$	$- 0.0118$ ( $- 0.0498$ , $- 0.0050$ )	$- 0.0120$ ( $- 0.0490$ , $- 0.0050$ )
$μ_{γ}$	0.0092 (0.0036, 0.0190)	0.0092 (0.0037, 0.0204)
Path coefficient of x ₁
$β_{10}$	0.2594 (0.1595, 0.4266)	0.2568 (0.1570, 0.4227)
$β_{11}$	0.2980 (0.1656, 0.5418)	0.2752 (0.1678, 0.4674)
$β_{12}$	0.2912 (0.1634, 0.5378)	0.2734 (0.1646, 0.4641)
$β_{1γ}$	— ^a	0.7689 (0.3162, NA ^b)
Path coefficient of x ₂
$β_{20}$	0.1747 (0.1063, 0.2865)	0.1722 (0.1059, 0.2830)
$β_{21}$	0.2234 (0.1104, 0.4097)	0.1860 (0.1114, 0.3068)
$β_{22}$	0.2234 (0.1091, 0.4122)	0.1834 (0.1104, 0.3061)
$β_{2γ}$	—	0.5625 (0.2179, NA)
Unexplained variance
$ψ_{00}$	0.0856 (0.0610, 0.1111)	0.0846 (0.0610, 0.1112)
$ψ_{11}$	0.1076 (0.0655, 0.1998)	0.0982 (0.0664, 0.1360)
$ψ_{22}$	0.1078 (0.0663, 0.2038)	0.0988 (0.0666, 0.1326)
$ψ_{γγ}$	—	0.6089 (0.2130, NA)

^a— indicates that the relative RMSEs are not available from the reduced model. ^b NA indicates that the upper bound of relative RMSE is not available. The model performance under the conditions with zero population value of difference in knot is of interest where the relative RMSE of those knot coefficients would go infinity.

Table 7 presents the median (range) of the CP for the full model and its reduced version under the conditions with 10 repeated measurements and the midway knot. Under the majority of conditions, the full model performed well in terms of empirical coverage as the median values of coverage probabilities of all parameters were around 0.95. We noticed that the coverage probabilities of knot variance and path coefficients to the knot could be conservative under the conditions with 0.6 knot SD.

Table 7.

Median (Range) of Coverage Probability of Each Parameter Across Conditions With 10 Repeated Measures and Midway Knot

Parameters	Reduced Model	Full Model
Parameters	Median (Range)	Median (Range)
Mean vector
$μ_{η_{0}}$	0.9470 (0.9280, 0.9680)	0.9459 (0.9203, 0.9677) ^a
$μ_{η_{1}}$	0.9430 (0.7800, 0.9660)	0.9418 (0.7660, 0.9744)
$μ_{η_{2}}$	0.9440 (0.7970, 0.9670)	0.9429 (0.7990, 0.9702)
$μ_{γ}$	0.9100 (0.7880, 0.9680)	0.9510 (0.9267, 0.9758)
Path coefficient of x ₁
$β_{10}$	0.9460 (0.9140, 0.9620)	0.9472 (0.8926, 0.9745)
$β_{11}$	0.9295 (0.6170, 0.9630)	0.9450 (0.9070, 0.9673)
$β_{12}$	0.9290 (0.6120, 0.9630)	0.9460 (0.9080, 0.9823)
$β_{1γ}$	— ^b	0.9463 (0.6640, 0.9751)
Path coefficient of x ₂
$β_{20}$	0.9440 (0.8950, 0.9590)	0.9479 (0.9220, 0.9685)
$β_{21}$	0.9090 (0.2780, 0.964)	0.9440 (0.8820, 0.9745)
$β_{22}$	0.9130 (0.2960, 0.9650)	0.9445 (0.8720, 0.9703)
$β_{2γ}$	—	0.9416 (0.3810, 0.9879)
Unexplained variance
$ψ_{00}$	0.9375 (0.8970, 0.9600)	0.9389 (0.9086, 0.9627)
$ψ_{11}$	0.9365 (0.4240, 0.9600)	0.9351 (0.8550, 0.9728)
$ψ_{22}$	0.9370 (0.4660, 0.9620)	0.9360 (0.8570, 0.9698)
$ψ_{γγ}$	—	0.9620 (0.0030, 0.9976)

^a For the full model, the reported coverage probabilities have four decimals since we calculated the coverage probabilities only based on the replications with proper solutions. ^b— indicates that the coverage probabilities are not available from the reduced model.

In summary, based on our simulation study, the full model performed better than the reduced model in terms of performance measures and generated unbiased estimates with small variance and proper target CP. As pointed out earlier, one necessary condition for applying the Taylor series approximation is that the knot variance approaches 0. Through the simulation study, we found that the full model was robust under the conditions with the medium level of knot SD (i.e., $s d (γ) = 0.3$ for the scaled measurement occasions) though it was less satisfactory when the knot SD was large, which is within our expectation to see. Factors, including sample size and measurement precision, also influenced model performance. Specifically, the larger sample size (i.e., $n = 500$ ) and more precise measurement ( $θ_{ε} = 1$ ) improved performance measures. The trajectory shape quantified by the standardized difference between two slopes affected the performance measures shown in Figure 2. We also examined the model performance beyond conditions summarized in Tables 4 –7 and found that shorter study duration (i.e., $J = 6$ ) and left- (right-) shifted knot locations only affected the performance measures slightly.

5. Application

This section demonstrates the use of the full model to approximate nonlinear trajectories and estimate a knot that varies across individuals. This application has two goals. Although researchers can employ the proposed model to estimate a knot and its variance of a developmental process directly, we realize that a considerable number of studies that examine the underlying change patterns might start from an exploratory stage, where researchers only have vague hypotheses in terms of the proper functional form of trajectories. One common approach for such exploratory studies is to fit a pool of candidate models and select the “best” one by prespecified criteria such as information criteria or the fit between the model implied trajectory and a smooth line of observed trajectories. We want to demonstrate this process in our application. Accordingly, we first fit the full and reduced LGC models with bilinear spline functional form with three common parametric change patterns: linear, quadratic, and Jenss-Bayley growth curve, demonstrate how to select the model, and then show how to interpret the coefficients obtained from the proposed model. Second, there are two possible time scales available in the data set: students’ age and their grade-in-school. We use children’s age (in months) to have individual measurement occasions in the main analysis and conduct a sensitivity analysis where students’ grade-in-school was used as the time scale. We hope to gain insights on how different time structures would affect estimation in practice through the sensitivity analysis. We randomly selected 400 students from the Early Childhood Longitudinal Study Kindergarten Cohort: 2010-11 (ECLS-K: 2011) with complete records of repeated mathematics item response theory (IRT) scaled scores, demographic information (sex, race, and age at each wave), school information (baseline school location and baseline school type), and socioeconomic status (baseline family income and the highest education level between parents).⁷

ECLS-K: 2011 is a nationally representative longitudinal sample of U.S. children enrolled in about 900 kindergarten programs starting from 2010−2011 school year. Children’s mathematics ability was assessed in nine waves: fall and spring of kindergarten (2010−2011), first (2011−2012) and second (2012−2013) grade, respectively, as well as spring of third (2014), fourth (2015), and fifth (2016), respectively. Only about 30% students were sampled in the fall of 2011 and 2012 (Lê et al., 2011). In the subsample, 51% of children were male, and 49% of children were female. Additionally, 44% of students were White, 5% were Black, 38% were Latinx, 8% were Asian, and 5% were others. We dichotomized the variable race to be White (44%) and others (56%). At the start of the study, 11% and 89% students were in private and public schools, respectively. We treated the variables school location (ranged between one and four), highest parents’ education (ranged between zero and eigth), and family income (ranged between one and 18) as continuous variables, and corresponding mean (SD) was $2.06 (1.00)$ , $5.18 (2.10)$ , and $11.56 (5.42)$ , respectively.

5.1. Main Analysis

In this section, we use children’s age (in months) as the time metric to have individual measurement occasions. We first fit the full and reduced linear-linear piecewise models, as well as three parametric growth curve functions, that is, linear, quadratic, and Jenss-Bayley. The graphical representations of the model implied curves together with the smooth line of the observed mathematics IRT scores are shown in Figure 3. From the figure, the nonlinear functional forms fit better than the linear function generally; and the models with bilinear spline change patterns fit better than the parametric nonlinear models to the first two or three measurement occasions. As shown in Figure 3, both the quadratic growth curve and the Jenss-Bayley tended to underestimate the mathematics achievement in the early stage.

Figure 3.

Model implied trajectory and smooth line of observed mathematics IRT scores.

Table 8 provides the obtained estimated likelihood, information criteria (including Akaike infromation criterion [AIC] and Bayesian Information Criterion [BIC]), and residuals of each LGC model. From the table, the model with the linear-linear functional form and a random knot has the largest estimated likelihood, smallest AIC, and smallest residual though its BIC is slightly larger than the LGC with the quadratic function. Moreover, from Table 8, the full model consistently performs better than its reduced version for this empirical example, indicating that we need to consider the variability of the knot, although the model implied trajectory of both captures the overall observed data well as shown in Figure 3. We then fit the full model to investigate the impact of covariates, such as sex, race, parents’ highest education (baseline), family income (baseline), school type (baseline), and school location (baseline), on the knot variance (and other growth factor variances). The model fit information is also provided in Table 8.

Table 8.

Summary of Model Fit Information for the Models

Model	−2ll	Akaike Infromation Criterion (AIC)	Bayesian Infromation Criterion (BIC)	# of Para.	Residual
Latent Growth Curve Models (Time Structure: Age in Months)
Linear	27,298	27,310	27,334	6	75.60
Quadratic	25,320	25,340	25,379	10	35.14
Jenss-Bayley	25,330	25,352	25,396	11	35.37
Linear-linear piecewise with a fixed knot	25,412	25,434	25,477	11	36.84
Linear-linear piecewise with a random knot	25,308	25,338	25,398	15	33.06
Latent Growth Curve Models With Time-Invariant Covariates (TICs; Time Structure: Age in Months)
Proposed model	31,582	31,714	31,977	66	33.00
Latent Growth Curve Models With TICs (Time Structure: Grade-in-School in Months)
Proposed model	31,580	31,912	31,975	66	34.28

Table 9 presents the estimates of parameters of interest. It is noticed that postknot development in mathematics skills slowed down substantially (on average 1.731 and 0.688 per month in the pre- and postknot stage, respectively). The knot was estimated at 102 months (i.e., 8.5 years old) on average. We noted that the initial mathematics performance (60-month old in our case) was positively associated with family income, private school, and parents’ highest education. Moreover, the variable family income also can explain individual differences in the development rate of mathematics IRT scores in the first stage. Further, boys improved more rapidly than girls in terms of mathematics ability in the first stage. Additionally, students coming from schools in more urban areas arrived at the change-point earlier.

Table 9.

Estimates of the Parameters of Mathematics Trajectories

Growth Factor	Intercept			First Slope
Parameter	Estimate	SE	p Value	Estimate	SE	p Value
Mean	26.855	0.589	<.0001*	1.731	0.019	<.0001*
Unexplained variance	106.321	9.822	<.0001*	0.083	0.010	<.0001*
Sex (0 = boy; 1 = girl)	0.377	0.591	.5235	−0.042	0.019	.0271*
Race (0 = White; 1 = Other)	−0.245	0.685	.7206	−0.001	0.022	.9637
School location	−0.585	0.634	.3562	−0.027	0.020	.1770
Income	1.743	0.838	.0375*	0.053	0.027	.0497*
School type (0 = public; 1 = private)	1.841	0.611	.0026*	−0.025	0.019	.1882
Parents’ highest education	2.863	0.787	.0003*	0.045	0.025	.0719
Growth Factor	Second Slope			Knot
Mean	0.688	0.017	<.0001*	101.741	0.523	<.0001*
Unexplained variance	0.014	0.008	.0801	42.477	8.308	<.0001*
Sex (0 = boy; 1 = girl)	0.021	0.017	.2167	−0.805	0.535	.1324
Race (0 = White; 1 = Other)	−0.043	0.019	.0236*	0.196	0.603	.7451
School location	−0.023	0.018	.2013	1.692	0.559	.0025*
Income	0.008	0.025	.7490	−0.529	0.799	.5079
School type (0 = public; 1 = private)	−0.028	0.018	.1198	0.097	0.554	.8610
Parents’ highest education	−0.026	0.023	.2583	0.337	0.739	.6484

* Statistical significance at .05 level.

5.2. Sensitivity Analysis

In this section, we use students’ grade-in-school as the time scale, where the time metric takes on discrete values (e.g., 0 is for the kindergarten fall semester and 0.5 in years or 6 in months is for the kindergarten spring semester). We list the estimated likelihood, AIC, BIC, and the residual variance of the linear-linear piecewise growth model in Table 8.⁸ Upon further examination, the mathematics ability slowed down between the sixth and the seventh wave (around the fall semester of Grade 3), matching the estimate in the main analysis: students in Grade 3 usually aged between 8- and 9-year old. The effect size of the estimates of the path coefficients to the knot increased slightly, but the direction of the estimates stayed the same.

6. Discussion

In this article, we extended an existing linear-linear piecewise growth model with an unknown random knot to the framework of individual measurement occasions to account for the heterogeneity of the measurement time. We explored the predictors to explain the between-individual differences in the within-individual change patterns. We also proposed using the multivariate delta method to derive the mean vector and variance-covariance matrix of growth factors as well as the path coefficients directly related to the underlying change patterns. More importantly, we examined the proposed method by extensive simulation studies to investigate whether the approximation introduced by the Taylor series expansion affects the model performance.

To demonstrate how well the introduced approximation works, we compared the full model (that requires the Taylor series approximation) with the reduced model and found that the full model performed better than the reduced one in terms of the performance measures under the conditions with varying knot locations. The simulation study also showed that the full model performed well under the conditions with a medium-level knot SD, which aligns with the essential condition of using the Taylor series expansion in our case. We also demonstrate the use of the proposed method on a subset with $n = 400$ from ECLS-K: 2011.

6.1. Practical Considerations

Based on the simulation study and the real-world data analysis, we provide a set of notes for empirical researchers. First, it has been shown that the bilinear latent growth model can be identified with at least five measurements with a prespecified knot at the halfway of study duration (i.e., there are two measures before and after the knot, respectively), although the information about model identification of linear-linear growth model with an unknown knot is scarce. Accordingly, we recommend visualizing the trajectories to examine the number of measurement occasions pre- and postknot. Second, for an empirical study, we still recommend examining the issue of improper solutions. As shown in the simulation study, greater than 80% of improper solutions were observed when we overspecified the model. Based on this result, an improper variance and/or correlation may simply suggest that the knot location is roughly similar across all individuals. We recommend employing the reduced model to estimate a fixed knot. Third, based on the simulation studies, the proposed model performed well in general. However, in cases with challenging conditions such as the large knot SD (i.e., 0.6), the small standardized difference between two slopes (i.e., $d_{z} = 1$ ), small sample size ( $n = 200$ ), and/or less precise measurements ( $θ_{ε} = 2$ ), the estimates of the path coefficients to the knot exhibited some bias greater than 10%.

It also suggests that, in practice, we need to interpret these knot-related estimates carefully. As shown in the simulation study, the unexplained knot variance was underestimated when the knot SD was 0.6 but could be overestimated when it was 0.3 when considering scaled measurement waves. In the Application section, the mathematics ability was evaluated every 6 months (the first six measurements) or every 12 months (the last three measurements). Given that the knot SD should be greater than 6.5 ( $\sqrt{42.477}$ , where 42.477 is the estimated unexplained knot variance) months, the true value of the unexplained knot variance might be larger than 42.477. Similarly, the path coefficients in the simulation study might be underestimated, suggesting that the covariate effects on the knot variance might be greater than the values (farther away from zero) presented in Table 9.

The proposed model is computationally complex, and thus, a set of proper starting values help improve the likelihood of convergence and lighten the computational burden. We first decide on an appropriate scale for each parameter. Summary statistics and visualization techniques are helpful in deciding the proper scales of growth factor means. In the Application section, we set the scales of the starting values of the growth factor variance and covariance to 1 and 0, respectively. The determination of covariates-related starting values is straightforward since we standardize those covariates: We set the means of covariates to 0 and the variance-covariance matrix of covariates to the observed values. Additionally, the scale of starting values of path coefficients was set as 0.5. For a complex model, we may need to further search starting values around the scales. One common approach for fine-tuning is the “grid search” method, where we create a list of starting values around the scale of each parameter, run the model with different combinations of the starting values across the parameter space, and see from which set we can obtain a converged solution. However, we do not suggest employing this method, given the size of the parameters. Instead, we recommend using the function mxTryHard() in the R package OpenMx that can help to obtain convergent solutions, as we did in the simulation study and the real-world data analysis. The function mxTryHard() allows for multiple attempts to fit a model until generating an acceptable solution. At each attempt, each starting value is multiplied by a random draw from a prespecified uniform distribution. With these random starts, OpenMx is often able to find a convergent solution in our experience.

Additionally, the proposed model can be employed for either theoretical or empirical consideration. Theoretically, as stated earlier, the proposed models allow for estimating the growth rate of each of two stages and the transition time from one stage to another. We recommend using either the full model or the reduced model if these growth coefficients are of research interest. When analyzing trajectories in an exploratory stage, where researchers only have a vague assumption about the underlying change patterns, the functional form can be decided by an empirical approach as we did in the Application section. To demonstrate the influence of the time scale on the estimates, we conducted a sensitivity analysis in the Application section. The fixed effects of growth factors obtained from the model with different time scales were almost the same, though the estimated effect size of path coefficients increased when we used the grade-in-school as the time scale. This difference is not surprising as students in the same grade-in-school could be of different ages. In practice, we recommend selecting the time scale by research questions instead of based on statistical criteria.

6.2. Methodological Considerations and Future Directions

This article focused on the growth model with a bilinear spline functional form because it is straightforward, intuitive, and useful. It can be generalized to a linear spline with multiple knots or a nonlinear spline (such as a linear-polynomial piecewise), where the corresponding (inverse-) transformation matrices can also be derived. Besides, it is feasible to extend the current work to address longitudinal studies with dropout under the assumption of missing at random thanks to the FIML technique. Additionally, we focused on TICs in the current work, which is helpful to identify baseline characteristics to explain the variability of growth factors. It can also be extended to investigate time-varying covariates to explain individual differences in the trajectories.

From the simulation studies, we observed that the improper solutions might also be generated under some challenging conditions such as the small standardized difference between two slopes (i.e., $d_{z} = 1$ ) and/or the large residual variance (i.e., $θ_{ε} = 2$ ). Accordingly, we propose to add constraints on the knot variance and its correlation with other growth factors and provide corresponding OpenMx code in the first author's GitHub website (https://github.com/Veronica0206/Dissertation_projects). The examination of the model’s performance with constraints through simulation studies is beyond the scope of the current study but can be examined in future work.

The (inverse-) transformation between growth factors in two parameter spaces is valuable for the reparameterized longitudinal model. Preacher and Hancock (2012, 2015) proposed to employ reparameterization to estimate parameters of the greatest research interest directly (e.g., a knot and its variance). One drawback of the reparameterized longitudinal model initially was that other reparameterized coefficients were no longer related to the developmental theory and therefore lacked meaningful, substantive interpretation, thus requiring a back transformation of the model parameters. In the linear-linear piecewise growth model with an unknown random knot, this transformation is at the individual level as each individual has a set of “personal” growth factors. We then proposed to use the multivariate delta method to simplify the calculation. The matrix form is especially convenient when we use the R package OpenMx that can carry out matrix multiplication automatically. The transformation functions and matrices for other reparameterized longitudinal models can also be proposed.

We construct the proposed model in the LGC modeling framework, which allows us to extend the model further. For example, we can extend the proposed model to examine a joint development of multiple repeated outcomes. There are at least two possible directions to realize this aim. One possible extension is a parallel process, and correlated growth model (McArdle, 1988), where each repeated outcome takes the linear-linear functional form with an unknown random knot. By utilizing this model, we can examine a joint development of reading ability and mathematics ability. The other possible extension is a second-order growth model (Hancock et al., 2001), where the factor-level change pattern is supposed to take a linear-linear functional form. To demonstrate our point, suppose we want to investigate the development of students’ overall academic performance measured by their reading, mathematics, and science ability. To construct the model, we consider the measurement model’s factors as first-order factors since they are immediately above the measured reading, mathematics, and science skills. Additionally, the growth factors (i.e., the intercept, two slopes, and the knot) can be viewed as second-order factors because they are two levels above the measured skills. In this model, the second-order growth factors are employed to account for the within-individual change and between-individual differences in the first-order factors.

6.3. Concluding Remarks

In this article, we propose transformation functions and matrices between coefficients in the original and reparameterized parameter spaces to derive interpretable coefficients related to developmental theory for an existing reparameterized longitudinal model in the LGC framework. The results of the simulation study and the real-world data analysis demonstrate the model’s valuable capabilities of estimating the knot and the path coefficients to the knot in the framework of individual measurement occasions and providing the interpretable coefficients in the original parameter space. Furthermore, as discussed above, the flexibility of the SEM framework and the ability to use latent variables in the LGC model allows for further extensions for the proposed model.

Supplemental Material

Supplemental Material, sj-docx-1-jeb-10.3102_10769986211052009 - Obtaining Interpretable Parameters From Reparameterized Longitudinal Models: Transformation Matrices Between Growth Factors in Two Parameter Spaces

Supplemental Material, sj-docx-1-jeb-10.3102_10769986211052009 for Obtaining Interpretable Parameters From Reparameterized Longitudinal Models: Transformation Matrices Between Growth Factors in Two Parameter Spaces by Jin Liu, Robert A. Perera, Le Kang, Roy T. Sabo and Robert M. Kirkpatrick in Journal of Educational and Behavioral Statistics

Supplemental Material

Supplemental Material, sj-pdf-1-jeb-10.3102_10769986211052009 - Obtaining Interpretable Parameters From Reparameterized Longitudinal Models: Transformation Matrices Between Growth Factors in Two Parameter Spaces

Supplemental Material, sj-pdf-1-jeb-10.3102_10769986211052009 for Obtaining Interpretable Parameters From Reparameterized Longitudinal Models: Transformation Matrices Between Growth Factors in Two Parameter Spaces by Jin Liu, Robert A. Perera, Le Kang, Roy T. Sabo and Robert M. Kirkpatrick in Journal of Educational and Behavioral Statistics

Footnotes

ORCID iD

Jin Liu

Notes

References

Bauer

D. J.

(2003). Estimating multilevel linear models as structural equation models. Journal of Educational and Behavioral Statistics, 28(2), 135–167. https://doi.org/10.3102/10769986028002135

Bauer

D. J.

Curran

P. J.

(2003). Distributional assumptions of growth mixture models: Implications for overextraction of latent trajectory classes. Psychological Methods, 8(3), 338–363. https://doi.org/10.1037/1082-989X.8.3.338

Biesanz

J. C.

Deeb-Sossa

Papadakis

A. A.

Bollen

K. A.

Curran

P. J.

(2004). The role of coding time in estimating and interpreting growth curve models. Psychological Methods, 9(1), 30–52. https://doi.org/10.1037/1082-989X.9.1.30

Boker

S. M.

Neale

M. C.

Maes

H. H.

Wilde

M. J.

Spiegel

Brick

T. R.

Estabrook

Bates

T. C.

Mehta

von Oertzen

Gore

R. J.

Hunter

M. D.

Hackett

D. C.

Karch

Brandmaier

A. M.

Pritikin

J. N.

Zahery

Kirkpatrick

R. M.

(2020). Openmx 2.17.2 user guide. https://vipbg.vcu.edu/vipbg/OpenMx2/docs/OpenMx/2.17.2/OpenMxUserGuide.pdf

Bollen

K. A.

Curran

P. J.

(2005). Latent curve models: A structural equation perspective. John Wiley & Sons, Inc.

Browne

M. W.

du Toit

S. H. C.

(1991). Models for learning data. In Collins

L. M.

Horn

J. L.

(Eds.), Best methods for the analysis of change: Recent advances, unanswered questions, future directions (pp. 47–68). American Psychological Association.

Chou

C.-P.

Yang

Pentz

M. A.

Hser

Y. I.

(2004). Piecewise growth curve modeling approach for longitudinal prevention study. Computational Statistics & Data Analysis, 46(2), 213–225. https://doi.org/10.1016/S0167-9473(03)00149-X

Cohen

(1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum Associates.

Cook

N. R.

Ware

J. H.

(1983). Design and analysis methods for longitudinal research. Annual Review of Public Health, 4(1), 1–23.

10.

Coulombe

Selig

J. P.

Delaney

H. D.

(2015). Ignoring individual differences in times of assessment in growth curve modeling. International Journal of Behavioral Development, 40(1), 76–86. https://doi.org/10.1177/0165025415577684

11.

Cudeck

du Toit

S. H. C.

(2003). Nonlinear multilevel models for repeated measures data. In Reise

S. P.

Duan

(Eds.), Multilevel modeling: Methodological advances, issues, and applications (pp. 1–24). Psychology Press.

12.

Cudeck

Harring

J. R.

(2010). Developing a random coefficient model for nonlinear repeated measures data. In Chow

S.-M.

Ferrer

Hsieh

(Eds.), The notre dame series on quantitative methodology. Statistical methods for modeling human dynamics: An interdisciplinary dialogue (pp. 289–318). Routledge/Taylor & Francis Group.

13.

Cudeck

Klebe

K. J.

(2002). Multiphase mixed-effects models for repeated measures data. Psychological Methods, 7(1), 41–63. https://doi.org/10.1037/1082-989X.7.1.41

14.

Curran

P. J.

(2003). Have multilevel models been structural equation models all along? Multivariate Behavioral Research, 38(4), 529–569. https://doi.org/10.1207/s15327906mbr3804_5

15.

Dominicus

Ripatti

Pedersen

N. L.

Palmgren

(2008). A random change point model for assessing variability in repeated measures of cognitive function. Statistics in Medicine, 27(27), 5786–5798. https://doi.org/10.1002/sim.3380

16.

Dumenci

Perera

R. A.

Keefe

F. J.

Ang

D. C.

Slover

Jensen

M. P.

Riddle

D. L.

(2019). Model-based pain and function outcome trajectory types for patients undergoing knee arthroplasty: A secondary analysis from a randomized clinical trial. Osteoarthritis and Cartilage, 27(6), 878–884. https://doi.org/10.1016/j.joca.2019.01.004

17.

Du Toit

S. H. C.

Cudeck

(2009). Estimation of the nonlinear random coefficient model when some random effects are separable. Psychometrika, 74, 65–82. https://doi.org/10.1007/s11336-008-9107-7

18.

Finkel

Reynolds

Mcardle

Gatz

Pedersen

N. L.

(2003). Latent growth curve analyses of accelerating decline in cognitive abilities in late adulthood. Developmental Psychology, 39, 535–550. https://doi.org/10.1037/0012-1649.39.3.535

19.

Flora

D. B.

(2008). Specifying piecewise latent trajectory models for longitudinal data. Structural Equation Modeling: A Multidisciplinary Journal, 15(3), 513–533. https://doi.org/10.1080/10705510802154349

20.

Grimm

K. J.

Ram

Estabrook

(2016). Growth modeling: Structural equation and multilevel modeling approaches. Guilford Press.

21.

Hancock

G. R.

Kuo

W. L.

Lawrence

F. R.

(2001). An illustration of second-order latent growth models. Structural Equation Modeling: A Multidisciplinary Journal, 8, 470–489.

22.

Harring

J. R.

Cudeck

du Toit

S. H. C.

(2006). Fitting partially nonlinear random coefficient models as sems. Multivariate Behavioral Research, 41(4), 579–596. https://doi.org/10.1207/s15327906mbr4104_7

23.

Hedeker

Mermelstein

R. J.

Flay

B. R.

(2006). Application of item response theory models for intensive longitudinal data. In Walls

T. A.

Schafer

J. L.

(Eds.), Models for intensive longitudinal data (pp. 84–108). Oxford University Press.

24.

Hunter

M. D.

(2018). State space modeling in an open source, modular, structural equation modeling environment. Structural Equation Modeling, 25(2), 307–324. https://doi.org/10.1080/10705511.2017.1369354

25.

Kohli

(2011). Estimating unknown knots in piecewise linear-linear latent growth mixture models [Doctoral dissertation, University of Maryland]. http://hdl.handle.net/1903/11973

26.

Kohli

Harring

J. R.

(2013). Modeling growth in latent variables using a piecewise function. Multivariate Behavioral Research, 48(3), 370–397. https://doi.org/10.1080/00273171.2013.778191

27.

Kohli

Harring

J. R.

Hancock

G. R.

(2013). Piecewise linear-linear latent growth mixture models with unknown knots. Educational and Psychological Measurement, 73(6), 935–955. https://doi.org/10.1177/0013164413496812

28.

Kohli

Harring

J. R.

Zopluoglu

(2016). A finite mixture of nonlinear random coefficient models for continuous repeated measures data. Psychometrika, 81(3), 851–880. https://doi.org/10.1007/s11336-015-9462-0

29.

Kohli

Hughes

Wang

Zopluoglu

Davison

M. L.

(2015). Fitting a linear-linear piecewise growth mixture model with unknown knots: A comparison of two common approaches to inference. Psychological Methods, 20(2), 259–275. https://doi.org/10.1037/met0000034

30.

Kohli

Sullivan

A. L.

Sadeh

Zopluoglu

(2015). Longitudinal mathematics development of students with learning disabilities and students without disabilities: A comparison of linear, quadratic, and piecewise linear mixed effects models. Journal of School Psychology, 53(2), 105–120. https://doi.org/10.1016/j.jsp.2014.12.002

31.

Kwok

Luo

West

S. G.

(2010). Using modification indexes to detect turning points in longitudinal data: A monte carlo study. Structural Equation Modeling: A Multidisciplinary Journal, 17(2), 216–240. https://doi.org/10.1080/10705511003659359

32.

Lê

Norman

Tourangeau

Brick

J. M.

Mulligan

(2011). Early childhood longitudinal study: Kindergarten class of 2010-2011—Sample design issues. JSM Proceedings 2011, 1629–1639. http://www.asasrms.org/Proceedings/y2011/Files/301090_66141.pdf

33.

Lehmann

E. L.

Casella

(1998). Theory of point estimation (2nd ed.). Springer-Verlag, Inc.

34.

Duncan

T. E.

Duncan

S. C.

Hops

(2001). Piecewise growth mixture modeling of adolescent alcohol use data. Structural Equation Modeling: A Multidisciplinary Journal, 8(2), 175–204. https://doi.org/10.1207/S15328007SEM0802_2

35.

Liu

(2019). Estimating knots in bilinear spline growth models with time-invariant covariates in the framework of individual measurement occasions [Doctoral dissertation, Virginia Commonwealth University]. https://doi.org/10.25772/9WDR-9R85

36.

Lock

E. F.

Kohli

Bose

(2018). Detecting multiple random changepoints in Bayesian piecewise growth mixture models. Psychometrika, 83(3), 733–750. https://doi.org/10.1007/s11336-017-9594-5

37.

Marcoulides

K. M.

(2018). Automated latent growth curve model fitting: A segmentation and knot selection approach. Structural Equation Modeling: A Multidisciplinary Journal, 25(5), 687–699. https://doi.org/10.1080/10705511.2018.1424548

38.

McArdle

J. J.

(1988). Dynamic but structural equation modeling of repeated measures data. In Nesselroade

Cattell

(Eds.), Handbook of multivariate experimental psychology (pp. 561–614). Springer.

39.

McArdle

J. J.

Wang

(2008). Modeling age-based turning points in longitudinal life-span growth curves of cognition. In Cohen

(Ed.), Applied data analytic techniques for turning points research (pp. 1–24). Routledge/Taylor & Francis Group.

40.

Mehta

P. D.

Neale

M. C.

(2005). People are variables too: Multilevel structural equations modeling. Psychological Methods, 10(3), 259–284. https://doi.org/10.1037/1082-989X.10.3.259

41.

Mehta

P. D.

West

S. G.

(2000). Putting the individual back into individual growth curves. Psychological Methods, 5(1), 23–43. https://doi.org/10.1037/1082-989x.5.1.23

42.

Morris

T. P.

White

I. R.

Crowther

M. J.

(2019). Using simulation studies to evaluate statistical methods. Statistics in Medicine, 38(11), 2074–2102. https://doi.org/10.1002/sim.8086

43.

Muniz Terrera

van den Hout

Matthews

F. E.

(2011). Random change point models: Investigating cognitive decline in the presence of missing data. Journal of Applied Statistics, 38(4), 705–716. https://doi.org/10.1080/02664760903563668

44.

Neale

M. C.

Hunter

M. D.

Pritikin

J. N.

Zahery

Brick

T. R.

Kirkpatrick

R. M.

Estabrook

Bates

T. C.

Maes

H. H.

Boker

S. M.

(2016). OpenMx 2.0: Extended structural equation and statistical modeling. Psychometrika, 81(2), 535–549. https://doi.org/10.1007/s11336-014-9435-8

45.

Preacher

K. J.

Hancock

G. R.

(2012). On interpretable reparameterizations of linear and nonlinear latent growth curve models. In Harring

J. R.

Hancock

G. R.

(Eds.), Cilvr series on latent variable methodology. Advances in longitudinal methods in the social and behavioral sciences (pp. 25–58). IAP Information Age Publishing.

46.

Preacher

K. J.

Hancock

G. R.

(2015). Meaningful aspects of change as novel random coefficients: A general method for reparameterizing longitudinal models. Psychological Methods, 20(1), 84–101. https://doi.org/10.1037/met0000028

47.

Pritikin

J. N.

Hunter

M. D.

Boker

S. M.

(2015). Modular open-source software for item factor analysis. Educational and Psychological Measurement, 75(3), 458–474. https://doi.org/10.1177/0013164414554615

48.

Riddle

D. L.

Perera

R. A.

Jiranek

W. A.

Dumenci

(2015). Using surgical appropriateness criteria to examine outcomes of total knee arthroplasty in a united states sample. Arthritis Care & Research, 67(3), 349–357. https://doi.org/10.1002/acr.22428

49.

Seber

G. A. F.

Wild

C. J.

(2003). Nonlinear regression. John Wiley & Sons, Inc.

50.

Sterba

S. K.

(2014). Fitting nonlinear latent growth curve models with individually varying time points. Structural Equation Modeling: A Multidisciplinary Journal, 21(4), 630–647. https://doi.org/10.1080/10705511.2014.919828

51.

Sullivan

A. L.

Kohli

Farnsworth

E. M.

Sadeh

Jones

(2017). Longitudinal models of reading achievement of students with learning disabilities and without disabilities. School Psychology Quarterly: The Official Journal of the Division of School Psychology, American Psychological Association, 32(3), 336–349. https://doi.org/10.1037/spq0000170

52.

Tishler

Zang

(1981). A maximum likelihood method for piecewise regression models with a continuous dependent variable. Journal of the Royal Statistical Society. Series C. Applied Statistics, 30, 116–124. https://doi.org/10.2307/2287599

53.

Venables

W. N.

Ripley

B. D.

(2002). Modern applied statistics with s (4th ed.). Springer.

54.

Wang

McArdle

J. J.

(2008). A simulation study comparison of Bayesian estimation with conventional methods for estimating unknown change points. Structural Equation Modeling: A Multidisciplinary Journal, 15(1), 52–74. https://doi.org/10.1080/10705510701758265

55.

Zopluoglu

Harring

J. R.

Kohli

(2014). Fitpmm: An r routine to fit finite mixtures of piecewise mixed-effects models with unknown random knots. Applied Psychological Measurement, 38, 583–584. https://doi.org/10.1177/0146621614540482

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.14 MB

0.12 MB