Using Bifactor Models to Examine the Predictive Validity of Hierarchical Constructs: Pros,Cons,and Solutions

Data generation and model estimation

We also manipulated the following factors to evaluate the performance of augmented bifactor predictive models: (a) loadings on the general factor (LG): .40 and .60; (b) loadings of the additional indicator on the general factor (LG1): .40 and .60; (c) loadings on specific factors (LS): .40 and .60; (d): whether the augmented bifactor model was used or not (augmented); (d) the number of indicators per facet (indicators): 4 and 8; and (f) sample size (SS): 200 to 800 with increments of 100. Finally, we fitted four different models to each generated data set to mimic strategies researchers are likely to adopt in realistic situations. For example, some researchers are only interested in the general factor. Therefore, they would regress outcome variables only on the general factor and omit all specific factors from the regression equation; this is Fitted Model 1 (FM1). Other researchers might be interested in the general factor and do not mind exploring the effects of specific factors. They might also be aware of the strategy suggested by Chen and colleagues (2006). Therefore, they would regress outcome variables on the general factor and some of the specific factors. FM2 using one specific factor and FM3 using two specific factors represent strategies that are likely to be adopted by these researchers. Many users of bifactor predictive models are simultaneously interested in the effects of the general factor and all specific factors and would include all factors in the model when predicting external variables (e.g., Chen et al., 2012, 2013; Debusscher, Hofmans, & De Fruyt, 2017; McAbee et al., 2014; Moshagen, Hilbig, & Zettler, 2018), which is FM4. FM1, FM2, and FM3 were labeled as partial models because only a subset of specific factors was included in the prediction, and FM4 was labeled as the full model because all the factors were used in the regression.

In total, there were 3,584 conditions (4 × 2 × 2 × 2 × 2 × 2 × 7 × 4 = 3,584). To achieve better representation of realistic situations and increase the generalizability of our findings, factor loadings within each replication were randomly drawn from a uniform distribution U (L – .10, L + .10), where L is the prespecified loading for that condition. Adopting .20 as the range for the uniform distribution is intended to represent realistically satisfactory bifactor structures where loadings on the same factor should have a similar size despite some variability. If loadings on the same factor vary greatly (e.g., from .20 to .80), one might question the appropriateness of the bifactor model. But we note that the choice of range may influence results given that a wider range should result in a lower probability of empirical underidentification. The outcome variable was set to be a single manifest variable. Population correlation matrices were generated using

∑ = Λ Φ Λ ' + Θ (2)

