Longitudinal Data Analyses Using Linear Mixed Models in SPSS: Concepts, Procedures and Illustrations

Although different methods are available for the analyses of longitudinal data, analyses based on generalized linear models (GLM) are criticized as violating the assumption of independence of observations. Alternatively, linear mixed models (LMM) are commonly used to understand changes in human behavior over time. In this paper, the basic concepts surrounding LMM (or hierarchical linear models) are outlined. Although SPSS is a statistical analyses package commonly used by researchers, documentation on LMM procedures in SPSS is not thorough or user friendly. With reference to this limitation, the related procedures for performing analyses based on LMM in SPSS are described. To demonstrate the application of LMM analyses in SPSS, findings based on six waves of data collected in the Project P.A.T.H.S. (Positive Adolescent Training through Holistic Social Programmes) in Hong Kong are presented.


INTRODUCTION
How can we analyze interindividual differences in intraindividual changes over time? Traditionally, researchers used generalized linear models (GLM), such as analysis of variance (ANOVA) and analysis of covariance (ANCOVA), to examine changes in behavior across time. However, these methods would only estimate the model accurately in a balanced, repeated-measures design (e.g., equal group sizes). Unfortunately, this condition is difficult to meet and the use of the traditional univariate and multivariate test statistics might increase Type I errors under the condition of an unbalanced repeated-measures design [1,2,3].
Furthermore, the assumption of independence of observations intrinsic to GLM is not easily met when longitudinal data are under examination. As longitudinal observations may not be truly independent because of a higher-level clustering unit (i.e., time), the data used for analysis will include data that are duplicated so that observations within the clustering unit are correlated. Although it is assumed that each observation contains unique information, this information will not be truly unique, which will eventually result in biased standard errors. While violation of independence of observations is not a must in longitudinal data and there are procedures to diagnose this problem, researchers must figure out ways to deal with this problem when it exists [2,4,5].
Against the above background, there is an increased interest to study the rate of change using individual growth curve (IGC) models. IGC is an advanced technique for modeling within-person systematic change and between-person differences in developmental outcomes across different measurement waves over time. By specifying different sets of models, researchers are able to examine change in the predictive effect when additional variables are added [6]. To determine individual growth profiles and to address the questions of stability over time, researchers call for the measurement of change using this strategy [2,7]. Although the term -individual growth curve‖ is commonly used, it is noteworthy that analyses are usually conducted to examine -aggregates‖ of individual curves, rather than separate analysis of each IGC. Discussion on the use of IGC models has been described by Singer and Willett [3].
Besides capturing developmental changes over time, many researchers advocated the use of IGC when examining the longitudinal pattern of treatment effects over time [1,8,9,10] and a number of advantages of using this method were identified [1,11]. First, it does not require balanced data across different waves of data. This provides researchers with a more flexible and powerful approach when handling unbalanced data (e.g., unequal sample size, inconsistent time interval, and missing data). For example, the number and spacing of measurement occasions may vary (i.e., different points in time for different individuals), instead of being fixed (i.e., regular spaced). This is important in longitudinal studies in which the problems of participant dropout and other forms of missing measurements within individuals are often encountered. This will overcome the limitation of other conventional statistical techniques (e.g., multivariate analysis of variance [MANOVA]) that do not allow for missing data.
Second, it allows researchers to study both intra-and interindividual differences in the growth parameters (e.g., slopes and intercepts). IGC retains all of the information and variability in the data when examining the rate of changes in the dependent variables [12]. This information is valuable in the field of developmental psychology as individuals vary not only in their initial status, but also their rates of changes. Most methods for repeated-measures designs (e.g., multiple regression analyses, ANOVA, MANOVA) only focus on group differences in patterns of change, but variations of growth curve parameters might also exist at the individual level. Understanding the patterns of change and the effects at both the individual and group levels would help researchers to analyze data appropriately and capture a comprehensive picture of developmental changes across time.
Third, IGC analyses estimate the change parameters with greater precision when the number of time waves is increased. This improves the reliability of the growth parameters by reducing standard errors of the within-subject change in the growth parameters estimates [11,13]. This is obviously an advantage when compared with traditional GLM.
Fourth, the effects of predictors at higher levels (e.g., family, classroom, community, etc.) and other predictors on individual growth can flexibly be added in the growth curve models [14]. IGC can be used to explore the causal links between the linkages of predictors and changes in outcome variables across time. In addition, it allows predictors of growth to be discrete or continuous as well as time variant or time invariant. Time-variant predictors refer to independent variables that change over time (e.g., age, weight, height). Time-invariant predictors refer to independent variables that remain constant over time (e.g., gender, ethnicity).
Lastly, IGC is more powerful than other methods (e.g., ANOVA, MANOVA, multiple regression analyses) in examining the effects associated with repeated measures as it models the covariance matrix (i.e., fitting the true covariance structure to the data [15]) rather than imposing a certain type of structure as commonly used in traditional univariate and multivariate approaches [16]. In particular, the error covariance structure of the repeated measurement can be specified in IGC models, and thus allow researchers to examine true change and possible determinants of this structure during hypothesis testing. By choosing an appropriate covariance structure for the growth curve model, error variance would be reduced and allow researchers to specify a correct model that conceptualizes the patterns of change over time.
When the search term -individual growth curve‖ was used in September 2010, there were 260 citations in PsycINFO and 11 citations in Social Work Abstracts. When the term -growth curve modeling‖ was used, there were 633 and 17 citations in PsycINFO and Social Work Abstracts, respectively. These figures clearly show that there is a strong need to conduct studies on IGC modeling in the social work research context. The paucity of this analytic tool research might be related to the lack of technical papers illustrating the practical use of growth curve modeling via SPSS. There are two reasons why we document the use of linear mixed methods (LMM) in SPSS. First, SPSS is popular software used by researchers in different disciplines. As such, many researchers would like to use SPSS to perform LMM instead of using additional software. Second, there are few publications illustrating how researchers can use SPSS to analyze longitudinal data in an experimental design. As such, an illustration of how to use SPSS to analyze longitudinal intervention research would be beneficial to researchers.
The purpose of this paper is to demonstrate the use of IGC in the analyses of longitudinal data using SPSS. The general strategy for model building, testing, and comparison are described. Previous studies have illustrated the application of IGC using PROC MIXED in SAS [16,17,18], HLM [19], R [20], and SPSS [21]. Nevertheless, the longitudinal analysis reported in Peugh and Enders [21] was only a simple example not conducted within an intervention context. Furthermore, as Francis et al. [1] pointed out, -more number of time points necessitated the use of polynomial models for the individual trajectories‖ (p. 36). As such, the pattern of change across six time points after participating in a positive youth development program in comparison to a control group was examined in the present study. By modeling the longitudinal data, the IGC method is described and SPSS commands and outputs are examined.

