Best Practices for Developing and Validating Scales for Health, Social, and Behavioral Research: A Primer

Scale development and validation are critical to much of the work in the health, social, and behavioral sciences. However, the constellation of techniques required for scale development and evaluation can be onerous, jargon-filled, unfamiliar, and resource-intensive. Further, it is often not a part of graduate training. Therefore, our goal was to concisely review the process of scale development in as straightforward a manner as possible, both to facilitate the development of new, valid, and reliable scales, and to help improve existing ones. To do this, we have created a primer for best practices for scale development in measuring complex phenomena. This is not a systematic review, but rather the amalgamation of technical literature and lessons learned from our experiences spent creating or adapting a number of scales over the past several decades. We identified three phases that span nine steps. In the first phase, items are generated and the validity of their content is assessed. In the second phase, the scale is constructed. Steps in scale construction include pre-testing the questions, administering the survey, reducing the number of items, and understanding how many factors the scale captures. In the third phase, scale evaluation, the number of dimensions is tested, reliability is tested, and validity is assessed. We have also added examples of best practices to each step. In sum, this primer will equip both scientists and practitioners to understand the ontology and methodology of scale development and validation, thereby facilitating the advancement of our understanding of a range of health, social, and behavioral outcomes.

section_content: 1 Introduction Muthén and Asparouhov (2012a) introduced a special type of Bayesian structural equation modeling that provides a more flexible confirmatory factor analysis (CFA) and structural equation model (SEM) analysis. Instead of translating hypotheses into fixed zero parameters which are characteristic of such analyses, the authors proposed the use of approximately zero parameters using zero-mean, small-variance informative priors via Bayesian analysis. The approach is referred to as BSEM and is mainly intended for cross-loadings and residual correlations in CFA and for direct effects representing measurement non-invariance in CFA with covariates, also referred to as MIMIC modeling. An application of BSEM to competing factor analysis models for the Wechsler Intelligence Scale is given in Golay et al. (2012). The current paper generalizes the BSEM approach to the analysis of measurement invariance across several groups or several time points, applying the zero-mean, small-variance prior idea to differences in measurement parameters, thereby introducing the concept of approximate measurement invariance. 

Measurement invariance CFA using maximum-likelihood estimation often takes as a baseline model the specification of fully-invariant measurement parameters and then uses modification indices as guidance to improve the model by freeing parameters that strongly violate invariance. This procedure is suitable for cases with only a few groups/time points, typically only two, and where only a small number of measurement parameters are non-invariant. Modification indices have not been used in item response analysis, presumable for computational reasons due to the numerical integration necessary for maximum-likelihood estimation with categorical variables. For categorical variables, maximumlikelihood estimation with likelihood ratio chi-square testing of item parameters across groups was discussed in Thissen et al. (1993). In this context, backward and forward maximum-likelihood approaches have been studied in (Kim & Yoon, 2011; Lee et al., 2010; Stark et al., 2006). The backward approach starts with a fully invariant model and then releases the invariance specification for one item at a time, producing a likelihood-ratio chi-square test of invariance for each item. The forward approach starts with a fully non-invariant model, comparing it to models where one item at a time is held invariant. Depending on the type of application, non-invariant items are either deleted or their invariance specification relaxed. Cai et al. (2011; pp. 239-243) discussed a multi-step Wald procedure with a designated set of invariant items. De Jong et al. (2007) and Fox (2010) proposed a novel random parameter approach to non-invariance where a measurement parameter for different groups is seen as randomly drawn from a population with shared mean and variance and where a significant variance corresponds to non-invariance. Bou and Satorra (2010) criticize the assumption of random draws in this approach in favor of a multiple-group approach. 

