From bc3e982409e66bed86ecdfa300873f3362375065 Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Sat, 15 Jun 2024 23:20:29 +0200 Subject: [PATCH 1/9] lav_test_print() skips $stat.group when NULL --- R/lav_test_print.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/lav_test_print.R b/R/lav_test_print.R index 4e141e6c..15d2b782 100644 --- a/R/lav_test_print.R +++ b/R/lav_test_print.R @@ -271,7 +271,7 @@ lav_test_print <- function(object, nd = 3L) { # multiple groups? ngroups <- ngroups - if (ngroups > 1L) { + if (ngroups > 1L && !is.null(TEST[[block]]$stat.group)) { c1 <- c2 <- c3 <- character(ngroups) for (g in 1:ngroups) { tmp <- sprintf(" %-40s", group.label[[g]]) From c64a0cd0e5d0193a207163170d52f134bc887f02 Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Sun, 16 Jun 2024 00:43:54 +0200 Subject: [PATCH 2/9] avoid lav_test_satorra_bentler() error when user h1.model's df==0 --- R/lav_test.R | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/R/lav_test.R b/R/lav_test.R index 3a797514..106c2135 100644 --- a/R/lav_test.R +++ b/R/lav_test.R @@ -723,6 +723,8 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { testNames1 <- names(lav_obj_h1@test) testNames <- intersect(testNames0, testNames1) + copyScaled <- FALSE # in case scaled test is NA (when standard == 0) + ## loop over those tests for (tn in testNames) { lrtCall <- lrtCallTemplate @@ -739,6 +741,16 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { } else if (tn %in% c("satorra.bentler", "yuan.bentler","yuan.bentler.mplus")) { lrtCall$test <- tn + + ## is LRT even necessary? + noScaledStat <- is.na(lav_obj_h1@test[[tn]]$stat) + noScalingFactor <- is.na(lav_obj_h1@test[[tn]]$scaling.factor) + perfectFit <- isTRUE( all.equal(lav_obj_h1@test$standard$stat, 0) ) + + if (perfectFit && (noScaledStat || noScalingFactor)) { + copyScaled <- TRUE + } + } else if (grepl(pattern = "browne", x = tn)) { lrtCall$type <- tn } else { @@ -748,7 +760,17 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { } ## get new test - if (lav_obj_h0@test[[1]]$df == lav_obj_h1@test[[1]]$df) { + if (tn %in% c("satorra.bentler", "yuan.bentler", "yuan.bentler.mplus") + && copyScaled) { + ## Don't run lavTestLRT(). Keep existing H0 test stat to avoid error. + ## But adjust df + newTEST[[tn]]$df <- newTEST[[tn]]$df - lav_obj_h1@test[[tn]]$df + ## for RMSEA, reverse-engineer $trace.UGamma from $scaling.factor & new df + newTEST[[tn]]$trace.UGamma <- newTEST[[tn]]$df * newTEST[[tn]]$scaling.factor + ## skip the rest of this loop + next + + } else if (lav_obj_h0@test[[1]]$df == lav_obj_h1@test[[1]]$df) { ## suppress warning about == df ANOVA <- suppressWarnings(eval(as.call(lrtCall))) } else { @@ -769,7 +791,7 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { newTEST[[tn]]$shift.parameter <- attr(ANOVA, "shift")[2] # first row is NA } else { ## unless scaled.shifted, RMSEA is calculated from $standard$stat and - ## df == sum($trace.UGamma). Reverse-engineer from $scaling factor: + ## df == sum($trace.UGamma). Reverse-engineer from $scaling.factor: newTEST[[tn]]$trace.UGamma <- newTEST[[tn]]$df * newTEST[[tn]]$scaling.factor } ## should not be necessary to replace $trace.UGamma2 From accab8f27b0f8917ba1e2a95ac68bf4289245081 Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Mon, 17 Jun 2024 13:09:24 +0200 Subject: [PATCH 3/9] fix order of models in lav_update_test_custom_h1() --- R/lav_test.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/lav_test.R b/R/lav_test.R index 106c2135..ee778869 100644 --- a/R/lav_test.R +++ b/R/lav_test.R @@ -715,8 +715,8 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { newTEST <- lav_obj_h0@test ## assemble a call to lavTestLRT() - lrtCallTemplate <- list(quote(lavTestLRT), object = quote(lav_obj_h1), - quote(lav_obj_h0)) # in ... + lrtCallTemplate <- list(quote(lavTestLRT), object = quote(lav_obj_h0), + quote(lav_obj_h1)) # in ... ## can only update tests available in both objects testNames0 <- names(lav_obj_h0@test) From cc4286831e0adc2cb0b45a7138783e7a1b48157e Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Mon, 17 Jun 2024 13:35:51 +0200 Subject: [PATCH 4/9] replace scaled.idx by block idx to print robust test --- R/lav_test_print.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/lav_test_print.R b/R/lav_test_print.R index 15d2b782..d167c534 100644 --- a/R/lav_test_print.R +++ b/R/lav_test_print.R @@ -283,7 +283,7 @@ lav_test_print <- function(object, nd = 3L) { justify = "right" ) } else { - tmp <- sprintf(num.format, TEST[[scaled.idx]]$stat.group[g]) + tmp <- sprintf(num.format, TEST[[block]]$stat.group[g]) c2[g] <- format(tmp, width = 8L + max(0, (nd - 3L)) * 4L, justify = "right" From 936bb2308da24adbd2c834b0ee9dcc5f5ee7d787 Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Tue, 18 Jun 2024 15:51:26 +0200 Subject: [PATCH 5/9] fix SatorraBentler2001 error when chisq=0 but df>0 --- R/lav_test_diff.R | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/R/lav_test_diff.R b/R/lav_test_diff.R index a3d57ad1..acb4fcf5 100644 --- a/R/lav_test_diff.R +++ b/R/lav_test_diff.R @@ -205,8 +205,15 @@ lav_test_diff_SatorraBentler2001 <- function(m1, m0, test = 2) { T1 <- m1@test[[1]]$stat r1 <- m1@test[[1]]$df c1 <- m1@test[[test]]$scaling.factor - if (r1 == 0) { # saturated model - c1 <- 1 + + ## check for situations when scaling.factor would be NA + if (r1 == 0) { + ## saturated model + c1 <- 1 # canceled out by 0 when calculating "cd" + + } else if (r1 > 0 && isTRUE(all.equal(T1, 0))) { + ## perfect fit + c1 <- 0 # cancels out r1 when calculating "cd" } T0 <- m0@test[[1]]$stat From b3a961e34b9ec8668fa75a4a65751961028c26e5 Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Wed, 19 Jun 2024 16:53:57 +0200 Subject: [PATCH 6/9] update_test_custom_h1() sets scaled.shifted=FALSE when test="mean.var.adjusted" --- R/lav_test.R | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/R/lav_test.R b/R/lav_test.R index ee778869..2712a029 100644 --- a/R/lav_test.R +++ b/R/lav_test.R @@ -723,7 +723,7 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { testNames1 <- names(lav_obj_h1@test) testNames <- intersect(testNames0, testNames1) - copyScaled <- FALSE # in case scaled test is NA (when standard == 0) + # copyScaled <- FALSE # in case scaled test is NA (when standard == 0) ## loop over those tests for (tn in testNames) { @@ -736,6 +736,7 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { lrtCall$method <- "mean.var.adjusted.PLRT" } else { lrtCall$method <- "satorra.2000" + if (tn == "mean.var.adjusted") lrtCall$scaled.shifted <- FALSE } lrtCall$scaled.shifted <- tn == "scaled.shifted" } else if (tn %in% c("satorra.bentler", @@ -743,13 +744,13 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { lrtCall$test <- tn ## is LRT even necessary? - noScaledStat <- is.na(lav_obj_h1@test[[tn]]$stat) - noScalingFactor <- is.na(lav_obj_h1@test[[tn]]$scaling.factor) - perfectFit <- isTRUE( all.equal(lav_obj_h1@test$standard$stat, 0) ) - - if (perfectFit && (noScaledStat || noScalingFactor)) { - copyScaled <- TRUE - } + # noScaledStat <- is.na(lav_obj_h0@test[[tn]]$stat) + # noScalingFactor <- is.na(lav_obj_h0@test[[tn]]$scaling.factor) + # perfectFit <- isTRUE( all.equal(lav_obj_h0@test$standard$stat, 0) ) + # + # if (perfectFit && (noScaledStat || noScalingFactor)) { + # copyScaled <- TRUE + # } } else if (grepl(pattern = "browne", x = tn)) { lrtCall$type <- tn @@ -760,17 +761,18 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { } ## get new test - if (tn %in% c("satorra.bentler", "yuan.bentler", "yuan.bentler.mplus") - && copyScaled) { - ## Don't run lavTestLRT(). Keep existing H0 test stat to avoid error. - ## But adjust df - newTEST[[tn]]$df <- newTEST[[tn]]$df - lav_obj_h1@test[[tn]]$df - ## for RMSEA, reverse-engineer $trace.UGamma from $scaling.factor & new df - newTEST[[tn]]$trace.UGamma <- newTEST[[tn]]$df * newTEST[[tn]]$scaling.factor - ## skip the rest of this loop - next - - } else if (lav_obj_h0@test[[1]]$df == lav_obj_h1@test[[1]]$df) { + # if (tn %in% c("satorra.bentler", "yuan.bentler", "yuan.bentler.mplus") + # && copyScaled) { + # ## Don't run lavTestLRT(). Keep existing H0 test stat to avoid error. + # ## But adjust df + # newTEST[[tn]]$df <- lav_obj_h0@test[[tn]]$df - lav_obj_h1@test[[tn]]$df + # ## for RMSEA, reverse-engineer $trace.UGamma from $scaling.factor & new df + # newTEST[[tn]]$trace.UGamma <- newTEST[[tn]]$df * newTEST[[tn]]$scaling.factor + # ## skip the rest of this loop + # next + # + # } else + if (lav_obj_h0@test[[1]]$df == lav_obj_h1@test[[1]]$df) { ## suppress warning about == df ANOVA <- suppressWarnings(eval(as.call(lrtCall))) } else { From 15460ad71dcdb5bb4ca696dee0266133126509ab Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Thu, 20 Jun 2024 15:40:25 +0200 Subject: [PATCH 7/9] remove "copyScaled" code, anova() fixed now --- R/lav_test.R | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/R/lav_test.R b/R/lav_test.R index 2712a029..24fe9af1 100644 --- a/R/lav_test.R +++ b/R/lav_test.R @@ -723,8 +723,6 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { testNames1 <- names(lav_obj_h1@test) testNames <- intersect(testNames0, testNames1) - # copyScaled <- FALSE # in case scaled test is NA (when standard == 0) - ## loop over those tests for (tn in testNames) { lrtCall <- lrtCallTemplate @@ -742,16 +740,6 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { } else if (tn %in% c("satorra.bentler", "yuan.bentler","yuan.bentler.mplus")) { lrtCall$test <- tn - - ## is LRT even necessary? - # noScaledStat <- is.na(lav_obj_h0@test[[tn]]$stat) - # noScalingFactor <- is.na(lav_obj_h0@test[[tn]]$scaling.factor) - # perfectFit <- isTRUE( all.equal(lav_obj_h0@test$standard$stat, 0) ) - # - # if (perfectFit && (noScaledStat || noScalingFactor)) { - # copyScaled <- TRUE - # } - } else if (grepl(pattern = "browne", x = tn)) { lrtCall$type <- tn } else { @@ -761,17 +749,6 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { } ## get new test - # if (tn %in% c("satorra.bentler", "yuan.bentler", "yuan.bentler.mplus") - # && copyScaled) { - # ## Don't run lavTestLRT(). Keep existing H0 test stat to avoid error. - # ## But adjust df - # newTEST[[tn]]$df <- lav_obj_h0@test[[tn]]$df - lav_obj_h1@test[[tn]]$df - # ## for RMSEA, reverse-engineer $trace.UGamma from $scaling.factor & new df - # newTEST[[tn]]$trace.UGamma <- newTEST[[tn]]$df * newTEST[[tn]]$scaling.factor - # ## skip the rest of this loop - # next - # - # } else if (lav_obj_h0@test[[1]]$df == lav_obj_h1@test[[1]]$df) { ## suppress warning about == df ANOVA <- suppressWarnings(eval(as.call(lrtCall))) From 415e976474ade06af1d51a59011514b59f213a09 Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Thu, 20 Jun 2024 16:01:09 +0200 Subject: [PATCH 8/9] pass args to lav_test_lrt_single_model() for >2 tests --- R/lav_test_LRT.R | 84 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 56 insertions(+), 28 deletions(-) diff --git a/R/lav_test_LRT.R b/R/lav_test_LRT.R index 750fa568..415fa30c 100644 --- a/R/lav_test_LRT.R +++ b/R/lav_test_LRT.R @@ -69,7 +69,7 @@ lavTestLRT <- function(object, ..., method = "default", test = "default", if (type == "cf") { lav_msg_warn(gettext("`type' argument is ignored for a single model")) } - return(lav_test_lrt_single_model(object)) + return(lav_test_lrt_single_model(object, method = method, test = test, type = type)) } # list of models @@ -474,7 +474,8 @@ lavTestLRT <- function(object, ..., method = "default", test = "default", # anova table for a single model -lav_test_lrt_single_model <- function(object) { +lav_test_lrt_single_model <- function(object, method = "default", + test = "default", type = "Chisq") { estimator <- object@Options$estimator aic <- bic <- c(NA, NA) @@ -482,35 +483,62 @@ lav_test_lrt_single_model <- function(object) { aic <- c(NA, AIC(object)) bic <- c(NA, BIC(object)) } - - if (length(object@test) > 1L) { - val <- data.frame( - Df = c(0, object@test[[2L]]$df), - AIC = aic, - BIC = bic, - Chisq = c(0, object@test[[2L]]$stat), - "Chisq diff" = c(NA, object@test[[2L]]$stat), - "Df diff" = c(NA, object@test[[2L]]$df), - "Pr(>Chisq)" = c(NA, object@test[[2L]]$pvalue), - row.names = c("Saturated", "Model"), - check.names = FALSE - ) - attr(val, "heading") <- "Chi-Squared Test Statistic (scaled)\n" + + ## determine which @test element + tn <- names(object@test) + if (length(tn) == 1L) { + TEST <- 1L # only choice + + ## More than 1. Cycle through possible user specifications: + } else if (method[1] == "standard") { + TEST <- 1L + } else if (grepl(pattern = "browne", x = type) && type %in% tn) { + TEST <- type + } else if (test %in% tn) { + TEST <- test } else { - val <- data.frame( - Df = c(0, object@test[[1L]]$df), - AIC = aic, - BIC = bic, - Chisq = c(0, object@test[[1L]]$stat), - "Chisq diff" = c(NA, object@test[[1L]]$stat), - "Df diff" = c(NA, object@test[[1L]]$df), - "Pr(>Chisq)" = c(NA, object@test[[1L]]$pvalue), - row.names = c("Saturated", "Model"), - check.names = FALSE - ) - attr(val, "heading") <- "Chi-Squared Test Statistic (unscaled)\n" + ## Nothing explicitly (or validly) requested. + ## But there is > 1 test, so take the second element (old default) + TEST <- 2L } + ## anova table + val <- data.frame( + Df = c(0, object@test[[TEST]]$df), + AIC = aic, + BIC = bic, + Chisq = c(0, object@test[[TEST]]$stat), + "Chisq diff" = c(NA, object@test[[TEST]]$stat), + "Df diff" = c(NA, object@test[[TEST]]$df), + "Pr(>Chisq)" = c(NA, object@test[[TEST]]$pvalue), + row.names = c("Saturated", "Model"), + check.names = FALSE + ) + ## scale/shift attributes + if (!is.null(object@test[[TEST]]$scaling.factor)) { + attr(val, "scale") <- c(NA, object@test[[TEST]]$scaling.factor) + } + if (!is.null(object@test[[TEST]]$shift.parameter)) { + attr(val, "shift") <- c(NA, object@test[[TEST]]$shift.parameter) + } + + ## heading + if (grepl(pattern = "browne", x = TEST)) { + attr(val, "heading") <- object@test[[TEST]]$label + + } else if (TEST == 1L) { + attr(val, "heading") <- "Chi-Squared Test Statistic (unscaled)\n" + + } else { + LABEL <- object@test[[TEST]]$label + attr(val, "heading") <- paste0("Chi-Squared Test Statistic (scaled", + ifelse(TEST == "scaled.shifted", + yes = " and shifted)", no = ")"), + ifelse(is.null(LABEL), + yes = "\n", no = paste("\n ", LABEL)), + "\n") + } + class(val) <- c("anova", class(val)) val From 5aae5f695b95610271b3ece54470dea2beb45453 Mon Sep 17 00:00:00 2001 From: TDJorgensen Date: Thu, 20 Jun 2024 16:56:26 +0200 Subject: [PATCH 9/9] remove redundant scaled.shifted setting from lav_update_test_custom_h1() --- R/lav_test.R | 1 - 1 file changed, 1 deletion(-) diff --git a/R/lav_test.R b/R/lav_test.R index 24fe9af1..bf9de45d 100644 --- a/R/lav_test.R +++ b/R/lav_test.R @@ -734,7 +734,6 @@ lav_update_test_custom_h1 <- function(lav_obj_h0, lav_obj_h1) { lrtCall$method <- "mean.var.adjusted.PLRT" } else { lrtCall$method <- "satorra.2000" - if (tn == "mean.var.adjusted") lrtCall$scaled.shifted <- FALSE } lrtCall$scaled.shifted <- tn == "scaled.shifted" } else if (tn %in% c("satorra.bentler",