LONGITUDINAL DATA SET
The data for this study were part of a multiyear positive youth development program. . The majority of missing data were the result of participant absence at the day of data collection rather than attrition from the study. The number of collected questionnaires was 7,846 in Wave 1; 7,388 in Wave 2; 6,939 in Wave 3; 6,697 in Wave 4; 6,876 in Wave 5; and 6,733 in Wave 6. The number of successfully matched responses of the overall sample was 98% in Wave 1, 96% in Wave 2, 97% in Wave 3, 98% in Wave 4, 99% in Wave 5, and 97% in Wave 6. Participants that completed all six waves were 4,712 (i.e., 60% of the sample). Details are shown in Table 1.
At different measurement points, participants were required to respond to an objective outcome questionnaire, which included the Chinese Positive Youth Development Scale (CPYDS) [22]. For the purpose of illustration, a variable based on the 36 most important items (i.e., KEY 36 indicator) of the scale was focused upon (Table 2). Previous studies [22,23] showed that the CPYDS is a valid and reliable instrument.

DATA ANALYTIC STRATEGY
The data were analyzed by using a mixed effect model with maximum likelihood (ML) estimation [24]. This method modeled individual change over time, determined the shape of the growth curves, explored systematic differences in change, and examined the effects of covariates (e.g., treatment) on group differences in the initial status and the rate of growth. It is an appropriate approach when studying individual change as it creates a two-level hierarchical model that nests time within individual [14,25].
The basic assumption of IGC is that the functional form of each individual trajectory is similar (e.g., linear growth over time in the whole sample, see Willett et al. [7]). To capture individual change over time,  each individual trajectory is summarized by fitting to a specific form of parametric model (i.e., regressing observed record into the average of the trajectories [26]). Generally speaking, ordinary least squares (OLS) regression is used to meet this exploratory purpose [3]. Interindividual differences in growth trajectories may be found in the individual growth parameters, such as intercepts (i.e., initial status) and slopes (i.e., steep or flat). Researchers are interested in whether this heterogeneity in change is systematically related to various contextual variables [3,7]. In other words, by fitting each person's OLS trajectory to a specific parametric model, the individual trajectory in a population was obtained and allowed us to further investigate whether the individual differences in growth parameters were related to other explanatory variables. There are two levels in IGC models. The Level 1 model refers to the within-person or intraindividual change model (i.e., repeated measurements over time). It focuses on the individual and describes the developmental changes for each individual (i.e., the variation within individual over time). The Level 1 model estimates the average within-person initial status and rate of change over time. No predictors are included in this model. The basic linear growth model is as shown below.
Level 1 model: In our example, β 0 is the initial status (i.e., Wave 1) of the KEY 36 indicator for individual i, β 1 is the linear rate of change for individual i, and r ij is the residual in the outcome variable for individual i at Time t. Y ij is the repeatedly measured KEY 36 indicator for an individual i at Time t. If the effect of linear growth (Time, β 1 ) is not statistically significant, there is no need to perform further growth curve modeling analysis. To test a nonlinear individual growth trajectory across time, other higher-order polynomial trends (i.e., quadratic and cubic slopes) can also be included for model testing. This is shown in Eq. 2, in which Time (i.e., the linear slope, β 1 ) remains, while Time 2 (i.e., quadratic slope, β 2 ) and Time 3 (i.e., cubic slope, β 3 ) are added in the model.
The linear slope suggests that the rate of growth remains constant across time (i.e., a straight line, see Fig. 1), whereas the higher-order polynomial trends indicate that the growth rates might not be the same over time. A quadratic (second-order polynomial) individual change trajectory has no constant common slope (i.e., accelerate/decelerate over time) and consists of a single stationary point (i.e., peak/trough) (i.e., a parabola shape, see Fig. 2). A cubic trajectory has two stationary points, with one peak and one trough (i.e., S-shaped, see      The growth parameters (i.e., the within-subjects intercepts and slope) of Level 1 are the outcome variables to be predicted by the between-subjects variables at Level 2. At this level (Eq. 3), an explanatory variable (W j ) is included to analyze the predictor's effect on interindividual variation on outcome variable. The errors are assumed to be independent and normally distributed, and the variance is equal across individuals [27]. The Level 2 model is: In our example, Y ij is the grand mean for the KEY 36 indicator for the whole sample at Time t. γ 0i is the initial status of the KEY 36 indicator for the whole sample at Time t. γ 1i is the linear slope of change relating to the KEY 36 indicator for the whole sample at Time t. γ 2i is the quadratic slope of change relating to the KEY 36 indicator for the whole sample at Time t. γ 3i is the cubic slope of change relating to the KEY 36 indicator for the whole sample at Time t. γ 4i is used to test whether the predictor (e.g., group) is associated with the growth parameters (i.e., initial status, linear growth, quadratic growth, and cubic growth). r ij refers to the random effects (i.e., amount of variance) that are unexplained by the predictor.
In this study, we tested whether treatment was predictive of students' initial status and different trajectory changes in positive youth development across time. A dummy/dichotomous variable was created (i.e., group-control vs. experimental groups) as a predictor. Participants in the control group were coded as -1 and those in the experimental group as 1. Two covariates (i.e., gender and initial age) were included when examining the predictive program effect on the outcome variable. Gender (k2) was coded as -1 = male and 1 = female. A similar coding method for a dichotomous variable was found in previous studies [14,24]. For the continuous variables, a grand mean centering method was generally recommended in order to simplify the interpretation of the results [2]. In our study, the mean age was 12. Initial age (k1) was then centered by subtracting the mean age (i.e., age = k1 -12) and, therefore, the centered initial age (age) was generated. To compute the centered initial age (age), the following syntax was used: Following the strategy suggested by Singer and Willet [3], several models were tested. These included (1) an unconditional model (Model 1) that was tested to examine any mean differences in the outcome variable across individuals, (2) an unconditional growth model (Model 2) that served as a baseline model to explore whether the growth curves are linear or curvilinear, (3) two higher-order polynomial models (Models 3 and 4) that were estimated to determine if the rate of change accelerated or decelerated across time, (4) a conditional model (Model 5) that was formed to investigate whether the predictor was related to the growth parameters (i.e., initial status, linear growth, quadratic growth, and cubic growth), and (5)  three different covariance structure models (Models 6, 7, and 8) that were generated to assess the error covariance structure of the longitudinal data. The intercept and linear slope were allowed to vary across individuals. Missing data were handled through pairwise/likewise deletion.
To select the best model, -2 log likelihood (i.e., likelihood ratio test/deviance test), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) were used. Generally, the smaller the statistical values, the better the model fit to the data. Analyses were performed using the mixed model procedure in SPSS 17.0 statistical software [28]. In SPSS, the restricted maximum likelihood method (REML) is the default option for model estimation. As we focused on the entire model (both fixed and random effects), the maximum likelihood (ML) method was used [3]. The KEY 36 indicator (Key) was the outcome variable. A high score of this variable suggested better positive youth development.