In many applications, measurement invariance needs to be studied for many groups or many time points where there is possibly a large number of non-invariant item parameters. The case of many groups/time points is the focus of this paper. Typical applications are country comparisons such as with the achievement studies of PISA (Program for International Student Assessment) conducted by the Organization for Economic Co-operation and Development, with cross-cultural studies such as the European Social Survey (see, e. g., Davidov et al., 2011), as well as with growth modeling with many time points. In these cases, an item parameter may be invariant across some groups/time points and non-invariant across some other groups/time points. Here, proposed backward and forward maximumlikelihood approaches and multi-step Wald procedures may be too cumbersome to be practical. Furthermore, although maximum-likelihood modification index procedures can be developed for categorical items, many applications involve highdimensional latent variable models where the computations would be too heavy. 

CFA is characterized by using strict parameter constraints. In single-group CFA these are fixed zero loadings corresponding to items hypothesized to not load on certain factors (so called cross-loadings). In multiple-group and multipletime point CFA these are equalities of loadings and intercepts across groups and time point. As illustrated in the top part of Figure 1, from a Bayesian perspective the zero cross-loading case can be seen as CFA using a very strong prior with mean zero and zero variance for certain factor loadings, here represented by a parameter λ. The bottom part of Figure 1 shows the more flexible BSEM approach of using a zero-mean, small-variance prior for the parameter as discussed in Muthén and Asparouhov (2012a). In multiple-group and multiple-time point modeling considered in this paper, the parameter in question concerns the difference in a certain measurement parameter across two groups or time points. This BSEM approach implies the specification of approximate measurement invariance. 

Appendix, Section 9. 1 presents the Mplus specification of such zero-mean, smallvariance priors for differences between parameters. 

In this paper it is shown that the BSEM approach provides a tool for detecting non-invariance that serves the same purpose as modification indices with maximum-likelihood estimation, but unlike maximum-likelihood is applicable also for high-dimensional latent variable models for categorical variables. Under certain forms of non-invariance, proper comparisons of factor means and variances are obtained using only approximate measurement invariance and without relaxing the invariance specifications or deleting non-invariant items. To ensure correct estimation, a two-step Bayesian analysis procedure is proposed, where step 1 uses BSEM to identify non-invariant parameters and step 2 frees those parameters. 

Section 2 discusses an alignment issue which reflects an indeterminacy in measurement invariance analysis. Section 3 presents the BSEM approach for detecting non-invariance. Section 4 discusses a simulation study of a multiplegroup factor analysis model. Section 5 present an application with binary items measuring an achievement factor in 40 countries. Section 6 discusses a simulation study of a multiple indicator growth model. Section 7 presents a binary multiple indicator growth model application for eight time points. Section 8 concludes.

section_content: As pointed out in Asparouhov and Muthén (2012b), the BSEM method for multiple-group, multiple-time point analysis, as well as other methods for such settings, faces a parameterization indeterminacy, or put differently an alignment issue. Briefly stated, the approximate invariance of BSEM essentially acts to find a solution where the variance across groups or time points for a measurement parameter is small. This is not the same simplicity criterion as for example seeking a solution that has many invariant parameters and few non-invariant parameters. 

The latter simplicity criterion may be preferred by researchers who hypothesize a large degree of invariance, while the BSEM criterion may be preferred by researchers who hypothesize that there is a large degree of minor non-invariance where different non-invariance may be in opposite directions and largely cancel each other out. 

As a hypothetical example of the BSEM approach applied to the case with few non-invariant parameters, consider an item where there is invariance of a measurement intercept or threshold across all time points except for one time point where there is a large positive deviation. The BSEM small-variance prior for the parameter differences tends to pull the deviating parameter towards the average of the parameters for all the time points. This means that the deviating parameter will be smaller and the invariant parameters larger than their true values. The same thing occurs for loadings. With intercepts/thresholds and loadings misestimated, the factor means and factor variances are misestimated. This is the alignment issue. As Asparouhov and Muthén (2012b) point out, it does not mean that the model does not fit the data, but that an equally well-fitting solution with a possibly simpler interpretation due to another simplicity criterion may be available. 

The alignment issue should be kept in mind in Monte Carlo simulation studies of the proposed approach. The BSEM analysis for multiple groups or time points is not expected to always recover parameter values used to simulate the data. Perfect recovery will only occur when the non-invariance is designed to be in line with BSEM, such as random, normally-distributed deviations from a parameter average over groups/time points. When non-invariant items have been detected and equalities relaxed, however, the parameter values will be recovered. 

