Error when doing measurement invariance

Roy Levy

unread,

Sep 10, 2025, 5:45:47 PM (7 days ago) Sep 10

to blavaan

We're working on measurement invariance analyses, and running into an error that has stumped us. I have a reproducible example below. The code simulates some data (we can't share private data) and then attempts to fit various versions of the model. The model has the following features:

- 3 groups

- 14 observed variables as indicators

- 3 correlated factors

- 17 loadings (i.e., with cross-loadings for some indicators)

- 14 intercepts and 14 error variances for the indicators

- 4 residual correlations among the indicators

We have run a model that constrains just the loadings to be equal across groups using group.equal = c("loadings"), and that runs. We have run a model that constrains just the intercepts to be equal across groups using group.equal = c("intercepts"), and that runs. However, when we run model that constrains both the loadings and the intercepts to be equal across groups using group.equal = c("loadings", "intercepts"), we encounter an error.

We have also tried to run the models using labels to impose constraints "manually." A model that imposes constraints for the 17 loadings works fine. A model that imposes constraints for the 14 intercepts works fine. However, a model that imposes constraints for the 17 loadings and the 14 intercepts yields an error. We have drilled down on how many parameters can be constrained, and it seems like it's 24. For example, we can run a model that constrains 17 loadings and 7 of the intercepts, or 16 loadings and 8 intercepts, but whenever we run it with 25 or more constraints we encounter the error.

Have you encountered this before, and is there a resolution? Ultimately we'd like to run a model that constrains all the loadings, intercepts, and possible error (co)variances to be equal across groups, and in principle we believe that should be doable.

Thanks,

Roy Levy

###########################################################

### bayesian measurement invariance analysis in blavaan ###

###########################################################

# load packages --------

library(simstandard)

library(blavaan)

# generate a data set with 3 groups for this analysis ------

# total N = 314, group1 n = 106, group2 n = 133, group3 n = 75

# 14 observed variables, 3 factors

group1n = 106

group2n = 133

group3n = 75

observed = 14

# data generating model

modDataGenGroup1 = '

f1 =~ 0.566*x1 + 0.583*x2 + 0.376*x3 + 0.418*x4

f2 =~ 0.523*x5 + 0.520*x6 + 0.591*x7 + 0.638*x8 + 0.462*x9 + 0.619*x10 + 0.248*x11 + -0.062*x12

f3 =~ 0.644*x13 + -0.173*x10 + 0.281*x11 + 0.262*x12 + 0.404*x14

'

modDataGenGroup2 = '

f1 =~ 0.548*x1 + 0.522*x2 + 0.425*x3 + 0.366*x4

f2 =~ 0.485*x5 + 0.530*x6 + 0.562*x7 + 0.521*x8 + 0.438*x9 + 0.496*x10 + 0.073*x11 + 0.169*x12

f3 =~ 0.455*x13 + -0.011*x10 + 0.531*x11 + 0.397*x12 + 0.279*x14

'

modDataGenGroup3 = '

f1 =~ 0.303*x1 + 0.088*x2 + 0.205*x3 + 0.467*x4

f2 =~ 0.639*x5 + 0.580*x6 + 0.659*x7 + 0.500*x8 + 0.456*x9 + 0.533*x10 + 0.077*x11 + 0.025*x12

f3 =~ 0.727*x13 + -0.169*x10 + 0.580*x11 + 0.501*x12 + 0.338*x14

'

# data generation

set.seed(117531)

dataGroup1 = as.data.frame(sim_standardized(modDataGenGroup1, n = group1n))

dataGroup1 = dataGroup1[1:observed]

dataGroup1$group = 1

dataGroup2 = as.data.frame(sim_standardized(modDataGenGroup2, n = group2n))

dataGroup2 = dataGroup2[1:observed]

dataGroup2$group = 2

dataGroup3 = as.data.frame(sim_standardized(modDataGenGroup3, n = group3n))

dataGroup3 = dataGroup3[1:observed]

dataGroup3$group = 3

data = rbind(dataGroup1, dataGroup2, dataGroup3)

# bayesian measurement invariance ---------

# settings ---------

nChains = 2

nWarmup = 500

nBurnin = 0

nItersPerChainAfterWarmupBurnin = 500

