diff --git a/R/test_likelihoodratio.R b/R/test_likelihoodratio.R index 6deadd348..6ab2d91bf 100644 --- a/R/test_likelihoodratio.R +++ b/R/test_likelihoodratio.R @@ -189,12 +189,14 @@ test_likelihoodratio.ListNestedRegressions <- function( } else { # lmtest::lrtest() lls <- sapply(objects, insight::get_loglikelihood, REML = REML, check_response = TRUE) - chi2 <- abs(c(NA, -2 * diff(lls))) + criterion <- -2 * lls + chi2 <- abs(c(NA, diff(criterion))) p <- stats::pchisq(chi2, abs(dfs_diff), lower.tail = FALSE) out <- data.frame( df = dfs, df_diff = dfs_diff, + Criterion = criterion, Chi2 = chi2, p = p, stringsAsFactors = FALSE @@ -228,10 +230,13 @@ test_likelihoodratio.ListNestedRegressions <- function( test_likelihoodratio_ListLavaan <- function(..., objects = NULL) { insight::check_if_installed("lavaan") - # Create data frame with info about model name and class + # Create data frame with info about model name, class, and criterion names_types <- data.frame( Model = names(objects), Type = sapply(objects, function(x) class(x)[1]), + Criterion = sapply(objects, function(x) { + -2 * as.numeric(lavaan::fitMeasures(x, "logl")) + }), stringsAsFactors = FALSE ) @@ -243,18 +248,20 @@ test_likelihoodratio_ListLavaan <- function(..., objects = NULL) { # Rename columns colnames(out)[names(out) == "Df"] <- "df" colnames(out)[names(out) == "Df diff"] <- "df_diff" - colnames(out)[names(out) == "Chisq"] <- "Chi2" + colnames(out)[names(out) == "Chisq diff"] <- "Chi2" colnames(out)[startsWith(names(out), "Pr(>")] <- "p" out$Model <- row.names(out) # Bind all data out <- merge(names_types, out[c("Model", "df", "df_diff", "Chi2", "p")], by = "Model") + out <- out[c("Model", "Type", "df", "df_diff", "Criterion", "Chi2", "p")] + out <- out[order(out$df), ] + class(out) <- c("test_likelihoodratio", "see_test_likelihoodratio", "data.frame") out } - # helper ---------------------- .is_lmer_reml <- function(x) { diff --git a/R/test_performance.R b/R/test_performance.R index af4a64ffc..293e06ef7 100644 --- a/R/test_performance.R +++ b/R/test_performance.R @@ -113,19 +113,18 @@ #' the Wald test for small sample sizes (under or about 30) or if the #' parameters are large. #' +#' The test statistic is calculated by comparing the -2 * log-likelihood +#' (-2LL) values for each model. In the output table, the `Criterion` column +#' represents this -2LL value. The difference in the criterion values +#' between the nested models corresponds to the `Chi2` statistic. This +#' Chi-square value is then used to compute the p-value based on the +#' difference in degrees of freedom (`df_diff`). +#' #' Note: for regression models, this is similar to #' `anova(..., test="LRT")` (on models) or `lmtest::lrtest(...)`, depending #' on the `estimator` argument. For **lavaan** models (SEM, CFA), the function #' calls `lavaan::lavTestLRT()`. #' -#' For models with transformed response variables (like `log(x)` or `sqrt(x)`), -#' `logLik()` returns a wrong log-likelihood. However, `test_likelihoodratio()` -#' calls `insight::get_loglikelihood()` with `check_response=TRUE`, which -#' returns a corrected log-likelihood value for models with transformed -#' response variables. Furthermore, since the LRT only accepts nested -#' models (i.e. models that differ in their fixed effects), the computed -#' log-likelihood is always based on the ML estimator, not on the REML fits. -#' #' - **Vuong's Test** - `test_vuong()`: Vuong's (1989) test can #' be used both for nested and non-nested models, and actually consists of two #' tests. diff --git a/man/test_performance.Rd b/man/test_performance.Rd index 6b2953ed6..2cc46dae4 100644 --- a/man/test_performance.Rd +++ b/man/test_performance.Rd @@ -150,18 +150,17 @@ efficient. Agresti (1990) suggests that you should use the LRT instead of the Wald test for small sample sizes (under or about 30) or if the parameters are large. +The test statistic is calculated by comparing the -2 * log-likelihood +(-2LL) values for each model. In the output table, the \code{Criterion} column +represents this -2LL value. The difference in the criterion values +between the nested models corresponds to the \code{Chi2} statistic. This +Chi-square value is then used to compute the p-value based on the +difference in degrees of freedom (\code{df_diff}). + Note: for regression models, this is similar to \code{anova(..., test="LRT")} (on models) or \code{lmtest::lrtest(...)}, depending on the \code{estimator} argument. For \strong{lavaan} models (SEM, CFA), the function calls \code{lavaan::lavTestLRT()}. - -For models with transformed response variables (like \code{log(x)} or \code{sqrt(x)}), -\code{logLik()} returns a wrong log-likelihood. However, \code{test_likelihoodratio()} -calls \code{insight::get_loglikelihood()} with \code{check_response=TRUE}, which -returns a corrected log-likelihood value for models with transformed -response variables. Furthermore, since the LRT only accepts nested -models (i.e. models that differ in their fixed effects), the computed -log-likelihood is always based on the ML estimator, not on the REML fits. \item \strong{Vuong's Test} - \code{test_vuong()}: Vuong's (1989) test can be used both for nested and non-nested models, and actually consists of two tests. diff --git a/tests/testthat/_snaps/model_performance.psych.md b/tests/testthat/_snaps/model_performance.psych.md index a557b2a12..ecba6645f 100644 --- a/tests/testthat/_snaps/model_performance.psych.md +++ b/tests/testthat/_snaps/model_performance.psych.md @@ -102,7 +102,7 @@ Chi2(1) | p (Chi2) | RMSR -------------------------- - 4.119 | 0.042 | 0.035 + 4.119 | 0.042 | 0.025 --- @@ -113,8 +113,8 @@ Model | Chi2 | df | p (Chi2) | RMSR | RMSR_corrected | TLI | RMSEA | RMSEA 90% CI | BIC | R2 | Correlation ----------------------------------------------------------------------------------------------------------------------------------- - 3-factor solution | 31.796 | 25 | 0.164 | 0.015 | 0.023 | | 0.087 | [0.000, 0.181] | -54.8 | | - g-model | 264.781 | 44 | < .001 | 0.393 | 0.440 | 0.195 | 0.395 | [0.356, 0.450] | 112.3 | 0.761 | 0.873 + 3-factor solution | 31.796 | 25 | 0.164 | 0.011 | 0.016 | | 0.087 | [0.000, 0.181] | -54.8 | | + g-model | 264.781 | 44 | < .001 | 0.278 | 0.311 | 0.195 | 0.395 | [0.356, 0.450] | 112.3 | 0.886 | 0.941 Compare the model fit of the 3-factor solution with the g-only model. If the g-model has smaller RMSR and RMSEA then your items are more @@ -131,8 +131,8 @@ Model | Chi2 | df | p (Chi2) | RMSR | RMSR_corrected | TLI | RMSEA | RMSEA 90% CI | BIC | R2 | Correlation ----------------------------------------------------------------------------------------------------------------------------------- - 3-factor solution | 31.796 | 25 | 0.164 | 0.015 | 0.023 | | 0.087 | [0.000, 0.181] | -54.8 | | - g-model | 264.781 | 44 | < .001 | 0.393 | 0.440 | 0.195 | 0.395 | [0.356, 0.450] | 112.3 | 0.761 | 0.873 + 3-factor solution | 31.796 | 25 | 0.164 | 0.011 | 0.016 | | 0.087 | [0.000, 0.181] | -54.8 | | + g-model | 264.781 | 44 | < .001 | 0.278 | 0.311 | 0.195 | 0.395 | [0.356, 0.450] | 112.3 | 0.886 | 0.941 Compare the model fit of the 3-factor solution with the g-only model. If the g-model has smaller RMSR and RMSEA then your items are more diff --git a/tests/testthat/_snaps/nestedLogit.md b/tests/testthat/_snaps/nestedLogit.md deleted file mode 100644 index 74188130e..000000000 --- a/tests/testthat/_snaps/nestedLogit.md +++ /dev/null @@ -1,12 +0,0 @@ -# model_performance - - Code - model_performance(mnl) - Output - # Indices of model performance - - Response | AIC | BIC | RMSE | Sigma | R2 - ------------------------------------------------ - work | 325.7 | 336.4 | 0.456 | 1 | 0.138 - full | 110.5 | 118.5 | 0.398 | 1 | 0.333 - diff --git a/tests/testthat/test-test_likelihoodratio.R b/tests/testthat/test-test_likelihoodratio.R index 6b1652193..59ac49f42 100644 --- a/tests/testthat/test-test_likelihoodratio.R +++ b/tests/testthat/test-test_likelihoodratio.R @@ -125,3 +125,74 @@ test_that("test_likelihoodratio - print p-digits", { expect_snapshot(test_likelihoodratio(m1, m2)) expect_snapshot(print_md(test_likelihoodratio(m1, m2), p_digits = 3)) }) + +test_that("test_likelihoodratio - Criterion values (lm)", { + m1 <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars) + m2 <- lm(mpg ~ wt + cyl + gear, data = mtcars) + + # Check ML estimator since OLS defaults to .test_wald() + rez <- test_likelihoodratio(m2, m1, estimator = "ML") + + # Check if column exists + expect_true("Criterion" %in% colnames(rez)) + + # Check exact values against manual -2LL calculation + ll1 <- as.numeric(stats::logLik(m1)) + ll2 <- as.numeric(stats::logLik(m2)) + expect_equal(rez$Criterion, -2 * c(ll2, ll1), tolerance = 1e-3) + + # Check that the difference in Criterion matches Chi2 exactly + expect_equal(rez$Chi2[2], abs(diff(rez$Criterion)), tolerance = 1e-3) +}) + +test_that("test_likelihoodratio - Criterion values (lme4)", { + skip_if_not_installed("lme4") + + m1 <- suppressMessages(lme4::lmer( + Sepal.Length ~ Petal.Width + (1 | Species), + data = iris, + REML = FALSE + )) + m2 <- suppressMessages(lme4::lmer( + Sepal.Length ~ Petal.Width + Petal.Length + (1 | Species), + data = iris, + REML = FALSE + )) + + rez <- test_likelihoodratio(m1, m2, estimator = "ML") + + expect_true("Criterion" %in% colnames(rez)) + + # Check values + ll1 <- as.numeric(stats::logLik(m1)) + ll2 <- as.numeric(stats::logLik(m2)) + expect_equal(rez$Criterion, -2 * c(ll1, ll2), tolerance = 1e-3) + + # Check math + expect_equal(rez$Chi2[2], abs(diff(rez$Criterion)), tolerance = 1e-3) +}) + +test_that("test_likelihoodratio - Criterion values (lavaan)", { + skip_if_not_installed("lavaan") + + structure1 <- " visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + visual ~~ textual + speed " + m1 <- suppressWarnings(lavaan::cfa(structure1, data = lavaan::HolzingerSwineford1939)) + + structure2 <- " visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + visual ~~ 0 * textual + speed " + m2 <- suppressWarnings(lavaan::cfa(structure2, data = lavaan::HolzingerSwineford1939)) + + rez <- test_likelihoodratio(m1, m2) + + expect_true("Criterion" %in% colnames(rez)) + + ll1 <- as.numeric(lavaan::fitMeasures(m1, "logl")) + ll2 <- as.numeric(lavaan::fitMeasures(m2, "logl")) + expect_equal(rez$Criterion, -2 * c(ll1, ll2), tolerance = 1e-3) + expect_equal(rez$Chi2[2], abs(rez$Criterion[2] - rez$Criterion[1]), tolerance = 1e-3) +})