{'fraction_non_alphanumeric': 0.042040807939671646, 'fraction_numerical': 0.015534140177122074, 'mean_word_length': 4.874939642684693, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Obtaining interpretable parameters from reparameterizing longitudinal models: transformation matrices between growth factors in two parameter spaces

The linear spline growth model (LSGM), which approximates complex patterns using at least two linear segments, is a popular tool for examining nonlinear change patterns. Among such models, the linear-linear piecewise change pattern is the most straightforward one. An earlier study has proved that other than the intercept and slopes, the knot (or change-point), at which two linear segments join together, can be estimated as a growth factor in a reparameterized longitudinal model in the latent growth curve modeling framework. However, the reparameterized coefficients were no longer directly related to the underlying developmental process and therefore lacked meaningful, substantive interpretation, although they were simple functions of the original parameters. This study proposes transformation matrices between parameters in the original and reparameterized models so that the interpretable coefficients directly related to the underlying change pattern can be derived from reparameterized ones. Additionally, the study extends the existing linear-linear piecewise model to allow for individual measurement occasions, and investigates predictors for the individual-differences in change patterns. We present the proposed methods with simulation studies and a real-world data analysis. Our simulation studies demonstrate that the proposed method can generally provide an unbiased and consistent estimation of model parameters of interest and confidence intervals with satisfactory coverage probabilities. An empirical example using longitudinal mathematics achievement scores shows that the model can estimate the growth factor coefficients and path coefficients directly related to the underlying developmental process, thereby providing meaningful interpretation. For easier implementation, we also provide the corresponding code for the proposed models.


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 pattern is a piecewise linear latent growth model (Chou et al., 2004;Harring et al., 2006;Cudeck and Harring, 2010;Kohli, 2011;Kohli et al., 2013;Kohli and Harring, 2013;Sterba, 2014), also known as the linear spline growth model (Grimm et al., 2016a). 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 a broad application of linear spline trajectory. For example, it has been employed to assess post-surgical rehabilitation (Riddle et al., 2015;Dumenci et al., 2019), the learning process of a specific task (Cudeck and 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, Flora (2008); Sterba (2014); Riddle et al. (2015); Dumenci et al. (2019), the knot is unknown in most cases. Marcoulides (2018) proposed to fit a pool of candidate models with a knot at different locations and select the 'best' one using the BIC criterion. By building a multi-phase 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 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) and Bayesian Dominicus et al. (2008); McArdle and Wang (2008); Wang and McArdle (2008) The most intuitive linear spline function is a bilinear spline growth model (or linear-linear piecewise model). This functional form can identify a process that has two stages with different rates of change (Riddle et al., 2015;Dumenci et al., 2019), but it can also help approximate more complex nonlinear change patterns (Kohli et al., 2015a,b;Sullivan et al., 2017). 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 re-expressing growth factors as linear combinations of their original forms. Kohli (2011), Kohli et al. (2013) and Grimm et al. (2016a) derived the growth coefficients that are directly related to the underlying change patterns from the estimated coefficients using Mplus. 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 latent growth curve (LGC) modeling framework. More importantly, they considered incorporating obesity status for predicting varying knot locations. By re-examining a data set containing longitudinal records of plasma phosphate, they demonstrated the linear-linear functional form with an unknown random knot is useful. However, one drawback of this model lies in that it is fitted in a reparameterized framework. The reparameterized coefficients were no longer directly related to the underlying developmental process and therefore lacked meaningful, substantive interpretation.
For these reasons, we propose inverse-transformation matrices to derive the growth factor coefficients in the original frame so that they are directly related to the underlying change pattern and interpretable. 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 the knot location. 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 and Ware, 1983;Mehta and West, 2000;Finkel et al., 2003). This can occur when time is measured precisely (for example, exact date and time of measurement), or if responses are self-initiated or randomly assigned. For instance, in an adolescent smoking study, participants were asked to complete a questionnaire on pocket computers immediately after smoking (Hedeker et al., 2006). 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 and West, 2000;Mehta and Neale, 2005). 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 time-invariant covariates 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 SEM 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 non-convergence 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 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 − 11 (ECLS-K: 2011), we give a set of feasible recommendations for real-world practice. Finally, a broad discussion is presented regarding the model's limitations as well as future directions.

Model Specification
This section briefly describes a linear-linear growth curve model with an unknown random knot and time-invariant covariates 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 LGC, 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 i th individual at the j th time point in the framework of individual measurement occasions as where y ij and t ij are the measurement and measurement time of the i th individual at time j. In Equation (1), η 0i , η 1i , η 2i 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 .

