From b12aef20db21a33c08a3fa2e820c83e4fc63a17b Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Fri, 3 May 2024 12:42:15 +0200
Subject: [PATCH 1/8] Add fmg test to lavTest.

---
 R/lav_test.R     |  47 +++++++----
 R/lav_test_fmg.R | 201 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 231 insertions(+), 17 deletions(-)
 create mode 100644 R/lav_test_fmg.R

diff --git a/R/lav_test.R b/R/lav_test.R
index e2ca95f0..5787a486 100644
--- a/R/lav_test.R
+++ b/R/lav_test.R
@@ -168,24 +168,33 @@ lav_test_rename <- function(test, check = FALSE) {
     test[target.idx] <- "browne.residual.nt.model"
   }
 
+  # report unknown values
+  bad.idx <- which(!test %in% c(
+    "standard", "none", "default",
+    "satorra.bentler",
+    "yuan.bentler",
+    "yuan.bentler.mplus",
+    "mean.adjusted",
+    "mean.var.adjusted",
+    "scaled.shifted",
+    "bollen.stine",
+    "browne.residual.nt",
+    "browne.residual.nt.model",
+    "browne.residual.adf",
+    "browne.residual.adf.model"
+  ))
+
+  fmgs = c()
+  for (j in seq_along(bad.idx)) {
+    fmg_parsed = lav_fmg_parse_input(test[bad.idx[j]])
+    if (!is.null(fmg_parsed)) {
+      fmgs = c(fmgs, test[bad.idx[j]])
+      bad.idx = bad.idx[-j]
+    }
+  }
 
   # check?
   if (check) {
-    # report unknown values
-    bad.idx <- which(!test %in% c(
-      "standard", "none", "default",
-      "satorra.bentler",
-      "yuan.bentler",
-      "yuan.bentler.mplus",
-      "mean.adjusted",
-      "mean.var.adjusted",
-      "scaled.shifted",
-      "bollen.stine",
-      "browne.residual.nt",
-      "browne.residual.nt.model",
-      "browne.residual.adf",
-      "browne.residual.adf.model"
-    ))
     if (length(bad.idx) > 0L) {
       lav_msg_stop(sprintf(
         ngettext(
@@ -212,7 +221,7 @@ lav_test_rename <- function(test, check = FALSE) {
     }
   }
 
-  # reorder: first nonscaled, then scaled
+  # reorder: first nonscaled, then scaled, then fmg.
   nonscaled.idx <- which(test %in% c(
     "standard", "none", "default",
     "bollen.stine",
@@ -229,7 +238,9 @@ lav_test_rename <- function(test, check = FALSE) {
     "mean.var.adjusted",
     "scaled.shifted"
   ))
-  test <- c(test[nonscaled.idx], test[scaled.idx])
+
+  scaled <- c(test[scaled.idx], sort(fmgs))
+  test <- c(test[nonscaled.idx], scaled)
 
   test
 }
@@ -684,6 +695,8 @@ lav_model_test <- function(lavobject = NULL,
         boot.larger = boot.larger,
         boot.length = boot.length
       )
+    } else if (!is.null(input <- lav_fmg_parse_input(this.test))) {
+      TEST[[this.test]] <- lav_test_fmg(lavobject, input)
     }
   } # additional tests
 
diff --git a/R/lav_test_fmg.R b/R/lav_test_fmg.R
new file mode 100644
index 00000000..2329b8e9
--- /dev/null
+++ b/R/lav_test_fmg.R
@@ -0,0 +1,201 @@
+lav_rls <- \(object) {
+  s_inv <- chol2inv(chol(lavaan::lavInspect(object, "sigma.hat")))
+  residuals <- lavaan::lavInspect(object, "residuals")
+  mm <- residuals[[1]] %*% s_inv
+  n <- lavaan::lavInspect(object, "nobs")
+  .5 * n * sum(mm * t(mm))
+}
+
+lav_test_fmg <- \(object, input) {
+
+  make_chisqs <- \(chisq) {
+    if (chisq == "ml") lavaan::fitmeasures(object, "chisq") else rls(object)
+  }
+
+  name <- input$name
+  degree <- input$degree
+  unbiased <- input$unbiased
+  chisq <- if (input$chisq == "ml") lavaan::fitmeasures(object, "chisq") else lav_rls(object)
+  df <- lavaan::fitmeasures(object, "df")
+  ug <- lav_ugamma_no_groups(object, unbiased)
+  lambdas_list <- Re(eigen(ug)$values)[seq(df)]
+
+  lambdas <- lambdas_list
+  list(
+    test = lav_fmg_reconstruct_input(input),
+    pvalue = if (name == "pEBA") {
+    lav_fmg_peba(chisq, lambdas, degree)
+  } else if (name == "EBA") {
+    lav_fmg_eba(chisq, lambdas, degree)
+  } else if (name == "pOLS") {
+    lav_fmg_pols(chisq, lambdas, degree)
+  },
+  label = lav_fmg_reconstruct_label(input))
+}
+
+lav_fmg_parse_input <- \(string) {
+  na_to_null <- \(x) if (is.na(x)) NULL else x
+  default <- \(x) if (x != "") as.numeric(x) else 2
+
+  tryCatch ({
+    string <- tolower(string)
+    splitted <- strsplit(string, "_")[[1]]
+    unbiased <- FALSE
+    chisq <- "ml"
+    type <- na_to_null(splitted[1])
+
+    if (length(splitted) == 3) {
+      unbiased <- if (splitted[2] == "ug") TRUE else FALSE
+      chisq <- splitted[3]
+    } else if (length(splitted) == 2) {
+      if (splitted[2] == "rls" || splitted[2] == "ml") {
+        chisq <- splitted[2]
+      } else if (splitted[2] == "ug") {
+        unbiased <- TRUE
+      }
+    }
+
+    out <- NULL
+    name <- NULL
+    if (startsWith(type, "peba")) {
+      degree <- default(substring(type, 5))
+      name <- "pEBA"
+    } else if (startsWith(type, "eba")) {
+      degree <- default(substring(type, 4))
+      name <- "EBA"
+    } else if (startsWith(type, "pols")) {
+      degree <- default(substring(type, 5))
+      name <- "pOLS"
+    }
+
+    if (is.null(degree) || (chisq != "ml" && chisq != "rls")) {
+      return(NULL)
+    }
+
+    list(
+      degree = degree,
+      name = name,
+      unbiased = unbiased,
+      chisq = chisq
+    )
+  }, error = \(e){
+    NULL
+  })
+}
+
+lav_fmg_reconstruct_input <- \(input) {
+  name <- input$name
+  degree <- input$degree
+  unbiased <- input$unbiased
+  chisq <- input$chisq
+
+  out <- paste0(name, degree)
+  if (unbiased) out <- paste0(out, "_ug")
+  if (chisq == "rls") out <- paste0(out, "_rls")
+  out
+}
+
+lav_fmg_reconstruct_label <- \(input) {
+  name <- tolower(input$name)
+  degree <- input$degree
+  unbiased <- input$unbiased
+  chisq <- input$chisq
+
+  out <- if (name == "peba") {
+    paste0("Penalized EBA with ", degree, " blocks, ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
+    chisq, ".")
+  } else if (name == "eba" && degree == 1) {
+    paste0("Satorra-Bentler with ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
+           chisq, ".")
+  } else if (name == "eba") {
+    paste0("EBA with ", degree, " blocks, ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
+           chisq, ".")
+  } else {
+    paste0("Penalized OLS with weighting parameter ", degree, ", ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
+           chisq, ".")
+  }
+  out
+}
+
+lav_fmg_eba <- \(chisq, lambdas, j) {
+  m <- length(lambdas)
+  k <- ceiling(m / j)
+  eig <- lambdas
+  eig <- c(eig, rep(NA, k * j - length(eig)))
+  dim(eig) <- c(k, j)
+  eig_means <- colMeans(eig, na.rm = TRUE)
+  repeated <- rep(eig_means, each = k)[seq(m)]
+  CompQuadForm::imhof(chisq, repeated)$Qq
+}
+
+lav_fmg_peba <- \(chisq, lambdas, j) {
+  m <- length(lambdas)
+  k <- ceiling(m / j)
+  eig <- lambdas
+  eig <- c(eig, rep(NA, k * j - length(eig)))
+  dim(eig) <- c(k, j)
+  eig_means <- colMeans(eig, na.rm = TRUE)
+  eig_mean <- mean(lambdas)
+  repeated <- rep(eig_means, each = k)[seq(m)]
+  CompQuadForm::imhof(chisq, (repeated + eig_mean) / 2)$Qq
+}
+
+lav_fmg_pols <- \(chisq, lambdas, gamma) {
+  x <- seq_along(lambdas)
+  beta1_hat <- 1 / gamma * stats::cov(x, lambdas) / stats::var(x)
+  beta0_hat <- mean(lambdas) - beta1_hat * mean(x)
+  lambda_hat <- pmax(beta0_hat + beta1_hat * x, 0)
+  CompQuadForm::imhof(chisq, lambda_hat)$Qq
+}
+
+lav_fmg_gamma_unbiased <- \(obj, gamma) {
+  gamma_est_nt <- \(sigma) {
+    lower_vec_indices <- \(n = 1L, diagonal = TRUE) {
+      rows <- matrix(seq_len(n), n, n)
+      cols <- matrix(seq_len(n), n, n, byrow = TRUE)
+      if (diagonal) which(rows >= cols) else which(rows > cols)
+    }
+
+    upper_vec_indices <- \(n = 1L, diagonal = TRUE) {
+      rows <- matrix(seq_len(n), n, n)
+      cols <- matrix(seq_len(n), n, n, byrow = TRUE)
+      tmp <- matrix(seq_len(n * n), n, n, byrow = TRUE)
+      if (diagonal) tmp[rows >= cols] else tmp[rows > cols]
+    }
+
+    n <- ncol(sigma)
+    lower <- lower_vec_indices(n)
+    upper <- upper_vec_indices(n)
+    y <- sigma %x% sigma
+    out <- (y[lower, , drop = FALSE] + y[upper, , drop = FALSE]) / 2
+    out[, lower, drop = FALSE] + out[, upper, drop = FALSE]
+  }
+
+  gamma_est_unbiased <- \(x, n = NULL, sigma = NULL, gamma_adf = NULL, gamma_nt = NULL) {
+    vech <- \(x) x[row(x) >= col(x)]
+    if (!missing(x)) n <- nrow(x)
+    sigma <- if (is.null(sigma)) stats::cov(x) * (n - 1) / n else sigma
+    gamma_adf <- if (is.null(gamma_adf)) gamma_est_adf(x) else gamma_adf
+    gamma_nt <- if (is.null(gamma_nt)) gamma_est_nt(sigma) else gamma_nt
+    gamma_rem <- tcrossprod(vech(sigma))
+    mult <- n / ((n - 2) * (n - 3))
+    mult * ((n - 1) * gamma_adf - (gamma_nt - 2 / (n - 1) * gamma_rem))
+  }
+
+  gamma_est_unbiased(
+    n = lavaan::lavInspect(obj, "nobs"),
+    sigma = obj@SampleStats@cov[[1]],
+    gamma_adf = gamma,
+    gamma_nt = NULL
+  )
+}
+
+lav_ugamma_no_groups <- \(object, unbiased) {
+  u <- lavaan::lavInspect(object, "U")
+  bias <- object@Options$gamma.unbiased
+  object@Options$gamma.unbiased <- FALSE
+  gamma <- lavaan::lavInspect(object, "gamma")
+  object@Options$gamma.unbiased <- bias
+  if (unbiased) gamma <- lav_fmg_gamma_unbiased(object, gamma)
+  u %*% gamma
+}

From acfd3987dd2dd3618a9000b4fb46654a267eef5e Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Thu, 27 Jun 2024 14:17:56 +0200
Subject: [PATCH 2/8] :hammer: Jobb jobb jobb

---
 R/lav_imhof.R    | 17 ++++++++++++++
 R/lav_test_fmg.R | 16 ++++++++++---
 lav_fmg_greier.r | 59 ++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 89 insertions(+), 3 deletions(-)
 create mode 100644 R/lav_imhof.R
 create mode 100644 lav_fmg_greier.r

diff --git a/R/lav_imhof.R b/R/lav_imhof.R
new file mode 100644
index 00000000..11b1284b
--- /dev/null
+++ b/R/lav_imhof.R
@@ -0,0 +1,17 @@
+imhof <- \(x, lambda) {
+  theta <- \(u, x, lambda) 0.5 * (sum(atan(lambda * u)) - x * u)
+  rho <- \(u, lambda) exp(1 / 4 * sum(log(1 + lambda^2 * u^2)))
+  integrand <- Vectorize(\(u) {
+    sin(theta(u, x, lambda)) / (u * rho(u, lambda))
+  })
+  z <- integrate(integrand, lower = 0, upper = Inf)$value
+  0.5 + z / pi
+}
+
+# k <- 50
+# lambda <- runif(k)
+# x <- k/2
+# imhof(x, lambda)
+# CompQuadForm::imhof(x, lambda)
+
+# microbenchmark::microbenchmark(imhof(x, lambda), CompQuadForm::imhof(x, lambda))
diff --git a/R/lav_test_fmg.R b/R/lav_test_fmg.R
index 2329b8e9..16d2128e 100644
--- a/R/lav_test_fmg.R
+++ b/R/lav_test_fmg.R
@@ -125,7 +125,7 @@ lav_fmg_eba <- \(chisq, lambdas, j) {
   dim(eig) <- c(k, j)
   eig_means <- colMeans(eig, na.rm = TRUE)
   repeated <- rep(eig_means, each = k)[seq(m)]
-  CompQuadForm::imhof(chisq, repeated)$Qq
+  lav_fmg_imhof(chisq, repeated)
 }
 
 lav_fmg_peba <- \(chisq, lambdas, j) {
@@ -137,7 +137,7 @@ lav_fmg_peba <- \(chisq, lambdas, j) {
   eig_means <- colMeans(eig, na.rm = TRUE)
   eig_mean <- mean(lambdas)
   repeated <- rep(eig_means, each = k)[seq(m)]
-  CompQuadForm::imhof(chisq, (repeated + eig_mean) / 2)$Qq
+  lav_fmg_imhof(chisq, (repeated + eig_mean) / 2)
 }
 
 lav_fmg_pols <- \(chisq, lambdas, gamma) {
@@ -145,7 +145,7 @@ lav_fmg_pols <- \(chisq, lambdas, gamma) {
   beta1_hat <- 1 / gamma * stats::cov(x, lambdas) / stats::var(x)
   beta0_hat <- mean(lambdas) - beta1_hat * mean(x)
   lambda_hat <- pmax(beta0_hat + beta1_hat * x, 0)
-  CompQuadForm::imhof(chisq, lambda_hat)$Qq
+  lav_fmg_imhof(chisq, lambda_hat)
 }
 
 lav_fmg_gamma_unbiased <- \(obj, gamma) {
@@ -199,3 +199,13 @@ lav_ugamma_no_groups <- \(object, unbiased) {
   if (unbiased) gamma <- lav_fmg_gamma_unbiased(object, gamma)
   u %*% gamma
 }
+
+lav_fmg_imhof <- \(x, lambda) {
+  theta <- \(u, x, lambda) 0.5 * (sum(atan(lambda * u)) - x * u)
+  rho <- \(u, lambda) exp(1 / 4 * sum(log(1 + lambda^2 * u^2)))
+  integrand <- Vectorize(\(u) {
+    sin(theta(u, x, lambda)) / (u * rho(u, lambda))
+  })
+  z <- integrate(integrand, lower = 0, upper = Inf)$value
+  0.5 + z / pi
+}
diff --git a/lav_fmg_greier.r b/lav_fmg_greier.r
new file mode 100644
index 00000000..af3fbaaa
--- /dev/null
+++ b/lav_fmg_greier.r
@@ -0,0 +1,59 @@
+n = 500
+data = psych::bfi
+
+model = "A =~ A1+A2+A3+A4+A5;
+         C =~ C1+C2+C3+C4+C5"
+object_chisq <- lavaan::sem(model, data[1:n, 1:10], test = "sb")
+
+
+#lavGamma(object_chisq, gamma.unbiased = TRUE)
+gamma_biased <- lav_samplestats_Gamma(data[1:n, 1:10])
+gamma_yves <- lav_samplestats_Gamma(data[1:n, 1:10], unbiased = TRUE)
+gamma_fmg <- lav_fmg_gamma_unbiased(object_chisq, gamma_biased)
+sum(abs(gamma_yves-gamma_fmg))
+
+cov_yves = object_chisq@SampleStats@cov[[1]]
+cov_fmg = cov(data[1:n, 1:10], use = "pairwise") * (n-1) / n
+sum(abs(cov_yves-cov_fmg))
+
+
+#object_rls <- lavaan::sem(model, psych::bfi[1:n, 1:10], test = "sb", scaled.test = "Browne.residual.nt")
+
+
+## UNBIASED
+
+GammaNT.cov <- 2 * lav_matrix_duplication_ginv_pre_post(COV %x% COV)
+
+Gamma.u <- (N * (N - 1) / (N - 2) / (N - 3) * Gamma.cov -
+              N / (N - 2) / (N - 3) * (GammaNT.cov -
+                                         2 / (N - 1) * tcrossprod(cov.vech)))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+lavaan::lavTest(object_chisq, test = "sb", scaled.test = "Browne.residual.nt")
+
+
+lavTest(object_chisq, scaled.test = "browne.residual.nt", test = "sb")
+
+lav_model_test(lavobject = object_chisq)
+lav_model_test(lavobject = object_rls)
+
+test.orig <- object_chisq@Options$test

From 6c0003a82f0c89a261ad981eb0e26a43c01e3290 Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Fri, 28 Jun 2024 12:27:05 +0200
Subject: [PATCH 3/8] :pear: Work.

---
 R/lav_options.R  |  12 ++-
 R/lav_test_fmg.R | 276 +++++++++++++++++++++++++----------------------
 lav_fmg_greier.r |  15 ++-
 3 files changed, 171 insertions(+), 132 deletions(-)

diff --git a/R/lav_options.R b/R/lav_options.R
index 55d496a1..04aee226 100644
--- a/R/lav_options.R
+++ b/R/lav_options.R
@@ -880,6 +880,16 @@ lav_options_set <- function(opt = NULL) {                     # nolint
     "browne.residual.adf.model",
     "bollen.stine"
   ))
+
+  fmgs = c()
+  for (j in seq_along(wrong.idx)) {
+    fmg_parsed = lav_fmg_parse_input(opt$test[wrong.idx[j]])
+    if (!is.null(fmg_parsed)) {
+      fmgs = c(fmgs, opt$test[wrong.idx[j]])
+      wrong.idx = wrong.idx[-j]
+    }
+  }
+
   if (length(wrong.idx) > 0L) {
     lav_msg_stop(gettextf(
     "invalid option(s) for test argument: %1$s. Possible options are: %2$s.",
@@ -889,7 +899,7 @@ lav_options_set <- function(opt = NULL) {                     # nolint
                     "browne.residual.nt.model", "satorra.bentler",
                     "yuan.bentler", "yuan.bentler.mplus",
                     "mean.var.adjusted", "scaled.shifted",
-                    "bollen.stine"), log.sep = "or")))
+                    "bollen.stine", "fmg", "fmgols"), log.sep = "or")))
   }
 
   # bounds