where Λ is the factor loading matrix (including the outcome variable whose loading on its factor was set to 1), Φ is the latent factor correlation matrix (those shown in Table 1), and Θ is the diagonal error variance-covariance matrix scaled to ensure the diagonal elements of Σ equal 1 (fixing all factor variance to 1). Individual-level sample data were generated from the population matrices using the psych R package (Revelle, 2018). Each simulated sample data set was fitted by models FM1 through FM4 using the R package lavaan (Rosseel, 2012). Maximum likelihood estimation was used. When software reported any anomalous result (e.g., nonconvergence, negative variance, or nonpositive definite residual correlation matrix) for a replication, it was rejected, and another draw was taken until 1,000 valid replications were obtained. The number of anomalous replications was recorded for each condition. All indices detailed in the following were calculated from valid replications. This strategy ensures that comparisons between conditions are based on the same number of valid replications so that the influence of sampling variability is held constant (Castro-Schilo et al., 2016) while simultaneously we can assess the likelihood of problematic solutions. The R package SimDesign was used to streamline the simulation process (Sigal & Chalmers, 2016). R script and results can be found on the Open Science Framework website (https://osf.io/kjt5m).

Analysis

Our primary interest in the present study was the four beta coefficients that represent the predictive validity of the general factor and the three specific factors. We focused on four criteria when evaluating the performance of traditional bifactor predictive models and augmented bifactor predictive models: (a) the frequency of anomalous results, (b) power/Type I error rate for hypothesis tests about the regression coefficients, (c) bias of regression coefficients, and (d) root mean square error (RMSE) of regression coefficients. Any replication reporting errors or warnings like nonconvergence, nonidentification, negative variance, and nonpositive definite residual matrix was considered an anomaly. We used the anomaly rate (AR) to index the frequency of estimation problems, calculated as follows:

A R = No . of Anomalies 1000+No . of Anomalies

Power was defined as the proportion of replications where a true nonzero effect was correctly detected, and the Type I error rate was defined as the proportion of replications where a true zero effect was falsely found to be significant. Bias was calculated as the difference between the average value of the estimate and the true parameter. Positive bias indicates overestimation, and negative bias indicates underestimation. Bias that falls between –10% and 10% of the true value is often considered acceptable (Muthén & Muthén, 2002). RMSE is given by the square root of the average of the squared difference between the estimated and true parameter. RMSE is an index of overall accuracy that combines bias and sampling variability. Note that although there are three specific factors in our models, when more than one specific factor was included in the regression, we only report power/Type I error rate for the first specific factor because the pattern was almost identical across specific factors.

Anomaly rate

Frequent anomalous results were observed in some conditions. For example, the AR was sometimes greater than .40 when the classic bifactor model (FM4) was fitted to the data in small samples. To provide a precise estimate of the influence of different simulation factors and their interactions, we computed an eight-factor analysis of variance (ANOVA) with the AR as the dependent variable. All main effects and two-way interactions were included in the model. This model explained 86.5% of the total variance in AR. We report predictors that explained at least 2% of the total variance. The strongest main effect was observed for sample size (η² = .205), followed by the fitted model (η² = .095), the type of bifactor model used (BiM vs. ABiM; η² = .073), number of indicators per facet (η² = .051), the magnitude of factor loadings on the general factor (η² = .040), and the magnitude of factor loadings on specific factors (η² = .032). The effect of data-generating model and the magnitude of factor loading of the extra indicator on the general factor had negligible effect on the AR (η²s < .01). Several interactions were also substantial. Specifically, sample size interacted significantly with the type of bifactor model used (η² = .041), fitted model (η² = .026), factor loadings on the general factor (η² = .047), number of indicators per facet (η² = .058), and factor loadings on specific factors (η² = .032). The fitted model also interacted significantly with the type of bifactor model used (η² = .022).

Figure 3 provides a visualization of the relationship between the AR and different factors. When both the general factor and all specific factors were used to predict the outcome, researchers are very likely to encounter anomalous results, especially when sample size is small and loadings on the general factor are relatively low. When loadings on the general factor were .40 and the number of indicators per facet was four, the anomaly rate exceeded 10% for all sample sizes less than 700 and was almost 35% for samples of 200. Augmenting traditional bifactor models with just one more indicator greatly reduced the rate of anomalous results. For example, when the sample size was 200, the number of indicators per facet was four, and average loading on the general factor was .40, adding one more indicator reduced the anomaly rate from nearly 35% to about 10%. When the sample size reached 300, the anomaly rate dropped below 5%. The same pattern held in other less extreme conditions.

OPEN IN VIEWER

Figure 3.

Anomaly rate in different conditions. LG = loadings on the general factor; indicators = number of indicators per facet; ABiM = augmented bifactor model; BiM = bifactor model; full = full predictive model (FM4); partial = partial predictive model (FM1-FM3).

Given the general statistic dictum that a greater number of parameters generally makes model estimation harder when sample size is held constant, it is surprising that increasing the number of indicators per facet ameliorated the anomaly problem (but see Idaszak, Bottom, & Drasgow, 1988, who found improved results for factor analysis when increasing the number of indicators per factor). For example, when we held the sample size to 200 and average loadings on the general factor to .40 and BiM was used, the anomaly rate was close to 35% when the number of indicators per facet was four. It was 17% when the number of indicators per facet was doubled. When the average loading on the general factor was .60, the same pattern was observed: The AR dropped from about 23% to 10% when number of indicators per facet was doubled. In the best condition where ABiM was used (average loading on the general factor was .60, and number of indictors per facet was eight), the anomaly rate was 2% even when sample size was only 200.

Power, bias, RMSE

General factor

The first column of Figure 4 shows results for the conditions where only the general factor predicted the outcome variable in the population (GM1). The three partial models had power greater than .80 even when the sample size was only 200. Furthermore, Figure 5 shows estimation of the regression coefficient for the general factor was virtually unbiased. When the full model FM4 was fitted, power patterns differed substantially from the three partial models. Specifically, when the classic bifactor model with four indicators per facet was fitted, power was below .50 for the sample size of 200 and remained below .80 even when the sample size was 800. For the augmented bifactor models, power reached .80 for samples of 400. For eight indicators per facet, samples of 300 sufficed to reach the .80 power level. More importantly, the additional power of the augmented model came without increased bias. In fact, the augmented model stabilized parameter estimates, as demonstrated by the shrinking RMSEs in the first column of Appendix C.

OPEN IN VIEWER

Figure 4.

Power or Type I error rate for the general factor. Columns 1 and 3 represent power under different conditions. Columns 2 and 4 stand for Type I error rate in different conditions. GM = generating model; FM = fitted model; ABiM = augmented bifactor model; BiM = bifactor model; full = full predictive model (FM4); partial = partial predictive model (FM1-FM3). The number following ABiM or BiM indicates number of indicators per facet. The three red dashed lines correspond to y values of .05, .10, and .8.

OPEN IN VIEWER

Figure 5.

Bias for the general factor. GM = generating model; FM = fitted model; ABiM = augmented bifactor model; BiM = bifactor model. The number following ABiM or BiM indicates number of indicators per facet. Full = full predictive model (FM4); partial = partial predictive model (FM1-FM3). The three red dashed lines correspond to y values of –.025 and .025, which reflected 10% relative bias.

When the general factor and all three specific factors were predictive of the outcome (GM3), the three partial models (FM1-FM3) had almost perfect power for the regression coefficient of the general factor. However, these high levels of power were obtained at the cost of severely overestimating this regression coefficient. As can be seen in the third column of Figure 5, regardless of sample size, even the lowest bias was higher than .125 (50% relative bias) for the three partial models. Apparently, the predictive validity of the omitted specific factors was shifted onto the general factor. The more variables omitted, the more serious the overestimation (as can be seen in the third column of Appendix B). Interestingly, the ABiM reduced the degree of overestimation by about 20%. The severe bias of estimates of the regression coefficient of the general factor was almost entirely eliminated when the fitted model FM4 matched the generating model GM3. But power was reduced: More than 800 respondents would be required to obtain power of .80 when the classic bifactor model with four indicators per facet was used. For the augmented model, 400 respondents sufficed to obtain power above .80. If the number of indicators per facet were doubled, only 300 respondents were needed.

Specific factors

As can be seen from the second column of Figure 6, when only one of the specific factors was predictive (GM2), fitting the partial models produced power higher than .80 as long as sample size exceeded 300, regardless of whether the classic or augmented bifactor model was used. Moreover, Figure 7 shows that bias was close to zero across conditions. However, when the full model FM4 was fitted, power was reduced; the classic model with four indicators per facet never reached a power of .80 when sample size was equal to or below 800. Moreover, Appendix F shows that RMSE increased substantially compared to when the partial models were fitted, indicating more volatile parameter estimates. Fortunately, the augmented model reduced the minimal sample size required for a power of .80 to 400. When the number of indicators per facet was doubled, only 300 respondents were sufficient. As displayed in Figure 7, bias stayed within the acceptable range, and the RMSEs of FM4 stayed close to those of the partial models when the augmented model was fitted.

OPEN IN VIEWER

Figure 6.

Power or Type I error rate for specific factor. Columns 2 and 3 represent power under different conditions. Columns 1 and 4 stand for Type I error rate in different conditions. GM = generating model; FM = fitted model; ABiM = augmented bifactor model; BiM = bifactor model; full = full predictive model (FM4); partial = partial predictive model (FM2-FM3). The number following ABiM or BiM indicates number of indicators per facet. The three red dashed lines correspond to y values of .05, .10, and .8.

OPEN IN VIEWER

Figure 7.

Bias for the first specific factor. GM = generating model; FM = fitted model; ABiM = augmented bifactor model; BiM = bifactor model; full = full predictive model (FM4); partial = partial predictive model (FM2-FM3). The number following ABiM or BiM indicates number of items per facet. The three red dashed lines correspond to y values of –.025, and .025, which reflected 10% relative bias.

When the general factor and all specific factors were nonzero predictors (GM3), fitting the two partial models greatly affected power, as shown in the third column of Figure 6. When the classic bifactor model with four indicators per facet was used, the power to detect the nonzero regression coefficient of the first specific factor was below 35% even in a condition (not shown) with N = 2,000. The augmented model increased power substantially over the classic model. Consistent with previous findings, doubling the number of indictors also increased power. The third column of Figure 7 shows severe bias in estimates of the regression coefficient of the specific factor when the two partial models (FM2 and FM3) were fitted. Bias was almost –.25 for the classic model with four indicators per facet. The reason why the augmented model and increased number of indicators per facet increase power is that they reduce the degree of underestimation and sampling variability simultaneously, as shown in the third columns of Figure 7 and Appendix F. When the full model FM4 was fitted, the pattern was nearly identical to what we reported in the previous paragraph. Augmenting one more indicator effectively reduced the minimal sample size required for power of .80 to 400. Further doubling the number of indicators reduced it to 300. Bias still stayed within the acceptable range. RMSE was also within the normal range.

Type I error rate, bias, and RMSE

Data generation and model estimation

When the number of facets is equal to three or two, HOM is nested within CFM. We wanted to avoid such special cases and therefore set the number of facets to four in this study. It also expands the breadth of conditions covered in Study 1. When setting up the population measurement model for predictors under different conceptualizations (ABiM, AHOM, and CFM), we kept item residuals constant across models. Specifically, when the generating model was ABiM, there were two types of items within a facet: One loaded more on the general factor ( λ g = .60 and λ s = .40), and one loaded more on the specific factor ( λ g = .40 and λ s = .60) to ensure that this model would not be statistically equivalent to HOM. When data were generated from AHOM, all items loaded on the first-order factor with λ = .50 ≈ .71 , and all first-order factors loaded on the higher-order factor λ = .50 ≈ .71 . The Schmid-Leiman transformation shows that such an AHOM is equivalent to a BiM with all loadings equal to .51. When CFM was the true model, all items loaded on their corresponding factors with λ = .52 ≈ .72 , and all first-order factors were correlated with each other with ρ = .50 . Item commonalities were kept at .52 across models. The additional indicator in ABiM, AHOM, and CFM had λ = .50 . When setting up the structural model, two types of models were generated for each conceptualization. In the first model (ABiM1, AHOM1, and CFM1), the regression coefficient associated with each predictor was set to .25. In the second model (ABiM2, AHOM2, and CFM2), the regression coefficients of the last two predictors were set to zero while the others remained at .25. Population parameters for measurement models, structural models, and population R ² are presented in Table 2. The outcome variable was always measured by four indicators with loadings equal to .70. Sample size was manipulated to range from 200 to 800 with increments of 100. The number of indicators per facet was set to four or eight. When data generated from ABiM or AHOM were fitted to a CFM, the additional indicator was forced to load on the last factor. ³ When CFM was used to generate the data, we simulated one more indicator for the last factor, and this item was used as the additional indicator that only loaded on to the general factor or higher-order factor in ABiM or AHOM.

OPEN IN VIEWER

Table 2.

Population Parameters Used to Simulate Data (Study 2).

	ABiM							AHOM							CFM
	Population factor loadings
	G	S1	S2	S3	S4	H		G	R1	R2	R3	R4	H		F1	F2	F3	F4	H
X1	.40	.60				.52		.51	.51				.52		.72				.52
X2	.40	.60				.52		.51	.51				.52		.72				.52
X3	.60	.40				.52		.51	.51				.52		.72				.52
X4	.60	.40				.52		.51	.51				.52		.72				.52
X5	.40		.60			.52		.51		.51			.52			.72			.52
X6	.40		.60			.52		.51		.51			.52			.72			.52
X7	.60		.40			.52		.51		.51			.52			.72			.52
X8	.60		.40			.52		.51		.51			.52			.72			.52
X9	.40			.60		.52		.51			.51		.52				.72		.52
X10	.40			.60		.52		.51			.51		.52				.72		.52
X11	.60			.40		.52		.51			.51		.52				.72		.52
X12	.60			.40		.52		.51			.51		.52				.72		.52
X13	.40				.60	.52		.51				.51	.52					.72	.52
X14	.40				.60	.52		.51				.51	.52					.72	.52
X15	.60				.40	.52		.51				.51	.52					.72	.52
X16	.60				.40	.52		.51				.51	.52					.72	.52
Additional	.50					.25		.50					.25					.50	.25
	Population correlations among latent factors
	G	S1	S2	S3	S4			G	R1	R2	R3	R4			F1	F2	F3	F4
G	1							1
S1/R1/F1	.00	1						.00	1						1
S2/R2/F2	.00	.00	1					.00	.00	1					.50	1
S3/R3/F3	.00	.00	.00	1				.00	.00	.00	1				.50	.50	1
S4/R4/F4	.00	.00	.00	.00	1			.00	.00	.00	.00	1			.50	.50	.50	1
	Population structural models
Regression	ABiM1: Y = .25G + .25S1 + .25S2 + .25S3 + .25S4, R 2 = .313							AHOM1: Y = .25G + .25R1 + .25R2 + .25R3 + .25R4, R 2 = .313							CFM1: Y = .25F1 + .25F2 + .25F3 + .25F4, R 2 = .625
	ABiM2: Y = .25G + .25S1 + .25S2 + 0S3 + 0S4, R 2 = .188							AHOM2: Y = .25G + .25R1 + .25R2 + 0R3 + 0R4, R 2 = .188							CFM2: Y = .25F1 + .25F2 + 0F3 + 0F4, R 2 = .186

Note: ABiM = augmented bifactor model; AHOM = augmented higher-order model; CFM = correlated-factors model; S = specific factors in BiM; R = first-order factor residuals in HOM; F = first-order factors in CFM; H = item commonality.

In total, there were 84 conditions (3 × 2 × 2 × 7 = 84). Sample data were generated using the same procedure as Study 1. Maximum likelihood estimation was used. Replications that produced anomalous results were rejected until 1,000 valid replications were obtained. In addition, we also selected a sample size of 10,000 and ran 500 replications for each of these conditions to study the asymptotic performance of each model. All analyses were performed on valid replications.

When reporting results for ABiM1, AHOM1, and CFM1, we only focus on the general factor/higher-order factor and the first specific factor/facet because results for the other specific factor/facets were the almost identical to the first one. When reporting results for ABiM2, AHOM2, and CFM2, we focus on the general factor/higher-order factor, the first and the third specific factor/facet. That is because results for the second and the fourth specific factor/facet were almost identical to the first and the third, respectively. Also, we only report results for conditions whose number of indicators per facet was four because the overall pattern was very similar to conditions that had eight indicators per facet. Detailed results for the eight-item conditions can be found in the Supplemental Material (available in the online version of the journal).

When the generation model was ABiM

As can be seen from Table 3, when the generation model was ABiM1 (βs = .25), fitting ABiM resulted in almost no bias in the regression coefficients and negligible bias in R ² once sample size reached 300. When sample size was 200, R ² was overestimated by 15%. Power for both the general factor and specific factors reached .85 when sample size reached 300. These findings are consistent with what we reported in Study 1. Fitting the AHOM led to somewhat different results. Although regression coefficients were still unbiased, the power to detect these effects decreased substantially compared to ABiM. Specifically, the power to detect the effect of the general factor and specific factor was only .42 and .30 when sample size was 200. Even when sample size reached 500, the power was still below 80%. Also, R ² was overestimated by 32% (bias = .100) when sample size was 200. Even though bias in R ² went down with increasing sample size and was asymptotically unbiased, it was overestimated by 11% (bias = .034) when sample size was 800. Fitting the CFM was more problematic. First, regression coefficient estimates of facets were underestimated by 44% (SOM = .138), and the power to detect the effect of these facets was far below .80 even for a sample size of 800. R ² was also constantly underestimated across sample sizes by 29% to approximately 36% (bias = .091∼.113).

OPEN IN VIEWER

Table 3.

Bias/Size of Misspecification and Power/Type I Error Rate When Data-Generating Model Was ABiM1 and ABiM2 (Study 2).

GM	Key Parameters	FM	Bias/Size of Misspecification									Power/Type I Error Rate
			200	300	400	500	600	700	800	10000		200	300	400	500	600	700	800	10000
ABiM1	R 2 = .313	ABiM	.048	.029	.024	.016	.017	.013	.011	–.003		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		AHOM	.101	.086	.063	.054	.050	.049	.035	.008		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		CFM	–.092	–.101	–.108	–.110	–.112	–.114	–.113	–.123		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	β G = .25	ABiM	.002	.003	.000	.002	–.001	–.001	.002	.000		.75	.89	.95	.99	.99	1.00	1.00	1.00
		AHOM	.012	.002	.001	.002	–.002	–.003	.003	.000		.42	.54	.65	.73	.81	.84	.90	1.00
		CFM	—	—	—	—	—	—	—	—		—	—	—	—	—	—	—	—
	βS1/R1/F1 = .25	ABiM	–.004	.001	–.001	–.003	–.001	–.002	–.001	.000		.70	.87	.92	.98	.99	1.00	1.00	1.00
		AHOM	–.012	.006	.000	.001	.002	.003	.001	.000		.30	.51	.58	.66	.75	.81	.87	1.00
		CFM	–.111	–.101	–.112	–.112	–.113	–.113	–.110	–.110		.23	.33	.39	.43	.50	.59	.65	1.00
ABiM2	R 2 = .188	ABiM	.058	.036	.027	.022	.015	.014	.015	.003		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		AHOM	.137	.093	.073	.056	.046	.042	.037	.003		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		CFM	.021	.008	.004	.000	–.003	–.003	–.002	–.008		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	β G = .25	ABiM	.006	–.001	.000	.003	–.001	.001	–.001	.000		.76	.88	.95	.98	1.00	1.00	1.00	1.00
		AHOM	.005	.001	–.001	.004	–.001	–.001	–.002	.000		.42	.53	.65	.74	.83	.86	.91	1.00
		CFM	—	—	—	—	—	—	—	—		—	—	—	—	—	—	—	—
	βS1/R1/F1 = .25	ABiM	–.005	–.002	–.001	–.003	–.005	–.002	.002	.000		.68	.82	.92	.96	.98	.99	1.00	1.00
		AHOM	–.001	.000	.003	–.002	–.002	.003	.007	.000		.33	.46	.60	.64	.73	.80	.86	1.00
		CFM	.045	.043	.043	.040	.039	.041	.046	.040		.67	.84	.91	.97	.99	1.00	1.00	1.00
	βS3/R3/F3 = 0	ABiM	–.007	.003	–.002	.002	–.002	–.002	–.001	.000		.08	.07	.07	.05	.06	.05	.05	.05
		AHOM	–.007	.000	–.001	.001	–.001	.001	.001	.000		.04	.05	.06	.04	.05	.05	.04	.05
		CFM	–.087	–.075	–.081	–.074	–.078	–.080	–.080	–.080		.12	.11	.17	.15	.19	.24	.27	.99

Note: GM = generation model; FM = fitted model; ABiM = augmented bifactor model; AHOM = augmented higher-order model; CFM = correlated-factors model.

When the ABiM2 was the generation model (β _G = β _{S 1} = .25, β _{S 3} = 0), fitting ABiM resulted in almost no bias in regression coefficients. The power to detect the effect of the general factor and the first specific factor exceeded .80 when sample size reached 300. Type I error rate associated with the third specific factor was close to .05 (.050∼.081) across sample sizes. However, R ² was overestimated by 31% (bias = .058) when sample size was 200. Overestimation of R ² only dropped below 10% when sample size reached 600. Fitting AHOM resulted in no bias in regression coefficients. However, the power to detect the effect of the general factor and the first specific factor was substantially lower, although the Type I error rate associated with the third specific factor was close to .05. Moreover, R ² was more overestimated with AHOM despite being asymptotically unbiased. Specifically, R ² was overestimated by 73% (bias = .139) at a sample size of 200 and by 30% (bias = .056) at a sample size of 500. Even when sample size reached 800, it was still overestimated by 20% (bias = .038). Fitting CFM yielded a different pattern of results. First, compared to the corresponding regression coefficients of specific factors in the population model, the regression coefficient of the first facet was overestimated by about 18% (SOM = .040∼.045), whereas those of the third facet were underestimated by about 33% relative to .25 (SOM = –.0875 to approximately –.075). Second, the power to detect the effect of the first facet was close to that in the bifactor model and higher than AHOM. But the Type I error rate associated with the third facet was .12 and increased with sample size. Surprisingly, the bias in R ² was negligible when sample size reached 300.

When the generation model was AHOM

As shown in Table 4, when the generation model was AHOM1 (βs = .25), fitting ABiM resulted in negligible bias in regression coefficients once sample size reached 300. However, the power to detect the effect of the general factor and specific factors did not reach .80 until sample size reached 600. R ² was also moderately overestimated by 12% to approximately 33% (bias = .038∼.119) until sample size was 800. Results for fitting AHOM were similar to ABiM: There was almost negligible bias in regression coefficients; the power was below .80 until sample size was 600; R ² was overestimated by 11% to approximately 26% (bias = .034∼.081) until sample size reached 800. Fitting CFM led to consistent underestimation of the regression coefficients of facets by 33% to approximately 44% (SOM = .083∼.110), which caused low power (the power was only .70 even at a sample size of 800). R ² was also consistently underestimated by 27% to approximately 35% (bias = .068∼.088).

OPEN IN VIEWER

Table 4.

Bias/Size of Misspecification and Power/Type I Error Rate When Data-Generating Model Was AHOM1 and AHOM2 (Study 2).

GM	Key Parameters	FM	Bias/Size of Misspecification									Power/Type I Error Rate
			200	300	400	500	600	700	800	1,0000		200	300	400	500	600	700	800	1,0000
AHOM1	R 2 = .313	ABiM	.104	.068	.062	.059	.052	.036	.028	.008		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		AHOM	.081	.053	.044	.050	.044	.034	.028	.008		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		CFM	–.084	–.094	–.102	–.105	–.106	–.108	–.109	–.113		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	β G = .25	ABiM	.042	.024	.009	–.002	–.003	.000	.002	.000		.62	.67	.71	.76	.82	.88	.92	1.00
		AHOM	.023	.018	.007	–.002	–.003	–.001	.000	.000		.48	.60	.68	.74	.79	.87	.92	1.00
		CFM	—	—	—	—	—	—	—	—		—	—	—	—	—	—	—	—
	βS1/R1/F1 = .25	ABiM	–.042	–.021	–.012	.001	.003	.000	–.004	.000		.43	.56	.67	.76	.83	.87	.89	1.00
		AHOM	–.021	–.015	–.008	.001	.003	.002	–.001	.000		.31	.47	.61	.70	.80	.85	.89	1.00
		CFM	–.108	–.105	–.108	–.110	–.107	–.106	–.109	–.110		.25	.36	.42	.52	.59	.67	.70	1.00
AHOM2	R 2 = .188	ABiM	.166	.107	.080	.060	.046	.041	.031	.003		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		AHOM	.111	.083	.061	.046	.037	.036	.026	.003		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		CFM	.014	.007	.003	–.003	–.004	–.005	–.007	–.008		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	β G = .25	ABiM	.015	.002	.004	.004	.001	–.002	.001	.000		.58	.62	.72	.77	.82	.87	.90	1.00
		AHOM	.006	–.001	.004	.003	–.001	–.003	–.001	.000		.43	.53	.69	.75	.80	.84	.90	1.00
		CFM	—	—	—	—	—	—	—	—		—	—	—	—	—	—	—	—
	βS1/R1/F1 = .25	ABiM	–.022	–.004	–.005	–.005	–.005	.000	–.004	.000		.48	.61	.67	.72	.79	.85	.88	1.00
		AHOM	–.010	–.001	–.004	–.003	–.004	.002	–.002	.000		.34	.52	.59	.70	.77	.84	.87	1.00
		CFM	.033	.037	.039	.037	.034	.036	.035	.030		.69	.86	.94	.98	.99	1.00	1.00	1.00
	βS3/R3/F3 = 0	ABiM	–.018	–.004	–.010	–.004	–.004	.001	–.004	.000		.16	.13	.12	.11	.11	.09	.07	.05
		AHOM	–.007	–.001	–.010	–.002	–.003	.002	–.003	.000		.03	.04	.05	.04	.05	.05	.04	.04
		CFM	–.070	–.068	–.075	–.066	–.069	–.068	–.070	–.070		.10	.10	.15	.14	.20	.20	.23	.99

Note: GM = generation model; FM = fitted model; ABiM = augmented bifactor model; AHOM = augmented higher-order model; CFM = correlated-factors model.

When the generation model was AHOM2 (β _G = β _{R 1} = .25, β _{R 3} = 0) and the fitted model was ABiM, there was little bias in regression coefficients when sample size reached 300. But the power to detect the effect of the general factor and the first specific factor did not reach .80 until sample size was 600, which was lower than when the generation model was ABiM2. The Type I error rate for the third specific factor was .16 when sample size was 200. It decreased to .07 when sample size reached 800. R ² was overestimated by 16% to approximately 88% (bias = .031∼.165), especially when sample size was small. The pattern for fitting AHOM was close to ABiM, except that AHOM led to somewhat lower power, better Type I error rate, and slightly lower overestimation of R ². Fitting CFM resulted in different results. First, the regression coefficient of the first factor was consistently overestimated by about 12% (SOM = .030), and the regression coefficient of the third factor was consistently underestimated by about 28% (SOM = .070). Correspondingly, the power to detect the effect of the first factor reached .80 at a sample size of 300. The Type I error rate for the third facet was .10 at a sample size of 200 and increased with sample size. As found in conditions where the generation model was ABiM2, there was negligible bias in R ² for CFM.

When the generation model was CFM

As shown in Table 5, when the generation model was CFM1 (βs = .25), fitting both ABiM and AHOM led to consistent underestimation of the regression coefficients by about 25% (SOM = –.06). Consequently, the power to detect the effect of the first facet did not reach .80 until sample size was 800. Surprisingly, bias in R ² was negligible once sample reached 400 for ABiM and 300 for AHOM. Fitting CFM resulted in no bias in regression coefficients and R ² across sample sizes. The power to detect the effect of the first facet was close to .80 even at a sample size of 200.

OPEN IN VIEWER

Table 5.

Bias/Size of Misspecification and Power/Type I Error Rate When Data-Generating Model Was CFM1 and CFM2 (Study 2).

GM	Key Parameters	FM	Bias/Size of Misspecification									Power/Type I Error Rate
			200	300	400	500	600	700	800	1,0000		200	300	400	500	600	700	800	1,0000
CFM1	R 2 = .625	ABiM	.094	.077	.060	.045	.039	.031	.029	–.005		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		AHOM	.076	.054	.044	.033	.030	.025	.023	–.005		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		CFM	.014	.011	.010	.004	.006	.004	.004	–.005		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	βS1/R1/F1 = .25	ABiM	–.053	–.068	–.063	–.060	–.060	–.063	–.066	–.060		.55	.57	.66	.71	.77	.80	.82	1.00
		AHOM	–.057	–.063	–.062	–.061	–.061	–.062	–.063	–.060		.36	.47	.58	.66	.73	.77	.81	1.00
		CFM	.003	–.001	.000	–.001	.001	–.001	.000	.000		.79	.91	.97	1.00	1.00	1.00	1.00	1.00
CFM2	R 2 = .186	ABiM	.188	.145	.100	.080	.060	.051	.048	.002		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		AHOM	.136	.101	.072	.056	.047	.038	.037	.002		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
		CFM	.028	.018	.014	.009	.010	.006	.006	.002		1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	βS1/R1/F1 = .25	ABiM	–.054	–.045	–.045	–.042	–.049	–.050	–.044	–.050		.43	.53	.58	.63	.68	.71	.77	1.00
		AHOM	–.046	–.041	–.047	–.045	–.049	–.050	–.045	–.050		.28	.40	.47	.54	.61	.67	.74	1.00
		CFM	.001	.003	.002	.000	.000	–.005	.002	.000		.60	.79	.88	.92	.96	.99	1.00	1.00
	βS3/R3/F3 = 0	ABiM	.013	.022	.026	.028	.019	.027	.025	.020		.16	.17	.14	.12	.11	.11	.13	.19
		AHOM	.020	.025	.025	.026	.020	.027	.024	.020		.03	.04	.04	.06	.05	.05	.05	.18
		CFM	–.005	–.005	.000	.000	–.003	.004	.000	.000		.06	.06	.05	.06	.05	.07	.05	.04

Note: GM = generation model; FM = fitted model; ABiM = augmented bifactor model; AHOM = augmented higher-order model; CFM = correlated-factors model.

When the generation model was CFM2 (β _{F 1} = .25, β _{F 3} = 0), fitting both ABiM and AHOM led to consistent underestimation of the regression coefficient of the first facet by about 18% (SOM = –.045) and overestimation of the regression coefficient of the third facet by about 10% relative to .25. The power to detect the effect of the first facet did not reach .80 even at a sample size of 800. The Type I error rate associated with the third facet of the ABiM ranged from about .10 to approximately .16, whereas that of the AHOM stayed closer to the nominal level across sample sizes. Regarding R ², even though it was asymptotically unbiased, it was severely overestimated for both ABiM and AHOM, especially when sample size was small. For example, at a sample size of 200, R ² of ABiM and AHOM was overestimated by 101% and 73% (bias = .188 and .136), respectively. Even when sample size reached 800, R ² was overestimated by 26% and 20% (bias = .048 and .037) for ABiM and AHOM. Fitting CFM resulted in virtually no bias in regression coefficients and R ² across sample sizes. The power to detect the effect of the first facet was close to .80 even at a sample size of 300. The Type I error rate associated with the third facet remained around .05 across sample sizes.

Bifactor models and the augmentation strategy

As typically implemented, BiM specifies general and specific factors as orthogonal components, thus allowing researchers to examine their psychometric properties independently and additively. These features are particularly relevant for studying the predictive validity of hierarchical domains; because all factors are orthogonal, regression coefficients can be directly interpreted as the independent contribution of each component. However, our first series of simulations showed that BiM is prone to problems (nonconvergence, etc.) and inadequate power when the number of indicators per facet is modest (e.g., four) and sample size is small to moderate (n < 600). Nonetheless, many constructs in organizational research are measured by a small number of indicators (Hinkin, 1995). Typical sample sizes used in organizational research are also often in the range we examined (Su et al., 2019). Furthermore, interpretation of factors in BiM is problematic in many cases (Eid et al., 2017). These concerns limit the potential usefulness of BiMs.

ABiM inherits the conceptual advantage of BiM but seems to mitigate the problem of nonidentification and construct ambiguity. The additional indicator largely eliminates anomalous results and thus improves the statistical performance of BiMs. ABiMs work reasonably well even in the settings where BiM fails. In addition, unlike BiM, whose factors have ambiguous construct validity, the augmenting indicator gives clear construct meanings to factors in ABiM. The general factor is defined as the augmenting indicator, and its construct meaning should not change as facets are added or removed. Similarly, specific factors are residualized facets relative to the augmenting indicator. Interested readers are encouraged to refer to Eid et al. (2017) for an excellent discussion. An additional advantage of ABiM is that it encourages researchers to think deeper about the construct’s meaning because an indicator of the general factor is needed. This requirement may, in fact, facilitate theory development. Moreover, ABiM seems more robust to some model misspecifications. Although AHOM has the same conceptual virtues and performs almost identically to ABiM asymptotically, its statistical performance at more modest sample levels is less optimal. CFM seems more limited. First, unambiguous interpretation of regression coefficients is not possible when the facets are moderately to highly correlated. Second, it may result in underestimation of the overall predictive validity of hierarchical constructs and lead to wrong conclusions. In sum, the augmentation strategy seems to solve the statistical and conceptual problems arising with BiM while maintaining its intuitive appeal. CFM is not recommended for studying hierarchical constructs from both a conceptual and pragmatic perspective. Whether a construct is hierarchical or not, however, is a topic that has been debated since the time of Spearman and Thurstone.

Considerations of adopting ABiMs and recommendations

Despite the advantages discussed previously, researchers should still be cautious when adopting ABiM because it still inherits some problems of BiM. First, BiMs have an overfitting tendency. In reality, researchers do not know the true model underlying the data but still need to choose one to answer research questions. The common practice is to fit different models to the data and use various model fit indices to select the best-fitting model. BiMs in past studies have often stood out as best fitting (Rodriguez, Reise, & Haviland, 2016). However, recent studies suggested that BiMs tend to beat other models (e.g., HOM and CFM) because they have an overfitting tendency. For example, BiMs can absorb cross-loadings, correlated residuals, and population heterogeneity through the general factor (Murray & Johnson, 2013; Raykov et al., 2019) and fit invalid or even random data (Bonifay & Cai, 2017; Reise et al., 2016). Therefore, researchers are recommended to not only examine model comparison results based on multiple fit indices (e.g., comparative fit index, Tucker-Lewis index, root mean square error of approximation, standardized root mean square residual, Akaike information criterion, Bayesian information criterion) but also to conduct additional robustness check before making their final decision. For example, researchers can use the iteratively reweighted least squares estimation method to ensure that ABiM fits better because it better describes the response process instead of just tolerating invalid responses (Reise et al., 2016). Moreover, researchers are strongly encouraged to evaluate ABiM solutions through statistics such as various omega reliability coefficients and factor determinacy. If these statistics are too low, the estimated ABiM solution may not be meaningful despite good fit statistics. Another less well-known but extremely useful statistic for robustness check is the construct reliability—or more recently termed as construct replicability—that quantifies “how well a latent variable is represented by a given set of observed indicators, and, thus, replicable across studies” (Rodriguez et al., 2016, p.143; Hancock & Mueller, 2001). A meaningful ABiM solution should have moderate to high construct replicability for both the general factor and specific factors. Researchers can refer to Rodriguez et al. (2016) for an excellent tutorial on the calculation and interpretation of these indices. In addition, replication using multiple samples is also encouraged.

Even when robustness checks and theories support the use of ABiMs, researchers should still be careful. Our first series of simulation showed that arbitrarily fixing some paths to zero can lead to overestimation in some cases. Therefore, it appears that all factors should be included in the prediction of external variables because this can keep estimation bias to a minimum and maintain power at an adequate level, even though the researcher is only interested in the general factor.

Our second study demonstrated that an R ² derived from ABiM is likely to be overestimated when sample size is small to moderate. A potential solution is to recalculate R ² by summing up the squares of only statistically significant standardized regression coefficients. Shrinkage formula may also be useful here. We sometimes noticed suspiciously large standard errors appear for regression coefficients (but not for other parameters) even though model estimation terminated normally. The appearance of very large standard errors is a sign of empirical underidentification (Rindskopf, 1984). Results obtained from such conditions are not trustworthy. Although there is no established guideline about how large is “suspiciously large,” researchers can compare standard errors of regression coefficients to the standard errors of other parameters. If the former is 3 or 4 times larger than the latter, there is cause for concern.

Sample size is another extremely important consideration for adopting ABiM. We recommend sample sizes of 500 or more in general. Such large sample sizes may be needed to keep the anomaly rate, power, Type I error rate, and estimation bias within acceptable ranges. This recommendation also coincides with Kretzschmar and Gignac’s (2019) finding that a sample size of 490 is required to achieve stable estimation of latent variable correlations in typical scenarios (ρ = .20, reliability = .70). Tailored simulations can be used to obtain a more precise estimate of sample size requirements because they vary substantially across contexts (Wolf et al., 2013). We provided R code on the Open Science Framework website that interested readers can use to assess their needs. Generally, we recommend against the use of ABiM if sample size is around or below 200.

How to find a good augmenting indicator is another important issue for researchers. We tentatively suggest two approaches. First, researchers can start from existing measures. Because most measures of hierarchical constructs were not developed using bifactor models, it is possible that some of these items may in fact only tap the general factor but are incorrectly assigned to a specific facet. Therefore, researchers can fit bifactor models to such measures. If some items consistently show high loadings on the general factor and trivial loadings on specific factors across multiple samples, such items may be considered as good candidates. Subsequently, researchers can ask a group of subject matter experts to make informed judgments about these candidate augmenting items. Items that survive statistical examinations and expert judgments can be used as augmenting items. Second, researchers can also write augmenting items from scratch. In some cases, the names of hierarchical constructs can be good candidates. For example, the item “I am an agreeable person” might be a good augmenting item for Agreeableness. In some other cases, items that do not refer to a specific target or domain also work (Eid et al., 2017). For example, Brunner and colleagues (2010) used items like “I’m good at most school subjects” to measure general academic self-concept. Fraley et al. (2015) also used items like “I find it easy to depend on people in close relationships” to measure global attachment styles. These items also need to go through statistical scrutiny before being used.

Last but not least, even though specific factors are constrained be orthogonal in the ABiM, it is not necessarily required for model identification. One of the anonymous reviewers recommended that researchers should always examine whether the orthogonal constraint is reasonable by comparing the fit of the ABiM and the more general bifactor-(S·I-1) model via a likelihood ratio test. If there is no statistically significant difference, then researchers can proceed with ABiM. If the ABiM fits significantly worse, researchers should adopt the more general bifactor-(S·I-1) model.

Considerations of choosing between models

A reviewer raised two interesting questions regarding model selection and model misfit that are not only pertinent to the current study but also may be important in broader contexts. While an elaborated discussion is beyond the scope of the current study, we decided to briefly address these questions in hope of raising awareness among researchers.

First, the reviewer asked, “What if a researcher’s specific data set is (1) supported by multiple models, (2) not supported by multiple models, (3) or there is mixed support?” (1) Although not widely recognized, some seemingly very different SEM models can indeed be equivalent in that they yield identical implied covariance or correlation matrices. Therefore, researchers are recommended to use the nesting and equivalence testing (NET) procedures to test for equivalence (Bentler & Satorra, 2010). If the NET procedure indicates equivalence, it is impossible to identify which model is the “true model” statistically. In such cases, researchers are encouraged to refer to their theories to decide which model is most appropriate. (2) If the data are not supported by multiple models, then the researcher might search for a better fitting model by looking for one that is very different from the one the researcher is using. Or one might just choose the model that fits the data the best. However, picking a model based on fit statistics, especially if they are not nested, can be very difficult: Some fit statistics may prefer one model, whereas other fit statistics may prefer another model (but see a recent description of the Vuong test for statistically comparing nonnested models by Merkle, You, & Preacher, 2016). (3) For mixed support, researchers should collect more or better data to help inform their decision. Researchers are encouraged to consider substantive theory while taking into consideration model parsimony.

The reviewer also asked when there is model misfit, how can a researcher probe the source of misfit (misfit from the measurement model or the structural model)? It is possible that the measurement model fits well but the structural model fits badly, or vice versa. However, traditional model fit evaluation focuses on global fit and does not differentiate misfit due to the measurement part and the structural part, which makes it hard to probe the source of misfit (if there is any). Psychometricians have considered these problems and proposed several solutions. For example, Anderson and Gerbing (1988) advocated a two-step strategy. A measurement model is first fitted to the data. If the measurement model shows acceptable fit, then researchers fit the full SEM model. A sequential chi-square test is then used to test specific hypotheses concerning the structural proportion of the model. There are also some fit indices for evaluating the fit of the structural model, including the target coefficients (TC2; Marsh, 1987), RMSEA-P (McDonald & Ho, 2002), and noncentrality structural covariance index, with a focus on the path model component of a composite model (NSCI-P; Williams & O'Boyle, 2011). Lance et al. (2016) provided a good summary of existing indices and developed six new indices. Recently, Thoemmes, Rosseel, and Textor (2018) also suggested the use of graphical criteria like d-separation and trek-separation for evaluating local fit. Modification indices and expected parameter change statistic can also aid with identifying model misfit. In sum, these methods can be used to evaluate fit of different components of an SEM model. Interested readers are encouraged to refer to these articles for more details.

Limitations and future directions

As with any simulation, the current studies are far from comprehensive. For example, all effect sizes were set to be .25, which is frequently seen and important, but other values should be investigated. We also did not consider conditions where the number of items per facet was unequal. Results may differ in these conditions. Another limitation is that we assumed a perfect structure with no cross-loadings and residual correlations. However, empirical data are likely to be messy with small cross-loadings (Asparouhov & Muthén, 2009). It is worth exploring how robust our conclusions are when nonzero cross-loadings are forced to zero rather than correctly modeled. Third, even though the augmentation strategy and increased sample size can reduce anomalous results, it is not always possible to add such an indicator (e.g., archival data analysis) or obtain more respondents (e.g., clinical population). In such cases, it might be interesting to examine whether Bayesian estimation can produce more accurate and stable parameter estimates by supplying informative priors. Fourth, we assumed orthogonality among BiM factors throughout this article because it allows straightforward explanations of regression coefficients of these factors. However, this is not required for the bifactor-(S·I -1) model (Eid et al., 2017). Future studies should examine the generalizability of our current findings in these conditions. Fifth, strictly speaking, adding new items to an existing scale is scale adaptation that might be problematic (Heggestad et al., 2019). Researchers should provide sufficient evidence to support validity of the adapted scale, particularly the construct validity of the augmenting item. Sixth, when comparing ABiM/AHOM and CFM, we compared regression coefficient estimates to the generating model parameters because researchers often interpret them similarly—as the unique contribution of a specific facet/factor. However, such comparisons are not fully statistically justified because predictors in these models may be scaled on somewhat different metrics. As suggested by one of the reviewers, a more justified comparison should be made between the expected regression coefficient for, say, CFM calculated based on the generating BiM model and estimated CFM regression coefficients. However, we are currently not aware of an established way to derive the expected value for CFM regression coefficients based on ABiM/AHOM generating models. In sum, small differences in estimation accuracy across misspecified and correctly specified models should not be overemphasized.