A note on estimating single-class piecewise mixed-effects models with unknown change points

Abstract

Piecewise mixed-effects models are useful for analyzing longitudinal educational and psychological data sets to model segmented change over time. These models offer an attractive alternative to commonly used quadratic and higher-order polynomial models because the coefficients obtained from fitting the model have meaningful substantive interpretation. The current study thus focuses on the estimation of piecewise mixed-effects model with unknown random change points using maximum likelihood (ML) as described in Du Toit and Cudeck (2009). Previous simulation work (Wang & McArdle, 2008) showed that Bayesian estimation produced reliable parameter estimates for the piecewise model in comparison to frequentist procedures (i.e., first-order Taylor expansion and the adaptive Gaussian quadrature) across all simulation conditions. In the current article a small Monte Carlo simulation study was conducted to assess the performance of the ML approach, a frequentist procedure, and the Bayesian approach for fitting linear–linear piecewise mixed-effects model. The obtained findings show that ML estimation approach produces reliable and accurate estimates under the conditions of small residual variance of the observed variables, and that the size of the residual variance had the most impact on the quality of model parameter estimates. Second, neither ML nor Bayesian estimation procedures performed well under all manipulated conditions with respect to the accuracy and precision of the estimated model parameters.

Keywords

Longitudinal data set piecewise function change points maximum likelihood

Introduction

Longitudinal data analytic methods have significantly advanced our understanding of various cognitive and behavioral developmental processes (e.g., change in aptitude or onset of alcohol abuse) across a particular span of time. The overall scientific objective is to effectively describe patterns of change for each individual as well as for a population as a whole. Measuring change over time necessarily requires a longitudinal perspective where repeated measurements are gathered for individual subjects. Such an understanding is especially valuable as it affords researchers and practitioners the opportunity to individualize treatment, and achieve more equity with respect to desirable outcomes. To be concrete, for example, developmental psychologists, educators, policy makers, and researchers can enrich their understanding of the development of students’ diverse learning trajectories, their understanding of whether the change process in cognitive and/or behavioral development of students can be predicted from covariates (e.g., from students’ sociodemographic characteristics, teacher characteristics, and the classroom environment), and their understanding of whether the teaching, generally, as well as specific interventions have desirable effects on students’ development over time (e.g., Grimm, 2008; Jordan, Hanich, & Kaplan, 2003; Miles & Miles, 1992; Shin, Davison, Long, Chan, & Heistad, 2013).

It is often seen in educational and psychological longitudinal data sets that the overall growth or development pattern over time in observed repeated measurements consists of distinct segments (phases) of growth (e.g., Kreisman, 2003; Paris, 2005; Silverman, Speece, Harring, & Ritchey, 2012). For example, in studies of interventions (e.g., evaluating the effectiveness of Head Start programs) the children usually do not improve much in the months before the intervention and that after the intervention they show rapid improvement; thus, the rate of growth is different pre- and post- intervention (Kreisman, 2003). Thus, the overall growth trajectory in this scenario incorporates two distinct segments of change, one segment that characterizes a slower acquisition trend that gives way to a later rapid acquisition trend. Because the rate of change is different in different segments, piecewise mixed-effects models (PMEs) offer a flexible framework to analyze this kind of segmented data set. These models allow the specification of separate growth profiles (functional forms) corresponding to multiple developmental stages of the overall change process from which repeated observations are made, such that each segment does not have to follow the same functional form (Chou, Yang, Pentz, & Hser, 2004).

One of the most interesting features of a PME model is the change point (knot), the value of the predictor—time at which the response function switches from one segment to another. For instance, in research related to the development of reading-related skills among elementary school students researchers have found that the ability of students to accurately and automatically decode fluently increases in the beginning of the second grade, but then it slows down somewhere in the middle of the third grade (Silverman et al., 2012). Hence, in this case the change point is located around the middle of the third grade. The location of change point can be known a priori or can be estimated. Additionally, the change point can be treated as a fixed parameter, indicating that each individual’s change point is assumed to be the same, or random, in which case the change point is an estimable subject-specific parameter. A PME model with unknown random change points is a type of nonlinear random-effects model (Cudeck & Klebe, 2002; Du Toit & Cudeck, 2009; Wang & McArdle, 2008). The estimation of the PME model with unknown random change points is computationally very challenging because of the nonlinear random parameter (i.e., change point) where we estimate both its mean and its variance (Kohli, Harring, & Zopluoglu, 2016). The two common approaches for estimation of longitudinal models are: maximum likelihood (ML) estimation and Markov chain Monte Carlo (MCMC) for Bayesian estimation.

