[SAM] handle exogenous covariates in the presence of interaction terms

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: lavaan
 Title: Latent Variable Analysis
-Version: 0.6-19.2169
+Version: 0.6-19.2170
 Authors@R: c(person(given   = "Yves", family = "Rosseel",
                     role    = c("aut", "cre"),
                     email   = "Yves.Rosseel@UGent.be",

R/ctr_pml_utils.R

@@ -98,6 +98,8 @@ uni_lik <- function(Y1, th.y1, eXo = NULL, PI.y1 = NULL) {
   uni_lik <- pnorm(th.y1.upper) - pnorm(th.y1.lower)
 
   uni_lik[is.na(uni_lik)] <- 0
+
+  uni_lik
 }
 
 #################################################################

R/lav_model_compute.R

@@ -735,7 +735,6 @@ computeNU <- function(lavmodel = NULL, GLIST = NULL,
   NU
 }
 
-
 # E(Y): expectation (mean) of observed variables
 # returns vector 1 x nvar
 computeEY <- function(lavmodel = NULL, GLIST = NULL, lavsamplestats = NULL,
@@ -778,6 +777,48 @@ computeEY <- function(lavmodel = NULL, GLIST = NULL, lavsamplestats = NULL,
   EY
 }
 
+# E(Y|x_i): conditional expectation (mean) of observed variables
+# returns matrix N x nvar
+computeEYx <- function(lavmodel = NULL, GLIST = NULL, lavsamplestats = NULL,
+                       eXo = NULL, delta = TRUE) {
+  # state or final?
+  if (is.null(GLIST)) GLIST <- lavmodel@GLIST
+
+  nblocks <- lavmodel@nblocks
+  nmat <- lavmodel@nmat
+  representation <- lavmodel@representation
+
+  # return a list
+  EYx <- vector("list", length = nblocks)
+
+  # compute E(Y) for each group
+  for (g in 1:nblocks) {
+    # which mm belong to group g?
+    mm.in.group <- 1:nmat[g] + cumsum(c(0, nmat))[g]
+    MLIST <- GLIST[mm.in.group]
+
+    if (representation == "LISREL") {
+      EYx.g <- computeEYx.LISREL(
+        MLIST = MLIST,
+		eXo = eXo[[g]],
+        sample.mean = lavsamplestats@mean[[g]],
+        ov.y.dummy.ov.idx = lavmodel@ov.y.dummy.ov.idx[[g]],
+        ov.x.dummy.ov.idx = lavmodel@ov.x.dummy.ov.idx[[g]],
+        ov.y.dummy.lv.idx = lavmodel@ov.y.dummy.lv.idx[[g]],
+        ov.x.dummy.lv.idx = lavmodel@ov.x.dummy.lv.idx[[g]],
+        delta = delta
+      )
+    } else {
+      lav_msg_stop(gettext(
+        "only representation LISREL has been implemented for now"))
+    }
+
+    EYx[[g]] <- EYx.g
+  }
+
+  EYx
+}
+
 
 # E(Y | ETA, x_i): conditional expectation (means) of observed variables
 # for a given value of x_i AND eta_i

R/lav_objective.R

@@ -552,7 +552,7 @@ estimator.PML <- function(Sigma.hat = NULL, # model-based var/cov/cor
         pstar.idx <- PSTAR[i, j]
         # cat("pstar.idx =", pstar.idx, "i = ", i, " j = ", j, "\n")
         if (ov.types[i] == "numeric" &&
-          ov.types[j] == "numeric") {
+            ov.types[j] == "numeric") {
           # ordinary pearson correlation
           LIK[, pstar.idx] <-
             lav_bvreg_lik(
@@ -565,7 +565,7 @@ estimator.PML <- function(Sigma.hat = NULL, # model-based var/cov/cor
               rho = Cor.hat[i, j]
             )
         } else if (ov.types[i] == "numeric" &&
-          ov.types[j] == "ordered") {
+                   ov.types[j] == "ordered") {
           # polyserial correlation
           ### FIXME: th.y2 should go into ps_lik!!!
           LIK[, pstar.idx] <-
@@ -579,7 +579,7 @@ estimator.PML <- function(Sigma.hat = NULL, # model-based var/cov/cor
               rho = Cor.hat[i, j]
             )
         } else if (ov.types[j] == "numeric" &&
-          ov.types[i] == "ordered") {
+                   ov.types[i] == "ordered") {
           # polyserial correlation
           ### FIXME: th.y1 should go into ps_lik!!!
           LIK[, pstar.idx] <-
@@ -593,7 +593,7 @@ estimator.PML <- function(Sigma.hat = NULL, # model-based var/cov/cor
               rho = Cor.hat[i, j]
             )
         } else if (ov.types[i] == "ordered" &&
-          ov.types[j] == "ordered") {
+                   ov.types[j] == "ordered") {
           LIK[, pstar.idx] <-
             pc_lik_PL_with_cov(
               Y1 = X[, i],

R/lav_sam_step1.R

@@ -482,7 +482,11 @@ lav_sam_step1_local <- function(STEP1 = NULL, FIT = NULL,
       FS[[b]] <- do.call("cbind", tmp)
       colnames(FS[[b]]) <- LABEL
 
-      # dummy lv's?
+      # dummy lv's? (both 'x' and 'y'!)
+	  dummy.ov.idx <- c(FIT@Model@ov.y.dummy.ov.idx[[b]],
+	                    FIT@Model@ov.x.dummy.ov.idx[[b]])
+      dummy.lv.idx <- c(FIT@Model@ov.y.dummy.lv.idx[[b]],
+	                    FIT@Model@ov.x.dummy.lv.idx[[b]])
       if (length(dummy.lv.idx) > 0L) {
         FS.obs <- FIT@Data@X[[b]][, dummy.ov.idx, drop = FALSE]
         colnames(FS.obs) <- FIT@Data@ov.names[[b]][dummy.ov.idx]