Data simulation in R for Measurement invariance I am trying to generate random data for an analysis of Measurement invariance. I will be looking at a one-factor model for two groups with 4 variables each.
I would need to generate a sample where I can control several parameters:

*

*Sample size

*number of categories (I will be using categorical variables)

*distribution of answers (skewness, kurtosis)

*Factor loadings and intercepts/thresholds (To induce measurement noninvariance)

While I can generate data based on the first three conditions I am struggling to find a way to control the fourth. So far, I have only found resources in Mplus that would cover this type of simulation or descriptive mentions in papers such as:

Measurement noninvariance was simulated by creating group differences in the factor loadings for some observed indicators. For the invariance condition, the difference in factor loadings between the groups was 0 (complete invariance). For the noninvariant cases, loadings were simulated to differ by 0.1, 0.2, 0.3, 0.4, or 0.5. The percent of indicators allowed to be noninvariant was 0, 10, 20, or 30%. As an example of how noninvariance was simulated, in the 10 indicators, 10% noninvariant, 0.1 noninvariance magnitude condition, the factor loading for indicator 2 was set to 0.9 in one group, and kept at 1.0 in the other group.
https://doi.org/10.1007/978-3-030-01310-3_21

I would be happy with any resources or even keywords that would guide me towards a better understanding of such simulations and how to practically implement them in R.

Simulating data with measurement invariance is easy with the lavaan package. See the R code below. To simulate measurement invariance, I simply specify two separate population models. Loadings for items 1-3 are the same in both population models. However, the loading for item 4 is different (0.5 for group 1 and 0.4 for group 2).
# Packages
library(lavaan)

# specify population model for group 1 and group 2
population_model_g1 <-'f1=~x1+.6*x2+.7*x3+0.5*x4'
population_model_g2 <-'f1=~x1+.6*x2+.7*x3+0.4*x4'

# Sumulate data
set.seed(1234)
model_data_g1<-lavaan::simulateData(population_model_g1,sample.nobs=500)
model_data_g2<-lavaan::simulateData(population_model_g2,sample.nobs=500)

# Add group label
model_data_g1$group <- "Group1"
model_data_g2$group <- "Group2"

# put data into one data frame
model_data_all <- cbind(model_data_g1,model_data_g2)