Given that in the above hypothetical example there is a majority of invariant time points, the underestimation of the non-invariant parameter will be smaller than the overestimation of the invariant parameters. Only the non-invariant parameter is therefore likely to be detected as significantly different from the average over time points. Such parameters can then be freed from their invariance restriction in a subsequent analysis that gives correct results. This procedure is discussed in the next section.
section_name: The alignment issue

section_content: The BSEM analysis is augmented with estimation of the difference between each measurement parameter and its average across the groups or time points. Assuming that the parameter is approximately invariant for most of the groups/time points, these differences can point to the groups/time points that have significant non-invariance. In a follow-up analysis, the equality constraint for the parameter can then be relaxed for those groups/time points. This follow-up analysis would use Bayesian analysis but not specify BSEM approximate measurement invariance for any of the parameters but instead hold all but the non-invariant parameters exactly equal across groups/time points. This two-step approach serves the same purpose as working with modification indices with maximum-likelihood estimation. Instead of the modification index approach of relaxing one equality restriction at a time, however, it is here proposed that all significantly mis-fitting equalities are relaxed at the same time. One consideration is the risk of relaxing too many equality restrictions. This should not, however, have much effect on the point estimates, but only increase standard errors to a small extent. With too many equality restrictions relaxed, a non-identification problem arises, but this should be clearly seen in the Bayes plots (for Bayes non-identification detection, see the rejoinder discussion in Muthén & Asparouhov, 2012a). 

The parameter differences across groups/time points are printed in an Mplus output section labeled DIFFERENCE OUTPUT. In Monte Carlo simulations, the output also contains information in parentheses about the proportion of replications for which the difference is significant, providing an estimate of the power to detect the non-invariance. With only two groups/time points, the difference relative to the average can be augmented by the difference across the two groups/timepoints which can be expressed in MODEL CONSTRAINT. 

The choice of prior variance for the differences between parameters is important for detecting non-invariance. As the prior variance is increased, the non-invariance of a parameter is allowed to be more freely estimated, that is, the estimate can escape from the invariance value to a larger degree. At the same time, the standard error of the parameter increases as the prior variance is increased (see also Muthén & Asparouhov, 2012a). The significance of a parameter difference considers a ratio where the numerator is the invariance value and the denominator is the standard error. The numerator increases with increased prior variance and the denominator increases with increased prior variance so that it is difficult to predict how significance is going to be affected. A balance needs to be found for the best prior variance choice. Simulation studies which include these choices are presented in the next section and in Section 6.
section_name: Detecting non-invariance

section_content: Consider a situation with continuous items y measuring a single factor in different groups, using the factor model for individual i, variable j, and group g,

where the residual is uncorrelated with the factor η. The measurement parameters are the intercepts ν jg and the loadings λ jg. The focus is on correctly estimating the factor means α g and factor variances ψ g. The factor metric is set as mean zero and variance one in one group. 

The Appendix Section 9. 2 shows the Mplus Monte Carlo simulation input using six items, 10 groups, and a small degree of non-invariance for seven of the 60 intercepts and for seven of the 60 loadings. Only item 1 is invariant across all groups. Groups 1 and 10 have no non-invariant items. The intercept noninvariance is an increase of 0. 2 in group 2 for item 6, in group 3 for item 5, in group 4 for item 3, in group 5 for item 2, in group 7 for item 6, in group 8 for item 5, and in group 9 for item 3. The loading non-invariance is an increase of 0. 2 in group 2 for item 3, in group 3 for item 2, in group 4 for item 6, in group 5 for item 5, in group 7 for item 3, in group 8 for item 2, and in group 9 for item 5. The intercept non-invariance corresponds to an average of about 0. 2 standard deviations of the items and the loading non-invariance about 0. 15 in a standardized metric. The small number of non-invariant items and the small magnitude of non-invariance correspond to a well-developed measurement instrument. 

As a devil's advocate approach, the settings for the simulation are chosen as the least favorable for the proposed two-step BSEM procedure in two regards. 

