Measurement invariance in CFA with ordinal data (converegence issues)

I think there are two problems. One is that the model specified is not identified. This can be deduced from the number of parameters being greater than the number of parameters for the saturated model. The main reason for under-identification is that when thresholds are estimated, it is not possible to simultaneously estimate residual variances. Mehta PD, Neale, M.C., Flay BR: Squeezing interval change from ordinal panel data: latent growth curves with ordinal outcomes. Psychol Methods 2004; 9:301-333 describes how to avoid the issue by fixing the first two thresholds to arbitrary values, which identifies means and variances. I have modified your script to fix the thresholds (0 and 1 are more commonly used for this purpose but the values you supplied work as well and it's a matter of a simple transform to get back to the original z-score scale for the thresholds if needed). You may wish to liberate parameters for item means as well - zero is probably a poor choice :) and likely contributes to the poor fit of the model.

The second issue is that the starting values were not that great, and with ordinal data it can be difficult to get to the solution precisely on the first try. Therefore I substituted mxTryHard() for mxRun(), which seemed to solve the optimization issue. There's an mxTryHardOrdinal() which might have worked better, but I forgot to use it here.


radloffdata <- read.table("https://www.guilford.com/add/kline/radloff-lavaan.txt", sep="\t", header = FALSE)
colnames(radloffdata) <- c("x1", "x2", "x3", "x4", "x5", "g")
radloffdata[,1] <- ordered(radloffdata[,1])
radloffdata[,2] <- ordered(radloffdata[,2])
radloffdata[,3] <- ordered(radloffdata[,3])
radloffdata[,4] <- ordered(radloffdata[,4])
radloffdata[,5] <- ordered(radloffdata[,5])
rawData1 <- subset(radloffdata, g == 1)
rawData2 <- subset(radloffdata, g == 2)
# common parameters
alphas <- mxPath(from = "one", to = c("x1", "x2", "x3", "x4", "x5"), arrows = 1,
free = FALSE, values = 0)
varDep <- mxPath(from = "depression", arrows = 2, free = TRUE, values = .5)
lambdas <- mxPath(from = "depression", to = c("x1", "x2", "x3", "x4", "x5"), arrows = 1,
free = c(FALSE, TRUE, TRUE, TRUE, TRUE), values = c(1, 1, 1, 1, 1))
# model for group 1
dataRaw1 <- mxData(observed = rawData1, type = "raw")
kappa1 <- mxPath(from = "one", to = "depression", arrows = 1, free = FALSE, values = 0, labels = "kappa1")
thetas1 <- mxPath(from = c("x1", "x2", "x3", "x4", "x5"), arrows = 2, free = FALSE, values = 1)
thresholds1 <- mxThreshold(vars = c("x1", "x2", "x3", "x4", "x5"), nThresh = c(3, 3, 3, 3, 3),
free = c(rep(c(F,F,T),5)), values = c(1, 1.5, 2.5,
1.5, 2, 2.5,
1, 2, 2.5,
.5, 1.5, 2,
1, 2, 2.5),
labels = c("t11", "t12", "t13", "t21", "t22", "t23", "t31", "t32", "t33",
"t41", "t42", "t43", "t51", "t52", "t53"))
m1 <- mxModel("m1", type = "RAM",
manifestVars = c("x1", "x2", "x3", "x4", "x5"), latentVars = "depression",
dataRaw1, alphas, kappa1, lambdas, thetas1, varDep, thresholds1)
# model for group 2
dataRaw2 <- mxData(observed = rawData2, type = "raw")
kappa2 <- mxPath(from = "one", to = "depression", arrows = 1, free = T, values = .5, labels = "kappa2")
thetas2 <- mxPath(from = c("x1", "x2", "x3", "x4", "x5"), arrows = 2, free = TRUE, values = .7)
thresholds2 <- mxThreshold(vars = c("x1", "x2", "x3", "x4", "x5"), nThresh = c(3, 3, 3, 3, 3),
free = c(rep(c(F,F,T),5)), values = c(1, 1.5, 2,
1.5, 2.2, 2.7,
1, 1.5, 2,
.5, .7, 1,
1, 1.5, 2),
labels = c("t11", "t12", "t2_13", "t21", "t2_22", "t2_23", "t31", "t2_32", "t2_33",
"t41", "t2_42", "t2_43", "t51", "t2_52", "t2_53"))
m2 <- mxModel("m2", type = "RAM",
manifestVars = c("x1", "x2", "x3", "x4", "x5"), latentVars = "depression",
dataRaw2, alphas, kappa2, lambdas, thetas2, varDep, thresholds2)
# combine and run models
MGOrd <- mxModel("MG Configural Invariance",
m1, m2,
mxFitFunctionMultigroup(c("m1.fitfunction", "m2.fitfunction")))
MGOrdFit <- mxTryHard(MGOrd)
MGOrdRef <- mxRefModels(MGOrdFit, run = TRUE)
summary(MGOrdFit, refModels = MGOrdRef)


Does that mean each model (m1 and m2) should be identified in its own in OpenMx (and not lavaan or Mplus??).

No. You can indeed use cross-group equality constraints to identify the overall fitted model, and it appears that what you've done actually is sufficient to identify the model. At least, the `mxCheckIdentification()` function says that your model is identified (at the start values and at the dubious solution). I also tried using `mxCheckIdentification()` to check your model using the v2.11.5 release, but it appears there was a bug in that function present in the release but repaired in the source repository since then. But, there appears to be a bug in `mxRefModels()` as well, which has NOT been repaired since v2.11.5. The saturated model in your case should have 50 free parameters, not 35: fifteen thresholds and 10 polychorics in each of 2 groups. But here are the free parameters in the "saturated" model:

free parameters:
name matrix row col Estimate lbound ubound
1 Saturated m1.satCov[1,2] Saturated m1.satCov x1 x2 0.4211160
2 Saturated m1.satCov[1,3] Saturated m1.satCov x1 x3 0.4658516
3 Saturated m1.satCov[2,3] Saturated m1.satCov x2 x3 0.4776376
4 Saturated m1.satCov[1,4] Saturated m1.satCov x1 x4 0.3975597
5 Saturated m1.satCov[2,4] Saturated m1.satCov x2 x4 0.4093033
6 Saturated m1.satCov[3,4] Saturated m1.satCov x3 x4 0.4575572
7 Saturated m1.satCov[1,5] Saturated m1.satCov x1 x5 0.4202486
8 Saturated m1.satCov[2,5] Saturated m1.satCov x2 x5 0.4876613
9 Saturated m1.satCov[3,5] Saturated m1.satCov x3 x5 0.6360375
10 Saturated m1.satCov[4,5] Saturated m1.satCov x4 x5 0.4631449
11 x1ThrDev1 Saturated m1.thresholdDeviations 1 x1 0.7544172 -Inf
12 x1ThrDev2 Saturated m1.thresholdDeviations 2 x1 0.6649991 0.01
13 x1ThrDev3 Saturated m1.thresholdDeviations 3 x1 0.4499984 0.01
14 x2ThrDev1 Saturated m1.thresholdDeviations 1 x2 1.0062369 -Inf
15 x2ThrDev2 Saturated m1.thresholdDeviations 2 x2 0.5285165 0.01
16 x2ThrDev3 Saturated m1.thresholdDeviations 3 x2 0.2991975 0.01
17 x3ThrDev1 Saturated m1.thresholdDeviations 1 x3 0.5143386 -Inf
18 x3ThrDev2 Saturated m1.thresholdDeviations 2 x3 0.5868987 0.01
19 x3ThrDev3 Saturated m1.thresholdDeviations 3 x3 0.3138043 0.01
20 x4ThrDev1 Saturated m1.thresholdDeviations 1 x4 0.2883979 -Inf
21 x4ThrDev2 Saturated m1.thresholdDeviations 2 x4 0.7056041 0.01
22 x4ThrDev3 Saturated m1.thresholdDeviations 3 x4 0.4837302 0.01
23 x5ThrDev1 Saturated m1.thresholdDeviations 1 x5 0.5507459 -Inf
24 x5ThrDev2 Saturated m1.thresholdDeviations 2 x5 0.6843228 0.01
25 x5ThrDev3 Saturated m1.thresholdDeviations 3 x5 0.4456533 0.01
26 Saturated m2.satCov[1,2] Saturated m2.satCov x1 x2 0.5204439
27 Saturated m2.satCov[1,3] Saturated m2.satCov x1 x3 0.3353616
28 Saturated m2.satCov[2,3] Saturated m2.satCov x2 x3 0.3702252
29 Saturated m2.satCov[1,4] Saturated m2.satCov x1 x4 0.3284878
30 Saturated m2.satCov[2,4] Saturated m2.satCov x2 x4 0.3170401
31 Saturated m2.satCov[3,4] Saturated m2.satCov x3 x4 0.3350630
32 Saturated m2.satCov[1,5] Saturated m2.satCov x1 x5 0.4858893
33 Saturated m2.satCov[2,5] Saturated m2.satCov x2 x5 0.3613383
34 Saturated m2.satCov[3,5] Saturated m2.satCov x3 x5 0.4695005
35 Saturated m2.satCov[4,5] Saturated m2.satCov x4 x5 0.4674975

What's happening in the "saturated" model is that all the thresholds are constrained equal across groups. I verified that by examining MGOrdRef$Saturated$'Saturated m1'$thresholdDeviations and MGOrdRef$Saturated$'Saturated m2'$thresholdDeviations. Both had the same set of labels.

So, now I don't think that identification is your problem. I think it's just a matter of getting an analysis of ordinal data to converge (which can be tricky). I have a few tips that I'll put in a subsequent post.

Here are some things you can do (they are not necessarily mutually exclusive):
1. SADMVN, the algorithm OpenMx uses for the multivariate-normal probability integral, can be given an upper limit ,'mvnRelEps', to how much error is acceptable in the probabilities it reports to OpenMx. You can reduce this limit (i.e., request less error from SADMVN) with something like

MGOrd <- mxOption(MGOrd, 'mvnRelEps', mxOption(MGOrd, 'mvnRelEps')/5)

before you run your MxModel. Reducing 'mvnRelEps' will make SADMVN, and therefore OpenMx, run more slowly. Be careful: it is possible to set 'mvnRelEps' so low as to make it mathematically impossible for SADMVN to satisfy it. In that case, OpenMx will warn you about it, but maybe after running for a very long time.
2. You can replace `mxRun()` with mxTryHardOrdinal().
3. You can run your model through mxAutoStart() before running it.