# Multiple regression in SEM when the covariance between predictors is fixed to zero

I did a multiple regression with $$x_1$$ and $$x_2$$ predicting $$y$$, and found that I could recreate in SEM software (lavaan or AMOS) the same results I got doing things the regular way in R with the lm function. I found I can obtain the standardized coefficients produced by the SEM just by calculating for example

$$\beta_1 = b_1 \times \frac{sd(x_1)}{sd(y)}$$.

Then I did a second SEM in which I forced the covariance between $$x_1$$ and $$x_2$$ to be zero.

I found that some things changed:

• The standardized coefficients changed, even though the unstandardized coefficients stayed the same. Therefore I can't calculate the standardized coefficients from the unstandardized coefficients using the equation given above.
• $$R^2$$ declines from 0.7609 to 0.721. I can calculate where these values come from in terms of the equation $$R^2 = \beta_1^2 + \beta_2^2 + 2 \times \beta_1 \times \beta_2 \times r(x_1,x_2)$$, where that final bit obviously cancels out if the correlation between the two predictors is fixed to 0. But I don't really understand why the correlation between $$x_1$$ and $$x_2$$ should be relevant to $$R^2$$ in the first place.
• lavaan now seems to be estimating variances and intercepts for $$x_1$$ and $$x_2$$, when it didn't before.

Some things I'm wondering:

• Why are the standardized coefficients different when I forced the covariance between $$x_1$$ and $$x_2$$ to be zero, when the unstandardized coefficients remained the same? How can I get to those standardized coefficients from the unstandardized ones?
• Why did $$R^2$$ decline when i forced the covariances between the predictor variables to be zero? Is it the case that in a multiple regression set-up of this sort deleting the covariances between predictors will always make $$R^2$$ decline?

R code for the analysis I did is below.

library(tidyverse)

# Get some data for the regression
data <- mtcars %>% dplyr::select(
x1=drat,
x2=wt,
y=mpg
)

# Doing things using the lm function
model <- lm(y ~ x1 + x2,data=data)
summary(model)

# First SEM model recreates a typical multiple regression

lavaan_m1 <- 'y ~ x1 + x2'

mlav_m1 <- sem(lavaan_m1,data=data,meanstructure=TRUE)
inspect(mlav_m1, 'r2')
param_estimates_m1 <- parameterEstimates(mlav_m1)
stand_estimates_m1 <- standardizedSolution(mlav_m1)

# In the first model I can easily generate the standardized estimates
# For example, here's the standardized coefficient for y ~ x1
param_estimates_m1[1,4] * (sd(data$$x1)/sd(data$$y))

# Second SEM model I fix it so correlate between predictors is zero

lavaan_m2 <- 'y ~ x1 + x2
x1 ~~ 0*x2'

mlav_m2 <- sem(lavaan_m2,data=data,meanstructure=TRUE)
# R-Squared has now changed
inspect(mlav_m2, 'r2')
# Unstandardized parameter estimates are all the same,
# except x1 ~ x2 which as expected is now 0.
# The standard errors are different, though.
param_estimates_m2 <- parameterEstimates(mlav_m2)
# Standardized parameter estimates have changed substantially
stand_estimates_m2 <- standardizedSolution(mlav_m2)

# Now my calculation to recreate the standardized coefficients
# doesn't work anymore
param_estimates_m2[1,4] * (sd(data$$x1)/sd(data$$y))

# The R-Squared formulas work out exactly correctly
stand_estimates_m1[2,4]^2 + stand_estimates_m1[3,4]^2 + (2 * stand_estimates_m1[2,4] * stand_estimates_m1[3,4] * cor(data$$x1,data$$x2))
print(inspect(mlav_m1, 'r2'),nd=7)
stand_estimates_m2[1,4]^2 + stand_estimates_m2[2,4]^2
print(inspect(mlav_m2, 'r2'),nd=7)

