Your model is not identified. That means that Lavaan cannot find a unique solution, and cannot compute standard errors. Because it can't find a unique solution, it's finishing in different places.
[Very simple example of identification:
$x + y = 5$
There are many (an infinite number) of solutions to this equation for x and y, and they're all equally good, so it's not identified.]
This has happened for a couple of reasons.
First:
x1 ~~ x2
Is not going to work. You need the latent variable to explain the covariance of these two items. You can't have a covariance between the two variables, because there's nothing left for the latent variable to explain. So let's get rid of that. Oh, still problematic.
Second, your sample size is small. Nothing we can do about that.
Third, the variances of your variables are quite different. This causes convergence problems for most SEM programs. A rule of thumb is to make sure your SDs are less than an order of magnitude different. Let's check:
> apply(dat, 2, na.rm=TRUE, sd)
s v x1 x2
2.686326 2.940680 9.397064 43.761456
OK, let's fix that.
dat$x1 <- dat$x1 / 10
dat$x2 <- dat$x2 / 10
Hmmmm... still doesn't work.
Using FIML probably doesn't make life any easier for Lavaan. Let's turn that off.
Nope, still not working.
Now we're getting a bit desperate. We need some additional restrictions in the model. You have two indicators of a latent variable - that's not many. Let's fix the loadings to be equal, that makes it easier to converge.
That worked, we have convergence. But we also have a negative error variance on x1 (a Heywood case). Is it just a little bit negative? If so, we'll ignore it. Nope, it's not. It's big. OK, we need to fix that.
Again, because we're desperate, we're going to fix the error variance of x1 and x2 to be equal.
x1 ~~ a*x1
x2 ~~ a*x2
That gives us convergence, but horrible fit. What about fixing it to zero?
x2 ~~ 0*x2
Nope, still horrible fit. Let's go back a step and see if we can work out what's going on (which is worth starting with). Let's look at the correlations.
> cor(dat, use="pairwise.complete.obs")
s v x1 x2
s 1.0000000 0.3246869 0.1811305 0.0808834
v 0.3246869 1.0000000 0.1479880 0.1192254
x1 0.1811305 0.1479880 1.0000000 0.9206230
x2 0.0808834 0.1192254 0.9206230 1.0000000
OK, x1 and x2 are very highly correlated, so they are going to have similar loadings on M, and if they have similar loadings on M, they need to have similar relationships with the other variables (because M is explaining their relationship with those other variables).
The correlation of S with x1 is about twice as high as the correlation of s with s2. This is the crux of the problem. Somehow we need to account for that. The two paths are:
s -> m -> x1
s -> m -> x2
The first path needs to get a value about twice as high as the second. It can't do it from the s->m path, because there's only one of that. So to account for the correlation difference it wants to have different loadings. But it can't have different loadings, because the correlation is so high between x1 and x2.
You have (as I see it) three possible solutions.
- Drop one of x1 or x2. Then you no longer have a latent variable. (But you hardly have a latent variable anyway, with only two indicators.
- Add a correlation between either x1 or x2 with S. That probably doesn't make theoretical sense.
Lose the latent, and have x1 and x2 as mediators of the relationship between s and v directly, so your model is:
v ~ s + x1 + x2
x1 ~ s
x2 ~ s
x1 ~~ x2
Which is a regular multiple mediator model.
Best Answer
Rules of thumb may say that multicollinearity is a problem only if two variables correlate above, say, .9 or even more. If two of your latent variables correlated that much, or even in the range of .7 / .8, then you have a problem before it comes to predicting the third variable: Your measurement model seems to be not very well defined. Maybe, for example, the two latents would be better modeled as only one? I would care much more about the measurement model then about multicollinearity.