correlation structures: composites support (ML and GLS routes)

Composite models (the <~ operator) now work with correlation = TRUE,
on both the D-augmented ML route (the default) and the GLS route.
Root cause of the old failure (the GLS combination never converged;
the ML route stopped with a clean error): the correlation completion
ran BEFORE lav_lisrel_comp_set_intresvar(), so the completed residual
variances were computed from stale composite quantities, and the
implied diagonal drifted away from the unit target during
optimization (objective/gradient mismatch). The composite completion
now runs FIRST in lav_model_set_parameters() and lav_model_x2glist();
plain composite models are unaffected (nothing runs between the two
blocks when no correlation structure is present).

Two more fixes in the same setting:
- lav_pt_ndat(): the fixed composite-indicator VARIANCES were removed
  from the statistic count on top of the correlation correction (fixed
  ~*~ rows) or the free ~*~ scales that fit them -- the df now match
  an ordinary composite ML fit (which has the same chi-square).
- lav_lisrel_dimplied_dx(): in the D-augmented mode the composite
  psi* chain-rule block must not reintroduce diagonal-row entries for
  the non-delta columns (the unit-diagonal completion makes
  diag(Sigma*) constant); the analytic gradient again matches
  numDeriv at and off the solution.

missing = "ml" works for composite correlation models as well (same
likelihood as an ordinary composite FIML fit).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: yrosseel

R/lav_lavaan_step02_options.R

@@ -262,11 +262,6 @@ lav_step02_options <- function(slot_options = NULL,
       opt$meanstructure <- TRUE
     }
 
-    # composites marker: the D-augmented ML correlation mode has no
-    # composite support (yet); lav_options_set() keeps the (defaulted)
-    # estimator at GLS in that case
-    opt$.flat.composites <- any(flat_model$op == "<~")
-
     # conditional.x
     if ((is.list(ov_names_x) &&
       sum(sapply(ov_names_x, FUN = length)) == 0L) ||
@@ -324,16 +319,6 @@ lav_step02_options <- function(slot_options = NULL,
          \"delta\" parameterization."))
     }
 