The primary objective of the current study is to extend the prior simulation-based literature related to the estimation of PME models with unknown random change points. Researchers have shown that Bayesian estimation produces more reliable parameter estimates for PME models with unknown random change points in comparison to frequentist procedures of first-order Taylor expansion and adaptive Gaussian quadrature across all the simulation conditions (Wang & McArdle, 2008). In our study we intend to employ a more recent version of modified ML estimation approach proposed by Du Toit and Cudeck (2009) which seems to provide better estimation of nonlinear mixed-effects models, including PME models. Thus, in this study we aim to demonstrate the performance of the ML estimation as described by Du Toit and Cudeck (2009) with respect to the accuracy and precision of estimated PME model parameters via a small Monte Carlo simulation study. Moreover, to provide a more complete and comprehensive discussion, we also present Bayesian estimation results using MCMC as implemented in the rjags package for R (Plummer, 2016). A secondary objective of this article is to promote the implementation of PME models with random change points. We provide an (online) annotated Supplemental Appendix B that describes the R code for the ML model estimation used in this study. The remainder of this paper is organized as follows. Earlier simulation studies on the estimation of PME models with unknown random change points are briefly reviewed and discussed. In a subsequent method section, the PME model equation and the design of the simulation study are described. In the final section, the results from the simulation study are presented.

Literature review

To date there are few prior studies that have used Monte Carlo simulation techniques to study the performance of different estimation techniques for estimating PME model parameters with unknown random change points. The earliest simulation study on this topic was by Morrell, Pearson, Carter, and Brant (1995). The authors demonstrated the development and the utility of the PME model with random change points for examining change in prostate-specific antigen (PSA) levels among men before their prostate cancers were detected clinically. The estimation algorithm employed by them was a combination of least squares estimators for nonlinear mixed-effects models and ML estimators for a linear mixed-effects model. Based on the results of the simulation study, they reported that the estimates of fixed-effects parameters were close to the true (population) values. The asymptotic standard errors of the fixed-effects parameter estimates, however, were underestimated resulting in shorter confidence intervals and hence small coverage rates. The estimated value of the residual variance was close to its true value, but the estimated variance-covariance matrix of the random-effects was underestimated.

In Dominicus, Ripatti, Pedersen, and Palmgren (2008), the authors used the PME model with random change points to analyze the change in cognitive function during old age. The estimation method used was MCMC simulation using Gibbs sampling. Additionally, the authors executed a small simulation study to investigate the performance of the MCMC Bayesian estimation approach. The simulation results showed that the estimated model parameters, including the estimated fixed-effects parameters, residual variance parameter and variance-covariance matrix of the random effects, were all close to the true parameter values.

In Wang and McArdle (2008), the authors compared the performance of Bayesian estimation to frequentist first-order Taylor series expansion and adaptive Gaussian quadrature for estimating PME models with unknown random change points. They reported that Bayesian estimation yields more accurate and precise parameter estimates as compared to the frequentist methods across all the manipulated simulation conditions. The model convergence rates were also reported to be low for the frequentist methods of estimation.

This leads to the purpose of the current study. That is, to implement the ML approach as described in Du Toit and Cudeck (2009) for estimating linear–linear PME model parameters with unknown random change points. The choice to incorporate the linear–linear function was made because it is a frequently applied function in many educational and psychological research settings. Moreover, its estimation is expected to be more challenging with respect to model convergence as well as the accuracy and precision of parameter estimates (see Appendix A). Using the estimation approach of Du Toit and Cudeck (2009), the parameters of a nonlinear mixed-effects model are separated into linear parameters and intrinsically nonlinear parameters. The dimension of the integral for computing the ML of the marginal distribution of y is reduced by integrating out the linear parameters. Thus, the dimension of the integral is equal to the number of intrinsically nonlinear parameters. To elaborate, in a two-phase linear–linear PME model there are three linear parameters (i.e., intercept and slope parameters of the first phase, and the slope parameter of the second phase) and one intrinsically nonlinear parameter, the change point. This estimation approach will reduce the four-dimensional integral to a one-dimensional integral. The reduction in the dimension of the integral significantly saves computational time and potentially improves the convergence rate as well as the accuracy of parameter estimation as opposed to using the adaptive Gaussian quadrature method. Additionally, in the current study we incorporated the use multiple random start value sets (Hipp & Bauer, 2006) with the ML estimation procedure proposed by Du Toit and Cudeck (2009) to help improve the convergence rate and the overall accuracy and precision of parameter estimates. Moreover, to provide a more complete comparative analysis, we also used MCMC for Bayesian estimation as implemented in the rjags package for R (Plummer, 2016) to provide a legitimate comparison to the results obtained from the ML estimation procedure.

Method

Data generation & simulation design

The following two-phase linear–linear piecewise function was used to generate longitudinal data sets for the purpose of the current study

y_{i} = f (t_{i}, β_{i}) + e_{i},

where y _i = (y_i ₁,…, y_ij ) is a vector of observations for the ith individual on j time points, t _i = (t_i ₁,…, t_ij ) denotes the vector of time points where the ith individual is observed, β _i = (β_i ₁,…, β_ip ) is a vector of individual specific regression parameters, e _i = (e_i ₁,…, e_ij ) is a vector of random errors for the ith individual on j time points_, and f(t _i ,β _i ) is the piecewise function defined as