• If I had to guess what's happening, it's because SEM is based on the covariance matrix implied by the model. You constrained the covariance between predictors to 0, this affects all other parameters in the covariance matrix implied by this model. Everything else follows from this covariance matrix: coefficients, R2, ... I could check this out in a few hours. It should be possible to find the model implied covariance matrix. – Heteroskedastic Jim Jun 19 at 15:37
• Thanks! Yes, I did consider that, and found that lavInspect(mlav_m2,"cov.ov") and lavInspect(mlav_m2,"cor.ov") produce estimates that mismatch the original data since they have to have a zero correlation/covariance between $x_1$ and $x_2$, whereas for the saturated mlav_m1 they reproduce the matrices from the original data. One thing I think I'm missing is why that would affect the standardized regression coefficients only, and not the unstandardized ones. – user1205901 - Reinstate Monica Jun 19 at 23:37

As I said in the comments, it's due to the covariance matrix implied by the model.

Running through your example:

(cov.1 <- lavInspect(mlav_m1, "cov.ov"))
#    y      x1     x2
# y  35.189
# x1  2.126  0.277
# x2 -4.957 -0.361  0.927


This is the standard situation. Once you constrain the covariance between x1 and x2 to 0, the variance of and covariances with y change:

(cov.2 <- lavInspect(mlav_m2, "cov.ov"))
#    y      x1     x2
# y  30.207
# x1  0.399  0.277
# x2 -4.436  0.000  0.927


You still get the same coefficients from both matrices:

# model 1:
solve(cov.1[2:3, 2:3], cov.1[1, 2:3])
#        x1        x2
#  1.442491 -4.782890
# model 2:
solve(cov.2[2:3, 2:3], cov.2[1, 2:3])
#        x1        x2
#  1.442491 -4.782890


I think what's happening is the SEM coefficients are still the ML coefficients for the regression problem. And the ML solution to the problem is unaffected by the covariance between x1 and x2. However, the resulting covariance matrix implied by the new ML solution is different.

When you check the correlation matrix implied by the model and solve for the coefficients:

(cor.1 <- cov2cor(cov.1))  # Same as: lavInspect(mlav_m1, "cor.ov")
#    y      x1     x2
# y   1.000
# x1  0.681  1.000
# x2 -0.868 -0.712  1.000

(cor.2 <- cov2cor(cov.2))  # Same as: lavInspect(mlav_m2, "cor.ov")
#    y      x1     x2
# y   1.000
# x1  0.138  1.000
# x2 -0.838  0.000  1.000


Different again, but these matrices give different coefficients:

# model 1
solve(cor.1[2:3, 2:3], cor.1[1, 2:3])
#         x1         x2
#  0.1279701 -0.7764883
# model 2
solve(cor.2[2:3, 2:3], cor.2[1, 2:3])
#         x1         x2
#  0.1381212 -0.8380823


I hope this answers why you get the different set of standardized coefficients. lavaan likely uses the standard formula of $$b_x \times \frac{sd(x)}{sd(y)}$$ but $$sd(y)$$ is the model-implied value not the observed value. Also, as you've seen, the $$R^2$$ also follows from the model-implied covariance matrix, and not the actual data. This is likely so because SEM models deal with latent variables all the time, so researchers typically use formulas for different metrics that follow from the model-implied covariance matrix.

• This makes lots of sense and you're right about standardized coefficients, though lavaan uses the model-implied value for $sd(x)$ as well as $sd(y)$. This can be shown as param_estimates_m2[1,4] * (sqrt(lavInspect(mlav_m2,"cov.ov")[2,2])/sqrt(lavInspect(mlav_m2,"cov.ov")[1,1])) equals print((stand_estimates_m2[2,4]),nd=7). Do you know whether $R^2$ can increase from fixing the covariance of predictors to $0$? That piece of the puzzle still eludes me. Intuitively it seems like it shouldn't increase in a multiple regression scenario, though maybe it could in some more complex path analysis? – user1205901 - Reinstate Monica Jun 20 at 2:13
• Yes, it's possible. I just tried with a 3-by-3 correlation matrix (2 predictors, one outcome) where all correlations are 0.3 but corr(x1, x2) = -0.3. Assuming this -0.3 to be 0 as you did above increased R^2 from .257 to .331. To guard against things like this, SEMers rely heavily on model fit indices, which in this example suggested model misspecification. – Heteroskedastic Jim Jun 20 at 3:25
• Also, I think you're correct about the model-implied SD of the predictors. But I think these should just be the same as the standard deviation of the predictors but without Bessel's correction since the estimation method is ML. This would change if one did: estimator="ULS", least squares. – Heteroskedastic Jim Jun 20 at 3:30