PREPARATION FOR SPSS ANALYSES
Before performing IGC analysis, a -person-period data, one record for each period‖ (univariate format) set is required [3]. Based on this dataset, each subject's temporary sequenced outcome values were recoded vertically. For example, Subject ID 1234 has six records (six waves) and Subject ID 10296 has two records (two waves). The VARSTOCASES statement restructured the dataset into a -multiple-record‖/stacked format [3]. The MAKE statement converted the values of a repeated-measurement variable (i.e., KEY36, SKEY, TKEY36, FKEY36, GKEY36, HKEY36) into a single variable (i.e., key). The INDEX statement created a new variable (i.e., Wave) to specify the time interval of the six measurement occasions (i.e., Wave 1 = KEY36; Wave 2 = SKEY; Wave 3 = TKEY36; Wave 4 = FKEY36; Wave 5 = GKEY36; Wave 6 = HKEY36). To convert the data into a stacked format, the following syntax was used: One of the strengths of IGC is that it allows the irregularity of number and spacing of waves. Singer and Willett [3] highlighted the use of a time-structured predictor (TIME) for analyzing irregularly spaced datasets. In the present study, the measurement occasion was used because every individual was assessed on the same occasions. The pretest (Wave 1) values of time were set at 0, and the number of month from pretest was calculated for each wave of subsequent data collection (i.e., Waves 2-6). The data collection was scheduled at 8, 12, 20, 24, and 32 months after the baseline data collection. Therefore  Table 3. To test a nonlinear developmental trend over the measurement period, higher-order parameters (i.e., Time 2 , Time 3 ) were included. Six waves of data are sufficient for a precise measurement of nonlinear individual trajectories change over time [5]. The quadratic time was formed by squaring the linear term (i.e., TIME 2 = 0 at Wave 1; TIME 2 = 0.45 at Wave 2; TIME 2 = 1 at Wave 3; TIME 2 = 2.79 at Wave 4; TIME 2 = 4 at Wave 5; TIME 2 = 7.13 at Wave 6). The cubic time was also calculated by powering the linear term to three (i.e., TIME 3 = 0 at Wave 1; TIME 3 = 0.30 at Wave 2; TIME 3 = 1 at Wave 3; TIME 3 = 4.66 at Wave 4; TIME 3 = 8 at Wave 5; TIME 3 = 19.03 at Wave 6).
To compute linear function of change Time variable, the following syntax was used: The quadratic function of time (i.e., -Time_sq‖) was computed by the following syntax: The cubic function of time (i.e., -Time_cub‖) was computed by the following syntax: COMPUTE Time_cub=Time*Time*Time. EXECUTE.

