Comparing lavaan::sem to probit regression output

Question

I am trying to improve my understanding of lavaan::sem models when using a probit link function by comparing the output to simple probit regressions.

Currently, I see that the coefficients for each approach are approximately the same, but the regression intercept is negative, and the sem value is positive - although they are also approximately the same value otherwise.

Why would each approach have a positive or negative value? Is there a different interpretation for the sem value?

df = structure(list(y = c(1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 
                        0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 
                        1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 
                        0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 
                        1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 
                        1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 
                        0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 
                        0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 
                        1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 
                        0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 
                        0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 
                        0L, 1L, 0L, 1L, 0L, 1L), 
               x = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                     1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 
                     0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 
                     1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
                     0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 
                     0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
                     0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 
                     1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 
                     0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 
                     1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
                     0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 
                     0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L)), class = "data.frame", row.names = c(NA, -176L))

fit.glm = glm(y ~ x, data = df, family = binomial(link = probit))
fit.sem = lavaan::sem("y ~ x", data = df, ordered = TRUE)

> coef(fit.glm)
(Intercept)           x 
 -0.7197189   2.0012704 
> coef(fit.sem)
  y~x  y|t1 
2.001 0.720 

Answer 1

Oh, this is a fun topic! First, in case you aren't already familiar with the latent-response interpretation of a probit model, you can read a bit about it on Wikipedia:

https://en.wikipedia.org/wiki/Probit_model#Conceptual_framework

In that page's formulas, notice that the intercept (in $\beta$, pre-multiplied by $X$) and threshold (represented only as the number $0$ when assigning scores to category 0 or 1) are distinct parameters.

Currently, I see that the coefficients for each approach are approximately the same, but the regression intercept is negative, and the sem value is positive

Actually, the intercept (y~1) is 0 in the lavaan model:

> coef(fit.sem, type = "user") # prints all parameters
  y~x  y|t1  y~~y  x~~x y~*~y   y~1   x~1 
2.001 0.720 1.000 0.241 1.000 0.000 0.398 

> # intercept
> coef(fit.sem, type = "user")["y~1"]
y~1 
  0 

> # threshold
> coef(fit.sem)["y|t1"]
     y|t1 
0.7197189 

It is the threshold (y|t1) that is the same absolute value as the intercept that is estimated by glm(). This is because the probit model fitted by glm() is identified by fixing the threshold to 0. Thus, the interpretation of the intercept (latent location parameter) is the same under both parameterizations:

• On the (arbitrary) latent-response scale, the intercept is 0.72 units below the threshold.

But you cannot simultaneously estimate both the intercept and the threshold (they are indeterminate). So you can only identify one estimate by fixing the other.

You can make lavaan's parameterization equivalent to the GLM by changing the default parameters.

## bare-bones lavaan() doesn't automatically free anything,
## so auto.th=FALSE (see ?lavOptions)
mod2 <- ' y ~ x 
          y | 0*t1  # fix threshold to 0
          y ~ NA*1  # free the intercept
'
fit2 <- sem(mod2, data = df, ordered = TRUE)
coef(fit2, type = "user")
#   y~x   y|t1    y~1   y~~y   x~~x  y~*~y    x~1 
# 2.001  0.000 -0.720  1.000  0.241  1.000  0.398 

The latent-response scale (y ~*~ y) is also arbitrarily fixed (to 1 by default) for identification, but it can be freely estimated for ordinal data if you fix 2 thresholds (e.g., to 0 and 1). We can even do so with binary data by arbitrarily choosing to set the distance between the threshold and intercept such that they are 1 "unit" apart (just as arbitrary, but bear with me :-) Since we already peaked at the data, we know the intercept is below the threshold, so let's set the latent intercept to 0 and the threshold to 1. Thus, the variance of the latent responses ($y^*$) underlying $y$ can be determined relative to that distance (which is 0.72 when we fix the latent scale to 1):

mod3 <- ' y ~ x 
          y ~ 0*1   # fix intercept to 0
          y | 1*t1  # fix threshold to 1
          y ~*~ y   # free latent scaling factor (1 / SD)
'
fit3 <- sem(mod3, data = df, ordered = TRUE)
coef(fit3, type = "user")
#   y~x   y~1  y|t1 y~*~y  y~~y  x~~x   x~1 
# 2.781 0.000 1.000 0.720 1.931 0.241 0.398 

> 1 / coef(fit3, type = "user")["y~*~y"] # SD of y*
   y~*~y 
1.389431 

Now, the interpretation is that the intercept is 1 "unit" below the threshold, but the units are arbitrarily "larger" than the models above (when the intercept is only 0.72 "units" below the threshold). But don't fret: although the units are arbitrary, the standardized solution (in units of SD) will be the same. So let's interpret the fit3 results in units of the latent-response's ($y^*$) $SD = \frac{1}{0.72} = 1.39$:

• The 1-unit difference between the intercept and threshold is a standardized difference of $\frac{1}{SD} = 0.72$ units.

This is the same result we saw from fit2 and your OP (when the latent scale was already $SD=1$).