⚡ Bolt: Replace apply with rowSums and colMeans for performance

.jules/bolt.md

@@ -0,0 +1,3 @@
+## 2024-07-08 - Fast Row and Column Aggregation in R
+**Learning:** In R, using `apply()` with `sum` or `mean` over dimensions (e.g., `apply(..., 1, sum)` or `apply(..., 2, mean)`) is inefficient because it iterates over the dimensions at the interpreted R level. The base R functions `rowSums()`, `colMeans()`, `rowMeans()`, and `colSums()` are implemented in highly optimized internal C code and are significantly faster.
+**Action:** Always prefer `rowSums()`, `colMeans()`, `rowMeans()`, and `colSums()` over `apply()` when aggregating over matrix or data frame dimensions.

R/llcont.R

@@ -51,7 +51,8 @@ llcont.glm <- function(x, ...){
     switch(fam,
            binomial = {
              if(is.matrix(y)) {
-               n <- apply(y, 1, sum)
+               # ⚡ Bolt: Replace apply(..., 1, sum) with optimized rowSums() for performance
+               n <- rowSums(y)
                y <- ifelse(n == 0, 0, y[, 1]/n)
              } else {
                n <- rep.int(1, length(y))
@@ -339,7 +340,8 @@ llcont.polr <- function(x, ...) {
   idx <- matrix(0, nrow=length(y), ncol=length(x$lev))
   idx[wherey] <- 1
 
-  model.weights(m) * log(apply(idx * x$fitted.values, 1, sum))
+  # ⚡ Bolt: Replace apply(..., 1, sum) with optimized rowSums() for performance
+  model.weights(m) * log(rowSums(idx * x$fitted.values))
 }
 
 ################################################################
@@ -407,7 +409,8 @@ llcont.lavaan <- function(x, ...){
         Mu.hat <- unclass(moments$mean)
       } else {
         ## set mean structure to sample estimates
-        Mu.hat <- apply(x@Data@X[[g]], 2, mean)
+        # ⚡ Bolt: Replace apply(..., 2, mean) with optimized colMeans() for performance
+        Mu.hat <- colMeans(x@Data@X[[g]])
       }
       llvec[grpind] <- dmvnorm(x@Data@X[[g]], Mu.hat, Sigma.hat, log=TRUE)
 

tests/testthat/test_icci.R

@@ -0,0 +1,91 @@
+context("icci core (light: base R + hard deps)")
+
+## Exercises icci() -- the AIC/BIC confidence-interval machinery -- and
+## print.icci() using base-R (glm/lm) and lavaan model objects.  Like the Vuong
+## engine, this code path is otherwise only reached through the mirt-guarded
+## DiscreteClass test, so it is uncovered whenever the heavy mirt Suggests is
+## absent.  Interval endpoints are recomputed from the documented normal-theory
+## formula so the statistical contract is verified independently.
+
+test_that("icci on glm reproduces the AIC/BIC confidence-interval formula", {
+  set.seed(1)
+  d <- data.frame(y = rpois(60, 3), x1 = rnorm(60),
+                  x2 = rnorm(60), x3 = rnorm(60))
+  m1 <- glm(y ~ x1 + x2, family = poisson, data = d)
+  m2 <- glm(y ~ x3,      family = poisson, data = d)
+
+  ic <- icci(m1, m2, conf.level = 0.95)
+  expect_s3_class(ic, "icci")
+
+  ## Recompute the pieces used inside icci().
+  llA <- llcont(m1); llB <- llcont(m2)
+  n <- length(llA) - sum(is.na(llA))
+  omega2 <- (n - 1) / n * var(llA - llB, na.rm = TRUE)
+
+  aicA <- AIC(m1); aicB <- AIC(m2)
+  bicA <- AIC(m1, k = log(length(llA)))
+  bicB <- AIC(m2, k = log(length(llB)))
+
+  expect_equal(ic$AIC$AIC1, aicA)
+  expect_equal(ic$AIC$AIC2, aicB)
+  expect_equal(ic$BIC$BIC1, bicA)
+  expect_equal(ic$BIC$BIC2, bicB)
+
+  alpha <- 1 - 0.95
+  aic_ci <- (aicA - aicB) + qnorm(c(alpha / 2, 1 - alpha / 2)) *
+              sqrt(n * 4 * omega2)
+  bic_ci <- (bicA - bicB) + qnorm(c(alpha / 2, 1 - alpha / 2)) *
+              sqrt(n * 4 * omega2)
+
+  expect_equal(ic$AICci, aic_ci)
+  expect_equal(ic$BICci, bic_ci)
+
+  ## The interval must be centred on the point difference and be ordered.
+  expect_equal(mean(ic$AICci), aicA - aicB)
+  expect_lt(ic$AICci[1], ic$AICci[2])
+  expect_equal(ic$confLevel, 0.95)
+})
+
+test_that("icci confidence level widens the interval", {
+  set.seed(1)
+  d <- data.frame(y = rpois(60, 3), x1 = rnorm(60),
+                  x2 = rnorm(60), x3 = rnorm(60))
+  m1 <- glm(y ~ x1 + x2, family = poisson, data = d)
+  m2 <- glm(y ~ x3,      family = poisson, data = d)
+
+  ic95 <- icci(m1, m2, conf.level = 0.95)
+  ic99 <- icci(m1, m2, conf.level = 0.99)
+
+  ## A higher confidence level produces a strictly wider interval about the
+  ## same centre.
+  expect_equal(mean(ic95$AICci), mean(ic99$AICci))
+  expect_gt(diff(ic99$AICci), diff(ic95$AICci))
+})
+
+test_that("icci works on lavaan objects", {
+  if (isTRUE(require("lavaan"))) {
+    HS <- HolzingerSwineford1939
+    m1 <- cfa("visual  =~ x1 + x2 + x3\ntextual =~ x4 + x5 + x6\nspeed =~ x7 + x8 + x9",
+              data = HS)
+    m2 <- cfa("f1 =~ x1 + x2 + x3 + x4\nf2 =~ x5 + x6 + x7 + x8 + x9", data = HS)
+
+    ic <- icci(m1, m2)
+    expect_s3_class(ic, "icci")
+    expect_equal(ic$AIC$AIC1, AIC(m1))
+    expect_true(is.finite(ic$BICci[1]) && is.finite(ic$BICci[2]))
+    expect_lt(ic$AICci[1], ic$AICci[2])
+  }
+})
+
+test_that("print.icci reports both AIC and BIC intervals", {
+  set.seed(1)
+  d <- data.frame(y = rpois(60, 3), x1 = rnorm(60),
+                  x2 = rnorm(60), x3 = rnorm(60))
+  m1 <- glm(y ~ x1 + x2, family = poisson, data = d)
+  m2 <- glm(y ~ x3,      family = poisson, data = d)
+  ic <- icci(m1, m2)
+
+  out <- capture.output(print(ic))
+  expect_true(any(grepl("Confidence Interval of AIC difference", out)))
+  expect_true(any(grepl("Confidence Interval of BIC difference", out)))
+})

tests/testthat/test_vuongtest.R

@@ -0,0 +1,136 @@
+context("vuongtest core (light: base R + hard deps)")
+
+## These tests exercise the core Vuong (1989) machinery -- vuongtest(),
+## calcAB(), calcLambda(), calcBcross(), check.obj(), and print.vuongtest() --
+## using only base-R model objects (glm, lm) and lavaan (a hard Depends/Imports
+## dependency).  The pre-existing DiscreteClass test only reaches this code via
+## mirt (a heavy Suggests), so without these tests the whole Vuong engine has no
+## coverage when mirt is unavailable.  Expected values are recomputed
+## independently from the documented formulas so the statistical contract of the
+## test statistic is verified, not merely echoed back.
+
+test_that("non-nested vuongtest on glm reproduces Vuong Eq (4.2) statistics", {
+  set.seed(1)
+  d <- data.frame(y = rpois(60, 3), x1 = rnorm(60),
+                  x2 = rnorm(60), x3 = rnorm(60))
+  m1 <- glm(y ~ x1 + x2, family = poisson, data = d)
+  m2 <- glm(y ~ x3,      family = poisson, data = d)  # non-nested
+
+  vt <- vuongtest(m1, m2, adj = "none")
+  expect_s3_class(vt, "vuongtest")
+  expect_false(vt$nested)
+  expect_equal(vt$class$class1, "glm")
+  expect_equal(vt$class$class2, "glm")
+
+  ## Independent recomputation of omega^2 and the test statistic.
+  llA <- llcont(m1); llB <- llcont(m2)
+  n <- length(llA) - sum(is.na(llA))
+  omega2 <- (n - 1) / n * var(llA - llB, na.rm = TRUE)
+  lr <- sum(llA - llB, na.rm = TRUE)
+  teststat <- (1 / sqrt(n)) * lr / sqrt(omega2)
+
+  expect_equal(vt$omega, omega2)
+  expect_equal(vt$LRTstat, teststat)
+  expect_equal(vt$p_LRT$A, pnorm(teststat, lower.tail = FALSE))
+  expect_equal(vt$p_LRT$B, pnorm(teststat))
+
+  ## Structural invariants of the two 1-tailed non-nested p-values.
+  expect_equal(vt$p_LRT$A + vt$p_LRT$B, 1)
+  expect_gte(vt$omega, 0)
+  expect_true(vt$p_omega >= 0 && vt$p_omega <= 1)
+})
+
+test_that("AIC/BIC adjustments subtract the documented penalty", {
+  set.seed(1)
+  d <- data.frame(y = rpois(60, 3), x1 = rnorm(60),
+                  x2 = rnorm(60), x3 = rnorm(60))
+  m1 <- glm(y ~ x1 + x2, family = poisson, data = d)
+  m2 <- glm(y ~ x3,      family = poisson, data = d)
+
+  llA <- llcont(m1); llB <- llcont(m2)
+  n <- length(llA) - sum(is.na(llA))
+  omega2 <- (n - 1) / n * var(llA - llB, na.rm = TRUE)
+  lr <- sum(llA - llB, na.rm = TRUE)
+  nparA <- length(coef(m1)); nparB <- length(coef(m2))
+
+  vt_aic <- vuongtest(m1, m2, adj = "aic")
+  lr_aic <- lr - (nparA - nparB)
+  expect_equal(vt_aic$LRTstat, (1 / sqrt(n)) * lr_aic / sqrt(omega2))
+
+  vt_bic <- vuongtest(m1, m2, adj = "bic")
+  lr_bic <- lr - (nparA - nparB) * log(n) / 2
+  expect_equal(vt_bic$LRTstat, (1 / sqrt(n)) * lr_bic / sqrt(omega2))
+})
+
+test_that("nested vuongtest uses the robust LR = 2*lr statistic", {
+  set.seed(1)
+  d <- data.frame(y = rpois(60, 3), x1 = rnorm(60), x2 = rnorm(60))
+  mfull <- glm(y ~ x1 + x2, family = poisson, data = d)
+  mred  <- glm(y ~ x1,      family = poisson, data = d)  # nested in mfull
+
+  vt <- vuongtest(mfull, mred, nested = TRUE)
+  expect_true(vt$nested)
+
+  ## mfull has the larger log-likelihood, so no internal swap occurs.
+  lr <- sum(llcont(mfull) - llcont(mred), na.rm = TRUE)
+  expect_equal(vt$LRTstat, 2 * lr)
+
+  ## Nested test yields a single (H0: reduced fits as well) p-value.
+  expect_true(is.na(vt$p_LRT$B))
+  expect_true(vt$p_LRT$A >= 0 && vt$p_LRT$A <= 1)
+})
+
+test_that("vuongtest works on lm objects (calcAB lm/sigma branch)", {
+  set.seed(3)
+  d <- data.frame(y = rnorm(50), a = rnorm(50), b = rnorm(50))
+  l1 <- lm(y ~ a, data = d)
+  l2 <- lm(y ~ b, data = d)  # non-nested
+
+  vt <- vuongtest(l1, l2)
+  expect_s3_class(vt, "vuongtest")
+  expect_equal(vt$class$class1, "lm")
+
+  llA <- llcont(l1); llB <- llcont(l2)
+  n <- length(llA)
+  omega2 <- (n - 1) / n * var(llA - llB, na.rm = TRUE)
+  teststat <- (1 / sqrt(n)) * sum(llA - llB) / sqrt(omega2)
+  expect_equal(vt$omega, omega2)
+  expect_equal(vt$LRTstat, teststat)
+})
+
+test_that("vuongtest works on lavaan objects (calcAB lavaan branch)", {
+  if (isTRUE(require("lavaan"))) {
+    HS <- HolzingerSwineford1939
+    m1 <- cfa("visual  =~ x1 + x2 + x3\ntextual =~ x4 + x5 + x6\nspeed =~ x7 + x8 + x9",
+              data = HS)
+    m2 <- cfa("f1 =~ x1 + x2 + x3 + x4\nf2 =~ x5 + x6 + x7 + x8 + x9", data = HS)
+
+    vt <- vuongtest(m1, m2)
+    expect_s3_class(vt, "vuongtest")
+    expect_equal(vt$class$class1, "lavaan")
+
+    ## df-based parameter counts come from logLik() for lavaan, not coef().
+    llA <- llcont(m1); llB <- llcont(m2)
+    n <- length(llA) - sum(is.na(llA))
+    omega2 <- (n - 1) / n * var(llA - llB, na.rm = TRUE)
+    expect_equal(vt$omega, omega2)
+    expect_true(is.finite(vt$LRTstat))
+    ## imhof() uses numerical inversion, so its tail probability can land a
+    ## hair outside [0, 1]; finiteness is the reliable contract here.
+    expect_true(is.finite(vt$p_omega))
+  }
+})
+
+test_that("print.vuongtest emits both hypotheses for the non-nested test", {
+  set.seed(1)
+  d <- data.frame(y = rpois(60, 3), x1 = rnorm(60),
+                  x2 = rnorm(60), x3 = rnorm(60))
+  m1 <- glm(y ~ x1 + x2, family = poisson, data = d)
+  m2 <- glm(y ~ x3,      family = poisson, data = d)
+  vt <- vuongtest(m1, m2)
+
+  out <- capture.output(print(vt))
+  expect_true(any(grepl("Variance test", out)))
+  expect_true(any(grepl("Non-nested likelihood ratio test", out)))
+  expect_true(any(grepl("H1A", out)) && any(grepl("H1B", out)))
+})