First, there is limited non-invariance so that an approach that simply ignores the non-invariance is a strong contender. Second, the form of non-invariance is not of the type that is assumed for BSEM approximate invariance. The non-invariance is in only one direction with higher values than the average over time points in each non-invariant instance for both loadings and intercepts. This implies that the alignment will not be perfect and the population values used to generate the data will not be perfectly recovered. It is of interest to see to which degree recovery occurs and what the coverage values are. A determinant of this is the power to detect non-invariance and the study of how well an analysis with relaxed invariance for non-invariant items performs. 

As a first step, three BSEM analyses are done to study the power to detect non-invariance. Prior variances for the differences are set at 0. 01, 0. 05, and 0. 10, respectively. Using the three different prior variances, Table 1, Table 3, and Table 4 show power estimates in terms of the percentage of replications that show significant differences from the mean across groups, that is, the proportion of replications where the 95% credibility interval for the difference does not include zero. Power estimates for the non-invariant parameters are bolded. In Table 1 the prior variance 0. 01 is used and all 14 non-invariant parameters have power estimates close to 100%. The largest power estimate for an invariant parameter is 52% for the loading of item 3 in group 3, showing that it is possible to make an incorrect diagnosis of non-invariance. Most power estimates for invariant parameter are, however, much smaller. It is therefore likely that the correct non-invariant parameters are identified. 

Table 2 shows the Mplus output in the DIFFERENCE section for the intercept of the second item using prior variance 0. 01. The fifth group has a difference estimate of 0. 140 with power estimate 1. 00. The other nine groups show smaller negative differences, so that all differences sum to zero. 

Loadings Item 1 00 00 00 00 00 00 00 00 00 00 Item 2 00 00 94 00 00 00 00 100 01 00 Item 3 02 94 12 09 00 00 98 00 01 17 Item 4 00 00 00 00 00 00 00 00 00 01 Item 5 01 00 00 00 78 00 00 02 78 00 Item 6 00 03 05 98 00 00 00 00 00 06
section_name: A multiple-group factor analysis simulation

section_content: Item 1 00 00 00 00 00 00 00 00 00 00 Item 2 00 00 00 00 86 00 00 01 00 00 Item 3 02 00 04 96 03 00 00 05 98 03 Item 4 00 00 00 00 00 00 00 00 00 00 Item 5 00 00 73 00 03 00 00 53 00 00 Item 6 00 89 09 02 03 01 92 02 02 00

Table 4 : Multiple-group factor analysis simulation: Power estimates with prior variance 0. 10

Loadings Item 1 00 00 00 00 00 00 00 00 00 Item 2 00 00 62 00 00 00 00 96 01 Item 3 00 75 05 01 00 00 80 00 00 Item 4 00 00 00 00 00 00 00 00 00 Item 5 00 00 00 00 59 00 00 00 56 Item 6 00 00 01 89 00 00 00 00 00
section_name: Intercepts

section_content: Item 1 00 00 00 00 00 00 00 00 00 Item 2 00 00 00 00 47 00 00 00 00 Item 3 01 00 01 79 00 00 00 00 80 Item 4 00 00 00 00 00 00 00 00 00 Item 5 00 00 30 00 00 00 00 15 00 Item 6 00 74 03 00 00 00 76 00 00 parameters free. This analysis is expected to perform optimally given that this is how the data are generated, assuming that the correct non-invariance is diagnosed in analysis 2

Analysis 5 is the proposed two-step procedure. A sample size of n = 500 for each group (total sample size 5,000) and 100 replications are used in the five simulation analyses. 

Table 5 shows the factor mean and variance estimates and coverage using the five approaches. The results for the 10th group are not shown because this is the metric-setting group with fixed values. The last column gives the average absolute bias over the first nine groups. 

Analysis 1 specifies exact invariance across all 10 groups for all loadings and intercepts. The biases are not large in this case. The 95% coverage for the factor variances is, however, too low in four of the groups. 

