A Spline Regression Model for Latent Variables

Abstract

Spline (or piecewise) regression models have been used in the past to account for patterns in observed data that exhibit distinct phases. The changepoint or knot marking the shift from one phase to the other, in many applications, is an unknown parameter to be estimated. As an extension of this framework, this research considers modeling the relation between endogenous and exogenous latent variables with a spline regression model, where each latent variable is measured by multiple observed indicators. Consequently, the spline regression model is modified to include a measurement model that explicitly expresses the relation of the observed variables to the latent constructs. Maximum likelihood estimation of the model is developed and executed on educational data to illustrate the utility of the model.

Keywords

piecewise regression latent variable regression spline regression models structural equation modeling nonlinear

Introduction

Spline (or piecewise)¹ regression models have been used in the past to account for nonlinearity in observed data that exhibit distinct phases. The general theory of two-phase, piecewise models is reviewed by Seber and Wild (1989, chap. 9). These models are flexible because each segment can be specified to conform to a particular aspect of the overall response–predictor relation. A typical scenario might break the observed pattern into piecewise linear components with linear rates of change differing across ranges of the predictor. For example, the data graphed in the scatterplot in Figure 1 consist of 72 measurements on the age in months (x) and the weight/height ratio (y) of preschool boys.² Note that age ranges from 0.5 to 71.5 months in 1-month intervals.

In the first phase of development, weight/height ratio improves rapidly until about the eighth or ninth month, after which the ratio follows a linear trend distinct from the first phase, in which the ratio continues to grow, but just not as rapidly. A piecewise model suggested by these data is a segmented, linear–linear process with distinct pieces for the early and later ages. An examination of the plot shows that the changepoint, or knot, marks the shift from one phase to the other.³ It is straightforward to consider different patterns such as quadratic–linear⁴ or linear–exponential to allow a nonlinear trajectory in one segment and a steady-state condition in the other.

In other contexts, the situation may arise when there are more than two observed phases over the range of the predictor (Cudeck & Codd, 2012). In these situations, multiple knots would be necessary to account for the various transitions. Although the number of options can appear overwhelming, it is important to underscore that the functional form of each segment can be tailored to capture the interesting facets of the relation (Cudeck & Harring, 2010). One of the most interesting features of a piecewise model is the knot or changepoint. The knot, a parameter to be estimated in many applications, is the value of the predictor where the segments meet. Often the knot has special significance scientifically (see, e.g., Bacon & Watts, 1971) as it may signify some watershed event, treatment, intervention, or other important aspect fundamental to the bivariate relation. For instance, in the forthcoming empirical example in which students’ verbal reasoning depends on their verbal acuity, there is a verbal acuity proficiency threshold (e.g., knot) beyond which verbal reasoning increases more rapidly than before reaching the threshold. It is just this type of insight that makes piecewise models invaluable and more appropriate perhaps than traditional functions such as higher order polynomials, whose model parameters provide less direct information about interesting facets fundamental to the underlying behavioral process.

In instances when several observed indicators measuring the same latent construct are collected, the piecewise relation between latent criterion and predictor can be investigated within a general latent variable regression modeling framework. Latent variable models corresponding to their manifest variable counterparts have the distinct advantage that latent variables are by definition error-free and thus the models are not tainted by unaccounted measurement error. Conceptually, the jump from an observed to latent variable framework for piecewise regression models is fairly clear-cut. However, a general methodological problem exists to estimate such a model in practice. Like piecewise regression models for observed variables, even if the regression coefficients enter the function linearly, the unknown knot makes the entire system intrinsically nonlinear. For observed variables, this would suggest using an iterative procedure such as Newton’s method for nonlinear least squares estimation or a Newton–Raphson algorithm for maximum likelihood. In the context of latent variable regression, the inherent nonlinearity poses similar estimation challenges.

Although general nonlinear structural equation models are a very broad class (Wall & Amemiya, 2007), a considerable amount of attention in recent years has been devoted to polynomial models, including cross-product terms among the latent variables representing interactions. Piecewise latent variable regression models with an unknown knot are subsumed within this class. Because these models have at least one nonlinear coefficient corresponding to a latent variable, perhaps more than one, maximum likelihood estimation assuming a normal distribution for the exogenous predictor(s) is more difficult than in standard structural equation models that specify only linear relations (Klein & Muthén, 2007; Wall & Amemiya, 2000).

