# Why do matrices of simulation data generated in R create strongly related latent variables?

When creating examples of matrices for CFA in R for students, I still have a problem with getting data that looks realistic - specifically, that between latent variables I often have covariance path coefficients greater than 1 (standardized!) which just seems unrealistic.

Let's consider an example:

library(lavaan)
library(simstudy)
library(lavaanPlot)

samples = 1000
number_of_items = 21
number_of_factors = 3
correlation = 0.5

# generate "equal" array of correlations
sigma <- array(rep(correlation, number_of_items*number_of_items),
dim = c(number_of_items,number_of_items))
diag(sigma) <- 1
# create the mean vector
mu <- rep(0, number_of_items)

# for uniform distribution params1 = min, params2 = max
# (which is the square of the standard deviation)
df = genCorGen(samples, nvars = number_of_items, params1 = mu,
params2 = (mu+1), dist = "uniform", rho = correlation,
corstr = "cs", wide = TRUE)

# prepare dataframe for work
df = as.data.frame(df[,-1])
colnames(df) <- do.call('paste0', expand.grid("q",
c(1:number_of_items)))

# let's make model with three latent variables, each with 7 items
model <- "f1 =~ q1 + q4 + q7 + q10 + q13 + q16 + q19\nf2 =~ q2 + q5 +
q8 + q11 + q14 + q17 + q20\nf3 =~ q3 + q6 + q9 + q12 + q15 +
q18 + q21"
# take note that is just one of every possibility,
# below is comment in text

fit <- cfa(model, data=df)
# aaaand we have strange standarized coefficients, just look:
lavaanPlot(model = fit, coefs=T, covs=TRUE,
stars = c("covs", "latent", "regress"), stand=TRUE)


And now we have latent variables related to each other at a standardized coefficient of 1 or greater:

Importantly, the selected model is only one of the possibilities of specifying the model, but the same problem will occur for most of the possible model arrangements - you can increase/decrease the number of items, you can increase/decrease the number of latent variables, you can swap them, you can change the number of cases, or even transpose variables to Likert scales (as for CFA standardized questionnaires) and the problem will still be the same. (Of course, the model can be specified so that not all latent variables are related to each other, but I am interested in the validation of measurement tools that usually assume that relationships between all latent variables should be checked.)

I've also checked these fit factors in SPSS AMOS and other programs and they are exactly the same, so I suspect it's not Lavaan's fault, just some kind of data problem.

So I have questions:

1. What does it mean for this specific data that the path coefficients are greater than or very close to 1? In many places on internet forums there are advice to "explore data", but here we all have simulation data ready and no matter how many times the process is repeated (draw a new data matrix) the problem is still the same.
2. Why does generated data behave this way and real data collected from research do not? I mean any database in which there is a matrix of the same number of variables with a similar average level of correlation. (I'm not attaching the data because it's not really possible here). So - at first glance, we have similar data, the same number of cases, the same variable ranges, the same average correlations between items, and in the case of real data, this problem practically does not occur to me.
3. Does this mean that the simulated data is somehow different from the real data and that the simulations performed on them are not realistic?
• My first blush thought was that you generated with uniform distribution, and the real world is a lot sparser and more "peaky". Uniform is a way of saying "they are all equal". This could be that you have enough samples that random correlations look very strong. Jun 29, 2023 at 13:39
• True, but you'll get the same with other data generation methods, regardless of what type of distribution is used to generate the data. When generating data using simulators in R (and so on for other packages as well), you really only have two options: normal distribution and homogeneous distribution. All others you have to formulate yourself in the form of some transformations. Jun 30, 2023 at 9:14