Analysis 2a shows the BSEM analysis with approximate measurement invariance and prior variance 0. 10. This analysis gives worse bias than analysis 1. due to the alignment issue. Analysis 2b uses BSEM with prior variance 0. 01 which was found to better match the generated non-invariance and to give better power to detect non-invariance. Analysis 2b improves on analysis 2a but still performs worse than analysis 1 in terms of bias. 

Analysis 3 uses item 1 as an anchor item with exact invariance and does perform better than analysis 2b. It performs better than analysis 1 in terms of bias for factor means although not for factor variances. 

Analysis 4 frees the non-invariant parameters and thereby improves on analysis 2b and analysis 3. It performs better than analysis 1. 

Analysis 5 frees the non-invariant parameters and holds other parameter exactly invariant. It is the second step in the proposed approach and performs better than all the other analyses. 

Although not considered here, a more thorough study of analysis 5 can be done where for each replication the parameters with significant difference from the mean in analysis 2 are freed in analysis 5 and the final estimates of this twostep procedure summarized. This takes into account that in some replications the wrong parameters will be diagnosed as non-invariant. 

Other simulation settings with different choices of parameter values are of interest but are not pursued here. For example, data may be generated with many non-invariant parameters and where the non-invariance is in both directions and more in line with BSEM specifications. In such settings, the first step of BSEM analysis may be sufficient and not requiring the freeing of non-invariant parameters. The analyses use Bayesian estimation with a probit link so that for subject i, item j, and group g the probability of the binary variable u ijg can be expressed in terms of the factor η using the standard normal distribution function Φ as

where for each group it is assumed conditional independence among the u's conditional on the factor and normality for the factor. The item response model can be expressed in terms of continuous-normal latent response variable u * underlying the binary outcomes,

with the usual probit standardization of unit variance for the residual. The factor metric can be set by fixing the factor mean in one group at 0 and by fixing either a factor loading at 1 for one item or by fixing the factor residual variance at 1 for one group. For these BSEM analyses, the choice is made to fix the factor mean and variance at 0 and 1 for the 40th group. 

A posterior predictive p-value (PPP) is obtained for the fit to the u * structure as discussed in Asparouhov and Muthén (2010) and Muthén and Asparouhov (2012a). As described in Muthén and Asparouhov (2012a), the PPP is a useful guide when studying the sensitivity of the results to the choice of prior variance. 

If the prior variance is small relative to the magnitude of non-invariance, PPP will be lower than if the prior variance corresponds better to the magnitude of non-invariance. 

For the PISA data a BSEM prior variance of 0. 01 for the measurement parameters results in PPP = 0. 136, a prior variance of 0. 05 results in PPP = 0. 449, a prior variance of 0. 10 results in PPP = 0. 489, and a prior variance of 0. 20 results in PPP = 0. 493. A prior variance 0. 10 implies a prior belief that 95% of the distribution of the non-invariance lies in the range [-0. 62, +0. 62]. 

Using this prior variance the estimated magnitude of non-invariance is less than 0. 55 in absolute value and is therefore within this range. Using prior variance 0. 10 in the BSEM analysis, Table 6 shows the PISA countries with significant differences relative to the average across countries. Relatively few of the loadings and thresholds are non-invariant, only 31 out of the 640 parameters. Several countries have no non-invariant parameters and only countries 18 and 27 have more than one non-invariant parameter. Item 8 shows the least degree of noninvariance. 

This analysis may be compared to that of Fox (2010) using the same data. Fox (2010) used a random parameter approach to non-invariance where a measurement parameter for different groups is seen as randomly drawn from a population with shared mean and variance and where a significant variance corresponds to noninvariance. This technique finds non-invariance for all item parameters except the threshold for item 8. The analysis does not, however, pinpoint which countries contribute to this non-invariance. For example, as Table 6 shows, there are only five countries who have any non-invariant loadings. 

The ordering of factor means across the countries can be compared between the different types of analysis using scatter plots. Figure 2 compares BSEM analysis (X axis) with analysis imposing exact invariance (Y axis). On the whole, the points in the figure align well along a line with only minor deviations, describing a similar ordering of countries. One notable exception is seen in the upper right part of the figure where the analysis under exact invariance ranks the country lower than using BSEM. 

