From ffcd4c33ddd12a0da408e04e12128c68c33afb03 Mon Sep 17 00:00:00 2001 From: kss2k Date: Sun, 25 Oct 2026 17:11:34 +0100 Subject: [PATCH 1/7] Add more detailed quadrature summary --- R/est_lms.R | 15 ++++++----- R/generics_modsem_da.R | 55 +++++++++++++++++++++++++++++---------- R/inspect_da.R | 4 ++- R/model_da.R | 2 +- R/quadrature.R | 58 ++++++++++++++++++++++++++++++------------ R/utils_da.R | 47 ++++++++++++++++++++++++++++++---- 6 files changed, 137 insertions(+), 44 deletions(-) diff --git a/R/est_lms.R b/R/est_lms.R index 5bec7fc0..5e5165fc 100644 --- a/R/est_lms.R +++ b/R/est_lms.R @@ -408,11 +408,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 +440,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 +619,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 +701,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 +733,7 @@ fischerScoring <- function(par, iterations = i )) } - + alpha <- max.alpha success <- FALSE @@ -740,7 +741,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..c76f80ed 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..3698311e 100644 --- a/R/quadrature.R +++ b/R/quadrature.R @@ -10,6 +10,7 @@ quadrature <- function(m, k, adaptive = FALSE, quad.range = Inf, adaptive.frequency = 3, + wtol = NULL, ...) { if (quad.range < 0) { mod_msg_warn("`quad.range` should be positive, using `-quad.range` instead!\n") @@ -45,6 +46,13 @@ quadrature <- function(m, k, nodes <- sqrt(2) * nodes weights <- weights * pi ^ (-k/2) + if (!is.null(wtol) && is.finite(wtol)) { + idx <- order(weights) + keep <- idx[cumsum(weights[idx]) >= wtol] + nodes <- nodes[keep, , drop=FALSE] + weights <- weights[keep] + } + list( n = nodes, w = weights, @@ -66,10 +74,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 +93,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 +113,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,6 +125,15 @@ adaptiveGaussQuadrature <- function(fun, quadn <- as.matrix(do.call(expand.grid, NODES)) quadw <- Reduce(kronecker, rev(WEIGHTS)) + + # prune by weigths first, avoiding expensive call to fun() + if (secondary.pruning && !is.null(wtol) && is.finite(wtol)) { + idx <- order(quadw) + keep <- idx[cumsum(quadw[idx]) >= wtol] + quadn <- quadn[keep, , drop=FALSE] + quadw <- quadw[keep] + } + quadf <- fun(quadn, ...) weighted <- sweep(quadf, 2, quadw, "*") integral <- quadratureLogLik(weighted) @@ -140,12 +160,15 @@ adaptiveGaussQuadrature <- function(fun, 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 +182,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,17 +193,16 @@ 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) - quad <- quadrature(m = m.ceil, k = 1) + m.ceil <- max(m, min(m.ceil, node.max)) + quad <- quadrature(m = m.ceil, k = 1, wtol = wtol) + 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 @@ -229,6 +252,7 @@ adaptiveGaussQuadratureK <- function(fun, iter.max = iter.max, node.max = node.max, mdiff.tol = mdiff.tol, + wtol = wtol, ... )) } @@ -248,13 +272,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 ) } 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 + ) +} From f87499eb4645950b50e886eba943890b75e6fdea Mon Sep 17 00:00:00 2001 From: kss2k Date: Sun, 25 Oct 2026 17:47:08 +0100 Subject: [PATCH 2/7] Improve quasi adaptive quadrature tuning --- R/quadrature.R | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/R/quadrature.R b/R/quadrature.R index 3698311e..e79913b9 100644 --- a/R/quadrature.R +++ b/R/quadrature.R @@ -223,18 +223,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) { From 0e69ce73e88fa6f1592a1d2696e6c48ed83dfda4 Mon Sep 17 00:00:00 2001 From: kss2k Date: Sun, 25 Oct 2026 17:50:42 +0100 Subject: [PATCH 3/7] Don't check negative ll delta when quadrature grid changes --- R/est_lms.R | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/R/est_lms.R b/R/est_lms.R index 5e5165fc..b4c18ec9 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-18) { if (verbose) cat("\n") mod_msg_warn_immediate( sprintf("Loglikelihood decreased by %.2g!", deltaLL), From 86c708e61eb353558381dbcf22241462ad8ba36d Mon Sep 17 00:00:00 2001 From: kss2k Date: Sun, 25 Oct 2026 17:54:10 +0100 Subject: [PATCH 4/7] Remove whitespace from formatting --- R/generics_modsem_da.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/generics_modsem_da.R b/R/generics_modsem_da.R index c76f80ed..df266d32 100644 --- a/R/generics_modsem_da.R +++ b/R/generics_modsem_da.R @@ -482,7 +482,7 @@ print.summary_da <- function(x, digits = 3, ...) { eff.values <- c( "Effective", - paste0(formatCount(quad$eff.nodes.per.dim), collapse = ", "), + paste0(formatCount(quad$eff.nodes.per.dim), collapse = ","), formatCount(quad$dimensions), formatCount(quad$eff.nodes.total) ) From 17c39a7abac4691db8f5562b50d85ada8aa46af0 Mon Sep 17 00:00:00 2001 From: kss2k Date: Sun, 25 Oct 2026 20:18:33 +0100 Subject: [PATCH 5/7] fix typo --- R/est_lms.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/est_lms.R b/R/est_lms.R index b4c18ec9..624519e7 100644 --- a/R/est_lms.R +++ b/R/est_lms.R @@ -391,7 +391,7 @@ emLms <- function(model, # 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-18) { + if (!quadChanged && deltaLL < -1e-8) { if (verbose) cat("\n") mod_msg_warn_immediate( sprintf("Loglikelihood decreased by %.2g!", deltaLL), From 84ad577ebe019aa6c15788a151f97b0aecdb14fa Mon Sep 17 00:00:00 2001 From: kss2k Date: Mon, 26 Oct 2026 09:28:08 +0100 Subject: [PATCH 6/7] Check for outliers in quadarture --- R/quadrature.R | 94 ++++++++++++++++-------- tests/testthat/test_quadrature_pruning.R | 76 +++++++++++++++++++ 2 files changed, 138 insertions(+), 32 deletions(-) create mode 100644 tests/testthat/test_quadrature_pruning.R diff --git a/R/quadrature.R b/R/quadrature.R index e79913b9..6789f00e 100644 --- a/R/quadrature.R +++ b/R/quadrature.R @@ -10,7 +10,6 @@ quadrature <- function(m, k, adaptive = FALSE, quad.range = Inf, adaptive.frequency = 3, - wtol = NULL, ...) { if (quad.range < 0) { mod_msg_warn("`quad.range` should be positive, using `-quad.range` instead!\n") @@ -46,13 +45,6 @@ quadrature <- function(m, k, nodes <- sqrt(2) * nodes weights <- weights * pi ^ (-k/2) - if (!is.null(wtol) && is.finite(wtol)) { - idx <- order(weights) - keep <- idx[cumsum(weights[idx]) >= wtol] - nodes <- nodes[keep, , drop=FALSE] - weights <- weights[keep] - } - list( n = nodes, w = weights, @@ -126,25 +118,10 @@ adaptiveGaussQuadrature <- function(fun, quadn <- as.matrix(do.call(expand.grid, NODES)) quadw <- Reduce(kronecker, rev(WEIGHTS)) - # prune by weigths first, avoiding expensive call to fun() - if (secondary.pruning && !is.null(wtol) && is.finite(wtol)) { - idx <- order(quadw) - keep <- idx[cumsum(quadw[idx]) >= wtol] - quadn <- quadn[keep, , drop=FALSE] - quadw <- quadw[keep] - } - - 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 @@ -154,6 +131,13 @@ 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( @@ -194,18 +178,17 @@ adaptiveGaussQuadratureK <- function(fun, m.ceil <- round(estMForNodesInRange(m, a = -5, b = 5)) m.ceil <- max(m, min(m.ceil, node.max)) - quad <- quadrature(m = m.ceil, k = 1, wtol = wtol) + quad <- quadrature(m = m.ceil, k = 1) quadw <- quad$w if (K > 1) { # more computationally efficient/scalable quadn <- matrix(0, nrow = length(quadw), ncol = K) quadn[,k] <- quad$n } else quadn <- quad$n - - quadf <- fun(quadn, ...) - 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 @@ -362,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/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)) From b9775d42f69ba94df77ea937bb5aaea6efe7fe7e Mon Sep 17 00:00:00 2001 From: kss2k Date: Mon, 26 Oct 2026 09:32:40 +0100 Subject: [PATCH 7/7] Update `vignettes/estimation_lms.Rmd` --- vignettes/estimation_lms.Rmd | 6 ++++++ 1 file changed, 6 insertions(+) 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")`.