STEPS IN SPSS LINEAR MIXED MODEL ANALYSIS
The six waves of data from the Project P.A.T.H.S. (Positive Adolescent Training through Holistic Social Programmes) were used. The Project P.A.T.H.S. is a large-scale positive youth development program designed for junior secondary school students (Secondary 1 to 3, i.e., Grades 7 to 9) in Hong Kong. The details of the program are reported elsewhere [28,29,30,31].

Step 1: Unconditional Mean Model (Model 1)
This is a one-way ANOVA model with a random effect. In this model, no predictor is included. It serves as a baseline model to examine individual variation in the outcome variable without regard to time [3]. This model assesses (1) the mean of the outcome variable and (2) the amount of outcome variation that exists in intra-and interindividual levels. This latter information is important as it helps determine which level (i.e., Level 1, time variant; or Level 2, time invariant) of predictors to add when fitting the subsequent models. If the variation is high, it suggests that certain amount of outcome variation could be explained by the predictors at that level.
One of the strengths of IGC is that it examines the proportion of total outcome variation that is related to interindividual differences (i.e., intraclass correlation coefficient [ICC]). ICC describes the amount of variance in the outcome that is attributed to differences between individuals. It evaluates the necessity of modeling the nested data structure (i.e., any significant variation in individual initial status of the outcome variable). It is also a measure of the average autocorrelation of the outcome variable over time [3]. The higher value indicates the estimated average stability of the dependent variable over time.
To test the unconditional mean model, the following syntax was used: mixed key /fixed intercept /random intercept | subject(id) covtype(un) /print solution testcov /method ml. Lists the random-effect variables (i.e., intercept). Specifies the classification variable (i.e., ID) and the error covariance structure type (i.e., UN).
Requests the printed output with specific results (i.e., fixedeffect estimates, its standard errors, a t-test for the parameter, significance tests for the estimated variance components). Specifies the use of estimation method (i.e., ML).
The MIXED statement requests the procedure. The FIXED and RANDOM statements list the fixed and random effect variables (i.e., intercept), respectively. The SUBJECT statement specifies that ID is a classification variable to indicate that the data represent multiple observations over time for individuals. The COVTYPE statement captures the error covariance structure for data analysis. The unstructured (UN) covariance matrix for the random effects was tested (note: other covariance structure models were tested in a later section). The PRINT SOLUTION statement requested the printed output of fixed-effect estimates, its standard errors, and a t-test for the parameter. The TESTCOV is used to conduct significance tests for the estimated variance components. Maximum likelihood (ML) was used to estimate the model.  44), suggesting that about 67% of the total variation in the KEY 36 indicator was due to interindividual differences. In other words, the estimated average stability of the KEY36 indicator was 0.67. ICC can be used to help researchers be aware of possible mediating/moderating effects on outcome variables. If ICC is low, IGC might not perform better than the traditional method (e.g., ANOVA) in estimating fixed effects [32]. Generally, IGC is required if ICC is 0.25 or above [33,34]. The full SPSS output can be seen in Appendix A.