Figure 3 compares BSEM results with those of analysis 5, freeing the noninvariant parameters and imposing exact equality for the other parameters. There is no longer a notable deviation from the line as in Figure 2. Also, the two countries which in Figure 2 have zero factor means according to exact invariance analysis and BSEM analysis are separated along the Figure 3 Y axis using the proposed approach so that the two analyses disagree in this regard. According to the simulations, the ordering according to analysis 5 shown on the Y axis in Figure 3 is the most trustworthy. 

Figure 2 : Estimated factor means for 40 countries: Comparing BSEM analysis (X axis) with analysis imposing exact invariance (Y axis) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1. Comparing BSEM analysis (X axis) with analysis that frees non-invariant parameters and imposes exact equality for other parameters (Y axis) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1. 5 -1. 0 -0. 5 0. 0 0. 5 1. 0 1. 5 The growth model parameterization holds the thresholds equal across time so that the intercept growth factor mean is zero. The focus is on correctly estimating the slope growth factor mean as well as the variances of the intercept and slope growth factors. 

The analyses use Bayesian estimation with a probit link so that for subject i, item j, and time point t, the probability of the binary variable u ijt can be expressed in terms of the factor η, the intercept growth factor η 0 and the slope growth factor η 1 as

with conditional independence among the u's conditional on the factors. The item response model can be expressed in terms of continuous-normal latent response variable u * underlying the binary outcomes,

with the usual probit standardization of unit variance for the residual. The factor metric can be set either by a factor loading fixed at 1 for one item or by the factor residual variance fixed at 1 for one time point. With binary items, it turns out that the former approach is preferable due to making the Bayes iterations converge faster. The variance of the prior for the non-invariance, that is the difference between each pair of loading or threshold parameters, is chosen as 0. 10. This implies a prior belief that 95% of the distribution of the (unstandardized) non-invariance lies in the range [-0. 62, +0. 62], a range which almost covers the 0. 7 magnitude noninvariance. Using a prior variance of 0. 15 corresponds to the range [-0. 76, +0. 76] and gives similar results. Using the much smaller prior variance 0. 01 corresponds to the narrower range [-0. 20, +0. 20] and is not large enough to let non-invariance be sufficiently expressed for this example and does not give good BSEM point estimates or coverage. The non-invariance detection is still reasonably good also in this case, however, although indicating too many invariant parameters. 

As in the multiple-group simulation, the form of non-invariance simulated in this example is not of the type that is assumed for BSEM approximate invariance. 

This implies that the alignment will not be perfect and the population values used to generate the data will not be perfectly recovered. A sample size of n = 1000 and 100 replications are used in the three simulation analyses. The approximate BSEM invariance is specified with prior variance 0. 10 for both loadings and thresholds. 

Appendix Section 9. 3 shows the input for the Monte Carlo simulation of the second analysis, the BSEM analysis. A total of 50 items are analyzed, 10 items for each of the 5 time points. The metric of the factor is determined by fixing the loading of the last item at the first time point to 1. This approach has been found to perform better than setting the metric in the factor variance at one time point. 

As shown in the table, MODEL PRIORS therefore excludes the loadings of the last item so that the DO loop for the DIFF statement only goes to 9 instead of 10. Instead, a MODEL CONSTRAINT section is added to express the differences between the last item's loadings and their average over time points. In MODEL PRIORS the approximate invariance for the last item is specified for time points 2-5 using a prior with mean 1, that is, the same value as the fixed loading of the last item at the first time point. 

Table 7 shows the results for the difference at each of the five time points for each of the 10 item's loading and threshold relative to their mean over time points. 

In parentheses are given the estimates of the power to detect non-invariance. All non-invariant parameters have a power estimate of at least 0. 80. These power estimates are bolded. No invariant parameters has a power estimate larger than 0. 28. This suggests that BSEM will most often diagnose the correct parameters as non-invariant. 