This article considers a piecewise regression model for investigating data that exhibit distinct phases. Most common applications of piecewise models center on a single observed response variable which changes in a systematic, but potentially different, manner across values of the observed predictor. The approach considered here treats two or more observed variables but with a latent variable structure. This is accomplished by augmenting the piecewise regression model with a measurement model to capture the relation between multiple indicators and their corresponding latent construct. Gaussian–Hermite quadrature is presented as a method for maximum likelihood estimation of the model. This approach has its liabilities, in particular, the slow speed to convergence as the number of dimensions of integration increases (Cudeck, Harring, & du Toit, 2009). In the current application, this deficit notwithstanding, because there are only two latent variables to integrate over, any impact to the computation time is thought to be negligible.

The remainder of the article adheres to the following outline. In the next section, the piecewise model will be explicated within a latent variable modeling framework. Estimation of the model will be subsequently discussed. Finally, an empirical data set will be used to demonstrate the utility of the model.

Piecewise Regression Model for Latent Variables

Given a p-dimensional vector of manifest variables y for the ith individual of the sample of size n and a q-dimensional vector of latent variables f such that f = (f₁, f₂)′, the measurement model and nonlinear structural model are specified as

y = τ + Λ f + e,

f_{1} = {\begin{matrix} β_{1} + β_{2} f_{2} + d & f_{2} \leq γ \\ β_{3} + β_{4} f_{2} + d & f_{2} > γ \end{matrix},

where in Equation (1), τ is the p× 1 vector of observed variable intercepts, Λ is the p× 2 matrix of factor loadings relating each latent variable in f to its measured variable indicators in y and e represents the p× 1 vector of unique factors independent of f with E(e) = 0 and var(e) = Ψ. For the structural model in Equation (2) it is assumed the residual d has E(d) = 0 and $var (d) = σ_{d}^{2}$ , and is independent of f₂ as well as e. From this point forward, the residual notation in the structural model will be suppressed to facilitate the formulation and development of the piecewise model. Assigning a scale to the latent variables is accomplished by setting the factor loading of a single indicator for each latent variable to 1 and its corresponding intercept to 0. Equation (1) can be any identified structure that is appropriate for the problem under investigation. In the example reviewed subsequently, p = 5 with a measurement model given by

\begin{matrix} (\begin{matrix} y_{1} \\ y_{2} \\ y_{3} \\ y_{4} \\ y_{5} \end{matrix}) = (\begin{matrix} 0 \\ τ_{2} \\ τ_{3} \\ 0 \\ τ_{5} \end{matrix}) + (\begin{matrix} 1 & 0 \\ λ_{21} & 0 \\ λ_{31} & 0 \\ 0 & 1 \\ 0 & λ_{52} \end{matrix}) (\begin{matrix} f_{1} \\ f_{2} \end{matrix}) + (\begin{matrix} e_{1} \\ e_{2} \\ e_{3} \\ e_{4} \\ e_{5} \end{matrix}) \\ y = τ + Λ f + e . \end{matrix}

From Equation (2), the intercept and slope of the first phase are β₁ and β₂, and those of the second phase are β₃ and β₄, respectively. Including the knot, γ, this structural model specifies five parameters. If it is presumed that the functions characterizing the two phases join at the changepoint, then the function values at γ are equal (i.e., $β_{1} + β_{2} γ = β_{3} + β_{4} γ$ ). This implies that one parameter is superfluous and can be eliminated. Of the four regression parameters, β₃ seems the least interesting as it corresponds to the value of f₁ at f₂ = 0—a point not pertinent to the second phase, although the choice is entirely arbitrary. By rearranging terms in the equality constraint, we see that $β_{3} = β_{1} + β_{2} γ - β_{4} γ$ . The modified function is now

f_{1} = {\begin{array}{l} β_{1} + β_{2} f_{2} & f_{2} \leq γ \\ β_{1} + β_{2} γ + β_{4} (f_{2} - γ) & f_{2} > γ \end{array},

resulting in a model with three linear coefficients, $ξ = (β_{1}, β_{2}, β_{4})'$ , and one nonlinear parameter—the knot, γ. A more compact expression is convenient and emphasizes that the two segments make up a single regression function. It is written as

f_{1} = β_{1} + β_{2} f_{2} + (β_{2} - β_{4}) {(γ - f_{2})}_{+}

The last term is the truncation operator, defined as zero or the positive value of the argument, depending on its sign:

v_{+} = {\begin{array}{l} v & v \geq 0 \\ 0 & otherwise \end{array}

As examples, (6.3 − 4)₊ = 2.3, (1 − 0.9)₊ = 0.1, but (−2)₊ = 0 and (1 − 5)₊ = 0.

Interestingly, the piecewise linear model in Equation (3) suggests a steady increasing rate in the first phase (β₂), yet changes instantly at f₂ = γ to potentially a very different rate in the second phase (β₄). This type of abrupt change in rate at γ does not seem very plausible if f₂ is viewed as a continuous variable,⁵ the implication being that values on the endogenous latent variable f₁ corresponding to a one-unit change in the latent exogenous predictor are associated with distinctly different slopes. A more realistic supposition, perhaps, is that the transition between phases is smooth and gradual. Several modifications to the linear–linear model have been suggested in the literature resulting in a less abrupt shift from one phase to the next (see, e.g., Bacon & Watts, 1971; Seber & Wild, 1989).⁶ Following Griffiths and Miller (1973), the form of the model adopted here is

f_{1} = α_{1} + α_{2} (f_{2} - γ) + α_{3} \sqrt{{(f_{2} - γ)}^{2} + δ} δ \geq 0,

where α₁ is the expected value of f₁ at the knot, γ (i.e., the knot is $(γ, α_{1})$ ); α₂ is the average slope, $(β_{2} + β_{4}) / 2$ ; and α₃ is half the difference in slopes, $(β_{4} - β_{2}) / 2$ . Harring et al. (2006) demonstrated how the value of the parameter δ regulated the transition between the two linear phases. As $δ \to 0$ , the function becomes the standard linear–linear model, and approaches the functional form expressed in Equation (4) from below. As the value of δ increases, the function rounds the corner in an increasingly gradual manner, approaching the linear function of the second segment. This smoothing parameter can either be estimated, or in cases where the data are sparse—such as individuals’ data in longitudinal studies—set equal to a value commensurate with the amount of smoothing warranted by the data. After some initial exploration, δ was set at 0.01 in the subsequent application.

It should be noted that the transformed parameters no longer correspond to change characteristics of the linear–linear processes like those parameters in the original formulation of the model (e.g., Equation 4). However, if there is empirical evidence that a smooth, rather than an abrupt, function is required, then the formulation in Equation (5) may be justified. Another intended benefit of the transformation will occur when carrying out hypothesis tests on the coefficients. Testing α₃ in Equation (5) is a test if the slopes between the two regimens are the same, that is, whether there is a significant phase shift or not. This can be carried out by comparing the estimate with its standard error. Furthermore, since the regression coefficients, α, in Equation (5) are linear functions of the parameters, β, in Equation (4), on convergence the multivariate delta method (Oehlert, 1992; see the appendix) could be used to compute estimates and standard errors of parameters of the original model if desired.

It should be noted that the values of the two functions representing the different phases are not restricted to be equal at the knot. That is, there may be a discontinuity in the endogenous–exogenous relation at this point of transition. For example, in the context of random coefficient models for repeated measures data, Cudeck and Codd (2012) examined a phenomenon called reminiscence, which described unrehearsed improvement from one time point to the next. The authors developed a piecewise model to account for the disjointed change across the span of the study period.

Model Estimation

Conditional Distributions

Nonlinear structural models such as Equation (5) have a combination of endogenous and exogenous latent variables with the caveat that the exogenous latent variable enters the model in a nonlinear fashion. This was also true with Equation (4). More generally, a broad class of nonlinear structural models follows this same formula—at least one exogenous latent variable enters the regression model in a nonlinear fashion (Wall, 2009). The primary reason for making the endogenous–exogenous distinction is to facilitate improved computational efficiency of estimation (du Toit & Cudeck, 2009). Equation (5) is reproduced here for convenience,

f_{1} = α_{1} + α_{2} (f_{2} - γ) + α_{3} \sqrt{{(f_{2} - γ)}^{2} + δ} + d

Following Cudeck et al. (2009), for a particular value of the nonlinear exogenous factor, $f_{2} = {\tilde{f}}_{2}$ the endogenous latent variable is now a linear function

\begin{matrix} f_{1} | {\tilde{f}}_{2} & = α_{1} + α_{2} ({\tilde{f}}_{2} - γ) + α_{3} \sqrt{{({\tilde{f}}_{2} - γ)}^{2} + 0.01} + d \\ = Γ ({\tilde{f}}_{2}, γ) α + d . \end{matrix}

where Γ is a design matrix in which elements are either scalar constants or functions of the nonlinear parameter, γ. A typical row of Γ could be expressed as,

Γ_{q •} = [\begin{matrix} 1 & ({\tilde{f}}_{2} - γ) & \sqrt{{({\tilde{f}}_{2} - γ)}^{2} + 0.01} \end{matrix}]

It is assumed that the distribution of the exogenous latent variables is normal with moments:

E [f_{2}] = κ_{2} var (f_{2}) = φ_{22}

The implication of conditioning on the nonlinear exogenous factor becomes clearer. Integrating the joint distribution of y and f over the latent factors now has a closed form solution (the marginal distribution of y can be derived algebraically). Of practical significance is that the integration problem has been reduced to a single dimension.

The latent variables, $f = (f_{1}, f_{2})'$ , have mean vector conditional on $f_{2} = {\tilde{f}}_{2}$ ,

κ = E (f | {\tilde{f}}_{2}) = (\begin{matrix} Γ ({\tilde{f}}_{2}, γ) α \\ {\tilde{f}}_{2} \end{matrix})

and corresponding covariance matrix

ϕ = cov (f | {\tilde{f}}_{2}) = (\begin{matrix} var (f_{1} | {\tilde{f}}_{2}) \\ 0 & var ({\tilde{f}}_{2}) \end{matrix}) = (\begin{matrix} φ_{11} \\ 0 & 0 \end{matrix})

and where

\begin{matrix} var (f_{1} | {\tilde{f}}_{2}) & = var (Γ ({\tilde{f}}_{2}, γ) α + d) \\ = var (d) \\ = σ_{d}^{2} \end{matrix} .

These results imply that the measurement model for y given ${\tilde{f}}_{2}$ is

y | {\tilde{f}}_{2} = τ + Λ (f | {\tilde{f}}_{2}) + e .

Although the distribution of y is not normal, the conditional distribution of $y | {\tilde{f}}_{2}$ is

y | {\tilde{f}}_{2} ~ N (µ, Σ),

where

µ = τ + Λ κ, Σ = Λ ϕ Λ' + Ψ .

The matrix of unique variable variances is $Ψ = Diag (ψ_{1}, \dots, ψ_{p})$ , with $Ψ_{k} = var (e_{k})$ .

The derivation leading to Equations (10) and (11) completes the description of the model. For this example, there are 18 parameters:

θ' = (α_{1}, α_{2}, α_{3}, γ, τ_{2}, τ_{3}, τ_{5}, λ_{2}, λ_{3}, λ_{5}, κ_{2}, φ_{22}, ψ_{1}, \dots, ψ_{5}, σ_{d}^{2}) .

Marginal Distribution

The nonlinear exogenous latent predictor has distribution $f_{2} ~ N (κ_{2}, φ_{22})$ with density $P (f_{2} | θ_{f_{2}})$ ,

P (f_{2}) = {(2 π)}^{- \frac{1}{2}} {| φ_{22} |}^{- \frac{1}{2}} \exp [- \frac{1}{2} {(f_{2} - κ_{2})}^{2} / φ_{22}] .

From Equation (10), the conditional distribution of $y | {\tilde{f}}_{2}$ has density

P (y | f_{2}) = {(2 π)}^{- \frac{p}{2}} {| Σ |}^{- \frac{1}{2}} \exp [- \frac{1}{2} (y - µ)' Σ^{- 1} (y - µ)] .

The joint density can then be specified,

P (y, f_{2}) = P (y | f_{2}) P (f_{2}) .

The marginal distribution of y can now be expressed by a one-dimension integral

\begin{matrix} P (y) = \int P (y, f_{2} | θ) \partial f_{2} \\ = \int P (y | f_{2}) P (f_{2}) \partial f_{2} \end{matrix} .

Among the many modern statistical techniques available to handle the integration suggested in Equation (13), techniques that are based on direct approximation, such as Gaussian–Hermite quadrature, may be preferred because of its straightforward implementation and accessibility to practitioners.

Gaussian–Hermite Quadrature

Gaussian–Hermite quadrature is used for numerical integration in statistically related applications because of its relation to Gaussian densities. Pinheiro and Bates (1995) and Davidian and Gallant (1993) established much of the theory surrounding the use of Gaussian–Hermite quadrature to approximate the likelihood for the nonlinear mixed-effects model, while Skrondal and Rabe-Hesketh (2004) demonstrated how this approach could be applied to a general class of latent variable models. Patefield (2002), and later Wall (2009), demonstrated how this same technique could be used for nonlinear structural equation models. Gaussian quadrature is a numerical integration method that approximates integrals of functions using a weighted average of the integrand evaluated at a predetermined set of abscissa (Davis & Rabinowitz, 1984).

With this rule, an integral over a function which can be expressed as $u (t) = f (t) \exp (- t^{2})$ is approximated by the sum

\int_{t} f (t) \exp (- t^{2}) dt \approx \sum_{k = 1}^{Q} w_{k} f (v_{k}),

where $w_{k}$ and $v_{k}$ are the weights and nodes of the Hermite polynomial of degree Q (Krommer & Ueberhuber, 1994, sec. 4.2.6). Define $t = {(2 φ_{22})}^{- \frac{1}{2}} (f_{2} - κ_{2})$ so that for a change of variables, $f_{2}^{(t)} = t \sqrt{2 φ_{22}} + κ_{2}$ with Jacobian $| \frac{\partial f_{2}^{t}}{\partial t} | = \sqrt{2 φ_{22}}$ . The marginal distribution in Equation (13) under the change of variables is,

P (y) = K \int_{t} {| Σ |}^{- \frac{1}{2}} \exp [- \frac{1}{2} (y - µ)' Σ^{- 1} (y - µ)] \exp (- t^{2}) \partial t

with $K = \sqrt{2} {(2 π)}^{- \frac{p + 1}{2}}$ . When the quadrature rule is applied, the distribution of y is then approximated as

P (y) \approx K \sum_{k = 1}^{Q} w_{k} {| Σ^{*} |}^{- \frac{1}{2}} \exp [- \frac{1}{2} (y - µ^{*})' {Σ^{*}}^{- 1} (y - µ^{*})] .

It should be noted that where the first two moments of f in Equation (12) were functions of ${\tilde{f}}_{2}$ , due to the change in variables the moments in Equation (16) are functions of the nodes $v$ (i.e., $µ^{*} = µ (f_{2}^{x})$ and $Σ^{*} = Σ (f_{2}^{x})$ , where $f_{2}^{x} = v_{k} \sqrt{2 φ_{22}} + κ_{2}$ ). By increasing the number of quadrature points to an arbitrarily large number, the approximation implied in Equations (15) and (16) becomes more refined. When integrating over a single dimension, like the current example, a large number of quadrature points could be used if warranted without distressing about potentially slow time to convergence. Yet there is generally a time/accuracy trade-off that governs most problems. This would more likely be relevant in situations that require integrating over multiple dimensions.

The log-likelihood function to be maximized for a sample of N observations, $y_{1}, \dots, y_{N}$ can be written as

\ln L (θ | y_{1}, \dots, y_{N}) = \sum_{i = 1}^{N} \ln P (y_{i}) .

For the data analysis in the next section, a Newton–Raphson optimization algorithm was used with Q = 30 quadrature points.⁷

Piecewise Model for Latent Verbal Reasoning

To illustrate the utility of the model and estimation method, data will be used that were obtained from a large educational monitoring study of high school students. Among the multitude of collected information, variables pertaining to verbal reasoning and verbal acuity were gathered.

Variables making up the measurement model in Equation (1) are

\begin{matrix} f_{1} : VerbalReasoning y_{1} : Verbal reasoning \\ y_{2} : 2 - D Spatial Perception \\ y_{3} : 3 - D Spatial Perception \end{matrix}

\begin{matrix} f_{2} : VerbalAcuity y_{4} : Classification of Word Pairs \\ y_{5} : Pattern Completion \end{matrix} .

Table 1 shows sample statistics for a sample of 10th-grader students (N = 250).

Table 1.

Sample Correlations, Means, and Standard Deviations for Five Verbal-Based Measures With a Sample of N = 250.

Variable	1	2	3	4	5
y₁: Verbal reasoning	1
y₂: 2-D Spatial perception	.569	1
y₃: 3-D Spatial perception	.436	.511	1
y₄: Classification of word pairs	.393	.367	.378	1
y₅: Pattern completion	.347	.385	.370	.480	1
Mean	19.14	20.82	15.03	15.69	12.89
Standard deviation	4.15	6.42	4.21	8.13	6.63

Fitting any nonlinear structural model requires a bit of finesse at the onset because convergence using methods such as Gaussian–Hermite quadrature to obtain maximum likelihood estimates are quite sensitive to, and actually depend on, the quality of starting values used in the algorithm. Initial values for the regression coefficients in the nonlinear structural model in Equation (5) were obtained by first fitting the manifest variables with a confirmatory factor analysis model allowing the factor means to be estimated and factors to correlate. Using the estimates of model parameters at convergence, predicted factor scores are then computed and submitted to a regression analysis using nonlinear least squares estimation. An approximate starting value for the variance of the regression error (i.e., var(d)]) in the structural model can also be obtained from this analysis. Starting values for the measurement model can be obtained vis-à-vis the confirmatory factor analysis, which can be handled by any structural equation modeling software program.

Maximum Likelihood Estimates

Table 2 contains the maximum likelihood estimates for the fitted structural model, while Table 3 contains estimates for the measurement portion of the model. Standard errors for the parameters were obtained as a by-product of the Newton–Raphson algorithm, which uses the inverse of the Hessian matrix (second partial derivates of the log-likelihood function with respect to the parameter vector, evaluated at $\hat{θ}$ after convergence):

SE ({\hat{θ}}_{j}) = {[- H {(\hat{θ})}^{- 1}]}_{jj}^{1 / 2} .

Table 2.

Parameter Estimates and Standard Errors for the Piecewise Latent Regression Model.

Parameter	Estimate	Standard Error
α ₁	12.52	0.58
α ₂	0.85	0.10
α ₃	0.35	0.04
γ	17.22	0.82
κ ₂	19.03	0.43
ϕ ₂₂	11.56	1.73
$σ_{d}^{2}$	0.43	1.43

Table 3.

Parameter Estimates and Standard Errors for the Measurement Model.

Parameters	Estimates	Standard Errors
λ ₂	1.09	0.02
λ ₃	1.06	0.03
λ ₅	0.87	0.03
ψ ₁	4.96	1.03
ψ ₂	21.69	2.11
ψ ₃	8.62	1.28
ψ ₄	46.86	6.04
ψ ₅	30.55	3.68

Each of the regression parameters are large compared with their standard errors, thus would be considered statistically significant at all nominal levels. At the knot ( $\hat{γ} = 17.2$ ), the value of verbal reasoning is ${\hat{α}}_{1} \approx 12.5$ . Recall that the knot is a point on the verbal acuity scale where the transition between the first- and second-phases occurs. Of particular interest, is the test of whether two linear functions are actually needed to represent the two phases ( $H_{0} : β_{4} - β_{2} = 0$ ). This effect, as was previously defined, is captured by the coefficient, α₃, which is statistically significant ( ${\hat{α}}_{3} = 0.35, se ({\hat{α}}_{3}) = 0.04$ ).

To see how well the model fits the data, Figure 2 displays the fitted structural model superimposed on a graph with individuals’ predicted latent variable scores on f₁ and f₂. Strikingly, the relation between verbal reasoning and verbal acuity in the first phase shows a slow, steady increase until verbal acuity is approximately 17 on the scale of the word pair classification measure. At that point verbal ability appears to take off—increasing in a more rapid manner than when there was less verbal acuity.

Of particular note in Figure 2 is the fact that some subjects’ bivariate factor scores in the first phase are separated from the mass of points showing more variability than those in the second phase. Consequently, the linear function in the first phase has been pulled toward these influential points. Although not acted on here, a sensible course of action may be to compute a diagnostic, case deletion statistic such as Cook’s distance to assess and better understand the effects of individual data on the analysis. Another reasonable approach is to fit the structural model with an estimation algorithm, such as generalized least-squares, that gives less weight to points in the regression that are further away from the fitted function.

Based on a visual inspection of Figure 2, the smooth linear–linear piecewise model fits the data well. One might wonder if other functional forms may fit slightly better. For completeness, three other functions were also fitted to the data: (a) a linear–linear piecewise model with abrupt transition at the changepoint (i.e., the function in Equation [4]), (b) a quadratic–linear piecewise model (Cudeck & Harring, 2010) with first- and second-order continuity constraints, and (c) a simple quadratic model of the form

f_{1} = β_{0} + β_{1} (f_{2} - κ_{2}) + β_{2} {(f_{2} - κ_{2})}^{2} + d .

The Bayesian information criterion (Schwarz, 1978) was used to adjudicate the fits of the models. These values are summarized in Table 4. Although all the competing models show similar data–model fit, the linear–linear model with smooth transition might well be preferred based on the theoretical notion that increase in verbal reasoning occurs more gradually at the changepoint in verbal acuity, but the rate of change in each phase tends to be constant.

Table 4.

Bayesian information criterion (BIC) Values for Four Competing Latent Regression Models Fitted to the Education Data.

Model	Parameters	BIC
Linear–linear smooth transition	15	7661.4
Linear–linear abrupt transition	15	7665.2
Quadratic–linear	15	7680.2
Conventional quadratic	14	7663.5

Discussion

This research study considers a piecewise regression model for describing the bivariate relation between latent variables where the latent constructs are measured by a set of observed variable indicators. To accommodate this added complexity, the nonlinear regression model is augmented with a measurement model that characterizes the relation between observed and latent variables.

When the knot is unknown (as is true in this current work), then the structural model becomes intrinsically nonlinear which complicates the estimation of model parameters via maximum likelihood. The proposed linear–linear piecewise function is nominally a five-parameter model; however, the parameters are constrained so that the segments join at the knot, resulting in four independent coefficients. Constraining the linear–linear model to make a smooth transition rather than an abrupt change may be reasonable given the continuous nature of the latent predictor. The model was subsequently transformed to incorporate this feature.

A version of maximum likelihood estimation was used here that reduced the number of dimensions of integration to one—a nontrivial savings in computation time. Whatever algorithm is used, however, it must be able to handle the intrinsic nonlinearity of the piecewise function. As Wall (2009) pointed out, conventional estimation methods appropriate for fitting linear structural equation models focus on minimizing a discrepancy function between the observed and modeled covariance matrix and, this cannot be extended in a straightforward way to handle nonlinear structural models of the type explicated in this article.

Spline regression is an attractive alternative to modeling nonlinear data that exhibit distinct phases. Frequently, the model can be tailored so that regression parameters correspond to interesting features of the underlying relation. The knot, an interesting parameter in its own right, is the value of the exogenous predictor at which the bivariate relation changes. In the current example, the knot represents a level of verbal acuity at which verbal reasoning increases, but at a much more rapid rate. This type of valuable information may be used by practitioners to suggest curricular modifications, timing of interventions, or simply broader programmatic changes.

Footnotes

Appendix

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) declared receipt of the following financial support for the research, authorship, and/or publication of this article: The research reported here was funded by a grant from the Institute of Education Sciences, U.S. Department of Education, to the University of Maryland (R305A090152). The opinions expressed are those of the author and do not represent views of the Institute or the U.S. Department of Education.

Notes

References

Bacon

D. W.

Watts

D. G.

(1971). Estimating the transition between two intersecting straight lines. Biometrika, 58, 525-534.

Bolger

Shrout

P. E.

(2007). Accounting for statistical dependency in longitudinal data on dyads. In Little

T. D.

Bovaird

J. A.

Card

N. A.

(Eds.), Modeling ecological and contextual effects in longitudinal studies of human development (pp. 285-298). Mahwah, NJ: Lawrence Erlbaum.

Bollen

K. A.

Curran

P. J.

(2006). Latent curve models: A structural equation approach. New York, NY: Wiley.

Cudeck

Codd

C. L.

(2012). A template for describing individual differences in longitudinal data with application to the connection between learning and ability. In Harring

J. R.

Hancock

G. R.

(Eds.), Advances in longitudinal methods in the social and behavioral sciences (pp. 3-24). Charlotte, NC: Information Age.

Cudeck

Harring

J. R.

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

S.-M.

Ferrer

Hsieh

(Eds.), Statistical methods for modeling human dynamics: An interdisciplinary dialogue (pp. 289-318). New York, NY: Routledge.

Cudeck

Harring

J. R.

du Toit

S. H. C.

(2009). Marginal maximum likelihood estimation of a latent variable model with interaction. Journal of Educational and Behavioral Statistics, 34, 131-144.

Cudeck

Klebe

K. J.

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

Davidian

Gallant

A. R.

(1993). The nonlinear mixed effects model with a smooth random effects density. Biometrika, 80, 475-488.

Davis

P. J.

Rabinowitz

(1984). Methods of numerical integration. New York, NY: Academic Press.

10.

du Toit

S. H. C.

Cudeck

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

11.

Eppright

E. S.

Fox

H. M.

Fryer

B. A.

Lamkin

G. H.

Vivian

V. M.

Fuller

E. S.

(1972). Nutrition of infants and preschool children in the north central region of the United States of America. World Review of Nutrition and Dietetics, 14, 269-332.

12.

Griffiths

D. A.

Miller

A. J.

(1973). Hyperbolic regression-a model based on two-phase piecewise linear regression with a smooth transition between regimes. Communications in Statistics, 2, 561-569.

13.

Harring

J. R.

Cudeck

du Toit

S. H. C.

(2006). Fitting partially nonlinear random coefficient models as SEMs. Multivariate Behavioral Research, 41, 579-596.

14.

Klein

Muthén

(2007). Quasi-maximum likelihood of structural equation modeling with multiple interaction and quadratic effects. Multivariate Behavioral Research, 42, 647-673.

15.

Krommer

A. R.

Ueberhuber

C. W.

Numerical integration on advanced computer systems (Lecture Notes in Computer Science 848). Berlin, Germany: Springer-Verlag; 1994.

16.

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 prostrate cancer. Journal of the American Statistical Association, 90, 44-53.

17.

Naumova

E. N.

Must

Laird

N. M.

(2001). Tutorial in biostatistics: Evaluating the impact of “critical” periods in longitudinal studies of growth using piecewise mixed-effects models. International Journal of Epidemiology, 30, 1332-1341.

18.

Oehlert

G. W.

(1992). A note on the delta method. American Statistician, 46, 27-29.

19.

Patefield

(2002). Fitting non-linear structural relationships using SAS procedure NLMIXED. Journal of the Royal Statistical Society: Series D (The Statistician), 51, 355-366.

20.

Pinheiro

J. C.

Bates

D. M.

(1995). Approximations to the loglikelihood function in the nonlinear mixed effects model. Journal of Computational and Graphical Statistics, 4, 12-35.

21.

Schwarz

G. E.

(1978). Estimating the dimension of a model. Annals of Statistics, 6, 461-464.

22.

Seber

G. A. F.

Wild

C. J.

(1989). Nonlinear regression. New York, NY: Wiley.

23.

Skrondal

Rabe-Hesketh

(2004). Generalized latent variable modeling: Multilevel, longitudinal and structural equation models. Boca Raton, FL: Chapman & Hall/CRC Press.

24.

Wall

M. M.

(2009). Maximum likelihood and Bayesian estimation of nonlinear structural equation models. In Millsap

(Ed.), The SAGE handbook of quantitative methods in psychology (pp. 540-567). Thousand Oaks, CA: Sage.

25.

Wall

M. M.

Amemiya

(2000). Estimation for polynomial structural equation models. Journal of the American Statistical Association, 95, 929-940.

26.

Wall

M. M.

Amemiya

(2007). A review of nonlinear factor analysis and nonlinear structural equation modeling. In Cudeck

MacCallum

R. C.

(Eds.), Factor analysis at 100: Historical developments and future directions (pp. 337-362). Mahwah, NJ: Lawrence Erlbaum.