nItersTotalPerChain = nWarmup + nBurnin + nItersPerChainAfterWarmupBurnin

# ********************************************************* -------

# try 1: use group.equal in the bcfa function -------

model1 = '

f1 =~ x1 + x2 + x3 + x4

f2 =~ x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12

f3 =~ x13 + x10 + x11 + x12 + x14

x1 ~~ x2

x1 ~~ x14

x3 ~~ x4

x6 ~~ x11

'

# this runs: ------

set.seed(19750122)

fit1 = bcfa(model = model1,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group",

group.equal = c("loadings"))

summary(fit1)

# this runs: ----------

set.seed(19750122)

fit2 = bcfa(model = model1,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group",

group.equal = c("intercepts"))

summary(fit2)

# this does not run: ------

set.seed(19750122)

fit3 = bcfa(model = model1,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group",

group.equal = c("loadings", "intercepts"))

summary(fit3)

# ********************************************************* -------

# try 2: manual constraints -----------

# this runs (manual constraints for all 17 loadings): ------

model2 = '

f1 =~ c(load1, load1, load1)*x1 + c(load2, load2, load2)*x2 + c(load3, load3, load3)*x3 + c(load4, load4, load4)*x4

f2 =~ c(load5, load5, load5)*x5 + c(load6, load6, load6)*x6 + c(load7, load7, load7)*x7 + c(load8, load8, load8)*x8 + c(load9, load9, load9)*x9 + c(load10, load10, load10)*x10 + c(load11, load11, load11)*x11 + c(load12, load12, load12)*x12

f3 =~ c(load13, load13, load13)*x13 + c(load14, load14, load14)*x10 + c(load15, load15, load15)*x11 + c(load16, load16, load16)*x12 + c(load17, load17, load17)*x14

x1 ~~ x2

x1 ~~ x14

x3 ~~ x4

x6 ~~ x11

'

set.seed(19750122)