=========================
Insert Figure 1 about here ========================= Note that we cannot fit the linear-linear piecewise change pattern specified in Equation (1) in the LGC framework for two reasons. First, the outcome variable y ij does not have the same model expression pre-and post-knot. Second, the model in Equation (1) specifies a nonlinear relationship between the y ij and the growth factor (or random coefficient in the mixed-effects model framework) γ i and then cannot be estimated in the structural equation modeling (SEM) framework directly (Grimm et al., 2016b). Following Tishler and Zang (1981); Seber and Wild (2003); Grimm et al. (2016a);Liu (2019), the expression of repeated outcome in Equation (1) before and after the knot can be unified by reparameterization (see Appendix A.1 for details of reparameterization). We then write the repeated outcome as a linear combination of all four growth factors by the Taylor series expansion (Browne and du Toit, 1991;Grimm et al., 2016a;Liu, 2019) (see Appendix A.2 for details of Taylor series expansion). Accordingly, the model specified in Equation (1) can be written as a standard LGC model with reparameterized growth factors where y i is a J × 1 vector of the repeated outcomes for the i th person (in which J is the number of measurements), η i is a 4 × 1 vector of individual-level reparameterized growth factors and Λ i , which is a function of time points and the unknown knot, is a J × 4 matrix of factor loadings 2 . Additionally, i is a J × 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 corresponding factor loadings are and respectively, where µ γ is the knot mean and δ i is the deviation from the knot mean of the i th individual. The reparameterized growth factors of the i th individual can be further regressed on TICs, where α is a 4 × 1 vector of reparameterized growth factor intercepts (which is the mean vector of reparameterized growth factors if the TICs are centered), B is a 4 × c matrix of regression coefficients (in which c is the number of TICs) from TICs to reparameterized growth factors, X i , which may either be continuous or be binary, is a c × 1 vector of covariates of the i th individual, and ζ i is a 4 × 1 vector of deviations of the i th subject from the reparameterized growth factor means. To simplify estimation, we assume that the (reparameterized) growth factors are normally distributed conditional on individual-level covariates, so ζ i ∼ MVN(0, Ψ η ), where Ψ η is a 4 × 4 unexplained variance-covariance matrix of reparameterized growth factors. We also assume that individual residuals, i in Equation (2), are identical and independent normal distribution over time, that is, i ∼ MVN(0, θ I), where I is a J × J identity matrix. Note that the intercepts and coefficients of reparameterized growth factors (α and B ) as well as the unexplained growth factor variances Ψ η are no longer directly related to the underlying developmental process and therefore lack meaningful, substantive interpretation. We demonstrate how to derive the coefficients in the original parameter-space so that all coefficients are interpretable.

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 initial values can improve the likelihood of convergence and then accelerate the computational process. Generally, descriptive statistics and visualization are tools for researchers to decide suitable initial 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 initial values. More importantly, all reparameterized coefficients are not directly related to the underlying change patterns and then need to be transformed back to be interpretable after estimation. So the inverse-transformation that help derive the coefficients related to the underlying trajectory directly from the estimates are also needed.
Suppose f : R 4 → R 4 is a function, which takes a point η i ∈ R 4 as input and produces the vector f (η i ) ∈ R 4 (i.e., η i ∈ R 4 ) as output. By the multivariate delta method (Lehmann and Casella, 1998), where µ η and Ψ η are the mean vector and variance-covariance matrix of original growth factors, respectively, and f is defined as Similarly, suppose h : R 4 → R 4 is a function, which takes a point η i ∈ R 4 as input and produces the vector h(η i ) ∈ R 4 (i.e., η i ∈ R 4 ) as output. By the multivariate delta method, where µ η and Ψ η are the mean vector and variance-covariance matrix of reparameterized growth factors, respectively, and h is defined as Based on Equations (6) and (7), we can make the transformation between the growth factor means of two parameter spaces by µ η ≈ f (µ η ) and µ η ≈ h(µ η ), respectively. We can also define the transformation matrix ∇ f (µ η ) and ∇ h (µ η ) between the variance-covariance matrix of two parameter spaces as 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 α = µ η ≈ f (µ η ) and α = µ η ≈ 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 ∇ f (µ η ) and ∇ h (µ η ), respectively. The detailed derivation is shown in Appendix A.3.
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., 2018), we only need to provide h(µ η ) and ∇ h (µ η ), and the package can derive the coefficients in the original setting and generate the point estimates with their standard errors of the growth factor coefficients in the original setting 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 NEW 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 since Mplus does not allow for matrix algebra. All these expressions are in Appendix A.4.