-    # D-augmented ML mode (correlation + ML-family estimator): no
-    # composites support (the composite completion and the unit-diagonal
-    # completion both reparameterize the same matrices)
-    if (isTRUE(lavoptions$.correlation.ml) && any(flat_model$op == "<~")) {
-      lav_msg_stop(gettext(
-        "estimator = \"ML\" with correlation structures is not supported
-        (yet) for models with composites (the \"<~\" operator); use
-        estimator = \"GLS\"."))
-    }
-
     # correlation = TRUE with observed product (interaction) terms: the
     # product of standardized variables does not have unit variance, so the
     # product terms must keep a free variance (github issue #427). We turn

R/lav_model_utils.R

@@ -100,6 +100,25 @@ lav_model_set_parameters <- function(lavmodel = NULL, x = NULL) {
   # parameters): the completion then targets diag(Sigma*) = 1
   correlation_ml <- correlation && lav_model_delta_free(lavmodel)
 
+  # composites first: the derived composite quantities (variances,
+  # intercepts) feed the correlation completion below
+  if (lavmodel@composites) {
+    # for package stdmod only! (vignette stdmod_lavaan uses old fit object)
+    # if (.hasSlot(lavmodel, "composites") && lavmodel@composites) {
+    nmat <- lavmodel@nmat
+    if (lavmodel@representation == "LISREL") {
+      for (g in 1:lavmodel@nblocks) {
+        # which mm belong to group g?
+        mm_in_group <- 1:nmat[g] + cumsum(c(0L, nmat))[g]
+
+        tmp[mm_in_group] <-
+          lav_lisrel_comp_set_intresvar(mlist = tmp[mm_in_group])
+      }
+    } else {
+      cat("FIXME: deal with Composites if representation = RAM")
+    }
+  }
+
   # categorical? set categorical theta elements (if any)
   if (lavmodel@categorical || correlation) {
     nmat <- lavmodel@nmat
@@ -140,23 +159,6 @@ lav_model_set_parameters <- function(lavmodel = NULL, x = NULL) {
     }
   }
 
-  if (lavmodel@composites) {
-    # for package stdmod only! (vignette stdmod_lavaan uses old fit object)
-    # if (.hasSlot(lavmodel, "composites") && lavmodel@composites) {
-    nmat <- lavmodel@nmat
-    if (lavmodel@representation == "LISREL") {
-      for (g in 1:lavmodel@nblocks) {
-        # which mm belong to group g?
-        mm_in_group <- 1:nmat[g] + cumsum(c(0L, nmat))[g]
-
-        tmp[mm_in_group] <-
-          lav_lisrel_comp_set_intresvar(mlist = tmp[mm_in_group])
-      }
-    } else {
-      cat("FIXME: deal with Composites if representation = RAM")
-    }
-  }
-
   lavmodel@GLIST <- tmp
 
   lavmodel
@@ -219,6 +221,23 @@ lav_model_x2glist <- function(lavmodel = NULL, x = NULL,
   #        }
   #    }
 
+  # composites first: the derived composite quantities feed the
+  # correlation completion below
+  if (lavmodel@composites) {
+    nmat <- lavmodel@nmat
+    if (lavmodel@representation == "LISREL") {
+      for (g in 1:lavmodel@nblocks) {
+        # which mm belong to group g?
+        mm_in_group <- 1:nmat[g] + cumsum(c(0L, nmat))[g]
+
+        glist[mm_in_group] <-
+          lav_lisrel_comp_set_intresvar(mlist = glist[mm_in_group])
+      }
+    } else {
+      cat("FIXME: deal with Composites when representation = RAM")
+    }
+  }
+
   # in 0.6-13: we always set theta/delta
   if ((lavmodel@categorical || correlation) && set_delta) {
     nmat <- lavmodel@nmat
@@ -250,21 +269,6 @@ lav_model_x2glist <- function(lavmodel = NULL, x = NULL,
     }
   }
 
-  if (lavmodel@composites) {
-    nmat <- lavmodel@nmat
-    if (lavmodel@representation == "LISREL") {
-      for (g in 1:lavmodel@nblocks) {
-        # which mm belong to group g?
-        mm_in_group <- 1:nmat[g] + cumsum(c(0L, nmat))[g]
-
-        glist[mm_in_group] <-
-          lav_lisrel_comp_set_intresvar(mlist = glist[mm_in_group])
-      }
-    } else {
-      cat("FIXME: deal with Composites when representation = RAM")
-    }
-  }
-
   glist
 }
 

R/lav_options.R

@@ -1653,7 +1653,7 @@ lav_options_set <- function(opt = NULL) {
     # do not depend on the scaling of the input). Since 0.7-1 this is
     # also the DEFAULT route: when the estimator is left unspecified,
     # we only fall back to GLS (quietly) for settings the
-    # D-augmentation does not support (yet): composites.
+    # D-augmentation does not support.
     # missing = "ml" (FIML) is supported on the ML route (since 0.7-2):
     # the casewise/pattern likelihood simply consumes the D-augmented
     # implied moments Sigma = Delta P(theta) Delta and Mu.
@@ -1664,10 +1664,6 @@ lav_options_set <- function(opt = NULL) {
     correlation_ml_ok <-
       (!opt$conditional.x && any(opt$missing == c("listwise", "ml"))) ||
       (opt$conditional.x && opt$missing == "listwise")
-    if (estimator_default) {
-      correlation_ml_ok <- correlation_ml_ok &&
-        !isTRUE(opt$.flat.composites)
-    }
     if (lav_options_estimatorgroup(opt$estimator) == "ML" &&
         correlation_ml_ok) {
       opt$.correlation.ml <- TRUE

R/lav_partable_utils.R

@@ -184,6 +184,21 @@ lav_pt_ndat <- function(partable) {
                          partable$lhs %in% ov_cind &
                          partable$rhs %in% ov_cind &
                          partable$free == 0L)
+      # correlation structures: the indicator VARIANCE moments are
+      # already accounted for -- a FIXED ~*~ row subtracts them via the
+      # correlation correction below, a FREE ~*~ row (D-augmented ML)
+      # fits them -- so only the fixed indicator COVARIANCES remove a
+      # sample statistic here
+      if (correlation && !categorical && length(covar_idx) > 0L) {
+        scaled_ov <- partable$lhs[partable$op == "~*~" &
+                                  partable$block == b]
+        drop_idx <- which(partable$lhs[covar_idx] ==
+                          partable$rhs[covar_idx] &
+                          partable$lhs[covar_idx] %in% scaled_ov)
+        if (length(drop_idx) > 0L) {
+          covar_idx <- covar_idx[-drop_idx]
+        }
+      }
       ndat[b] <- ndat[b] - length(covar_idx)
     }
 

R/lav_representation_lisrel.R

@@ -4199,6 +4199,21 @@ lav_lisrel_dimplied_dx <- function(mlist           = NULL,
       sigma_chain <- sigma_chain * delta_weight
     }
 
+    # D-augmented ML mode (FREE ~*~ scales, full vech layout): the
+    # unit-diagonal completion makes diag(Sigma*) constant, so the
+    # composite chain must not reintroduce diagonal-row entries for the
+    # non-delta columns (the num_idx variables keep a free variance)
+    if (!categorical && correlation && n_del > 0L) {
+      diag_pos_c <- lav_mat_diagh_idx(nvar)
+      if (length(num_idx) > 0L) {
+        diag_pos_c <- diag_pos_c[-num_idx]
+      }
+      nondelta_cols <- setdiff(seq_len(ncol(sigma_chain)), x_delta_idx)
+      if (length(diag_pos_c) > 0L && length(nondelta_cols) > 0L) {
+        sigma_chain[diag_pos_c, nondelta_cols] <- 0
+      }
+    }
+
     # categorical / correlation layouts subset and reorder the sigma rows;
     # bring the (full-vech) chain term into the same layout
     if (!is.null(sigma_row_map)) {

man/lavOptions.Rd

@@ -111,12 +111,15 @@ Model features:
      scales are free \code{~*~} parameters, and the unit-variance
      constraint applies to the marginal variances; the result is
      identical (same conditional likelihood) to the joint
-     \code{fixed.x = TRUE} analysis. The ML route requires
-     single-level data and no composites (and \code{missing = "ml"}
-     requires \code{conditional.x = FALSE}); in all other settings, a
-     defaulted estimator quietly falls back to GLS. Estimator WLS (=
-     ADF weights in the correlation metric) covers the non-normal
-     case.
+     \code{fixed.x = TRUE} analysis. Composite models (the \code{<~}
+     operator) are supported as well (since 0.7-2, on both the ML and
+     the GLS route): the composite-indicator covariances are fixed in
+     the correlation metric, and the indicator variances are absorbed
+     by the \code{~*~} scales. The ML route requires single-level data
+     (and \code{missing = "ml"} requires \code{conditional.x =
+     FALSE}); in all other settings, a defaulted estimator quietly
+     falls back to GLS. Estimator WLS (= ADF weights in the
+     correlation metric) covers the non-normal case.
      Both
      \code{fixed.x = FALSE} and \code{fixed.x = TRUE} are supported (the
      latter affects the way the standard errors are computed), and (since