diff --git a/R/est_lms.R b/R/est_lms.R
index 5bec7fc0..624519e7 100644
--- a/R/est_lms.R
+++ b/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,7 +740,7 @@ fischerScoring <- function(par,
         iterations = i
       ))
     }
-        
+
     alpha   <- max.alpha
     success <- FALSE
 
@@ -740,7 +748,7 @@ fischerScoring <- function(par,
       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
     }
diff --git a/R/generics_modsem_da.R b/R/generics_modsem_da.R
index 937c93b6..df266d32 100644
--- a/R/generics_modsem_da.R
+++ b/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 ------------------------------------------------------------
diff --git a/R/inspect_da.R b/R/inspect_da.R
index 16162c1f..2f4c6463 100644
--- a/R/inspect_da.R
+++ b/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]]
     )
diff --git a/R/model_da.R b/R/model_da.R
index 0364ab3c..0748a945 100644
--- a/R/model_da.R
+++ b/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,
diff --git a/R/quadrature.R b/R/quadrature.R
index ff83fb6a..6789f00e 100644
--- a/R/quadrature.R
+++ b/R/quadrature.R
@@ -66,10 +66,11 @@ adaptiveGaussQuadrature <- function(fun,
                                     m.ceil = m + m / 2,
                                     k = 1,
                                     iter.max = 10,
-                                    node.max = 500,
+                                    node.max = max(m, 750),
                                     tol = 1e-12,
                                     mdiff.tol = 2,
                                     secondary.pruning = TRUE,
+                                    wtol = .Machine$longdouble.eps,
                                     ...) {
   if (k == 0 || m == 0)
     return(list(n = matrix(0), w = 1, f = NA, m = 1, k = 1))
@@ -84,7 +85,8 @@ adaptiveGaussQuadrature <- function(fun,
     out <- adaptiveGaussQuadratureK(
       fun = fun, a = a, b = b, m = m, m.ceil = m.ceil,
       k = 1, K = k, iter = 1, iter.max = iter.max,
-      node.max = node.max, tol = tol, mdiff.tol = mdiff.tol, ...
+      node.max = node.max, tol = tol, mdiff.tol = mdiff.tol,
+      wtol = wtol, ...
     )
 
     return(out)
@@ -103,7 +105,8 @@ adaptiveGaussQuadrature <- function(fun,
     QUAD <- adaptiveGaussQuadratureK(
       fun = fun, a = a[i], b = b[i], m = m, m.ceil = m.ceil[i],
       k = i, K = k, iter = 1, iter.max = iter.max,
-      node.max = node.max, tol = tol, mdiff.tol = mdiff.tol, ...
+      node.max = node.max, tol = tol, mdiff.tol = mdiff.tol,
+      wtol = wtol, ...
     )
 
     NODES[[i]]     <- QUAD$n[, i]
@@ -114,17 +117,11 @@ adaptiveGaussQuadrature <- function(fun,
 
   quadn <- as.matrix(do.call(expand.grid, NODES))
   quadw <- Reduce(kronecker, rev(WEIGHTS))
-  quadf <- fun(quadn, ...)
-  weighted <- sweep(quadf, 2, quadw, "*")
-  integral <- quadratureLogLik(weighted)
-  error <- 0
-  error.abs <- 0
-  error.rel <- 0
 
   if (secondary.pruning) {
-    pruned <- pruneQuadratureNodes(
-      quadw = quadw, quadn = quadn, quadf = quadf,
-      a = a, b = b, tol = tol
+    pruned <- prunedQuadrature(
+      fun = fun, quadw = quadw, quadn = quadn,
+      a = a, b = b, tol = tol, wtol = wtol, ...
     )
 
     quadw     <- pruned$quadw
@@ -134,18 +131,28 @@ adaptiveGaussQuadrature <- function(fun,
     error     <- pruned$I.err
     error.abs <- pruned$I.err.abs
     error.rel <- pruned$I.err.rel
+
+  } else {
+    quadf     <- fun(quadn, ...)
+    integral  <- quadratureLogLik(sweep(quadf, 2, quadw, "*"))
+    error     <- 0
+    error.abs <- 0
+    error.rel <- 0
   }
 
   list(
     n = quadn,
     w = quadw,
     F = quadf,
+    m.target = m,
+    k = k,
     m.ceil = new.ceils,
     capped = any(capped),
     error = error,
     error.abs = error.abs,
     error.rel = error.rel,
-    integral = integral
+    integral = integral,
+    adaptive = TRUE
   )
 }
 
@@ -159,9 +166,10 @@ adaptiveGaussQuadratureK <- function(fun,
                                      K = 1,
                                      iter = 1,
                                      iter.max = 10,
-                                     node.max = 500,
+                                     node.max = max(750, m),
                                      tol = 1e-12,
                                      mdiff.tol = 2,
+                                     wtol = .Machine$longdouble.eps,
                                      ...) {
   if (k == 0 || m == 0)
     return(list(n = matrix(0), w = 1, f = NA, m = 1, k = 1))
@@ -169,20 +177,18 @@ adaptiveGaussQuadratureK <- function(fun,
   if (is.null(m.ceil) || is.na(m.ceil) || m.ceil <= 0)
     m.ceil <- round(estMForNodesInRange(m, a = -5, b = 5))
 
-  m.ceil <- min(m.ceil, node.max)
+  m.ceil <- max(m, min(m.ceil, node.max))
   quad   <- quadrature(m = m.ceil, k = 1)
+  quadw  <- quad$w
 
   if (K > 1) { # more computationally efficient/scalable
-    quadn <- matrix(0, nrow = m.ceil, ncol = K)
+    quadn <- matrix(0, nrow = length(quadw), ncol = K)
     quadn[,k] <- quad$n
   } else quadn <- quad$n
 
-  quadf <- fun(quadn, ...)
-  quadw <- quad$w
-
-  pruned <- pruneQuadratureNodes(
-    quadw = quadw, quadn = quadn, quadf = quadf,
-    a = a, b = b, tol = tol
+  pruned <- prunedQuadrature(
+    fun = fun, quadw = quadw, quadn = quadn,
+    a = a, b = b, tol = tol, wtol = wtol, ...
   )
 
   quadw     <- pruned$quadw
@@ -200,18 +206,22 @@ adaptiveGaussQuadratureK <- function(fun,
   diff.m <- NROW(quadn) - m
 
   # Calculate next ceiling
-  if (abs(diff.m) > 2 * mdiff.tol) {
+  if (abs(diff.m) > 4 * mdiff.tol) {
     # If the difference is large, use a rough estimate
     new.ceil <- round(estMForNodesInRange(k = m, a = lower, b = upper))
 
     # else do a step change
-  } else new.ceil <- m.ceil - diff.m
+  } else {
+    pct.diff <- abs(diff.m) / m # relative error
+    rel.diff <- sign(diff.m) * max(round(pct.diff * m.ceil), 1L)
+    new.ceil <- m.ceil - rel.diff
+  }
 
   # Never request more than `node.max` nodes. If the target can't be reached
   # within `node.max`, use `node.max` instead of keeping the current ceiling
   new.ceil   <- min(new.ceil, node.max)
   OKNextCeil <- new.ceil != m.ceil
-  converged <- abs(diff.m) <= mdiff.tol
+  converged <- abs(diff.m) <= 0 || (abs(diff.m) <= mdiff.tol && iter > 1L)
   OKIter <- iter < iter.max
 
   if (!converged && OKIter && OKNextCeil) {
@@ -229,6 +239,7 @@ adaptiveGaussQuadratureK <- function(fun,
       iter.max = iter.max,
       node.max = node.max,
       mdiff.tol = mdiff.tol,
+      wtol = wtol,
       ...
     ))
   }
@@ -248,13 +259,15 @@ adaptiveGaussQuadratureK <- function(fun,
     k = k,
     m = nrow(quadn) ^ (1 / k),
 
-    m.ceil   = m.ceil,
-    capped   = m.ceil >= node.max && NROW(quadn) < m - mdiff.tol,
-    iter     = iter,
-    error    = I.err,
+    m.target  = m,
+    m.ceil    = m.ceil,
+    capped    = m.ceil >= node.max && NROW(quadn) < m - mdiff.tol,
+    iter      = iter,
+    error     = I.err,
     error.abs = I.err.abs,
     error.rel = I.err.rel,
-    integral = I.cur
+    integral  = I.cur,
+    adaptive  = TRUE
   )
 }
 
@@ -332,6 +345,53 @@ estMForNodesInRange <- function(k, a, b,
 }
 
 
+prunedQuadrature <- function(fun, quadw, quadn, a, b, tol, wtol = NULL,
+                             edge.ratio = 1e4, edge.ratio.eval = 100, ...) {
+  keep <- seq_along(quadw)
+  drop <- integer(0L)
+
+  # pre pruning based on weights only
+  if (!is.null(wtol) && is.finite(wtol) && wtol > 0) {
+    idx  <- order(quadw)
+    drop <- idx[cumsum(quadw[idx]) < wtol]
+
+    if (length(drop) < length(quadw))
+      keep <- setdiff(keep, drop) # keeps the original node order
+    else
+      drop <- integer(0L) # never remove all nodes...
+  }
+
+  quadf  <- fun(quadn[keep, , drop = FALSE], ...)
+  pruned <- pruneQuadratureNodes(
+    quadw = quadw[keep], quadn = quadn[keep, , drop = FALSE], quadf = quadf,
+    a = a, b = b, tol = tol
+  )
+
+  if (!length(drop)) # no nodes were (pre) pruned
+    return(pruned)
+
+  # Check for non-trimmed tails. If we don't trim a tail, we might have an
+  # outlier, and we should evaluate the pre-pruned weights/nodes
+  w.edge <- max(
+    edge.ratio * max(quadw[drop]), # largest removed weight
+    edge.ratio.eval * min(quadw[keep]) # smallest weight evaluated
+  )
+
+  if (all(pruned$quadw > w.edge))
+    return(pruned) # no (kept) weights are close to the cut off
+
+  # If the tails weren't trimmed, we just evaluate everything
+  quadf.full <- matrix(NA_real_, nrow = NROW(quadf), ncol = length(quadw))
+  quadf.full[, keep] <- quadf
+  quadf.full[, drop] <- fun(quadn[drop, , drop = FALSE], ...)
+
+  pruneQuadratureNodes(
+    quadw = quadw, quadn = quadn, quadf = quadf.full,
+    a = a, b = b, tol = tol
+  )
+}
+
+
 quadratureLogLik <- function(weighted) {
   rs <- rowSums(weighted)
 
diff --git a/R/utils_da.R b/R/utils_da.R
index f8cee84e..d9a0b47d 100644
--- a/R/utils_da.R
+++ b/R/utils_da.R
@@ -431,11 +431,6 @@ nFreeInterceptsDA <- function(model) {
 }
 
 
-getInfoQuad <- function(quad) {
-  list(dim = quad$k, nodes.dim = quad$m, nodes.total = quad$m ^ quad$k)
-}
-
-
 getFixedInterceptSyntax <- function(indicator, syntax, parTable) {
   if (is.null(indicator) || is.null(syntax) ||
       NROW(parTable[parTable$lhs == indicator &
@@ -1469,3 +1464,45 @@ removeTempOV_RowsParTable <- function(parTable) {
   tmp <- startsWith(parTable$lhs, TEMP_OV_PREFIX) | startsWith(parTable$rhs, TEMP_OV_PREFIX)
   parTable[!tmp, , drop = FALSE]
 }
+
+
+getQuadInfo <- function(quad, method) {
+  if (is.null(quad) || is.null(method) || tolower(method) != "lms")
+    return(list())
+
+  # some basic info can be retrieved from the first quadrature spec alone
+  adaptive <- quad[[1L]]$adaptive && !is.null(quad[[1L]]$m.target)
+
+  if (adaptive) {
+    nodes.per.dim <- quad[[1L]]$m.target
+    dimensions <- quad[[1L]]$k
+    nodes.total <- nodes.per.dim ^ dimensions
+
+    # effective number of nodes
+    M <- do.call(rbind, lapply(quad, FUN = \(Q) Q$m.ceil))
+
+    eff.nodes.per.dim <- apply(
+      X = M,
+      MARGIN = 2L,
+      FUN = mean
+    )
+
+    eff.nodes.total <- prod(eff.nodes.per.dim)
+
+  } else {
+    nodes.per.dim <- quad[[1L]]$m
+    dimensions <- quad[[1L]]$k
+    nodes.total <- nodes.per.dim ^ dimensions
+    eff.nodes.per.dim <- nodes.per.dim
+    eff.nodes.total <- nodes.total
+  }
+
+  list(
+    adaptive          = adaptive,
+    nodes.per.dim     = nodes.per.dim,
+    dimensions        = dimensions,
+    nodes.total       = nodes.total,
+    eff.nodes.per.dim = eff.nodes.per.dim,
+    eff.nodes.total   = eff.nodes.total
+  )
+}
diff --git a/tests/testthat/test_quadrature_pruning.R b/tests/testthat/test_quadrature_pruning.R
new file mode 100644
index 00000000..41ef9bf2
--- /dev/null
+++ b/tests/testthat/test_quadrature_pruning.R
@@ -0,0 +1,76 @@
+devtools::load_all()
+
+# Pruning nodes by their weights before evaluating the integrand should not
+# change the integral, also when an outlier has its posterior mass beyond the
+# weight cutoff (see `prunedQuadrature()`)
+
+set.seed(2938)
+s <- 0.2
+y <- rnorm(200, mean = 0, sd = sqrt(1 + s^2))
+
+counter <- new.env()
+densityFun <- function(y) {
+  function(z, ...) {
+    counter$n <- counter$n + NROW(z)
+    outer(y, z[, 1], FUN = \(yy, zz) dnorm(yy, mean = zz, sd = s))
+  }
+}
+
+integral1D <- function(fun, m, wtol) {
+  q <- quadrature(m = m, k = 1)
+  counter$n <- 0
+  pruned <- prunedQuadrature(fun = fun, quadw = q$w, quadn = q$n,
+                             a = -Inf, b = Inf, tol = 1e-12, wtol = wtol)
+  list(I = pruned$I.cur, evaluated = counter$n, total = length(q$w))
+}
+
+wtol <- .Machine$longdouble.eps
+
+for (m in c(32, 100, 300)) {
+  # No outlier: tails are trimmed well inside the cutoff, so the dropped
+  # nodes should never be evaluated
+  fun  <- densityFun(y)
+  full <- integral1D(fun, m = m, wtol = NULL)
+  wcut <- integral1D(fun, m = m, wtol = wtol)
+
+  testthat::expect_equal(wcut$I, full$I, tolerance = 1e-10)
+  testthat::expect_lt(wcut$evaluated, wcut$total)
+
+  # Outlier beyond the weight cutoff (posterior mode at ~11.5, cutoff at
+  # ~9-10). GH(32) is a coarse rule, where neighbouring weights in the tails
+  # differ by more than `edge.ratio`
+  fun  <- densityFun(c(y, 12))
+  full <- integral1D(fun, m = m, wtol = NULL)
+  wcut <- integral1D(fun, m = m, wtol = wtol)
+
+  testthat::expect_equal(wcut$I, full$I, tolerance = 1e-10)
+  testthat::expect_equal(wcut$evaluated, wcut$total)
+}
+
+# Adaptive quadrature, single dimension
+for (y.out in c(10, 12)) {
+  fun  <- densityFun(c(y, y.out))
+  full <- suppressWarnings(adaptiveGaussQuadrature(fun, a = -Inf, b = Inf, m = 32, k = 1,
+                                                   m.ceil = NULL, wtol = NULL))
+  wcut <- suppressWarnings(adaptiveGaussQuadrature(fun, a = -Inf, b = Inf, m = 32, k = 1,
+                                                   m.ceil = NULL, wtol = wtol))
+
+  testthat::expect_equal(wcut$integral, full$integral, tolerance = 1e-10)
+  testthat::expect_equal(wcut$m.ceil, full$m.ceil)
+}
+
+# Adaptive quadrature, product grid (two dimensions)
+s2 <- 0.3
+Y  <- rbind(matrix(rnorm(400, mean = 0, sd = sqrt(1 + s2^2)), ncol = 2), c(11, 0.5))
+fun2D <- function(z, ...) {
+  vapply(seq_len(NROW(z)), FUN.VALUE = numeric(NROW(Y)), FUN = \(j)
+         dnorm(Y[, 1], mean = z[j, 1], sd = s2) * dnorm(Y[, 2], mean = z[j, 2], sd = s2))
+}
+
+full <- suppressWarnings(adaptiveGaussQuadrature(fun2D, a = -Inf, b = Inf, m = 16, k = 2,
+                                                 m.ceil = NULL, wtol = NULL))
+wcut <- suppressWarnings(adaptiveGaussQuadrature(fun2D, a = -Inf, b = Inf, m = 16, k = 2,
+                                                 m.ceil = NULL, wtol = wtol))
+
+testthat::expect_equal(wcut$integral, full$integral, tolerance = 1e-10)
+testthat::expect_equal(NROW(wcut$n), NROW(full$n))
diff --git a/vignettes/estimation_lms.Rmd b/vignettes/estimation_lms.Rmd
index c3a35b47..33c1141d 100644
--- a/vignettes/estimation_lms.Rmd
+++ b/vignettes/estimation_lms.Rmd
@@ -79,3 +79,9 @@ fit <- modsem(tpb_uk,
               max.step = 1) # Maximum number of steps in the maximization step
 summary(fit)
 ```
+
+When using the quasi-adaptive quadrature, the `Numerical Integration` section of the summary
+reports two sets of node counts. `Standard` is the requested number of nodes (`nodes`), which
+is roughly the number of nodes kept after pruning. `Effective` is the number of nodes you would
+need with a fixed Gauss-Hermite quadrature (`adaptive.quad = FALSE`) to get effectively the same
+results. The same information can be retrieved using `modsem_inspect(fit, "quad")`.