Small improvement to adaptive quadrature

R/est_lms.R

@@ -307,6 +307,12 @@ emLms <- function(model,
         testSimpleGradient <- FALSE
       }
 
+      # Did the recalculated quadrature end up on a different grid?
+      quadChanged <- recalcQuad && !is.null(lastQuad) && !identical(
+        lapply(lastQuad, FUN = \(q) q$m.ceil),
+        lapply(P$quad,   FUN = \(q) q$m.ceil)
+      )
+
       lastQuad <- P$quad
 
       logLikNew  <- P$obsLL
@@ -382,9 +388,10 @@ emLms <- function(model,
         }
       }
 
-      # If quadrature has just been recalculated, we should be less strict
-      quad.err <- as.integer(recalcQuad) * P$quad.err
-      if (deltaLL < -max(1e-8, quad.err, na.rm = TRUE)) {
+      # If the quadrature grid itself changed `logLikOld` and `logLikNew`
+      # are computed on different grids, and a decrease says nothing about the
+      # EM step, so we skip the check
+      if (!quadChanged && deltaLL < -1e-8) {
         if (verbose) cat("\n")
         mod_msg_warn_immediate(
           sprintf("Loglikelihood decreased by %.2g!", deltaLL),
@@ -408,11 +415,11 @@ emLms <- function(model,
       } else if (algorithm == "EMA" && n.em.iter > ema.min.iter) {
         # Check if expected number of iterations is increasing
         change <- diff(utils::tail(history$expected.remaining.iterations, ema.min.iter))
-       
+
         if (all(is.finite(change)) && mean(change) > ema.trend.min) {
           remaining <- utils::tail(history$expected.remaining.iterations, 1L)[[1L]]
 
-          # Increase the max step in case remaining > 500 or 
+          # Increase the max step in case remaining > 500 or
           # QN/FS fails, and we fall back to an EM step
           tmp.max.step <- max(fallback.max.step.min,
                             fallback.max.step.mult * max.step,
@@ -440,7 +447,7 @@ emLms <- function(model,
 
       .obsgr <- function(x) {
         grad <- gradientObsLogLikLms(theta = x, model = model, P = P, sign = -1, epsilon = epsilon)
-   
+
         non.finite <- !is.finite(grad)
         if (any(non.finite)) {
           # This might happen during a bad line search in optim()
@@ -619,6 +626,7 @@ emLms <- function(model,
       includeStartModel = TRUE,
       startModel        = model
     )
+
     fit$iteration.history <- history
     fit$convergence.grad  <- list(max.abs.grad = grad.max, tol = convergence.grad,
                                   steps = grad.steps, status = grad.status)
@@ -700,12 +708,12 @@ fischerScoring <- function(par,
 
   truncate <- function(x) {
     underflow <- x < lower
-    overflow  <- x > upper 
+    overflow  <- x > upper
 
     x[underflow] <- lower[underflow]
     x[overflow]  <- upper[overflow]
 
-    x 
+    x
   }
 
   f0 <- fn(par, ...)
@@ -732,15 +740,15 @@ fischerScoring <- function(par,
         iterations = i
       ))
     }
-        
+
     alpha   <- max.alpha
     success <- FALSE
 
     while (alpha > min.alpha) {
       parTrial  <- truncate(par + alpha * direction)
       f1  <- suppressWarnings(fn(parTrial, ...))
 
-      ok <- is.finite(f1) && f1 <= f0 
+      ok <- is.finite(f1) && f1 <= f0
       if (ok) { success <- TRUE; break }
       alpha <- alpha / 2
     }

R/generics_modsem_da.R

@@ -451,22 +451,49 @@ print.summary_da <- function(x, digits = 3, ...) {
                    width.out = width.out), "\n")
 
   # Intergration ---------------------------------------------------------------
-  if (!is.null(x$quad)) {
-    nodesDim   <- x$quad$nodes.dim
-    dimVal     <- x$quad$dim
-    nodesTotal <- x$quad$nodes.total
-
-    if (length(nodesDim) || length(dimVal) || length(nodesTotal)) {
-      cat("Numerical Integration:\n")
-      names <- c("Points of integration (per dim)", "Dimensions",
-                 "Total points of integration")
-      values <- c(formatCount(nodesDim),
-                  formatCount(dimVal),
-                  formatCount(nodesTotal))
-      cat(allignLhsRhs(lhs = names, rhs = values, pad = "  ",
-                       width.out = width.out), "\n", sep = "")
+  if (length(x$quad)) {
+    quad <- x$quad
+    adaptive <- quad$adaptive
+
+    cat("Numerical Integration:\n")
+    names <- c(
+      "Points of integration (per dim)",
+      "Dimensions",
+      "Total points of integration"
+    )
+
+    values <- c(
+      formatCount(quad$nodes.per.dim),
+      formatCount(quad$dimensions),
+      formatCount(quad$nodes.total)
+    )
+
+    eff.values <- NULL
+    if (adaptive) {
+      names <- c(
+        "",
+        names
+      )
+
+      values <- c(
+        "Standard",
+        values
+      )
+
+      eff.values <- c(
+        "Effective",
+        paste0(formatCount(quad$eff.nodes.per.dim), collapse = ","),
+        formatCount(quad$dimensions),
+        formatCount(quad$eff.nodes.total)
+      )
     }
 
+    cat(allignLhsRhs(
+      lhs = names, rhs = values, pad = "  ",
+      rhs.scaled = eff.values,
+      width.out = width.out),
+      "\n", sep = ""
+    )
   }
 
   # Comparative fit ------------------------------------------------------------

R/inspect_da.R

@@ -3,7 +3,7 @@ inspectDA_Matrices <- c("lambda", "theta", "wmat", "tmat",
 
 
 inspectDA_Optim <- c("coefficients.free", "vcov.free", "information",
-                     "loglik", "iterations", "convergence")
+                     "loglik", "iterations", "convergence", "quadrature")
 
 
 modsem_inspect_da <- function(model,
@@ -206,6 +206,7 @@ modsem_inspect_da <- function(model,
      iterations        = model$iterations,
      convergence       = model$convergence,
      ovs               = ovs,
+     quadrature        = model$info.quad,
 
      ngroups     = model$model$info$n.groups,
      group       = model$args$group,
@@ -269,6 +270,7 @@ modsem_inspect_da <- function(model,
       all       = info,
       matrices  = info[inspectDA_Matrices],
       optim     = info[optim.fields],
+      quad      = info[["quadrature"]], # quad -> info[[what]]
       fit       = FIT(),
       info[[what]]
     )

R/model_da.R

@@ -966,7 +966,7 @@ finalizeModelEstimatesDA <- function(model,
     convergence.msg  = getConvergenceMessage(converged, iterations),
     type.se          = typeSE,
     type.estimates   = "unstandardized",
-    info.quad        = if (method == "lms") getInfoQuad(model$models[[1L]]$quad) else NULL,
+    info.quad        = getQuadInfo(quad = P$quad, method = method),
     FIM              = FIMo$FIM,
     vcov.all         = FIMo$vcov.all,
     vcov.free        = FIMo$vcov.free,