Update sem.R

It seems there is a scaling issue wrt to the parameters of the second factor. If the first loading is scaled to 1, using the same scaling on the rest of the parameters yield the "correct" result. Maybe there's an extra parameter not accounted for somewhere.
haziqj · Jun 21, 2024 · 1ad4f71 · 1ad4f71
1 parent a8374fe
commit 1ad4f71
Showing 1 changed file with 47 additions and 27 deletions.
diff --git a/inst/sem.R b/inst/sem.R
@@ -9,8 +9,8 @@ sd_y <- 0.01
 
 eta1 <- rnorm(n, sd = 1)
 eta2 <- 0.3 * eta1 + rnorm(n, sd = 1)
-Lambda <- matrix(c(1, 1.2, 1.5, 0, 0, 0,
-                   0, 0, 0, 1, 1.2, 1.5), ncol = 2)
+Lambda <- matrix(c(1, 1.2, 1.5, 0,   0,   0,
+                   0,   0,   0, 1, 1.2, 1.5), ncol = 2)
 dat <- t(Lambda %*% rbind(eta1, eta2)) + rmvnorm(n, sigma = diag(rep(sd_y, p)))
 
 # Prepare data for INLA --------------------------------------------------------
@@ -37,34 +37,34 @@ IDX <- bind_rows(
 
 # This version has fixed precision for latent variables and fixed loadings for
 # first item (lambda1 := 1 and lambda4 := 1).
-form <- Y ~ -1 +
-  f(idy1, hyper = list(prec = list(initial = 0, fixed = TRUE))) +
-  f(idy2, copy = "idy1", hyper = list(beta = list(fixed = FALSE))) +
-  f(idy3, copy = "idy1", hyper = list(beta = list(fixed = FALSE))) +
-
-  f(idy4, hyper = list(prec = list(initial = 0, fixed = TRUE))) +
-  f(idy5, copy = "idy4", hyper = list(beta = list(fixed = FALSE))) +
-  f(idy6, copy = "idy4", hyper = list(beta = list(fixed = FALSE))) +
-
-  f(idzero1, neg_one, copy = "idy1", hyper = list(beta = list(initial = 0.3, fixed = FALSE))) +
-  f(idzero2, copy = "idy4")
-
-# This version has "anchor" latent variables with fixed precision (:= 1)
 # form <- Y ~ -1 +
-#   f(idxNA1, model = "iid", hyper = list(prec = list(initial = 0, fixed = TRUE)),
-#     values = 1:n, constr = TRUE) +
-#   f(idy1, copy = "idxNA1", hyper = list(beta = list(fixed = FALSE))) +
-#   f(idy2, copy = "idxNA1", hyper = list(beta = list(fixed = FALSE))) +
-#   f(idy3, copy = "idxNA1", hyper = list(beta = list(fixed = FALSE))) +
+#   f(idy1, hyper = list(prec = list(initial = 0, fixed = TRUE))) +
+#   f(idy2, copy = "idy1", hyper = list(beta = list(fixed = FALSE))) +
+#   f(idy3, copy = "idy1", hyper = list(beta = list(fixed = FALSE))) +
 #
-#   f(idxNA2, model = "iid", hyper = list(prec = list(initial = 0, fixed = TRUE)),
-#     values = 1:n, constr = TRUE) +
-#   f(idy4, copy = "idxNA2", hyper = list(beta = list(fixed = FALSE))) +
-#   f(idy5, copy = "idxNA2", hyper = list(beta = list(fixed = FALSE))) +
-#   f(idy6, copy = "idxNA2", hyper = list(beta = list(fixed = FALSE))) +
+#   f(idy4, hyper = list(prec = list(initial = 0, fixed = TRUE))) +
+#   f(idy5, copy = "idy4", hyper = list(beta = list(fixed = FALSE))) +
+#   f(idy6, copy = "idy4", hyper = list(beta = list(fixed = FALSE))) +
 #
-#   f(idzero1, neg_one, copy = "idxNA1", hyper = list(beta = list(initial = 0.3, fixed = FALSE))) +
-#   f(idzero2, copy = "idxNA2")
+#   f(idzero1, neg_one, copy = "idy1", hyper = list(beta = list(initial = 0.3, fixed = FALSE))) +
+#   f(idzero2, copy = "idy4")
+
+# This version has "anchor" latent variables with fixed precision (:= 1)
+form <- Y ~ -1 +
+  f(idxNA1, model = "iid", hyper = list(prec = list(initial = 0, fixed = TRUE)),
+    values = 1:n, constr = TRUE) +
+  f(idy1, copy = "idxNA1", hyper = list(beta = list(fixed = FALSE))) +
+  f(idy2, copy = "idxNA1", hyper = list(beta = list(fixed = FALSE))) +
+  f(idy3, copy = "idxNA1", hyper = list(beta = list(fixed = FALSE))) +
+
+  f(idxNA2, model = "iid", hyper = list(prec = list(initial = 0, fixed = TRUE)),
+    values = 1:n, constr = TRUE) +
+  f(idy4, copy = "idxNA2", hyper = list(beta = list(fixed = FALSE))) +
+  f(idy5, copy = "idxNA2", hyper = list(beta = list(fixed = FALSE))) +
+  f(idy6, copy = "idxNA2", hyper = list(beta = list(fixed = FALSE))) +
+
+  f(idzero1, neg_one, copy = "idxNA1", hyper = list(beta = list(initial = 0.3, fixed = FALSE))) +
+  f(idzero2, copy = "idxNA2")
 
 # INLA call --------------------------------------------------------------------
 fit_inla <- inla(
@@ -75,8 +75,28 @@ fit_inla <- inla(
     rep(list(list()), p),
     list(list(control.sem = list(B = matrix(c("", "idzero1", "", ""), 2, 2), idx = 2)))
   ),
+  control.inla = list(int.strategy = "eb"),
   safe = FALSE,
   verbose = TRUE
 )
 summary(fit_inla)
 
+plot(fit_inla$summary.random$idxNA2$mean, eta2)
+abline(0, 1)
+
+plot(fit_inla$summary.random$idxNA1$mean, eta1)
+abline(0, 1)
+
+library(lavaan)
+dat_lav <- as.data.frame(dat)
+names(dat_lav) <- paste0("y", 1:p)
+fit_lav <- sem(
+  mod = "
+  eta1 =~ y1 + y2 + y3
+  eta2 =~ y4 + y5 + y6
+  eta2 ~ eta1
+  ",
+  data = dat_lav,
+  std.lv = TRUE
+)
+