diff --git a/R/lav_test_fmg.R b/R/lav_test_fmg.R
index 16d2128e..f4df097a 100644
--- a/R/lav_test_fmg.R
+++ b/R/lav_test_fmg.R
@@ -13,7 +13,7 @@ lav_test_fmg <- \(object, input) {
   }
 
   name <- input$name
-  degree <- input$degree
+  param <- input$param
   unbiased <- input$unbiased
   chisq <- if (input$chisq == "ml") lavaan::fitmeasures(object, "chisq") else lav_rls(object)
   df <- lavaan::fitmeasures(object, "df")
@@ -24,111 +24,57 @@ lav_test_fmg <- \(object, input) {
   list(
     test = lav_fmg_reconstruct_input(input),
     pvalue = if (name == "pEBA") {
-    lav_fmg_peba(chisq, lambdas, degree)
+    lav_fmg_peba(chisq, lambdas, param)
   } else if (name == "EBA") {
-    lav_fmg_eba(chisq, lambdas, degree)
+    lav_fmg_eba(chisq, lambdas, param)
   } else if (name == "pOLS") {
-    lav_fmg_pols(chisq, lambdas, degree)
+    lav_fmg_pols(chisq, lambdas, param)
   },
   label = lav_fmg_reconstruct_label(input))
 }
 
 lav_fmg_parse_input <- \(string) {
   na_to_null <- \(x) if (is.na(x)) NULL else x
-  default <- \(x) if (x != "") as.numeric(x) else 2
+  default <- \(x) if (x != "") as.numeric(x) else 4
 
   tryCatch ({
     string <- tolower(string)
     splitted <- strsplit(string, "_")[[1]]
-    unbiased <- FALSE
-    chisq <- "ml"
-    type <- na_to_null(splitted[1])
-
-    if (length(splitted) == 3) {
-      unbiased <- if (splitted[2] == "ug") TRUE else FALSE
-      chisq <- splitted[3]
-    } else if (length(splitted) == 2) {
-      if (splitted[2] == "rls" || splitted[2] == "ml") {
-        chisq <- splitted[2]
-      } else if (splitted[2] == "ug") {
-        unbiased <- TRUE
-      }
-    }
-
+    name <- na_to_null(splitted[1])
     out <- NULL
-    name <- NULL
-    if (startsWith(type, "peba")) {
-      degree <- default(substring(type, 5))
-      name <- "pEBA"
-    } else if (startsWith(type, "eba")) {
-      degree <- default(substring(type, 4))
-      name <- "EBA"
-    } else if (startsWith(type, "pols")) {
-      degree <- default(substring(type, 5))
-      name <- "pOLS"
-    }
 
-    if (is.null(degree) || (chisq != "ml" && chisq != "rls")) {
+    if (name != "fmg" && name != "fmgols") {
       return(NULL)
     }
 
+    param <- if (length(splitted) == 2) {
+      as.numeric(splitted[2])
+    } else if (name == "fmg") {
+      4
+    } else {
+      0.5
+    }
+
     list(
-      degree = degree,
-      name = name,
-      unbiased = unbiased,
-      chisq = chisq
+      param = param,
+      name = name
     )
   }, error = \(e){
-    NULL
+    cat(e)
   })
 }
 
-lav_fmg_reconstruct_input <- \(input) {
-  name <- input$name
-  degree <- input$degree
-  unbiased <- input$unbiased
-  chisq <- input$chisq
+lav_fmg_reconstruct_input <- \(input) paste0(input$name, "_", input$param)
 
-  out <- paste0(name, degree)
-  if (unbiased) out <- paste0(out, "_ug")
-  if (chisq == "rls") out <- paste0(out, "_rls")
-  out
-}
-
-lav_fmg_reconstruct_label <- \(input) {
-  name <- tolower(input$name)
-  degree <- input$degree
-  unbiased <- input$unbiased
-  chisq <- input$chisq
-
-  out <- if (name == "peba") {
-    paste0("Penalized EBA with ", degree, " blocks, ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
-    chisq, ".")
-  } else if (name == "eba" && degree == 1) {
-    paste0("Satorra-Bentler with ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
-           chisq, ".")
-  } else if (name == "eba") {
-    paste0("EBA with ", degree, " blocks, ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
-           chisq, ".")
+lav_fmg_construct_label <- \(input) {
+  if (input$name == "fmg") {
+    paste0("FMG with ", input$param, " blocks.")
   } else {
-    paste0("Penalized OLS with weighting parameter ", degree, ", ", if (unbiased) "unbiased" else "biased", " gamma, based on ",
-           chisq, ".")
+    paste0("FMGOLS with weighting parameter ", input$param, ".")
   }
-  out
 }
 
-lav_fmg_eba <- \(chisq, lambdas, j) {
-  m <- length(lambdas)
-  k <- ceiling(m / j)
-  eig <- lambdas
-  eig <- c(eig, rep(NA, k * j - length(eig)))
-  dim(eig) <- c(k, j)
-  eig_means <- colMeans(eig, na.rm = TRUE)
-  repeated <- rep(eig_means, each = k)[seq(m)]
-  lav_fmg_imhof(chisq, repeated)
-}
-
-lav_fmg_peba <- \(chisq, lambdas, j) {
+lav_fmg <- \(chisq, lambdas, j) {
   m <- length(lambdas)
   k <- ceiling(m / j)
   eig <- lambdas
@@ -140,7 +86,7 @@ lav_fmg_peba <- \(chisq, lambdas, j) {
   lav_fmg_imhof(chisq, (repeated + eig_mean) / 2)
 }
 
-lav_fmg_pols <- \(chisq, lambdas, gamma) {
+lav_fmgols <- \(chisq, lambdas, gamma) {
   x <- seq_along(lambdas)
   beta1_hat <- 1 / gamma * stats::cov(x, lambdas) / stats::var(x)
   beta0_hat <- mean(lambdas) - beta1_hat * mean(x)
@@ -148,64 +94,136 @@ lav_fmg_pols <- \(chisq, lambdas, gamma) {
   lav_fmg_imhof(chisq, lambda_hat)
 }
 
-lav_fmg_gamma_unbiased <- \(obj, gamma) {
-  gamma_est_nt <- \(sigma) {
-    lower_vec_indices <- \(n = 1L, diagonal = TRUE) {
-      rows <- matrix(seq_len(n), n, n)
-      cols <- matrix(seq_len(n), n, n, byrow = TRUE)
-      if (diagonal) which(rows >= cols) else which(rows > cols)
-    }
+lav_fmg_imhof <- \(x, lambda) {
+  theta <- \(u, x, lambda) 0.5 * (sum(atan(lambda * u)) - x * u)
+  rho <- \(u, lambda) exp(1 / 4 * sum(log(1 + lambda^2 * u^2)))
+  integrand <- Vectorize(\(u) {
+    sin(theta(u, x, lambda)) / (u * rho(u, lambda))
+  })
+  z <- integrate(integrand, lower = 0, upper = Inf)$value
+  0.5 + z / pi
+}
 
-    upper_vec_indices <- \(n = 1L, diagonal = TRUE) {
-      rows <- matrix(seq_len(n), n, n)
-      cols <- matrix(seq_len(n), n, n, byrow = TRUE)
-      tmp <- matrix(seq_len(n * n), n, n, byrow = TRUE)
-      if (diagonal) tmp[rows >= cols] else tmp[rows > cols]
-    }
 
-    n <- ncol(sigma)
-    lower <- lower_vec_indices(n)
-    upper <- upper_vec_indices(n)
-    y <- sigma %x% sigma
-    out <- (y[lower, , drop = FALSE] + y[upper, , drop = FALSE]) / 2
-    out[, lower, drop = FALSE] + out[, upper, drop = FALSE]
+#' Calculate non-nested ugamma for multiple groups.
+#' @param object A `lavaan` object.
+#' @param unbiased If `TRUE`, uses the unbiased gamma estimate.
+#' @keywords internal
+#' @return Ugamma for non-nested object.
+ugamma_non_nested <- \(object) {
+  lavmodel <- object@Model
+
+  ceq_idx <- attr(lavmodel@con.jac, "ceq.idx")
+  if (length(ceq_idx) > 0L) {
+    stop("Testing of models with groups and equality constraints not supported.")
   }
 
-  gamma_est_unbiased <- \(x, n = NULL, sigma = NULL, gamma_adf = NULL, gamma_nt = NULL) {
-    vech <- \(x) x[row(x) >= col(x)]
-    if (!missing(x)) n <- nrow(x)
-    sigma <- if (is.null(sigma)) stats::cov(x) * (n - 1) / n else sigma
-    gamma_adf <- if (is.null(gamma_adf)) gamma_est_adf(x) else gamma_adf
-    gamma_nt <- if (is.null(gamma_nt)) gamma_est_nt(sigma) else gamma_nt
-    gamma_rem <- tcrossprod(vech(sigma))
-    mult <- n / ((n - 2) * (n - 3))
-    mult * ((n - 1) * gamma_adf - (gamma_nt - 2 / (n - 1) * gamma_rem))
+  test <- list()
+  lavsamplestats <- object@SampleStats
+  lavmodel <- object@Model
+  lavoptions <- object@Options
+  lavimplied <- object@implied
+  lavdata <- object@Data
+  test$standard <- object@test[[1]]
+
+  if (test$standard$df == 0L || test$standard$df < 0) {
+    stop("Df must be > 0.")
   }
 
-  gamma_est_unbiased(
-    n = lavaan::lavInspect(obj, "nobs"),
-    sigma = obj@SampleStats@cov[[1]],
-    gamma_adf = gamma,
-    gamma_nt = NULL
+  e <- lavaan:::lav_model_information(
+    lavmodel = lavmodel,
+    lavimplied = lavimplied,
+    lavsamplestats = lavsamplestats,
+    lavdata = lavdata,
+    lavoptions = lavoptions,
+    extra = TRUE
   )
-}
 
-lav_ugamma_no_groups <- \(object, unbiased) {
-  u <- lavaan::lavInspect(object, "U")
-  bias <- object@Options$gamma.unbiased
-  object@Options$gamma.unbiased <- FALSE
-  gamma <- lavaan::lavInspect(object, "gamma")
-  object@Options$gamma.unbiased <- bias
-  if (unbiased) gamma <- lav_fmg_gamma_unbiased(object, gamma)
-  u %*% gamma
+  delta <- attr(e, "Delta")
+  wls_v <- attr(e, "WLS.V")
+
+  gamma <- lavaan:::lav_object_gamma(object)
+  if (is.null(gamma[[1]])) {
+    gamma <- lapply(lavaan::lavInspect(object, "gamma"), \(x) {
+      class(x) <- "matrix"
+      x
+    })
+  }
+
+  gamma_global <- as.matrix(Matrix::bdiag(gamma))
+  delta_global <- do.call(rbind, delta)
+  v_global <- as.matrix(Matrix::bdiag(wls_v))
+  x <- v_global %*% delta_global
+  u_global <- v_global - crossprod(t(x), solve(t(delta_global) %*% x, t(x)))
+  u_global %*% gamma_global
 }
 
-lav_fmg_imhof <- \(x, lambda) {
-  theta <- \(u, x, lambda) 0.5 * (sum(atan(lambda * u)) - x * u)
-  rho <- \(u, lambda) exp(1 / 4 * sum(log(1 + lambda^2 * u^2)))
-  integrand <- Vectorize(\(u) {
-    sin(theta(u, x, lambda)) / (u * rho(u, lambda))
-  })
-  z <- integrate(integrand, lower = 0, upper = Inf)$value
-  0.5 + z / pi
+#' Calculate nested ugamma.
+#'
+#' This can also be used with restrictions.
+#'
+#' @param m0,m1 Two nested `lavaan` objects.
+#' @param a The `A` matrix. If if `NULL`, gets calculated by
+#'    `lavaan:::lav_test_diff_A` with `method = method`.
+#' @param method Method passed to `lavaan:::lav_test_diff_A`.
+#' @param unbiased If `TRUE`, uses the unbiased gamma estimate.
+#' @keywords internal
+#' @return Ugamma for nested object.
+lav_ugamma_nested <- \(m0, m1, a = NULL, method = "delta", unbiased = FALSE) {
+  m0@Options$gamma.unbiased <- unbiased
+  m1@Options$gamma.unbiased <- unbiased
+
+  # extract information from m1 and m2
+  t1 <- m1@test[[1]]$stat
+  r1 <- m1@test[[1]]$df
+
+  t0 <- m0@test[[1]]$stat
+  r0 <- m0@test[[1]]$df
+
+  # m = difference between the df's
+  m <- r0 - r1
+
+  # check for identical df setting
+  if (m == 0L) {
+    return(list(
+      T.delta = (t0 - t1), scaling.factor = as.numeric(NA),
+      df.delta = m, a = as.numeric(NA), b = as.numeric(NA)
+    ))
+  }
+
+  gamma <- lavaan:::lav_object_gamma(m0) # the same for m1 and m0
+
+  if (is.null(gamma)) {
+    stop("lavaan error: Can not compute gamma matrix; perhaps missing \"ml\"?")
+  }
+
+  wls_v <- lavaan::lavTech(m1, "WLS.V")
+  pi <- lavaan::lavInspect(m1, "delta")
+
+  p_inv <- lavaan::lavInspect(m1, what = "inverted.information")
+
+  if (is.null(a)) {
+    a <- lavaan:::lav_test_diff_A(m1, m0, method = method, reference = "H1")
+    if (m1@Model@eq.constraints) {
+      a <- a %*% t(m1@Model@eq.constraints.K)
+    }
+  }
+
+  paapaap <- p_inv %*% t(a) %*% MASS::ginv(a %*% p_inv %*% t(a)) %*% a %*% p_inv
+
+  # compute scaling factor
+  fg <- unlist(m1@SampleStats@nobs) / m1@SampleStats@ntotal
+
+  # We need the global gamma, cf. eq.~(10) in Satorra (2000).
+  gamma_rescaled <- gamma
+  for (i in (seq_along(gamma))) {
+    gamma_rescaled[[i]] <- fg[i] * gamma_rescaled[[i]]
+  }
+  gamma_global <- as.matrix(Matrix::bdiag(gamma_rescaled))
+  # Also the global V:
+  v_global <- as.matrix(Matrix::bdiag(wls_v))
+  pi_global <- do.call(rbind, pi)
+  # U global version, eq.~(22) in Satorra (2000).
+  u_global <- v_global %*% pi_global %*% paapaap %*% t(pi_global) %*% v_global
+  return(u_global %*% gamma_global)
 }
diff --git a/lav_fmg_greier.r b/lav_fmg_greier.r
index af3fbaaa..d5c9fb3d 100644
--- a/lav_fmg_greier.r
+++ b/lav_fmg_greier.r
@@ -1,10 +1,21 @@
 n = 500
 data = psych::bfi
 
+model = "A =~ A1+b*A2+A3+A4+A5;
+         C =~ C1+b*C2+C3+C4+C5
+         "
+object <- sem(model, data[1:n, ], test = "sb")
+sum(diag(ugamma_non_nested(object))) / 34
+
+dim(lavaan::lavInspect(object, "delta"))
+
 model = "A =~ A1+A2+A3+A4+A5;
-         C =~ C1+C2+C3+C4+C5"
-object_chisq <- lavaan::sem(model, data[1:n, 1:10], test = "sb")
+         C =~ C1+C2+C3+C4+C5
+         "
+object <- sem(model, data[1:n, ], test = "sb")
+dim(lavaan::lavInspect(object, "delta"))
 
+object_chisq <- sem(model, data[1:n, ], test = "sb")
 
 #lavGamma(object_chisq, gamma.unbiased = TRUE)
 gamma_biased <- lav_samplestats_Gamma(data[1:n, 1:10])

From cb99f559e7ac11b4c9bc0ecd17a05229d9eb54d2 Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Mon, 12 Aug 2024 14:04:07 +0200
Subject: [PATCH 4/8] :hammer: FMG / FMGols works in SOME circumstances.

* Works for no groups only.
* Only works with semtest.
* Probably no support for unbiased gamma.
* Probably no support for rls.
---
 R/lav_test_fmg.R | 24 +++++++++++++-----------
 lav_fmg_greier.r |  4 ++--
 2 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/R/lav_test_fmg.R b/R/lav_test_fmg.R
index f4df097a..df0585b9 100644
--- a/R/lav_test_fmg.R
+++ b/R/lav_test_fmg.R
@@ -7,28 +7,30 @@ lav_rls <- \(object) {
 }
 
 lav_test_fmg <- \(object, input) {
-
   make_chisqs <- \(chisq) {
     if (chisq == "ml") lavaan::fitmeasures(object, "chisq") else rls(object)
   }
 
   name <- input$name
   param <- input$param
-  unbiased <- input$unbiased
-  chisq <- if (input$chisq == "ml") lavaan::fitmeasures(object, "chisq") else lav_rls(object)
-  df <- lavaan::fitmeasures(object, "df")
-  ug <- lav_ugamma_no_groups(object, unbiased)
+  #unbiased <- input$unbiased
+  #chisq <- if (input$chisq == "ml") lavaan::fitmeasures(object, "chisq") else lav_rls(object)
+  chisq <- fitmeasures(object, "chisq")
+  df <- fitmeasures(object, "df")
+  ug <- lav_object_gamma(object)[[1]] # Only for one group.
+  #ug <- lav_ugamma_no_groups(object, unbiased)
+
   lambdas_list <- Re(eigen(ug)$values)[seq(df)]
 
   lambdas <- lambdas_list
   list(
-    test = lav_fmg_reconstruct_input(input),
-    pvalue = if (name == "pEBA") {
-    lav_fmg_peba(chisq, lambdas, param)
+    test = lav_fmg_reconstruct_label(input),
+    pvalue = if (name == "fmg") {
+    lav_fmg(chisq, lambdas, param)
   } else if (name == "EBA") {
     lav_fmg_eba(chisq, lambdas, param)
-  } else if (name == "pOLS") {
-    lav_fmg_pols(chisq, lambdas, param)
+  } else if (name == "fmgols") {
+    lav_fmgols(chisq, lambdas, param)
   },
   label = lav_fmg_reconstruct_label(input))
 }
@@ -64,7 +66,7 @@ lav_fmg_parse_input <- \(string) {
   })
 }
 
-lav_fmg_reconstruct_input <- \(input) paste0(input$name, "_", input$param)
+lav_fmg_reconstruct_label <- \(input) paste0(input$name, "_", input$param)
 
 lav_fmg_construct_label <- \(input) {
   if (input$name == "fmg") {
diff --git a/lav_fmg_greier.r b/lav_fmg_greier.r
index d5c9fb3d..56ffe3b0 100644
--- a/lav_fmg_greier.r
+++ b/lav_fmg_greier.r
@@ -1,8 +1,8 @@
 n = 500
 data = psych::bfi
 
-model = "A =~ A1+b*A2+A3+A4+A5;
-         C =~ C1+b*C2+C3+C4+C5
+model = "A =~ A1+A2+A3+A4+A5;
+         C =~ C1+C2+C3+C4+C5
          "
 object <- sem(model, data[1:n, ], test = "sb")
 sum(diag(ugamma_non_nested(object))) / 34

From d6e745ae168a97919e8b6674f775032653e35d40 Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Thu, 15 Aug 2024 10:27:28 +0200
Subject: [PATCH 5/8] Update lav_test_fmg.R

---
 R/lav_test_fmg.R | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/R/lav_test_fmg.R b/R/lav_test_fmg.R
index df0585b9..972bebba 100644
--- a/R/lav_test_fmg.R
+++ b/R/lav_test_fmg.R
@@ -17,12 +17,13 @@ lav_test_fmg <- \(object, input) {
   #chisq <- if (input$chisq == "ml") lavaan::fitmeasures(object, "chisq") else lav_rls(object)
   chisq <- fitmeasures(object, "chisq")
   df <- fitmeasures(object, "df")
-  ug <- lav_object_gamma(object)[[1]] # Only for one group.
+  ug <- lav_object_inspect_UGamma(object) # Only for one group.
   #ug <- lav_ugamma_no_groups(object, unbiased)
+  lambdas <- Re(eigen(ug)$values)[seq(df)]
 
-  lambdas_list <- Re(eigen(ug)$values)[seq(df)]
-
-  lambdas <- lambdas_list
+  print(chisq)
+  print(lambdas)
+  print(param)
   list(
     test = lav_fmg_reconstruct_label(input),
     pvalue = if (name == "fmg") {

From 166a4a7103f49e0736ba0f4ec2e5c426d656a47b Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Tue, 20 Aug 2024 07:36:47 +0200
Subject: [PATCH 6/8] Imhof upper = 1000 when maximum subdivision reached.

---
 R/lav_imhof.R    | 17 -----------------
 R/lav_test_fmg.R |  7 ++++++-
 2 files changed, 6 insertions(+), 18 deletions(-)
 delete mode 100644 R/lav_imhof.R

diff --git a/R/lav_imhof.R b/R/lav_imhof.R
deleted file mode 100644
index 11b1284b..00000000
--- a/R/lav_imhof.R
+++ /dev/null
@@ -1,17 +0,0 @@
-imhof <- \(x, lambda) {
-  theta <- \(u, x, lambda) 0.5 * (sum(atan(lambda * u)) - x * u)
-  rho <- \(u, lambda) exp(1 / 4 * sum(log(1 + lambda^2 * u^2)))
-  integrand <- Vectorize(\(u) {
-    sin(theta(u, x, lambda)) / (u * rho(u, lambda))
-  })
-  z <- integrate(integrand, lower = 0, upper = Inf)$value
-  0.5 + z / pi
-}
-
-# k <- 50
-# lambda <- runif(k)
-# x <- k/2
-# imhof(x, lambda)
-# CompQuadForm::imhof(x, lambda)
-
-# microbenchmark::microbenchmark(imhof(x, lambda), CompQuadForm::imhof(x, lambda))
diff --git a/R/lav_test_fmg.R b/R/lav_test_fmg.R
index 972bebba..722b4b54 100644
--- a/R/lav_test_fmg.R
+++ b/R/lav_test_fmg.R
@@ -103,7 +103,12 @@ lav_fmg_imhof <- \(x, lambda) {
   integrand <- Vectorize(\(u) {
     sin(theta(u, x, lambda)) / (u * rho(u, lambda))
   })
-  z <- integrate(integrand, lower = 0, upper = Inf)$value
+  z <- tryCatch({
+    integrate(integrand, lower = 0, upper = Inf)$value
+  },
+  error = \(e) {
+    integrate(integrand, lower = 0, upper = 1000)$value
+  })
   0.5 + z / pi
 }
 

From 754b74f6c011d495a1ff1f1d95a3331ed4bd98fb Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Thu, 22 Aug 2024 13:43:33 +0200
Subject: [PATCH 7/8] :bug: remove eba option.

---
 R/lav_test_fmg.R | 2 --
 1 file changed, 2 deletions(-)

diff --git a/R/lav_test_fmg.R b/R/lav_test_fmg.R
index 722b4b54..05edd6f8 100644
--- a/R/lav_test_fmg.R
+++ b/R/lav_test_fmg.R
@@ -28,8 +28,6 @@ lav_test_fmg <- \(object, input) {
     test = lav_fmg_reconstruct_label(input),
     pvalue = if (name == "fmg") {
     lav_fmg(chisq, lambdas, param)
-  } else if (name == "EBA") {
-    lav_fmg_eba(chisq, lambdas, param)
   } else if (name == "fmgols") {
     lav_fmgols(chisq, lambdas, param)
   },

From 8ce096b6976aa5a6a3c41429522053f49085b11a Mon Sep 17 00:00:00 2001
From: Jonas Moss <jonas.gjertsen@gmail.com>
Date: Thu, 29 Aug 2024 13:20:00 +0200
Subject: [PATCH 8/8] :sparkles: Implement multiple groups + NO TESTING

---
 R/lav_test.R     |  3 +++
 lav_fmg_greier.r | 15 +++++++++++----
 steps.txt        |  2 ++
 3 files changed, 16 insertions(+), 4 deletions(-)
 create mode 100644 steps.txt

diff --git a/R/lav_test.R b/R/lav_test.R
index 5787a486..f3a37569 100644
--- a/R/lav_test.R
+++ b/R/lav_test.R
@@ -566,6 +566,9 @@ lav_model_test <- function(lavobject = NULL,
       # which test statistic shall we scale?
       unscaled.TEST <- TEST[[1]]
       if (lavoptions$scaled.test != "standard") {
+        print(test.orig)
+        print(lavoptions$scaled.test)
+        print(TEST)
         idx <- which(test.orig == lavoptions$scaled.test)
         if (length(idx) > 0L) {
           unscaled.TEST <- TEST[[idx[1]]]
diff --git a/lav_fmg_greier.r b/lav_fmg_greier.r
index 56ffe3b0..7f8f74ec 100644
--- a/lav_fmg_greier.r
+++ b/lav_fmg_greier.r
@@ -1,14 +1,21 @@
-n = 500
+n = 50
 data = psych::bfi
 
-model = "A =~ A1+A2+A3+A4+A5;
-         C =~ C1+C2+C3+C4+C5
+model = "A =~ A1+a*A2+A3+A4+A5;
+         C =~ C1+a*C2+C3+C4+C5
          "
-object <- sem(model, data[1:n, ], test = "sb")
+
+object <- sem(model, data[1:n, ], group = "gender", test = "sb")
+lavTest(object, test = "fmg_1")
+lavTest(object, test = "sb")
+
+
 sum(diag(ugamma_non_nested(object))) / 34
 
 dim(lavaan::lavInspect(object, "delta"))
 
+
+
 model = "A =~ A1+A2+A3+A4+A5;
          C =~ C1+C2+C3+C4+C5
          "
diff --git a/steps.txt b/steps.txt
new file mode 100644
index 00000000..2b281bc2
--- /dev/null
+++ b/steps.txt
@@ -0,0 +1,2 @@
+1. [ ] Sjekk ugamma for nested.
+2. [ ] Verifiser ugamma for groups med constraints.