Table 8 shows the results of the five analyses to be discussed. Five key parameters are considered. First the mean of the slope growth factor (the mean of the intercept growth factor is zero in this parameterization where the thresholds are held equal or approximately equal across time) and the variances of the intercept and slope growth factors. Next, the standardized slope mean is considered to include a scale-free quantity. A final quantity considered is the product of the last item's factor loading and the slope mean which relates to how quickly the last item probability changes over time. 

Analysis 1 with exact invariance gives unacceptable results in that the bias and coverage is poor for the slope mean and the loading × slope mean product. 

Analysis 2 using BSEM with small-variance prior variance 0. 10 gives improved results but due to the alignment issue it still gives a large bias for the slope mean and poor coverage for the loading × slope mean product. Analysis 3 with exact invariance for invariant items based on analysis 2 while maintaining BSEM approximate invariance for other items does not give uniform improvement over analysis 2. Analysis 4 with free non-invariant parameters based on analysis 2 while maintaining BSEM approximate invariance for other parameters also does not give uniform improvement over analysis 2. Analysis 5 that frees the noninvariant parameters while having exact invariance for the other parameters is the preferred analysis and gives good bias results and coverage for all five quantities. 

The remaining bias seen in analysis 5 is due to sampling error and improves with a larger sample size. 

As for the multiple-group simulation study, the outcome of the growth simulation suggests that the proposed two-step approach works well where the BSEM analysis (analysis 2) is followed by analysis 5. The strong power results

shown in Table 7 are key to the good performance. Table 9 shows the time points with significant differences relative to the average across time points for both loadings and thresholds. Among the 144 possible instances of non-invariance, 50 are detected as significant in this BSEM analysis. 

No item is invariant across all time point with respect to both loadings and thresholds. Time points one and eight show the largest number of non-invariant parameters which may be natural given that these time points span the time range. The table is informative about the invariance properties of the items. Of the nine items, the items stubborn and take property show the largest degree of non-invariance, while the item lies shows the least degree of non-invariance. 

Table 10 shows the growth factor estimates using the three analyses: Exact invariance, BSEM, and exact invariance except freeing non-invariants (analyses 1, 2, and 5 in Section 6). The results are quite different among the three analyses with a much steeper linear increase using the third analysis. The third analysis was found to be the better performing in the simulation studies of Section 6.
section_name: Intercepts

section_content: The multiple-group and multiple-time point simulations show that the BSEM approach with approximate measurement invariance is capable of detecting noninvariant parameters. It is a convenient way to identify the specific combination of item and group/time point that gives rise to the non-invariance. The simulations also show that the proposed two-step procedure of BSEM followed by a regular Bayes analysis that frees the non-invariant parameters works well. The procedure easily scales to many groups and time points, as well as many latent variable The DO DIFF approach is used in Table 11 to analyze multiple groups, which for Bayes is carried out using KNOWNCLASS (not shown) and TYPE=MIXTURE. 