Model Estimation
For the i th 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 and where µ X and Φ are the mean vector (c × 1) and the variance-covariance matrix (c × 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 re-expressed path coefficients as well as the means and the variancecovariance structure of covariates. Then we can obtain the parameters in the original setting, as shown in Section 2.2. Θ 1 and Θ 1 shown below list the parameters in the original and reparameterized frames, respectively.
We then use full information maximum likelihood (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 and respectively, where C is a constant, n is the number of individuals, µ i and Σ i are the mean vector and the variancecovariance 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 (Neale et al., 2016;Pritikin et al., 2015;Hunter, 2018;Boker et al., 2018), with which we are able to fit the model and obtain interpretable coefficients efficiently, as shown in Section 2.2.

Reduced Model
With an assumption that the change-point is roughly similar for each individual, 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 where η i is a 3 × 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 × 3 matrix of factor loadings. We can further express the reparameterized growth factors and corresponding factor loadings as The growth factor intercepts α , path coefficients B and the deviation ζ i of the i th individual from the reparameterized growth factor means also reduce to be a 3 × 1 vector, 3 × c matrix and 3 × 1 vector, respectively. The transformation and inverse-transformation functions 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 ∇ f (µ η ) and ∇ h (µ η ), since only three growth factors need to be reparameterized.
For the i th individual, the expected mean vector and the variance-covariance matrix of the repeated outcomes of this model are respectively. For this reduced model, Θ 2 and Θ 2 are defined as and they list the parameters in the original and reparameterized setting, 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 the FIML technique to estimate the parameters. We provide the OpenMx syntax for the full model and its reduced version as well as a demonstration in the online appendix (https://github.com/ Veronica0206/Dissertation_projects). For the researchers who are interested in using Mplus, we also provide Mplus 8 syntax for both models in the online appendix.

Model Evaluation
The proposed method is evaluated using a Monte Carlo simulation study with two goals. The first goal is to evaluate how the approximation based on the Taylor series expansion affects 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. We provide the definitions and estimates of these four performance metrics in Table 1. The second goal is to examine how the reduced model performed sufficiently well as compared to the full model.

Design of Simulation Study
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 time-invariant covariates. 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. As shown in Table 2, accordingly, 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 a tenth scale, guided by previous simulation studies for growth models (Bauer and Curran, 2003;Kohli, 2011;Kohli et al., 2015b). Further, the growth factors were set to be positively correlated to a moderate degree (ρ = 0.3).

=========================
Insert Table 2 about here ========================= We list all conditions that we considered in the simulation design in Table 2. 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: 6 and 10. We selected 6 as the minimum number of repeated outcomes to make the full model fully identified 4 . For the conditions with 6 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 Appendix A.2, one necessary condition to apply this approximation is that the knot variance approaches to 0. We then investigated 3 levels of magnitude of between-person differences in the knot to evaluate how the introduced approximation affects the model performance and to compare the full model versus the reduced model. We set the knot standard deviation 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 standard deviation is relative to the scale of measurement occasions.
Additionally, we examined how the shape of change patterns, which was quantified by the standardized difference between two slopes, affects the full model. It is noted that we considered both positive and negative differences (η 1i < η 2i and η 1i > η 2i , respectively) to assess model performance in the two cases as shown in Figure A.1. For each case, we considered 3 levels of standardized difference between two slopes to evaluate the models. We also considered 2 levels of effects of exogenous variables on endogenous growth factors 5 , 2 levels of sample size, and 2 levels of measurement precision as shown in Table 2. 3 Monte Carlo SE(Bias) = V ar(θ)/S (Morris et al., 2019). 4 Although the information of model identification of linear-linear growth model with an unknown knot is scarce, it has been proved that the bilinear latent growth model can be identified with at least five waves with a specified knot at the halfway of study duration (Bollen and Curran, 2005). 5 We set up TICs to explain moderate (13%) and substantial (26%) variances of such factors (Cohen, 1988).

Data Generation and Simulation Step
For each condition which is listed in Table 2, the general steps of the simulation study for the models were outlined as follows: 1. Generated data for growth factors and TICs simultaneously using the R package MASS (Venables and Ripley, 2002) (details are provided in Appendix A.5), 2. Generated the time structure with J scaled and equally-spaced waves t j and obtained individual measurement occasions: t ij ∼ U (t j − ∆, t j + ∆) by allowing disturbances around each wave, 3. Calculated factor loadings, which are functions of individual measurement occasions and the knot, for each individual, 4. Calculated values of the repeated measurements based on growth factors, factor loadings, and residual variances.
5. Implemented the full model and the reduced model on simulated data, estimated the parameters, constructed corresponding 95% Wald CIs, 6. Repeated the above steps until after obtaining 1, 000 convergent solutions.

Model Convergence and Proper Solution
In this section, we first examined the convergence rate and component fit measures of each condition. The convergence 6 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%). Out 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 non-convergence rate across all conditions was 46/1046, 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 standard deviation (i.e., sd(γ) = 0.6), and the smallest positive standardized difference between two slopes (i.e., d z = 1).
We also conducted diagnostics to investigate improper solutions, such as estimates of growth factor variances less than 0, and/or estimates of correlations between growth factors beyond [−1, 1]. Table 3

=========================
Insert Table 3 about here ========================= 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 standard deviation. 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.

Performance Measures
In this section, we present simulation results in terms of performance measures, which include the relative bias, empirical SE, relative RMSE and empirical coverage probability for a nominal 95% confidence interval for each parameter of interest. In general, the full model estimates coefficients unbiasedly, precisely, with satisfactory confidence interval coverage. 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 the relative bias/empirical SE of each parameter across 1, 000 replications under each condition and then summarized the relative biases/empirical SEs of each parameter over all aforementioned conditions 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 0.12, and that of unexplained variances of the intercept and two slopes was under 0.07. From Table 4, the full model may produce biased estimates for the 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.

=========================
Insert Table 5 about here ========================= 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 knot variance and the path coefficients to the knot. From the figure, we observed how the conditions we set up in the simulation design affected these estimates. First, under the conditions where the knot standard deviation 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 standard deviation while could be overestimated under the conditions with 0.3 knot standard deviation, 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.

=========================
Insert Figure 2 about here ========================= From Table 5, estimates from the full model and its reduced version were precise: the magnitude of empirical standard errors of slope or knot parameters were around 0.10, although those values of intercept parameters were relatively large: the empirical standard error 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.

=========================
Insert Table 6 about here ========================= Table 7 presents the median (range) of the coverage probability (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 standard deviation.

=========================
Insert Table 7 about here ========================= 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 coverage probability. 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 standard deviation (i.e., sd(γ) = 0.3 for the scaled measurement occasions) though was less satisfactory when the knot standard deviation 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, 5, 6 and 7 and found that shorter study duration (i.e., J = 6) and left-(right-) shifted knot locations only affected the performance measures slightly.

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. The first goal is to compare the fitness of bilinear spline functional form to three common parametric change patterns: linear, quadratic and Jenss-Bayley growth curve. The second goal is to show how to select the model between the full model and the reduced version. We randomly selected 400 students from the Early Childhood Longitudinal Study Kindergarten Cohort: 2010-11 (ECLS-K: 2011) with complete records of repeated mathematics IRT scaled scores, demographic information (sex, race, and age at each wave), school information (baseline school location and baseline school type), and social-economic status (baseline family income and the highest education level between parents) 7 . . Noting that two possible time scales are available in the data set: students' age and their grade-in-school, we decided to 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. 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 Hispanic, 8% were Asian, and 5% were others. For this analysis, 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 school, respectively. We treated the variables school location (ranged between 1 and 4), highest parents' education (ranged between 0 and 8) and family income (ranged between 1 and 18) as continuous variables, and corresponding mean (SD) was 2.06 (1.00), 5.18 (2.10) and 11.56 (5.42), respectively.

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, i.e., 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.

=========================
Insert Figure 3 about here ========================= Table 8 provides the obtained estimated likelihood, information criteria (including AIC and BIC), residuals of each latent growth curve 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). Its model fit information is also provided in Table 8.

=========================
Insert Table 8 about here ========================= Table 9 presents the estimates of parameters of interest. It is noticed that post-knot development in mathematics skills slowed down substantially (on average 1.731 and 0.688 per month in the pre-and post-knot stage, respectively). The knot was estimated at 102 month on average (8.5-year old). 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 family income also can explain the individual difference 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.

Sensitivity Analysis
In this section, we use students' grade-in-school as the time scale, where the time metric takes on discrete values (for example, 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 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 8to 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.

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 standard deviation, which aligns with the essential condition using the Taylor series expansion. Based on the result of simulation studies, we recommend using the bilinear spline growth model with an unknown random knot if the assumption of the same knot location is indefensible. In the scenario where a priori knowledge regarding the knot variability is scarce, we recommend fitting the full model and the reduced model and then select the 'best' one as we did in the Application section. Additionally, based on the simulation studies, the reduced model, which only estimates the knot itself, is a good backup if the full model fails to converge or generates improper solutions or the knot variance is out of research interest, although sometimes at a little cost of a small increase in bias.
The (inverse-) transformation between growth factors in two parameter-spaces are valuable for the reparameterized longitudinal model. One drawback of the reparameterized longitudinal model initially was that the reparameterized coefficients were no longer related to the underlying developmental process 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 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 demonstrated the use of the proposed method on a subset with n = 400 from ECLS-K: 2011 and provided a set of recommendations for possible issues that researchers may face in a real-world data analysis. First of all, we recommend selecting a proper functional form to depict trajectories before identifying the possible time-invariant covariates to explain the between-individual difference in the within-individual change, though this selection is not straightforward. Many factors, including but not limited to the research question, the fit between the model implied and the observed trajectory, and the information criteria, drive this selection process. It is worth considering either the full or the reduced model to estimate the knot if it is part of the research interest. If the interest lies in the fitting of nonlinear trajectories, it may be appropriate to fit the change patterns with several functional forms and then select the 'best' one. It is noted that the model fit information, sometimes, fails to lead unequivocal selection as in our example. Although either the full model or the quadratic functional forms can be viewed as an 'acceptable' model, the full model is considered better if it is important to capture children's mathematics achievement in the early stage of the study. The selection between the full and the reduced model can only rely on the information criteria such as the BIC as they usually fit the trajectories equally well, as shown in Figure 3. To demonstrate the influence of the time scale has 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.
To summarize, in this article, we focused on the growth model with a bilinear spline functional form because it is the most 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 time-invariant covariates 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 the individual-difference in the trajectories.

Appendix A Formula Derivation
A.1 The Reparameterizing Procedure for Growth Factors Tishler and Zang (1981) and Seber and Wild (2003) showed that the bilinear spline regression model can be written as either the minimum or maximum response value of two trajectories. By extending such expressions to the LGC framework, two forms of bilinear spline for the i th individual are shown in Figure A.1. In the left panel (η 1i > η 2i ), the measurement y ij should always be the minimum value of two lines; that is, y ij = min (η 0i + η 1i t ij , η 02i + η 2i t ij ). To unify the expression of measurements pre-and post-knot, we have the following equation where η 0i , η 1i , η 2i and γ i are the measurement at the knot, mean of two slopes, the half-difference between two slopes and knot for the i th individual. Through straightforward algebra, the measurement y ij of the bilinear spline in the right panel, in which the measurement y ij should always be the maximum value of two lines, has the identical final form as Equation (A.1).

A.2 Taylor Series Expansion
For the i th individual, suppose we define a function f (γ i ) and its first derivative with respect to γ i 9 , shown below respectively. Then the Taylor series expansion of f (γ i ) can be expressed as from which we then have the reparameterized growth factors and the corresponding factor loadings for the i th individual.

A.3 Details of the transformations between intercept coefficients and path coefficients
A.4 Expression of each cell of the re-reparameterized mean vector and variance-covariance matrix

A.5 Generate Exogenous Variables and Growth Factors Simultaneously
To generate growth factors and exogenous variables simultaneously as a multivariate normal distribution, we need to specify the mean vector and the variance-covariance matrix of the distribution. In our setting, the mean vector can be represented in equation where α and (0 0) are the mean vector of the growth factors and that of the TICs, respectively. Moreover, through straightforward linear algebra, the underlying variance-covariance matrix of the distribution is where B is the matrix of path coefficients, Ψ η and Φ are the variance-covariance matrices of the growth factors and the TICs, respectively. Accordingly, we can generate the growth factors and TICs simultaneously, which follows a multivariate normal distribution N 6 (µ, Σ).    Time (t j ) 6 scaled and equally spaced (j = 0 · · · , J − 1, J = 6) 10 scaled and equally spaced (j = 0, · · · , J − 1, J = 10) (Mean of) the knot µ γ at t = 2.5 for J = 6 µ γ at t = 3.5 or 4.5 or 5.5 for J = 10 Individual t ij t ij ∼ U (t j − ∆, t j + ∆)(j = 0, · · · , J − 1; ∆ = 0.25) Sample Size n = 200 n = 500 Full Factorial Design Variables Conditions