Step 2: Unconditional Linear Growth Curve Model (Model 2)
This is a baseline growth curve model that examines individual variation of the growth rates (i.e., any significant variations in individual trajectory changes over time). Unlike the unconditional mean model, which only assesses the outcome variation across individuals (i.e., the differences between the observed mean value of each person and the true mean from the population), this model also examines individual changes over time (i.e., how each person's rate of change deviates from the true rate of change of the population) [3]. If there is no interindividual difference in trajectory change over time (i.e., Time is not statistically significant), further model testing would not be performed.
To test the unconditional growth curve model, the following syntax was used: mixed key with Time /fixed intercept Time /random intercept Time | subject(id) covtype(un) /print solution testcov /method ml. Lists the random-effect variables (i.e., intercept, Time). Specifies the classification variable (i.e., ID) and the error covariance structure type (i.e., UN). 4 /print solution testcov /method ml.
Requests the printed output with specific results (i.e., fixed-effect estimates, its standard errors, a t-test for the parameter, significance tests for the estimated variance components). Specifies the use of estimation method (i.e., ML).
In this model, WITH was specified in the MIXED statement, and Time was added in the MIXED and FIXED statements to test the linear growth of the KEY 36 indicator over time. Furthermore, linear slopes were allowed to randomly vary across individuals by listing Time in the RANDOM statement. The significant values in both the intercept and linear slope parameters indicate that the initial status and linear growth rate were not constant over time. There was a significant linear decrease in the KEY 36 indicator scores (β = -0.33, SE = 0.12, p < 0.01). The mean estimated initial status and linear growth rate for the sample were 155.91 and -0.33, respectively. This suggested that the mean KEY 36 indicator was 155.91 and decreased with time. The random error terms associated with the intercept and linear effect were significant (p < 0.01), suggesting that the variability in these parameters could be explained by between-individual predictors.

Unconditional linear growth model (degrees of freedom=6)
Comparing within-individual variation in initial status between Model 1 and Model 2, there was a decline in the residual variance of 36.63 (228.74 to 192.11). This suggested that about 37% of the withinindividual variation in the KEY 36 indicator was associated with linear rate of change.
The correlation (β = -55.62, SE = 3.64, p < 0.01) between the intercept and the linear growth parameter was negative. This suggests that students with high KEY 36 indicator scores had a slower linear decrease, whereas students with low KEY 36 indicator scores had a faster decrease in linear growth over time.

Step 3: Quadratic Growth Curve Model (Higher-Order Change Trajectories) (Model 3)
Individual growth trajectories are usually nonlinear over time as shown in previous developmental studies [35,36]. As such, two higher-order polynomial models were tested. The analyses examined whether the rate of growth accelerated or decelerated over time. To test the quadratic rate of change, a model with quadratic time (Time_sq) was examined by adding quadratic parameter in the previous model.
To test the quadratic effect, the following syntax was used: mixed key with Time Time_sq /fixed intercept Time Time_sq /random intercept time | subject(id) covtype(un) /print solution testcov /method ml. Lists the random-effect variables (i.e., intercept, Time). Specifies the classification variable (i.e., ID) and the error covariance structure type (i.e., UN). 4 /print solution testcov /method ml. Requests the printed output with specific results (i.e., fixed-effect estimates, its standard errors, a t-test for the parameter, significance tests for the estimated variance components). Specifies the use of estimation method (i.e., ML). Results showed that all growth parameters were significant (p < 0.01), indicating that there were significant between-subjects variations in the initial status, and linear and quadratic time trajectories (i.e., reliably different from zero). The initial status (grand mean score at Wave 1) of the KEY

Step 4: Cubic Growth Curve Model (Higher-Order Change Trajectories) (Model 4)
Researchers noted that more number of time points was required when testing different types of polynomial models for the individual trajectories, even though it is difficult to interpret such a complex model [3]. A cubic model was also tested to summarize individual change for the whole sample. The syntax of this model was similar to the previous one, except with the inclusion of Time_cub in the FIXED statement. The purpose of this model is to test any cubic changes in individual trajectories over time (i.e., examine whether another nonlinear growth model fits the data better).
The following syntax was used: mixed key with Time Time_sq Time_cub /fixed intercept Time Time_sq Time_cub /random intercept time | subject(id) covtype(un) /print solution testcov /method ml.

mixed key with Time Time_sq Time_cub
Requests the mixed-level analysis procedure.
Requests the printed output with specific results (i.e., fixed-effect estimates, its standard errors, a t-test for the parameter, significance tests for the estimated variance components). Specifies the use of estimation method (i.e., ML). Time, Time_sq, and Time_cub had a significant contribution in the model (p < 0.01). The negative effect of linear growth (β = -9.46, SE = 0.63, p < 0.01) suggested that the KEY 36 indicator decreased at the beginning. The positive effect of quadratic growth (β = 6.26, SE = 0.57, p < 0.01) indicated a deceleration in the rate of change (i.e., initially decreased and then began to increase). However, the negative effect of cubic growth (β = -1.12, SE = 0.14, p < 0.01) revealed that such deceleration gradually diminished over time. Given that the cubic model improved model fit over the previous model (χ 2

Step 5: Adding Predictors (Model 5)
To test the predictor effect on the shape of individual growth trajectories, a dichotomous variable -group‖ was examined as a time-invariant covariate to explore any group differences in change over time (i.e., interaction with time). It examined whether group was a predictor of the intercept, linear, quadratic, and cubic parameters. In particular, the relationships between the KEY 36 indicator and group were estimated after controlling the effect of gender -k2‖ and initial age -age‖.
The following syntax was used: mixed key with Time Time_sq Time_cub group k2 age /fixed Time Time_sq Time_cub group k2 age group*Time group*Time_sq group*Time_cub k2*Time k2*Time_sq k2*Time_cub age*Time age*Time_sq age*Time_cub /random intercept time | subject(id) covtype(un) /print solution testcov /method ml. Group was a significant predictor of the linear, quadratic, and cubic changes in the KEY36 indicator (p < 0.01), but not associated with the initial status (β = 0.06, SE = 0.36, p < 0.01). Regarding the linear slope of the KEY 36 indicator, the control group showed a faster rate of change as compared with the experimental group (β = 2.92, t = 4.16, p < 0.01). In terms of quadratic growth, the control group had a slower rate of change in the KEY 36 indicator when compared with the experimental group (β = -2.38, t = -3.65, p < 0.01). Lastly, the control group had a faster rate of cubic change than the experimental group (β = 0.53, t = 3.31, p < 0.01). In other words, a stable trajectory of the KEY 36 indicator was found in the experimental group, but not in the control group. The predictor (group) accounted for 7% [(192.11 -179.36) / 192.11 = 0.066] of the within-individual variations in the KEY 36 indicator. This shows that only 7% of the overall variability in the KEY36 indicator is explained by Group.
Singer and Willett [3] proposed using prototypical values to demonstrate the effect of treatment on initial status and the rate of change across time. The step in creating prototypical plots is generally identical to the method of plotting graphs in regression [37]. We can obtain the fitted trajectories by . This method was also used in previous studies [38,39]. The trajectories of the control and experimental groups are shown in Fig. 4. In general, the experimental group had a steady growth rate of change in the KEY 36 indicator compared to the control group.

Step 6: Examining Covariance Structure
One of the advantages of IGC is the availability to specify the within-individual error covariance structure that best fits the data. The purpose of testing different error covariance matrices is to describe how the error is distributed [18]. It examines whether the properties imposed on the error covariance structure of the parametric model fit well to the data [3]. This is very important when we examine unequally spaced and unbalanced data, which are commonly found in longitudinal studies.
In fact, studies showed that the estimated variances of the parameter estimates are likely to be biased and inconsistent when repeated measurements are taken on the same individual across time (i.e., failure to take account for heteroscedasticity) [40,41] and consequently affect the precision of estimating the appropriate model [42]. Researchers advocated the use of this variance-covariance testing approach as it improves model predictions and statistic inferences, especially when examining random effects models [43,44]. In the present study, three types of covariance structures (i.e., unstructured, compound symmetric, and first-order autoregressive) that were commonly examined in previous studies were tested [18,45,46,47].

Unstructured (UN) Covariance Structure (Model 6)
The unstructured covariance structure model often offers the best fit and is most commonly found in longitudinal data as it is the most parsimonious, which requires no assumption in the error structure [18]. In this model, the syntax is the same as the previous one, except for the inclusion of the REPEATED statement, which substitutes the RANDOM statement. The REPEATED statement lists -wave‖, which specified the number of each repeated measurement (e.g., 1, 2, 3….

Compound Symmetry (CS) Covariance Structure (Model 7)
To examine whether the variance and correlation between each pair of observations are constant across time points, a compound symmetry (CS) covariance structure was tested. Therefore, the specification of CS was added in the COVTYPE statement.

First-Order Autoregressive (AR1) Covariance Structure (Model 8)
In this model, the variance is assumed to be heterogeneous and the correlations between the two adjacent time points decline across measurement occasions. The AR1 was specified in the COVTYPE statement.
The following syntax was used: mixed key with Time Time_sq Time_cub group k2 age /fixed Time Time_sq Time_cub group k2 age group*Time group*Time_sq group*Time_cub k2*Time k2*Time_sq k2*Time_cub age*Time age*Time_sq age*Time_cub /repeated wave | subject(id) covtype(AR1) /print solution testcov /method ml. Based on Table 4, it was observed that the smallest values in the three fit criterion (i.e., -2 times the log-likelihood [-2LL], Akaike's Information Criteriaon [AIC], and the Bayesian Information Criterion [BIC]) were found in the UN model. These differences were statistically significant when compared to the results of a chi-square distribution with degrees of freedom (i.e., 37 parameters-18 parameters, Δdf = 19, p < 0.01). This suggested that the UN model was the best model in fitting the data, although it required many parameters. The correlated error terms and heterogeneous variances might be the result of unequally spaced time points of measurement. If the time points were closely spaced, the possibility of modeling correlated errors might be higher than those that were scheduled far apart [2]. The use of this variance-covariance approach would improve model predictions. Readers interested in different error covariance structures can read the work of others for details [47,48,49,50].

DISCUSSION AND CONCLUSIONS
The present study demonstrates the application of IGC analyses using SPSS. We first describe the basic growth curve modeling framework and demonstrate how various growth curve models fit to empirical multiwave data via SPSS. To explore nonlinear changes with longitudinal time-structured data, we used IGC to examine the intra-and interindividual differences in longitudinal trajectories over time. Lastly, we examined the effectiveness of the program by comparing differences in the longitudinal patterns of positive youth development between the control and experimental groups. As the SPSS manuals on IGC do not have many examples in the intervention context, the present paper is a significant contribution to the literature. With reference to the lack of IGC studies in the social work literature, the present paper is a pioneer contribution to the field. Furthermore, it is noteworthy that very few longitudinal intervention studies have been conducted by employing this advanced technique in different Chinese contexts.
One of the strengths of IGC analyses is that the numbers of observations collected on each individual can maximize the degree to estimate a complex nonlinear growth curve model. The high number of time points allows us to model different types of polynomial growth curve models [51]. Furthermore, the power of the test is greatly improved by adding only a few additional waves of data collection [26,52]. Many researchers have commented on the need for conducting more developmental studies examining individual growth by using appropriate statistical methods, and using nonlinear growth curves to describe between-and within-individual change over time [1,9,10,16]. In particular, using a developmental approach in understanding the dynamic process of psychosocial development and risk-taking behaviors during adolescence may be useful for designing successful prevention programs [8,53,54,55]. Cleary, this study is a positive response and attempts to fill this research gap. We hope that this method will be more accessible to researchers in the field of social work and other allied professions.