The ANALYSIS specification MODEL=ALLFREE lets parameter arrays be different across groups. The auto labeling statement f BY y1-y4* (lam# 1 -lam# 4); saves an increasingly large amount of typing with an increasing number of groups. 

9. 2 Mplus input for multiple-group simulation
section_name: Conclusions

section_name: A growth model application

section_content: As a growth model example, consider a teacher-rated measurement instrument capturing aggressive-disruptive behavior among a sample of U. S. students in Baltimore public schools (Ialongo et al., 1999). A total of 1174 boys are observed in 41 classrooms at eight different time points, from the Fall of Grade 1 to the Spring of Grade 7. Nine items scored as 1 (Almost Never) through 6 (Almost Always) are considered. The item distribution is very skewed with a high percentage in the Almost Never category. The items are therefore dichotomized into Almost Never versus the other categories combined. A multipleindicator growth model is considered for the eight time points. With maximumlikelihood estimation, this requires eight-dimensional numerical integration which is feasible with Monte Carlo integration, for instance using 5000 integration points. Maximum-likelihood, however, gives heavy computations with this many latent variable dimensions, has difficulty in obtaining precise likelihood values (Asparouhov & Muthén, 2012a), does not automatically produce indices of noninvariance, and does not scale to more time points as easily as Bayesian analysis. 

The multilevel aspect of students within classrooms is ignored in this analysis, but leads to even more latent variable dimensions; see Muthén and Asparouhov (2012b). 

Appendix Section 9. 4 shows the input for BSEM analysis with approximate measurement invariance across the eight time points. A quadratic growth model is used with zero quadratic growth factor variance. The aim of the BSEM analysis is to detect non-invariance in the measurement parameters. In line with the simulations of the previous section, the metric is set in the factor loading for dimensions. 

The proposed BSEM approach has the advantage that it can also be combined with the types of BSEM applications discussed in Muthén and Asparouhov (2012a). For example, small-variance priors can be used to accommodate small cross-loadings in multifactorial confirmatory factor analysis models. This flexibility may be necessary in for example cross-cultural studies where a confirmatory factor analysis model with exactly zero cross-loadings is often unrealistic. The COVARIANCE option can also be used to specify small differences between parameters. Note that The Mplus DIFF option is used with MODEL PRIORS to simplify specifying differences between parameters. The above example can be specified as: DIFF(a, b) ∼ N(0, 0. 01); DO DIFF is used to express parameter differences between large sets of parameters and groups/timepoints. Consider the example of group differences for the loadings on a factor measured by 4 variables in 3 groups. Let lamjk denote a factor loading for group j and variable k. Table 11 shows the Mplus input for this analysis. Here the MODEL PRIORS statement DO(1,4) DIFF(lam1# -lam3#) ∼ N(0,0. 01); results in approximate invariance across the 3 groups, so that for the first variable lam11 ≈ lam21, lam11 ≈ lam31, lam21 ≈ lam31, etc. for the lam parameters for variables 2-4. f1*2; y1-y6*. 5; %c#1% [u$1@15 u$2@16 u$3@17 u$4@18 u$5@19 u$6@20 u$7@21 u$8@22 u$9@23]; f1 BY y1*1 y2*. 7 y3*. 5 y4*1 y5*. 7 y6*. 5 (lam1 1-lam1 6); [y1*0 y2*. 5 y3*1 y4*1. 5 y5*2 y6*2. 5] (nu1 1-nu1 6); [f1*1]; f1*2; y1-y6*. 5; %c#2% [u$1@-15 u$2@16 u$3@17 u$4@18 u$5@19 u$6@20 u$7@21 u$8@22 u$9@23]; f1 BY y1*1 y2*. 7 y3*. 7 y4*1 y5*. 7 y6*. 5 (lam2 1-lam2 6); [y1*0 y2*. 5 y3*1 y4*1. 5 y5*2 y6*2. 7] (nu2 1-nu2 6); [f1*1]; f1*1. 8; y1-y6*. 5; f1 BY u11*1. 70 u12-u19*1 (lam1 1-lam1 9) u20@1; f2 BY u21*1 u22*1. 70 u23-u30*1 (lam2 1-lam2 10); f3 BY u31-u32*1 u33*1. 35 u34-u40*1 (lam3 1-lam3 10); f4 BY u41-u43*1 u44*1. 35 u45-u50*1 (lam4 1-lam4 10); f5 BY u51-u54*1 u55*1. 35 u56-u60*1 (lam5 1-lam5 10); [u11$1-u15$1*0 u16$1-u17$1*. 7 u18$1-u20$1*0] (tau1 1-tau1 10); [u21$1-u30$1*0] (tau2 1-tau2 10); [u31$1-u40$1*0] (tau3 1-tau3 10); [u41$1-u50$1*0] (tau4 1-tau4 10); [u51$1-u55$1*0 u56$1-u57$1*-. 7 u58$1-u60$1*0] (tau5 1-tau5 10); [f1-f5@0]; f1*3 f2*2. 5 f3*2 f4*1. 5 f5*1; i s | f1@0 f2@. 5 f3@1 f4@1. 5 f5@2; [i@0 s*. 5]; i*. 5 s*.
section_name: A growth model application

section_name: Appendix