Calculating path coefficients manually from correlations

Question

As I understand it, in a path model, we can use Wright's approach to manually calculate (standardized) path coefficients from the correlations between the variables. In the saturated case, we perfectly reproduce the correlations, while in the non-saturated one, we get a potential difference between model-implied correlations and true correlations (if the paths set to zero are not zero in reality). This is the basis for the calculation of "model fit".

When I run a path model e.g. with lavaan in R, the program goes through several "iterations" in order to fit the best model. However, if we can manually calculate the coefficients from the correlations, why do we need an optimization procedure that goes through several iterations? I'd be grateful if someone could explain.

Answer 1

Indeed manual calculation using Wright’s method is limited to simple saturated models where there're enough paths (free parameters) to perfectly reproduce all correlations among the variables. But in most real-world applications your model is not saturated like your example where the paths set to zero are not zero in reality or where the specified latent constructs cannot fully explain the observed correlations. This means the model cannot perfectly reproduce all observed correlations and the model-implied correlations will differ from the observed ones.

Therefore iterative optimization is necessary to find the parameters that minimize the difference between the observed correlations and the model-implied correlations. This process often uses a maximum likelihood or similar estimation method.

Answer 2

Lavaan doesn't know that the path model is saturated - so it fits the model iteratively. (Even if it knew, I don't think it's been programmed to do that).

There's nothing wrong with fitting it iteratively - all your regression models could be fit this way, and they'd get the same answer.

For example, consider a mediation model, with two predictors, x1 and x2. Two mediators (m1 and m2) are regressed on x1 and x2, and y is regressed on all four variables. This model is saturated, and so we can estimate it based with the regular regression equations.

set.seed(42)
d <- data.frame(
  x1 = rnorm(1000),
  x2 = rnorm(1000),
  m1 = rnorm(1000),
  m2 = rnorm(1000),
  y  = rnorm(1000)
)

summary(lm(y ~ x1 + x2 + m1 + m2, data = d))
summary(lm(m1 ~ x1 + x2, data = d))
summary(lm(m2 ~ x1 + x2, data = d))

Or we could estimate it using lavaan.

library(lavaan)

sem_model = '
  y ~ m1 + m2 + x1 + x2
  m1 ~ x1 + x2
  m2 ~ x1 + x2
'

sem_fit <- lavaan::sem(
  model,
  data = d
  )
summary(sem_fit)

Lavaan took 11 iterations to get to a solution. But compare the lavaan and lm() solutions:

Here are the lm() solutions:

Call:
lm(formula = y ~ x1 + x2 + m1 + m2, data = d)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0152039  0.0322527  -0.471    0.637
x1          -0.0007105  0.0321810  -0.022    0.982
x2           0.0124892  0.0327703   0.381    0.703
m1          -0.0204112  0.0313420  -0.651    0.515
m2           0.0294935  0.0326551   0.903    0.367

Call:
lm(formula = m1 ~ x1 + x2, data = d)
Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.003867   0.032589  -0.119    0.906
x1          -0.021799   0.032514  -0.670    0.503
x2          -0.018288   0.033057  -0.553    0.580

Call:
lm(formula = m2 ~ x1 + x2, data = d)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.02110    0.03128  -0.675   0.5000  
x1           0.01413    0.03121   0.453   0.6509  
x2           0.05984    0.03173   1.886   0.0596 .
---

And here's lavaan:

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  y ~                                                 
    m1               -0.020    0.031   -0.653    0.514
    m2                0.029    0.033    0.906    0.365
    x1               -0.001    0.032   -0.022    0.982
    x2                0.012    0.033    0.382    0.702
  m1 ~                                                
    x1               -0.022    0.032   -0.671    0.502
    x2               -0.018    0.033   -0.554    0.580
  m2 ~                                                
    x1                0.014    0.031    0.453    0.650
    x2                0.060    0.032    1.889    0.059

The results are pretty close to identical. We could have coded the mediation model to estimate the parameters using a closed form equation, which would take zero iterations, and therefore presumable be a little faster. But the amount would be minimal.

Sometimes there's a tool that's good enough, and if you don't have the ideal tool, you can use the one that's good enough.