The easiest way to avoid this problem is probably to simulate data directly from a pre-specified population CFA model in which all parameters are "proper" (e.g., all factor correlations are smaller than |1|). Have you considered using the simsem package (https://github.com/simsem/simsem/wiki/Example-1:-Getting-Started)? It allows you to simulate data directly from a population model where you can set the parameters to plausible/meaningful/real-world values and then simulate data from the model. The improper-value problem that you encountered should not occur too frequently that way, as long as you simulate large enough samples from population models in which the parameters are not too close to their theoretical boundary values.

• It's an interesting approach. Within simulation research, however, this is the opposite approach to what I'm interested in - with simsem we create a model that simulate the data instead of simulating the data to discover the model. I am curious, however, why simulation models in R use this algorithm (simstudy or MASS) and rather do not pay attention to it? And most importantly - how to understand this problem? What exactly does it come from? I don't want to avoid the problem, but to understand what causes it. Jun 30, 2023 at 20:27

I'm not familiar with GenCorGen(), but you specify the number of factors at the start of your script and you never use it again.

It looks like you are generating data from a single factor model, so your three latent variables would be correlated 1.0. 1.01 is sampling variability.

Edit: The data are samples from a population where the latent correlations are all equal to 1.0. Real data don't do this, because you would never find a population with three latent variables correlated 1.0.

Here's how I would do it (in case it helps) to generate 2 latent variables with 4 items loading on each, with loadings of 0.7 (lambda) and latent variables correlating 0.5:

library(dplyr)
library(MASS)

lambda <-
matrix(
c(0, 0.7,
0, 0.7,
0, 0.7,
0, 0.7,
0.7, 0,
0.7, 0,
0.7, 0,
0.7, 0),
nrow = 8, ncol = 2, byrow = TRUE)

phi <- matrix(
c(1, 0.5, 0.5, 1), nrow = 2, byrow = TRUE
)

s <- lambda %*% phi %*% t(lambda)
diag(s) <- 1
df <- MASS::mvrnorm(n = 1000, mu = c(rep(0, 8)), Sigma = s, empirical = FALSE) %>%
as.data.frame()
cor(df)



Which gives me:

          V1        V2        V3        V4        V5        V6        V7
V1 1.0000000 0.5125808 0.5372649 0.5159997 0.3191535 0.3142949 0.3152547
V2 0.5125808 1.0000000 0.4805850 0.4994852 0.2458347 0.2425625 0.2631177
V3 0.5372649 0.4805850 1.0000000 0.4914618 0.2038897 0.2055747 0.2326917
V4 0.5159997 0.4994852 0.4914618 1.0000000 0.2448801 0.2397230 0.2487003
V5 0.3191535 0.2458347 0.2038897 0.2448801 1.0000000 0.4963736 0.4996688
V6 0.3142949 0.2425625 0.2055747 0.2397230 0.4963736 1.0000000 0.5200427
V7 0.3152547 0.2631177 0.2326917 0.2487003 0.4996688 0.5200427 1.0000000
V8 0.3090881 0.2620999 0.2403555 0.2712580 0.4830959 0.4939030 0.5149979
V8
V1 0.3090881
V2 0.2620999
V3 0.2403555
V4 0.2712580
V5 0.4830959
V6 0.4939030
V7 0.5149979
V8 1.0000000

• 1/2 This line only informs reader about the model specification and is meaningless. I do not generate a model based on one factor, but a matrix with given correlations. Within the package, there is no distinction between the number of latent factors, only the correlation matrix is established. Your suspicion would mean that at the level of the algorithm, there is simply a projection of the latent variable onto the data matrix. Jun 30, 2023 at 20:34
• 2/2 But I'm more interested in whether it is possible to distinguish a one-dimensional matrix of randomly generated numbers with such parameters from a multidimensional matrix (as far as I know, this is impossible). Also, the latent variable problem is not solvable, is it? What do you mean by "1.01 is sampling variability"? Jun 30, 2023 at 20:34
• Sorry, I don't quite get what you're trying to do. You generate a matrix of equal correlations. That fits a one factor model. Then you fit a three factor model - the population values of the latent correlations is 1.0. However, because samples vary, sometimes these are too high, and sometimes too low. I'll edit my answer too. Jun 30, 2023 at 21:05

This is not exactly a simstudy issue, because GenCorGen is generating a set of correlated data as you wanted. Here I am just generating 9 variables to simplify, and you can see that the desired correlation structure is more or less recovered:

library(simstudy)

samples = 1000
number_of_items = 9
correlation = 0.5

mu <- rep(0, number_of_items)

# for uniform distribution params1 = min, params2 = max
# (which is the square of the standard deviation)
df = genCorGen(samples, nvars = number_of_items, params1 = mu,
params2 = (mu+1), dist = "uniform", rho = correlation, corstr = "cs", wide = TRUE)

df = as.data.frame(df[,-1])
round(cor(as.matrix(df)), 2)
#>      V1   V2   V3   V4   V5   V6   V7   V8   V9
#> V1 1.00 0.50 0.51 0.48 0.47 0.48 0.48 0.50 0.48
#> V2 0.50 1.00 0.49 0.52 0.49 0.48 0.50 0.51 0.49
#> V3 0.51 0.49 1.00 0.48 0.47 0.48 0.47 0.52 0.52
#> V4 0.48 0.52 0.48 1.00 0.47 0.49 0.50 0.50 0.49
#> V5 0.47 0.49 0.47 0.47 1.00 0.50 0.49 0.50 0.46
#> V6 0.48 0.48 0.48 0.49 0.50 1.00 0.48 0.48 0.48
#> V7 0.48 0.50 0.47 0.50 0.49 0.48 1.00 0.48 0.46
#> V8 0.50 0.51 0.52 0.50 0.50 0.48 0.48 1.00 0.51
#> V9 0.48 0.49 0.52 0.49 0.46 0.48 0.46 0.51 1.00



I agree with one of the responses above that it might be more direct to simulate the latent factors with some correlation structure, and then generate the observed data based on a second correlation structure. It is a little tricky to generate uniform data this way, so I understand why you opted to use GenCorGen, but by doing this, you completely lose the ability to generate correlations between latent factors that are less than 1. (Incidentally, when I use lavaan to fit the model with your example, I mostly get warning messages that the covariance matrix of latent variables is not positive definite; in those cases, the estimated correlations can exceed 1. However, without that warning, the correlations do not exceed 1.)

Here is a simple example that generates the latent data and then the observed data.

library(lavaan)
#> This is lavaan 0.6-15
#> lavaan is FREE software! Please report any bugs.

options(digits = 3)

d1 <- defDataAdd(varname = "q1", formula = "f1", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q2", formula = "f1", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q3", formula = "f1", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q4", formula = "f2", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q5", formula = "f2", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q6", formula = "f2", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q7", formula = "f3", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q8", formula = "f3", variance = "..v", dist = "normal")
d1 <- defDataAdd(d1, varname = "q9", formula = "f3", variance = "..v", dist = "normal")

samples = 100000
number_of_factors = 3

sigma <- matrix(c(1, .4, .2, .4, 1, .6, .2, .6, 1), nrow = 3)

mu <- rep(0, number_of_factors)

df <- genCorGen(samples, nvars = number_of_factors, params1 = mu,
params2 = 1, dist = "normal", corMatrix = sigma, wide = TRUE,
cnames = c("f1", "f2", "f3"))
df
#>             id      f1     f2      f3
#>      1:      1 -0.0235  1.024  0.6585
#>      2:      2  1.3679  1.398  1.2690
#>      3:      3 -0.9396 -0.574 -1.3512
#>      4:      4 -0.4710  1.848 -0.0812
#>      5:      5  0.8711  0.282 -0.4162
#>     ---
#>  99996:  99996 -0.3487 -0.243 -0.5527
#>  99997:  99997  0.0659  1.103  1.5921
#>  99998:  99998  0.1843 -0.409  0.8698
#>  99999:  99999 -0.3679  1.248  0.9249
#> 100000: 100000 -0.9233 -0.281  0.4592

v <- 1.5
df
#>             id      f1     f2      f3     q1     q2      q3     q4      q5
#>      1:      1 -0.0235  1.024  0.6585 -0.260 -2.398 -0.9547  1.763 -0.0966
#>      2:      2  1.3679  1.398  1.2690  2.223  2.648  0.6692  1.876  0.9024
#>      3:      3 -0.9396 -0.574 -1.3512 -2.763 -1.736 -1.3478  0.764 -0.1768
#>      4:      4 -0.4710  1.848 -0.0812 -2.290 -1.828 -0.4864  0.176  2.5500
#>      5:      5  0.8711  0.282 -0.4162  0.296  2.203 -1.0675 -3.297  1.0174
#>     ---
#>  99996:  99996 -0.3487 -0.243 -0.5527 -1.585  1.334  0.0122  2.295  1.0039
#>  99997:  99997  0.0659  1.103  1.5921 -1.884  2.402  0.8528 -1.747  1.0097
#>  99998:  99998  0.1843 -0.409  0.8698  1.248  0.744  1.1510 -1.342  1.5825
#>  99999:  99999 -0.3679  1.248  0.9249 -2.404  0.838  0.7844 -1.075  0.6437
#> 100000: 100000 -0.9233 -0.281  0.4592 -1.927 -2.694  1.3322  0.585  0.6662
#>             q6     q7     q8      q9
#>      1:  0.879  1.804 -2.132  0.7788
#>      2:  2.675  0.811  1.468  0.0527
#>      3: -1.771 -1.300 -0.213 -2.1955
#>      4:  2.673 -0.588  0.499  0.2840
#>      5: -0.644 -1.921  0.933 -3.0063
#>     ---
#>  99996:  0.280 -0.366 -1.085 -1.5904
#>  99997:  0.293  2.886 -0.173  0.8924
#>  99998: -0.992  4.012  0.893 -1.0576
#>  99999:  3.087 -0.104 -0.503 -0.0809
#> 100000: -2.272  1.247  0.427 -1.1628

df <- df[,-(1:4)]
cor(as.matrix(df))
#>        q1     q2     q3    q4    q5    q6     q7     q8     q9
#> q1 1.0000 0.4007 0.3994 0.163 0.165 0.164 0.0805 0.0805 0.0812
#> q2 0.4007 1.0000 0.4006 0.163 0.159 0.162 0.0749 0.0773 0.0721
#> q3 0.3994 0.4006 1.0000 0.165 0.159 0.164 0.0820 0.0801 0.0781
#> q4 0.1626 0.1635 0.1652 1.000 0.401 0.401 0.2478 0.2438 0.2406
#> q5 0.1654 0.1594 0.1590 0.401 1.000 0.402 0.2432 0.2465 0.2403
#> q6 0.1638 0.1616 0.1642 0.401 0.402 1.000 0.2370 0.2394 0.2434
#> q7 0.0805 0.0749 0.0820 0.248 0.243 0.237 1.0000 0.4020 0.4004
#> q8 0.0805 0.0773 0.0801 0.244 0.247 0.239 0.4020 1.0000 0.3987
#> q9 0.0812 0.0721 0.0781 0.241 0.240 0.243 0.4004 0.3987 1.0000

model <- "f1 =~ q1 + q2 + q3\nf2 =~ q4 + q5 + q6\nf3 =~ q7 + q8 + q9"

fit <- cfa(model, data=df)


And we can see that the between-factor correlations were recovered. You can control the correlation of observed values and the factors with v.

This may not be exactly what you are aiming for, but it might be closer. And yes, the simsem package might be a better fit.