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.
It's not really a hard and fast rule, it's just that if you only have two first order factors, your model will be equivalent to a model where they are correlated. Your point 2 is correct.
It's just like if you have two measured variables, and want to fit a single factor model. If you have two variables (latent or measured), you have one covariance / correlation to account for - so you have only got 1 degree of freedom.
Say that these two variables are correlated 0.64. You want to estimate two loadings, but you only have 1 df, so you can only estimate one parameter. The way to fix that is to constrain the two loadings to be equal. If you do that, the loadings will be sqrt(0.64) = 0.8.
So (to address point 3) you can fit the model, but it won't tell you anything you didn't already know, and you won't learn anything from it.
(I'd argue that three first order factors doesn't tell you something - the factor loadings are a transformation of the factor covariance matrix, and are completely determined - you need four factors to over-identify the second order factor).
Best Answer
The same reason a simple 1-factor model with 3 indicators is not identified if you try to estimate all loadings and residual variances. Its df=0 even without residual covariances, so trying to estimate even 1 underidentifies the model.
Then there is no reason to expect their residuals to correlate, after the covariances are accounted for by g. SEM parameters are residual, not marginal, (co)variances, unless the factor(s) happen to be exogenous. In your model, the 3 dimensions are endogenous, explained by g.