fit4 = bcfa(model = model2,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group")

summary(fit4)

# this runs (manual constraints for all 14 intercepts): --------

model3 = '

f1 =~ x1 + x2 + x3 + x4

f2 =~ x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12

f3 =~ x13 + x10 + x11 + x12 + x14

x1 ~ c(int1, int1, int1) * 1

x2 ~ c(int2, int2, int2) * 1

x3 ~ c(int3, int3, int3) * 1

x4 ~ c(int4, int4, int4) * 1

x5 ~ c(int5, int5, int5) * 1

x6 ~ c(int6, int6, int6) * 1

x7 ~ c(int7, int7, int7) * 1

x8 ~ c(int8, int8, int8) * 1

x9 ~ c(int9, int9, int9) * 1

x10 ~ c(int10, int10, int10) * 1

x11 ~ c(int11, int11, int11) * 1

x12 ~ c(int12, int12, int12) * 1

x13 ~ c(int13, int13, int13) * 1

x14 ~ c(int14, int14, int14) * 1

x1 ~~ x2

x1 ~~ x14

x3 ~~ x4

x6 ~~ x11

'

set.seed(19750122)

fit5 = bcfa(model = model3,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group")

summary(fit5)

# this runs (manual constraints for all 17 loadings and the first 7 intercepts; 24 constraints total): ----------

model4 = '

f1 =~ c(load1, load1, load1)*x1 + c(load2, load2, load2)*x2 + c(load3, load3, load3)*x3 + c(load4, load4, load4)*x4

f2 =~ c(load5, load5, load5)*x5 + c(load6, load6, load6)*x6 + c(load7, load7, load7)*x7 + c(load8, load8, load8)*x8 + c(load9, load9, load9)*x9 + c(load10, load10, load10)*x10 + c(load11, load11, load11)*x11 + c(load12, load12, load12)*x12

f3 =~ c(load13, load13, load13)*x13 + c(load14, load14, load14)*x10 + c(load15, load15, load15)*x11 + c(load16, load16, load16)*x12 + c(load17, load17, load17)*x14

x1 ~ c(int1, int1, int1) * 1

x2 ~ c(int2, int2, int2) * 1

x3 ~ c(int3, int3, int3) * 1

x4 ~ c(int4, int4, int4) * 1

x5 ~ c(int5, int5, int5) * 1

x6 ~ c(int6, int6, int6) * 1

x7 ~ c(int7, int7, int7) * 1

x1 ~~ x2

x1 ~~ x14

x3 ~~ x4

x6 ~~ x11

'

set.seed(19750122)

fit6 = bcfa(model = model4,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group")

summary(fit6)

# this runs (manual constraints for all 17 loadings and the 2nd to 8th intercepts; 24 constraints total): ----------

model5 = '

f1 =~ c(load1, load1, load1)*x1 + c(load2, load2, load2)*x2 + c(load3, load3, load3)*x3 + c(load4, load4, load4)*x4

f2 =~ c(load5, load5, load5)*x5 + c(load6, load6, load6)*x6 + c(load7, load7, load7)*x7 + c(load8, load8, load8)*x8 + c(load9, load9, load9)*x9 + c(load10, load10, load10)*x10 + c(load11, load11, load11)*x11 + c(load12, load12, load12)*x12

f3 =~ c(load13, load13, load13)*x13 + c(load14, load14, load14)*x10 + c(load15, load15, load15)*x11 + c(load16, load16, load16)*x12 + c(load17, load17, load17)*x14

x2 ~ c(int2, int2, int2) * 1

x3 ~ c(int3, int3, int3) * 1

x4 ~ c(int4, int4, int4) * 1

x5 ~ c(int5, int5, int5) * 1

x6 ~ c(int6, int6, int6) * 1

x7 ~ c(int7, int7, int7) * 1

x8 ~ c(int8, int8, int8) * 1

x1 ~~ x2

x1 ~~ x14

x3 ~~ x4

x6 ~~ x11

'

set.seed(19750122)

fit7 = bcfa(model = model5,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group")

summary(fit7)

# this runs (manual constraints for the first 16 loadings and the first 8 intercepts; 24 constraints total): --------------

model6 = '

f1 =~ c(load1, load1, load1)*x1 + c(load2, load2, load2)*x2 + c(load3, load3, load3)*x3 + c(load4, load4, load4)*x4

f2 =~ c(load5, load5, load5)*x5 + c(load6, load6, load6)*x6 + c(load7, load7, load7)*x7 + c(load8, load8, load8)*x8 + c(load9, load9, load9)*x9 + c(load10, load10, load10)*x10 + c(load11, load11, load11)*x11 + c(load12, load12, load12)*x12

f3 =~ c(load13, load13, load13)*x13 + c(load14, load14, load14)*x10 + c(load15, load15, load15)*x11 + c(load16, load16, load16)*x12 + x14

x1 ~ c(int1, int1, int1) * 1

x2 ~ c(int2, int2, int2) * 1

x3 ~ c(int3, int3, int3) * 1

x4 ~ c(int4, int4, int4) * 1

x5 ~ c(int5, int5, int5) * 1

x6 ~ c(int6, int6, int6) * 1

x7 ~ c(int7, int7, int7) * 1

x8 ~ c(int8, int8, int8) * 1

x1 ~~ x2

x1 ~~ x14

x3 ~~ x4

x6 ~~ x11

'

set.seed(19750122)

fit8 = bcfa(model = model6,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group")

summary(fit8)

# this does not run (manual constraints for all 17 loadings and all 14 intercepts; 31 constraints total): ----------

model7 = '

f1 =~ c(load1, load1, load1)*x1 + c(load2, load2, load2)*x2 + c(load3, load3, load3)*x3 + c(load4, load4, load4)*x4

f2 =~ c(load5, load5, load5)*x5 + c(load6, load6, load6)*x6 + c(load7, load7, load7)*x7 + c(load8, load8, load8)*x8 + c(load9, load9, load9)*x9 + c(load10, load10, load10)*x10 + c(load11, load11, load11)*x11 + c(load12, load12, load12)*x12

f3 =~ c(load13, load13, load13)*x13 + c(load14, load14, load14)*x10 + c(load15, load15, load15)*x11 + c(load16, load16, load16)*x12 + c(load17, load17, load17)*x14

x1 ~ c(int1, int1, int1) * 1

x2 ~ c(int2, int2, int2) * 1

x3 ~ c(int3, int3, int3) * 1

x4 ~ c(int4, int4, int4) * 1

x5 ~ c(int5, int5, int5) * 1

x6 ~ c(int6, int6, int6) * 1

x7 ~ c(int7, int7, int7) * 1

x8 ~ c(int8, int8, int8) * 1

x9 ~ c(int9, int9, int9) * 1

x10 ~ c(int10, int10, int10) * 1

x11 ~ c(int11, int11, int11) * 1

x12 ~ c(int12, int12, int12) * 1

x13 ~ c(int13, int13, int13) * 1

x14 ~ c(int14, int14, int14) * 1

x1 ~~ x2

x1 ~~ x14

x3 ~~ x4

x6 ~~ x11

'

set.seed(19750122)

fit9 = bcfa(model = model7,

data = data,

meanstructure = TRUE,

n.chains = nChains,

burnin = nWarmup,

sample = nItersPerChainAfterWarmupBurnin,

mcmcfile = TRUE,

save.lvs = FALSE,

bcontrol = list(cores = nChains),

seed = 20210117,

group = "group")

summary(fit9)

Ed Merkle

unread,

Sep 11, 2025, 2:55:40 PM (6 days ago) Sep 11

to Roy Levy, blavaan

Hi Roy,

Thanks for the report. This looks like a problem caused by some restructuring of how blavaan handles covariance matrices. I am currently looking into it.

I can't remember whether we have discussed it before, but constraining error covariances can get messy. blavaan does the "SRS" thing where the covariance matrix is decomposed into a diagonal matrix of standard deviations and a correlation matrix. So constraining a covariance to be equal actually involves constraining two parameters to be equal: the residual variance (or SD) and the residual correlation. I am not yet sure whether this is also impacted by the error you found.

Ed

--
You received this message because you are subscribed to the Google Groups "blavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to blavaan+u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/blavaan/ebc44944-0800-42a4-8d89-704ebedfc3f5n%40googlegroups.com.

Roy Levy

unread,

Sep 12, 2025, 2:33:24 PM (5 days ago) Sep 12

to blavaan

Hi Ed,

Thanks for attending to this. I imagine the covariance matrix decomposition might affect lots of things, the imposition of constraints being one. In the models we've run, we hadn't (yet) tried constraining covariances. And I've just tried running a version of the code I sent (model1) that doesn't include the residual covariances at all, so the only covariances are among the 3 latent variables, and they are not constrained across group. I still ran into the error when constraining all the loadings and intercepts (fit3 in the code). And I realize I didn't put the exact error in my initial post. Here is the specific error (though if the code I posted runs for you I take it you saw it):

Error in eval(mc, parent.frame()) : Exception: model_stanmarg_namespace::model_stanmarg: v10[2][4] is 0, but must be greater than or equal to 1.000000 (in 'string', line 629, column 2 to column 38) failed to create the sampler; sampling not done

Again, thanks for thinking about and looking into this!

Roy

Ed Merkle

unread,

Sep 13, 2025, 4:42:54 PM (4 days ago) Sep 13

to Roy Levy, blavaan

Hi Roy,

This error should now be fixed on github. Some floating point thing was happening where a covariance fixed to 0 was actually equal to something like 1e-17, which led to unexpected behavior.

And if you move toward equality constraints on covariances, we might talk again. In addition to decomposing into correlations and SDs, I have been working on a method where you try to make the covariance matrix block diagonal so that a prior can be assigned to each individual block. It is described in a bit of detail here:

https://ecmerkle.github.io/cs/blockmat.html

But this adds an additional layer of complication for equality constraints in covariance matrices, and I believe you will see an error if you try to specify the equality constraints. One idea here is to try target = "jags", which does not use this new functionality.

Ed

To view this discussion visit https://groups.google.com/d/msgid/blavaan/bc1931c4-e307-4e3d-8813-5fc81b1ce6e7n%40googlegroups.com.

Roy Levy

unread,

Sep 16, 2025, 5:37:49 PM (23 hours ago) Sep 16

to blavaan

Hi Ed,

Yep, it seems fixed on my end too. Thank you for looking into it!

Re: fixing covariances, the work we're doing here might move (slowly) in that direction, and I will reach out to discuss as warranted. And thanks for pointing me to your work and discussion on making covariance matrices into block diagonal. It seems like that would make for a step forward in allowing for more flexibility to specify prior on covariance structures.