f (t_{i}, β_{i}) = {\begin{matrix} β_{1 i} + β_{2 i} t_{i j} t_{i j} < γ_{i} \\ β_{1 i} + β_{2 i} γ_{i} + β_{3 i} (t_{i j} - γ_{i}) t_{i j} \geq γ_{i} \end{matrix},

where β _1i and β _2i are the subject-specific intercept and slope parameters of the first phase, β _3i is the subject-specific slope parameter of the second phase, and γ_i is the subject-specific change point for the ith individual. The distribution of random errors is assumed to be multivariate normal with a zero mean vector and a covariance structure of Θ _i , which is typically assumed to follow the independent structure (i.e., diagonal matrix with equal error variances (σ_e ²) are on the diagonal). The distribution of random effects is assumed to be multivariate normal with a zero mean vector and a covariance matrix of Φ.

The simulation design (including the population values for parameters chosen) for this study is close to reported values of previous simulation design used in the literature (Wang & McArdle, 2008). The total number of measurement occasions was eight where individuals were assumed to be measured at time points of 0, 10, 20, 30, 40, 50, 60, and 70. The manipulated conditions in the simulation design were sample size, residual variance, and location of change point. The two levels of the sample size condition were: small (N = 100) and medium (N = 500); the two levels of the change point location were: early occurrence of change point (γ = 20) and change point occurring somewhere in the middle (γ = 35); and for the residual variance the two levels were: small (σ_e ² = 1) and large (σ_e ² = 9). The manipulated conditions were fully crossed, yielding a total of eight conditions in the simulation. Each simulated condition contained 100 replications. Across all simulated conditions, the set of fixed regression parameters were specified as β ₁ = 80, β ₂ = 0.5, and β ₃ = 0.1 and the covariance structure of random effects was specified as a diagonal matrix.

Φ = (\begin{matrix} σ_{β_{1}}^{2} = 36 \\ 0 & σ_{β_{2}}^{2} = 0.0225 \\ 0 & 0 & σ_{β_{3}}^{2} = 0.0025 \\ 0 & 0 & 0 & σ_{γ}^{2} = 1 \end{matrix}) .

Model estimation

The parameters of the PME model under the ML framework were estimated using an R routine¹ (Zopluoglu, Harring, & Kohli, 2014) built upon the ML estimation procedure described in Du Toit and Cudeck (2009). The main advantage of this procedure is that it integrates out the linear parameters (β _1i, β _2i, β _3i) and reduces the dimension of the integral to the number of intrinsically nonlinear parameter(s), which is one in this particular model (i.e., γ), when computing the likelihood of the marginal distribution. In contrast, the adaptive Gaussian quadrature method as implemented in SAS PROC NLMIXED approximates a four-dimensional integral. The R routine based on Du Toit and Cudeck’s (2009) procedure also combines the idea of using multiple sets of random starting values as suggested by Hipp and Bauer (2006) in order to avoid local convergence or non-convergence and thereby, increasing the chance of reaching the global maxima. The number of random starting value sets for the initial stage iterations were fixed to 75.

As a first step, each set of random starting values was generated by fitting a piecewise regression model with unknown change point location to each individual’s data using the segmented package in R (Muggeo, 2008). The individual estimates for β ₁, β ₂, β ₃, γ, and σ_e obtained from the first step subsequently gave us the distribution of each fixed parameter in the model. In the second step, the first and the third quartiles of each fixed parameter distribution were computed. These two values were used as the lower and upper boundaries of the uniform distribution corresponding to respective fixed model parameters, and then a random number was drawn from this uniform distribution as the starting value for these parameters. The starting value for the residual variance was the average of the estimated residual variance across all the individuals (i.e., the rows of the generated data sets). To get the starting values for the variances of β ₁, β ₂, and β ₃ following Hipp and Bauer (2006), a random number was drawn from a uniform distribution whose minimum value was 0.05 and the maximum value was the corresponding variance estimate for the respective parameter distribution obtained from the first step. The starting values for growth factor covariances were randomly drawn from a uniform distribution such that the correlations were between −0.75 and 0.75. Following these aforementioned steps we generated multiple 75 random initial starting value sets.

The initial stage iterations were run for each set of the random starting value sets. The initial stage iterations were stopped when either the number of iterations exceeded 15 or the difference in the log-likelihood values between two successive iterations was smaller than 0.001. For the final stage iterations, the best set obtained from the initial stage iterations was re-estimated until the convergence criterion was met. That is, the final stage iterations were stopped when either the number of iterations exceeded 2000 or the difference in the log-likelihood values between two successive iterations was smaller than 0.0001. Furthermore, for both initial and final stage iterations, 20 Gauss–Hermite quadrature points were used to approximate the integration.

For the Bayesian estimation framework, non-informative prior distributions were considered for all model parameters, including the fixed mean parameters, random error, and the variance-covariance matrix of random effects. The prior distribution for the fixed linear parameters was normal with mean zero and variance defined by a diagonal matrix with the value of 100 on the diagonal. The prior distribution for the random error was inverse gamma with shape and rate parameters equal to 0.001. The prior distribution for the inverse of the linear random effects covariance matrix was a Wishart distribution with a 3×3 identity matrix as the scale matrix and 3 degrees of freedom. The prior distribution for the change point was a uniform distribution in the interval [10, 60] to avoid model identifiability issues that may occur at the boundaries of the time points (i.e., 0 or 70). And the prior distribution for the standard deviation of the random effect of the change point was uniform distribution in the interval [0, 35/2] to warrant the range of time points as feasible values of the random change point. To elaborate, the time points range between 0 and 70 and we have two simulated change point values (γ = 20, 35) whose variation must lie within such range. Taking the middle point of the range and dividing it by two (35/2) yields a reasonable value for the upper bound interval of the uniform distribution to assure that the values of the random change point lay within [0, 70]. In addition, given previous findings related to Bayesian estimation being not sensitive to starting values (Wang & McArdle, 2008), initial values for MCMC were randomly generated by rjags.

Outcome measures

Two outcome measures were computed for both estimation frameworks: average parameter bias and coverage rate. The average parameter bias is equal to the average difference between the respective estimated parameter and the corresponding true value. It was computed using the formula

Average Bias = \frac{\sum_{r = 1}^{R} ({\hat{θ}}_{q r} - θ_{T r u e})}{R},

where ${\hat{θ}}_{q r}$ is the estimated value for the qth parameter in the rth replication (in Bayesian estimation framework this corresponds to the posterior means), θ_TRUE is the corresponding population value, and R is the number of converged replications. Note that average bias and coverage rate were computed for only successfully converged replications in each cell.

In ML, assuming estimation errors are normally distributed, the coverage rate for a 95% confidence interval was computed. One would like (1 − α)100 of one hundred confidence intervals for a given parameter to contain the true value of that parameter, i.e. one desire a coverage rate of (1 − α)100% (Kohli, Hughes, Wang, Zopluoglu, & Davison, 2015). In Bayesian framework coverage rate is computed based on credible intervals, which are analogous to confidence intervals in ML framework. Thus, for the Bayesian approach coverage rate using credible intervals according to the frequentist convention was computed. Note that the coverage rate was computed using only successfully converged replications in each cell.

Results of simulation study

Model convergence rate

Table 1 shows the model convergence rate across the different manipulated conditions for both the ML and the Bayesian estimation procedures. The rate of converged replications for the parameters of the PME model across all the manipulated conditions for the ML estimation (total number of simulation conditions was 8) was found to be between 65% and 96%. At first, it seems that the lower bound of the range is small (i.e., 65%), however, upon closer inspection there is a reason why the convergence rate is small for some of the manipulated conditions. For conditions with large values of residual variance (i.e., σ_e ² = 9) the convergence rate is between 65% and 76%, and for the conditions with small values of residual variance (i.e., σ_e ² = 1) the convergence rate is between 83% and 96%. Thus, we can conclude that overall large residual variance has an adverse effect on the model convergence rate, which is consistent with the results from prior simulation studies related to the estimation of PME models (see, e.g., Kohli & Harring, 2013; Wang & McArdle, 2008).

Table 1.

Model convergence rate across all conditions.

Sample size	Residual variance	Change point location	Convergence rate
Sample size	Residual variance	Change point location	ML	Bayes
100	1	20	0.83	1.00
100	1	35	0.91	1.00
100	9	20	0.71	0.98
100	9	35	0.75	0.98
500	1	20	0.91	1.00
500	1	35	0.96	1.00
500	9	20	0.65	0.97
500	9	35	0.76	0.99

Model convergence rate for the Bayesian estimation approach was, in general, higher than the ML estimation approach. It was found to be between 97% and 100%. The lower bound of convergence rate corresponds to high residual variance and early occurrence of change point (i.e., σ_e ² = 9 and γ =20). However, these results (related to high convergence rate) should be taken with caution because as we will see in the following section, higher convergence rate for the Bayesian estimation framework did not necessarily translate to more accuracy and precision with respect to the estimation of some of the model parameters.

Average parameter bias

Table 2 presents the average bias for each parameter across all the manipulated conditions for both the ML and the Bayesian estimation procedures. As shown in Table 2, the size of the average bias for the mean parameter estimates (i.e., ${\hat{β}}_{1}, {\hat{β}}_{2}, {\hat{β}}_{3} and \hat{γ}$ ) for ML estimation was within the absolute range of 0.00 to 0.22, indicating that the mean parameter estimates, on average, were close to the true values (i.e., β ₁ = 80, β ₂ = 0.5, β ₃ = 0.1, and γ = 20 or 35). The size of the average bias for the residual variance (i.e., ${\hat{σ}}_{e}^{2}$ ) was within the absolute range of 0.00 to 0.40, indicating that overall the recovery of residual variance was decent (true values, $σ_{e}^{2} = 1 or 9$ ). The size of the average bias for the variance of the intercept parameter estimate (i.e., ${\hat{σ}}_{β_{1}}^{2}$ ) was within the absolute range of 0.21 to 2.35, indicating that the intercept variance, on average, was recovered well (true value, $σ_{β_{1}}^{2} = 36$ ). The size of the average bias for the variance of slope parameter estimates (i.e., ${\hat{σ}}_{β_{2}}^{2} and {\hat{σ}}_{β_{3}}^{2}$ ) was within the absolute range of 0.00 to 0.04. Furthermore, from Table 2 it appears that the true value of the variance of slope 2 parameter recovered really well (true value, $σ_{β_{3}}^{2} = 0.0025$ ); however, with respect to the true value of slope 1 parameter the recovery was uneven (true value, $σ_{β_{2}}^{2} = 0.0225$ ). That is, when the value of residual variance was large ( $σ_{e}^{2} = 9$ ), the sample size was small (n = 100) and the change point occurred early (γ =20) we observed more bias in the estimation of slope 1 variance parameter than otherwise. Finally, the size of the average bias for the variance of the change point parameter estimate (i.e., $σ_{γ}^{2}$ ) was within the absolute range of 0.06 to 51.42, indicating poor recovery of the true value (i.e., ${\hat{σ}}_{γ}^{2} = 1$ ).

Table 2.

Average parameter bias across all conditions.

Sample size	Residual variance	Change point location	ML									Bayes
Sample size	Residual variance	Change point location	β ₁	β ₂	β ₃	γ	$σ_{e}^{2}$	$σ_{β_{1}}^{2}$	$σ_{β_{2}}^{2}$	$σ_{β_{3}}^{2}$	$σ_{γ}^{2}$	β ₁	β ₂	β ₃	γ	$σ_{e}^{2}$	$σ_{β_{1}}^{2}$	$σ_{β_{2}}^{2}$	$σ_{β_{3}}^{2}$	$σ_{γ}^{2}$
100	1	20	−0.13	0.00	0.00	−0.03	0.00	−0.93	0.00	0.00	0.69	−0.27	0.01	0.00	−0.18	−0.01	0.84	0.01	0.01	0.66
100	1	35	0.08	0.00	0.00	−0.02	0.00	−0.48	0.00	0.00	0.49	−0.31	0.00	0.00	0.04	−0.01	0.12	0.01	0.01	0.53
100	9	20	−0.15	0.05	−0.01	0.05	−0.40	−1.41	0.04	0.00	51.42	−0.38	0.03	0.00	−0.83	−0.29	0.47	0.03	0.02	2.50
100	9	35	−0.14	0.01	0.00	−0.07	−0.32	−2.35	0.00	0.00	24.83	−0.27	0.00	0.00	0.13	−0.39	0.92	0.02	0.02	4.59
500	1	20	−0.02	0.00	0.00	0.02	0.00	−0.24	0.00	0.00	−0.06	−0.04	0.00	0.00	−0.03	−0.01	−0.22	0.00	0.00	0.17
500	1	35	−0.04	0.00	0.00	−0.01	0.00	0.21	0.00	0.00	0.02	−0.07	0.00	0.00	0.03	−0.02	0.46	0.00	0.00	−0.21
500	9	20	−0.09	0.02	0.00	−0.22	0.01	0.37	0.01	0.00	23.24	−0.04	0.00	0.00	−0.10	−0.30	1.00	0.01	0.01	0.33
500	9	35	0.01	0.00	0.00	−0.02	−0.06	−0.25	0.00	0.00	7.04	−0.07	0.00	0.00	0.24	−0.32	0.19	0.00	0.01	0.20

Upon closer inspection, we found that the estimation of change point variance was poor for the manipulated conditions where the value of residual variance was large (i.e., σ_e ² = 9). There may be a possible explanation for this finding. The estimation of intrinsically nonlinear models (e.g., PME models) can run into the issue of a flat likelihood function, especially when signal to noise ratio is small (e.g., σ_e ² = 9), which may result in a large average bias. Thus, it may explain the poor recovery of the true value of the change point variance parameter.

Average parameter bias for the Bayesian estimation approach was decent for all parameter estimates, except for the variances of the two linear slopes. The average bias for variances of slope parameters ranges from 0.01 to 0.03. This means that relative to true parameter values, a relative bias of up to 33% and 700% for variances of slope parameters of the first and second phases, respectively, was found using a Bayesian estimation approach. The average bias for the variance of the change point (i.e., ${\hat{σ}}_{γ}^{2}$ ) ranged from −0.21 to 4.59 for the Bayesian procedure. This range is small compared to the ML framework, however the absolute value of average bias for change point variance for the manipulated conditions of small residual variance (i.e., ${\hat{σ}}_{γ}^{2} = 1$ ) and large sample size (N = 500) seems a bit higher for the Bayesian approach when compared to the ML approach. Similar pattern was found for the manipulated condition of small residual variance (i.e., ${\hat{σ}}_{γ}^{2} = 1$ ) with small sample size and middle change point (N = 100 and γ =35). In short, the Bayesian estimation approach, on average, produced higher average bias than the ML framework for the mean intercept and the change point parameters (i.e., ${\hat{β}}_{1} and \hat{γ}$ ), the residual variance (i.e., ${\hat{σ}}_{e}^{2}$ ), and the variances of slope parameters (i.e., ${\hat{σ}}_{β_{2}}^{2} and {\hat{σ}}_{β_{3}}^{2}$ ). On the contrary, the Bayesian approach, on average, produced lower bias than the ML framework for the mean of the first slope (i.e., ${\hat{β}}_{2}$ ), the variance of the intercept (i.e., ${\hat{σ}}_{β_{1}}^{2}$ ), and the variance of the change point (i.e., ${\hat{σ}}_{γ}^{2}$ ).

In sum, both ML and Bayesian estimation approaches had challenges with respect to the recovery of some of the model parameters.

Parameter coverage rate

Table 3 shows the average coverage rates of the 95% confidence interval and credible intervals across all the model parameters. The ML procedure provided average coverage rates for the mean parameter estimates (i.e., ${\hat{β}}_{1}, {\hat{β}}_{2}, {\hat{β}}_{3} and \hat{γ}$ ) between 85% and 100%, 86% and 100%, 79% and 100%, and 84% and 100%, respectively. The average coverage rate for the residual variance of the observed variables was between 91% and 100%. The average coverage rates for the estimates of variance parameters of the linear coefficients (i.e., ${\hat{σ}}_{β_{1}}^{2}, {\hat{σ}}_{β_{2}}^{2}, and {\hat{σ}}_{β_{3}}^{2}$ ) was between 75% and 98%, 88% and 100%, and 82% and 100%, respectively. Last, the average coverage rate for the variance of change point parameter estimate was between 44% and 100%.

Table 3.

Coverage rate across all conditions.

Sample size	Residual variance	Change point location	ML									Bayes
Sample size	Residual variance	Change point location	β ₁	β ₂	β ₃	γ	$σ_{e}^{2}$	$σ_{β_{1}}^{2}$	$σ_{β_{2}}^{2}$	$σ_{β_{3}}^{2}$	$σ_{γ}^{2}$	β ₁	β ₂	β ₃	γ	$σ_{e}^{2}$	$σ_{β_{1}}^{2}$	$σ_{β_{2}}^{2}$	$σ_{β_{3}}^{2}$	$σ_{γ}^{2}$
100	1	20	0.93	0.90	0.95	0.88	0.93	0.82	0.91	0.85	0.91	0.96	0.98	1.00	0.93	0.95	0.99	0.06	0.00	0.97
100	1	35	0.93	0.90	0.92	0.92	0.93	0.92	0.92	0.88	0.90	0.95	0.99	1.00	0.97	0.93	0.95	0.11	0.00	0.99
100	9	20	0.85	0.86	0.79	0.84	0.91	0.75	0.88	0.94	0.45	0.91	0.93	1.00	0.90	0.92	0.97	0.00	0.00	1.00
100	9	35	1.00	1.00	1.00	1.00	1.00	0.87	1.00	1.00	0.44	0.94	0.97	1.00	0.95	0.90	1.00	0.02	0.00	1.00
500	1	20	1.00	0.93	0.96	0.87	0.98	0.90	0.96	0.82	0.84	0.89	0.93	0.99	0.91	0.94	0.93	0.74	0.00	0.97
500	1	35	1.00	1.00	1.00	1.00	1.00	0.98	1.00	1.00	1.00	0.92	0.97	0.98	0.94	0.91	0.92	0.70	0.00	0.94
500	9	20	1.00	1.00	1.00	0.94	1.00	0.86	1.00	1.00	0.74	0.97	0.97	0.99	0.94	0.75	0.95	0.14	0.00	1.00
500	9	35	1.00	1.00	1.00	1.00	1.00	0.92	1.00	1.00	0.88	0.97	0.93	0.98	0.86	0.75	0.93	0.59	0.00	1.00

The Bayesian estimation approach provided a high coverage rate of 90–100% in cells pertaining to the estimation of parameters: ${\hat{β}}_{2}, {\hat{β}}_{3}, {\hat{σ}}_{β_{1}}^{2},$ and ${\hat{σ}}_{γ}^{2}$ but not for the rest of the model parameter estimates. The coverage rate was especially poor (between 0% and 74%) for the variances of slope parameters (i.e., ${\hat{σ}}_{β_{2}}^{2} and {\hat{σ}}_{β_{3}}^{2}$ ), which was expected given the high values of relative bias associated with them. Coverage rate for the variance of change point parameter estimates was higher using the Bayesian procedure. It was found to be between 94% and 100%.

Overall, we can conclude that the ML estimation approach in this study did a good job of estimating model parameters under the condition of small residual variance of the observed variables on the outcome measures—convergence rate, average parameter bias, and coverage rate. Additionally, we found that, in general, Bayesian estimation performed poorly with respect to the recovery of the variances of slope parameters across all simulated conditions, but it performed better with respect to the estimation of the variance of the change point parameter.

One of the main findings of this study is that the ML estimation approach performs satisfactorily under the condition of small residual variance (i.e., σ_e ² = 1). Upon closer inspection of the simulated data, we found some general guidelines for when a researcher may take ML estimation results with caution. We computed the mean ratio of the residual variance to the total variance of the response variable across time points for each data set and averaged this index over the 100 replications within each simulated condition. That is, we computed the mean error variance to overall variance ratio for each data set using the formula

Mean ratio = \frac{1}{J} \sum_{j = 1}^{J} \frac{σ_{e}^{2}}{Var (y_{j})}

where $y_{j} = (y_{1 j}, \dots, y_{N j})$ is a vector of observations for all individuals at time point j. Table 4 shows this ratio averaged over all 100 replicated data sets within each simulated scenario, as well as its minimum and maximum values. We found that the condition of large residual variance (i.e., σ_e ² =9) presents ratios between 0.12 and 0.22, while the ratios for the condition of small residual variance (i.e., σ_e ² = 1) ranged from 0.01 to 0.03. Therefore, if the mean error variance and overall variance ratio for a given data set is 0.12 or bigger and the estimation method is ML, model results should be viewed with caution.

Table 4.

Error variance and overall variance ratio.

Sample size	Residual variance	Change point location	Error variance to overall variance ratio
Sample size	Residual variance	Change point location	Mean	Min	Max
100	1	20	0.02	0.02	0.03
100	1	35	0.02	0.01	0.03
100	9	20	0.17	0.13	0.22
100	9	35	0.15	0.12	0.19
500	1	20	0.02	0.02	0.03
500	1	35	0.02	0.02	0.02
500	9	20	0.17	0.15	0.19
500	9	35	0.15	0.13	0.17

Discussion

In summary, the current study expanded on prior simulation-based literature related to the estimation of PME models with unknown random change points. We used and evaluated the performance of the ML estimation using the approach proposed by Du Toit and Cudeck (2009) and the MCMC for Bayesian estimation with respect to the accuracy and precision of estimated linear–linear PME model parameters via a small Monte Carlo simulation study. The key findings from the simulation study were twofold. First, neither ML nor Bayesian estimation performed well under all manipulated conditions with respect to the accuracy and precision of the estimated model parameters, which is inconsistent with previous findings (Wang & McArdle, 2008). Wang and McArdle (2008) concluded that across all the manipulated simulation conditions frequentist first-order Taylor series expansion and adaptive Gaussian quadrature performed poorly in comparison to the Bayesian approach with respect to the quality of parameter estimates as well as the model convergence rates. In our study, results were not conclusive in favor of either estimation procedures. Second, we learned that under the condition of small residual variance of the observed variables the ML estimation approach did recover reliable and precise model parameter estimates, and that the size of the residual variance had the most impact on the quality of parameter estimates. To put this finding in the context of real data setting, often it is not possible to know the reliability of data collected beforehand. Some research settings by default produce data that has lots of noise, thus has large residual variance in the observed variable. For example, in the research area of Curriculum-Based Measurement of Oral Reading or CBM-R data collected weekly from second grade students across several weeks tend to have small signal to noise ratio, in general. But in other research settings, such as longitudinal controlled clinical trials, the data may have less noise and thus be more reliable. The decision to use ML estimation procedure for fitting PMEs should be based on several factors, such as sample size, residual variance, complexity of the models, and so on. An aspect to consider when deciding on the quality of the estimation results using the ML framework is the ratio of the mean error variance to overall variance. If this ratio is equal to 0.12 or higher, estimation results should be viewed cautiously. Last, the potential occurrence of flat likelihood function in the event of small signal to noise ratio when estimating intrinsically nonlinear parameter(s) is a well-known issue.

Overall, PME models provide a flexible framework to researchers and practitioners that enables them to characterize individual behavior that exhibits distinct phases of development in each segment, and thereby allow researchers to address key questions such as developmental studies or studies seeking to measure the effect of treatment/ intervention, and so forth such as when individuals may need to seek professional services at the timing when mental ability decreases. The estimation of different variations of PME models is an active area of research, where both frequentist and Bayesian procedures have been used for model fitting purposes.

Footnotes

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Note

Appendix A

The estimation of linear–linear piecewise function is more challenging relative to the estimation of quadratic-linear piecewise function used in the Wang and McArdle (2008) study. To elaborate on this point, let’s take a quadratic-linear piecewise function expressed as

where y_ij is the observation for the ith individual at the jth time point, x_ij refers to the time point at which the ith individual is observed (the time scale includes eight time points—0, 10, 20, 30, 40, 50, 60 and 70), τ_i is the subject-specific nonlinear change point, and $β_{i} = (b_{0 i}, b_{1 i}, b_{2 i})^{'}$ is a vector of subject-specific linear regression parameters of intercept, linear and quadratic term, respectively. Because of the presence of quadratic term in the first segment of the estimated and population model, the model estimation benefitted greatly with respect to model convergence, and accuracy and precision of the estimated model parameters. The reason for making this conclusion is attributed to a smaller ratio of the residual variance to the total variance of the response variable, y, as a result of the inclusion of the quadratic term. In addition, in a quadratic-linear PME model, any variation in the location of the change point produces larger effects on y than an analogous linear–linear model would. For example, let’s assume the quadratic coefficient is 1 and consider time point 0. A difference in the change point from 16 to 15 would imply a difference of 31 in y (i.e., $[1 * {(16 - 0)}^{2}] - [1 * {(15 - 0)}^{2}]$ ), a difference in the change point from 21 to 20 would imply a difference of 41 in y, etc. This allows for precise estimation of the change point and the quadratic slope (occurring before the change point) under the simulation scenarios considered. The advantage to the model estimation by incorporating the quadratic term in the first segment would not exist if the first segment was modeled by a linear function. That is, if the function is linear–linear piecewise, then the component $b_{2 i} {(τ_{i} - x_{i j})}^{2}$ would be replaced by $b_{2 i} (τ_{i} - x_{i j})$ which will not produce as big of an effect on the scale of y. Thus, one would expect the estimation of linear–linear piecewise function, in general, to be more challenging with respect to model convergence as well as the accuracy and precision of estimated model parameters, irrespective of the estimation approach used.

References

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, 213–225.

Cudeck

Klebe

K. J.

(2002). Multiphase mixed-effects models for repeated measures data. Psychological Methods, 7, 41–63.

Dominicus

Ripatti

Pedersen

N. L

Palmgren

. (2008). A random change point model for assessing the variability in repeated measures of cognitive function. Statistics in Medicine, 27, 5786–5798.

Du Toit

S. H. C.

Cudeck

(2009). Estimation of the nonlinear random coefficient model when some random effects are separable. Psychometrika, 74, 65–82.

Grimm

K. J.

(2008). Longitudinal associations between reading and mathematics achievement. Developmental Neuropsychology, 33, 410–426.

Hipp

J. R.

Bauer

D. J.

(2006). Local solutions in the estimation of growth mixture models. Psychological Methods, 11, 36–53.

Jordan

N. C.

Hanich

Kaplan

(2003). A longitudinal study of mathematical competencies in children with specific mathematics difficulties versus children with co-morbid mathematics and reading difficulties. Child Development, 74, 834–850.

Kohli

Harring

J. R.

(2013). Modeling growth in latent variables using a piecewise function. Multivariate Behavioral Research, 48, 370–397.

Kohli

Harring

J. R.

Zopluoglu

. (2016). A finite mixture of nonlinear random coefficient models for continuous repeated measures data. Psychometrika, 81, 851–880.

10.

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, 259–275.

11.

Kreisman

M. B.

(2003). Evaluating academic outcomes of Head Start: An application of general growth mixture modeling. Early Childhood Research Quarterly, 18, 238–254.

12.

Miles

T. R.

Miles

(1992). Dyslexia and mathematics. London, UK: Routledge.

13.

Morrell

C. H.

Pearson

J. D.

Carter

H. B.

Brant

L. J.

(1995). Estimating unknown transition times using a piecewise nonlinear mixed-effects model in men with prostate cancer. Journal of the American Statistical Association, 90, 45–53.

14.

Muggeo

V. M. R.

(2008). Segmented: an R package to fit regression models with broken-line relationships. R News, 8/1, 20–25. Retrieved from http://cran.r-project.org/doc/Rnews/.

15.

Paris

S. G.

(2005). Reinterpreting the development of reading skills. Reading Research Quarterly, 40, 184–202.

16.

Plummer

. (2016). rjags: Bayesian graphical models using MCMC. R package version 4-6.

17.

Shin

Davison

M. L.

Long

Chan

C.-K.

Heistad

(2013). Exploring gains in reading and mathematics among regular and exceptional students using growth curve modeling. Learning and Individual Differences, 23, 92–100.

18.

Silverman

R. D.

Speece

D. L.

Harring

J. R.

Ritchey

K. D.

(2012). Fluency has a role in the simple view of reading. Scientific Studies of Reading, 16, 1–26.

19.

Wang

McArdle

J. J.

(2008). A simulation study comparison of Bayesian estimation with conventional methods for estimating unknown change points. Structural Equation Modeling, 15, 52–74.

20.

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.

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.45 MB