Bifactor always fits best

Recently, I reviewed a paper where the authors wanted to test whether a second-order factor or bifactor model fitted their data best. They had about 580 observations for 6 items, measuring two factors. As this amounts to three item indicators per factor, it is like two just-identified latent variable models in one. But with two factors, the model is clearly overidentied. However, I wondered whether their tests would be able to find the true model. In general, I expect bifactor models to often fit better, as they are less restricted versions of second-order factor models. I wondered whether bifactor models even fit better when the true model is the second-order factor model. I performed a small simulation to see.

I generated 100 datasets with 6 items, measuring two second-order factors, and one first-order factor:

library(lavaan)
set.seed(642)
rdata <- list()
for (i in 1:100) {
  # generate first-order factor:
  FO_factor <- rnorm(580)
  # generate second-order factors:
  SO_factor1 <- FO_factor + rnorm(580, sd = 1)
  SO_factor2 <- FO_factor + rnorm(580, sd = 1)
  # generate the observed data:
  rdata[[i]] <- data.frame(X1 = SO_factor1 + rnorm(580, sd = .75), X2 = SO_factor1 + rnorm(580, sd = .75), 
                           X3 = SO_factor1 + rnorm(580, sd = .75), X4 = SO_factor2 + rnorm(580, sd = .75),
                           X5 = SO_factor2 + rnorm(580, sd = .75), X6 = SO_factor2 + rnorm(580, sd = .75))
}

Then I created a unidimensional model, a model with two correlated factors, a second-order factor model, and two bifactor models to fit on the datasets:

library(semPlot)
F1.mod <- '
  F =~ X1 + X2 + X3 + X4 + X5 + X6
'
semPaths(F1.mod)

F2.mod <- '
  F1 =~ X1 + X2 + X3
  F2 =~ X4 + X5 + X6
'
semPaths(F2.mod)

SO.mod <- '
  F1 =~ X1 + X2 + X3 
  F2 =~ X4 + X5 + X6
  F =~ a*F1 + a*F2
'
semPaths(SO.mod)

BF.mod <- '
  G =~ X1 + X2 + X3 + X4 + X5 + X6
  F1 =~ X1 + X2 + X3
  F2 =~ X4 + X5 + X6
'
semPaths(BF.mod)

BF2.mod <- '
  G =~ X1 + X2 + X3 + X4 + X5 + X6
  F1 =~ X1 + X2 + X3
'
semPaths(BF2.mod)

Then I fitted the models to each of the datasets and calculated p-values (of the chi-square statistic), CFIs and RMSEAs:

F1.fit <- list()
F2.fit <- list()
SO.fit <- list()
BF.fit <- list()
BF2.fit <- list()
cfis <- list()
rmseas <- list()
pvals <- list()
results <- data.frame(dataset = rep(1:100, each = 5), pval = NA, cfi = NA, 
                      rmsea = NA, mod = rep(1:5, times = 100))
fit.ind <- c("pvalue", "cfi", "rmsea")
for (i in 1:100) {
  F1.fit[[i]] <- cfa(F1.mod, data = rdata[[i]], std.lv = TRUE)
  F2.fit[[i]] <- cfa(F2.mod, data = rdata[[i]], std.lv = TRUE)
  SO.fit[[i]] <- cfa(SO.mod, data = rdata[[i]], std.lv = TRUE)
  BF.fit[[i]] <- cfa(BF.mod, data = rdata[[i]], orthogonal = TRUE, std.lv = TRUE)
  BF2.fit[[i]] <- cfa(BF2.mod, data = rdata[[i]], orthogonal = TRUE, std.lv = TRUE)
  results[(i-1)*5+1, 2:4] <- fitmeasures(F1.fit[[i]], fit.ind)
  results[(i-1)*5+2, 2:4] <- fitmeasures(F2.fit[[i]], fit.ind)
  results[(i-1)*5+3, 2:4] <- fitmeasures(SO.fit[[i]], fit.ind)
  results[(i-1)*5+4, 2:4] <- fitmeasures(BF.fit[[i]], fit.ind)
  results[(i-1)*5+5, 2:4] <- fitmeasures(BF2.fit[[i]], fit.ind)
}
results$mod = rep(c("one-factor", "two-factor", "second-order", "bifactor", "bifactor2"), times = 100)
head(results)
aggregate(results[2:4], by = list(results$mod), FUN = mean, na.rm = TRUE)

Not surprisingly, there are a lot of warnings for the first bifactor model, because it may not be identified. We should take the fit indices for that model with a grain of salt.

On average, p-values, RMSEA and CFI indicate bad fit for the one-factor model, as they should. The second-order and two correlated-factor models fit the data well, on average, according to the p-value, CFI and RMSEA. Note that these models are equivalent models, so should fit the data identical.

The most striking findings is that the second bifactor model always seems to fit the data better than the second-order factor model, which was the model that generated the data!

How about the variances of the fit indices, do they differ much between models?

aggregate(results[2:4], by = list(results$mod), FUN = sd, na.rm = TRUE)

