Double check separate vs. joint model in lavaan

lavaan_separate_joint.R

library(lavaan)

# sim_dat function from my blog post
# https://erikjanvankesteren.nl/blog/multivariate_regression
sim_dat <- function(N = 100, P = 5, Q = 2, varprop = 0.5, rho_X = 0, rho_B = 0, rho_E = 0, exact = TRUE) {
  # Create X matrix
  X <- matrix(rnorm(N*P), N)
  RX <- matrix(rho_X, P, P)
  diag(RX) <- 1
  CRX <- chol(RX)
  X <- if (exact) X %*% solve(chol(var(X)), CRX) else X %*% CRX
  
  # Create B matrix
  B <- matrix(rnorm(P*Q), P)
  RB <- matrix(rho_B, Q, Q)
  diag(RB) <- 1
  CRB <- chol(RB)
  B <- if (exact) B %*% solve(chol(var(B)), CRB) else B %*% CRB
  
  # Create E matrix
  E <- matrix(rnorm(N*Q), N)
  XB <- X %*% B
  vXB <- diag(var(XB))
  sE <- sqrt(vXB*varprop/(1-varprop))
  RE <- matrix(rho_E, Q, Q)
  diag(RE) <- 1
  CVE <- chol(RE) %*% diag(sE)
  E <- if (exact) E %*% solve(chol(var(E)), CVE) else E %*% CVE
  
  # Create Y matrix
  Y <- XB + E
  
  return(list(Y = Y, X = X, B = B, E = E))
}

# generate data with residual correlations
dat <- sim_dat(rho_E = 0.4)

# turn it into a data frame for lavaan
df <- data.frame(dat$Y) |> setNames(c("Y1", "Y2")) |> cbind(data.frame(dat$X))


# model without including residual correlation
model_separate <- "
  Y1 ~ X1 + X2 + X3 + X4 + X5
  Y2 ~ X1 + X2 + X3 + X4 + X5
  Y1 ~~ Y1
  Y2 ~~ Y2
  Y1 ~~ 0*Y2 # do not allow residual covariance!
"
fit_separate <- sem(model_separate, df)

# model that includes residual correlation
model_joint <- "
  Y1 ~ X1 + X2 + X3 + X4 + X5
  Y2 ~ X1 + X2 + X3 + X4 + X5
  Y1 ~~ Y1
  Y2 ~~ Y2
  Y1 ~~ Y2 # allow residual covariance!
"
fit_joint <- sem(model_joint, df)

# regression parameter estimates are the same
partable(fit_separate)[1:13,]
partable(fit_joint)[1:13,]

# also the same as lm
fit_lm <- lm(Y ~ X, dat)
summary(fit_lm)

# double check that we indeed have residual covariance
cor(residuals(fit_lm))