They don't seem to differ much. If anything, the variance of the p-values seems comparable between bifactor and second-order factor models. The variance of the CFI is lower for bifactor models, and the variance of RMSEA is lower for the second-order factor model.


Conclusion 

There is not much sense in comparing the fit between second-order and bifactor models if we want to know which model generated the data. A bifactor model will often fit the data better, even when the data were generated using a second-order factor model. I am curious whether this is also the case with more indicators per factor.

On calculation of factor means and item intercepts in LISREL, using ML estimation

Options for identification

Among the options for identifying factor means are:
- Fixing the first items intercept to a constant (in most cases: zero)
- Fixing the latent mean to a constant (in most cases: zero)
- Forcing the intercepts to sum to a constant (in most cases: zero)

The second one is not suitable if you expect latent means to differ, and are interested in latent mean or intercept differences. The first and third one are more suitable in such a case, but offer some challenges as well. The third offers a challenge, because I don't know how to do this in LISREL. The first offers a challenge, because item intercepts have to be interpreted as some kind of deviation from the item with intercept fixed to a constant.



Identification by fixing the first item's intercept to zero

In this case, factor means and item intercepts are very dependent on the item chosen for identification. This is an important consideration if one wants to test an hypothesis concerning factor means, for example in multiple group comparisons. It almost amounts to a test concerning observed means for the items used for identification.

Formulae

Let i denote an item used for identification
Let j denote any other item
Let k denote a factor

Let y_ denote an observed item mean
Let alpha denote a factor mean
Let lambda denote a factor loading
Let tau denote an item intercept

alpha(k) = lambda(i) * y_(i)
tau(j) = y_(j) - lambda(j) * alpha(k)

When the model is identified by constraining the first items intercept to 0, and the first items loading to 1 is used, this simplifies to:

alpha(k) = y_(i)
tau(j) = y_(j) - lambda(j) * y_(i)

In words

The factor mean equals the observed mean of the identification item. Other items intercepts are determined by the items factor loading, times the identification items observed mean. An illustration is provided below, for a model of 3 factors fitted to 2 groups.


Examples

Values are taken from LISREL in- and output provided below. Small discrepancies arise from rounding.

1)
The first item is used as an identification item for the first factor.
In the first group, the item intercept for the second item (item 2) of the first factor is calculated as follows:
y_(item1) = 0.425 = alpha(factor1)
tau(item2) = 0.10 = y_(item2) - lambda(item2)*y_(item1) = 0.511 - 0.98*0.425

2)
The fourth item is used as an identification item for the second factor. In the second group, the item intercept for the second item of factor 2 is calculated as follows:
y_(item1) = 0.673 = alpha(factor2)
tau(item11) = -0.01 = y_(item11) - lambda(item11)*y_(item4) = 0.511 - 0.71*0.673


LISREL input

Observed means (y_; item 1-21):

G1: 0.425 0.511 0.488 0.618 0.380 0.320 0.611 0.582 0.182 0.222 0.406 0.389 0.387 0.602 0.530 0.570 0.588 0.092 0.042 0.206 0.375

G2: 0.566 0.459 0.475 0.673 0.289 0.275 0.590 0.559 0.247 0.288 0.474 0.463 0.343 0.553 0.422 0.603 0.494 0.202 0.023 0.217 0.642


Factor pattern (item 1-21):

item factor
1 1
2 1
3 1
4 2
5 1
6 1
7 1
8 1
9 1
10 1
11 2
12 2
13 2
14 1
15 2
16 3
17 2
18 3
19 3
20 2
21 3


LISREL output

Estimated factor means (alpha, factor 1-3):
G1: 0.425 0.618 0.570
G2: 0.566 0.673 0.603


Estimated item intercepts (tau):

item G1 G2

factor 1
1 0.00 0.00
2 0.10 -0.14
3 0.01 -0.23
5 0.03 -0.10
6 0.01 -0.34
7 0.11 -0.26
8 0.19 -0.40
9 0.05 -0.20
10 0.01 -0.13
14 0.16 -0.14

factor 2
4 0.00 0.00
11 0.05 -0.01
12 0.06 0.04
13 0.11 -0.21
15 0.09 -0.09
17 0.00 0.10
20 0.11 0.14

factor 3
16 0.00 0.00
18 0.23 -0.21
19 0.21 -0.04
21 0.16 -0.15


estimated factor loadings (lambda):

item G1 G2

factor 1
1 1.00 1.00
2 0.98 1.07
3 1.12 1.24
5 0.83 0.68
6 0.78 1.09
7 1.18 1.50
8 0.92 1.70
9 0.54 0.79
10 0.51 0.74
14 1.03 1.23

factor 2
4 1.00 1.00
11 0.57 0.71
12 0.72 0.63
13 0.81 0.82
15 1.00 0.75
17 0.96 0.58
20 0.15 0.12

factor 3
16 1.00 1.00
18 0.56 0.69
19 0.45 0.11
21 0.94 1.31