From 1b5756c82fe4266c7a45fd8aa536064c6f3d8c7d Mon Sep 17 00:00:00 2001 From: jmgirard Date: Wed, 26 Mar 2025 12:58:10 -0500 Subject: [PATCH 1/3] Add Deviances to LRT output --- R/test_likelihoodratio.R | 7 +- tests/testthat/_snaps/check_dag.md | 275 ------------------ tests/testthat/_snaps/mclogit.md | 22 -- .../testthat/_snaps/model_performance.rma.md | 7 - tests/testthat/_snaps/nestedLogit.md | 12 - 5 files changed, 6 insertions(+), 317 deletions(-) delete mode 100644 tests/testthat/_snaps/check_dag.md delete mode 100644 tests/testthat/_snaps/mclogit.md delete mode 100644 tests/testthat/_snaps/model_performance.rma.md delete mode 100644 tests/testthat/_snaps/nestedLogit.md diff --git a/R/test_likelihoodratio.R b/R/test_likelihoodratio.R index 20a58f6fb..8fd65682e 100644 --- a/R/test_likelihoodratio.R +++ b/R/test_likelihoodratio.R @@ -71,6 +71,9 @@ print.test_likelihoodratio <- function(x, digits = 2, ...) { if ("LogLik" %in% names(x)) { best <- which.max(x$LogLik) footer <- c(sprintf("\nModel '%s' seems to have the best model fit.\n", x$Model[best]), "yellow") + } else if ("Dev" %in% names(x)) { + best <- which.min(x$Dev) + footer <- c(sprintf("\nModel '%s' seems to have the best model fit.\n", x$Name[best]), "yellow") } else { footer <- NULL } @@ -118,10 +121,12 @@ test_likelihoodratio.ListNestedRegressions <- function(objects, estimator = "ML" } else { # lmtest::lrtest() lls <- sapply(objects, insight::get_loglikelihood, REML = REML, check_response = TRUE) - chi2 <- abs(c(NA, -2 * diff(lls))) + devs <- -2 * lls + chi2 <- abs(c(NA, diff(devs))) p <- stats::pchisq(chi2, abs(dfs_diff), lower.tail = FALSE) out <- data.frame( + Dev = devs, df = dfs, df_diff = dfs_diff, Chi2 = chi2, diff --git a/tests/testthat/_snaps/check_dag.md b/tests/testthat/_snaps/check_dag.md deleted file mode 100644 index 779ec5a71..000000000 --- a/tests/testthat/_snaps/check_dag.md +++ /dev/null @@ -1,275 +0,0 @@ -# check_dag - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: y - - Exposure: x - - Identification of direct and total effects - - Model is correctly specified. - No adjustment needed to estimate the direct and total effect of `x` on `y`. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: y - - Exposure: x - - Adjustment: b - - Identification of direct and total effects - - Model is correctly specified. - All minimal sufficient adjustments to estimate the direct and total effect were done. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: y - - Exposure: x - - Identification of direct and total effects - - Incorrectly adjusted! - To estimate the direct and total effect, at least adjust for `b`. Currently, the model does not adjust for any variables. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: y - - Exposure: x - - Adjustment: c - - Identification of direct and total effects - - Incorrectly adjusted! - To estimate the direct and total effect, at least adjust for `b` and `c`. Currently, the model only adjusts for `c`. You possibly also need to adjust for `b` to block biasing paths. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: y - - Exposure: x - - Adjustment: c - - Identification of direct and total effects - - Incorrectly adjusted! - To estimate the direct and total effect, at least adjust for `b` and `c`. Currently, the model only adjusts for `c`. You possibly also need to adjust for `b` to block biasing paths. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: mpg - - Exposure: wt - - Adjustments: cyl, disp and gear - - Identification of direct and total effects - - Model is correctly specified. - All minimal sufficient adjustments to estimate the direct and total effect were done. - - -# check_dag, multiple adjustment sets - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: exam - - Exposure: podcast - - Identification of direct and total effects - - Incorrectly adjusted! - To estimate the direct and total effect, at least adjust for one of the following sets: - - alertness, prepared - - alertness, skills_course - - mood, prepared - - mood, skills_course. - Currently, the model does not adjust for any variables. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: exam - - Exposure: podcast - - Adjustments: alertness and prepared - - Identification of direct and total effects - - Model is correctly specified. - All minimal sufficient adjustments to estimate the direct and total effect were done. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: exam - - Exposure: podcast - - Adjustment: alertness - - Identification of direct and total effects - - Incorrectly adjusted! - To estimate the direct and total effect, at least adjust for one of the following sets: - - alertness, prepared - - alertness, skills_course. - Currently, the model only adjusts for `alertness`. You possibly also need to adjust for `prepared` and `skills_course` to block biasing paths. - - -# check_dag, different adjustements for total and direct - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: outcome - - Exposure: exposure - - Identification of direct effects - - Incorrectly adjusted! - To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model does not adjust for any variables. - - Identification of total effects - - Incorrectly adjusted! - To estimate the total effect, at least adjust for `x1`. Currently, the model does not adjust for any variables. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: outcome - - Exposure: exposure - - Adjustment: x1 - - Identification of direct effects - - Incorrectly adjusted! - To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model only adjusts for `x1`. You possibly also need to adjust for `x2` to block biasing paths. - - Identification of total effects - - Model is correctly specified. - All minimal sufficient adjustments to estimate the total effect were done. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: outcome - - Exposure: exposure - - Adjustment: x2 - - Identification of direct effects - - Incorrectly adjusted! - To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model only adjusts for `x2`. You possibly also need to adjust for `x1` to block biasing paths. - - Identification of total effects - - Incorrectly adjusted! - To estimate the total effect, do not adjust for `x2`. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: outcome - - Exposure: exposure - - Adjustments: x1 and x2 - - Identification of direct effects - - Model is correctly specified. - All minimal sufficient adjustments to estimate the direct effect were done. - - Identification of total effects - - Incorrectly adjusted! - To estimate the total effect, do not adjust for some or all of `x1` and `x2`. - - -# check_dag, collider bias - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: SMD_ICD11 - - Exposure: agegroup - - Adjustments: edgroup3, gender_kid, pss4_kid_sum_2sd and residence - - Identification of direct effects - - Incorrectly adjusted! - To estimate the direct effect, at least adjust for `edgroup3`, `gender_kid`, `pss4_kid_sum_2sd`, `residence` and `sm_h_total_kid`. Currently, the model only adjusts for `edgroup3`, `gender_kid`, `pss4_kid_sum_2sd` and `residence`. You possibly also need to adjust for `sm_h_total_kid` to block biasing paths. - - Identification of total effects - - Model is correctly specified. - All minimal sufficient adjustments to estimate the total effect were done. - - ---- - - Code - print(dag) - Output - # Check for correct adjustment sets - - Outcome: SMD_ICD11 - - Exposure: agegroup - - Adjustments: edgroup3, gender_kid, pss4_kid_sum_2sd, residence and sm_h_total_kid - - Collider: sm_h_total_kid - - Identification of direct effects - - Incorrectly adjusted! - Your model adjusts for a potential collider. To estimate the direct effect, do not adjust for `sm_h_total_kid` to avoid collider-bias. It is recommended to double-check for the collider-bias on the dagitty-website. - - Identification of total effects - - Incorrectly adjusted! - Your model adjusts for a potential collider. To estimate the total effect, do not adjust for `sm_h_total_kid` to avoid collider-bias. It is recommended to double-check for the collider-bias on the dagitty-website. - - diff --git a/tests/testthat/_snaps/mclogit.md b/tests/testthat/_snaps/mclogit.md deleted file mode 100644 index 6a69f5db3..000000000 --- a/tests/testthat/_snaps/mclogit.md +++ /dev/null @@ -1,22 +0,0 @@ -# model_performance - - Code - model_performance(mod_mb) - Output - # Indices of model performance - - AIC | BIC | Nagelkerke's R2 | RMSE | Sigma - ------------------------------------------------- - 38.823 | 47.618 | 0.836 | 0.298 | 1.016 - ---- - - Code - model_performance(mod_mc) - Output - # Indices of model performance - - AIC | BIC | Nagelkerke's R2 | RMSE | Sigma - ------------------------------------------------- - 13.228 | 24.424 | 0.998 | 0.009 | 0.068 - diff --git a/tests/testthat/_snaps/model_performance.rma.md b/tests/testthat/_snaps/model_performance.rma.md deleted file mode 100644 index f7087b31b..000000000 --- a/tests/testthat/_snaps/model_performance.rma.md +++ /dev/null @@ -1,7 +0,0 @@ -# check_outliers.rma - - Code - out - Output - 2 outliers detected: studies 4 (Hart & Sutherland) and 8 (TPT Madras). - diff --git a/tests/testthat/_snaps/nestedLogit.md b/tests/testthat/_snaps/nestedLogit.md deleted file mode 100644 index f5c9a0bdd..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.733 | 336.449 | 0.456 | 1.000 | 0.138 - full | 110.495 | 118.541 | 0.398 | 1.000 | 0.333 - From 9b402782d250755973dc76c4b6ec8e55d184d86b Mon Sep 17 00:00:00 2001 From: jmgirard Date: Wed, 26 Mar 2025 14:18:35 -0500 Subject: [PATCH 2/3] Change Dev to Deviance --- R/test_likelihoodratio.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/test_likelihoodratio.R b/R/test_likelihoodratio.R index 8fd65682e..403c8dea2 100644 --- a/R/test_likelihoodratio.R +++ b/R/test_likelihoodratio.R @@ -126,7 +126,7 @@ test_likelihoodratio.ListNestedRegressions <- function(objects, estimator = "ML" p <- stats::pchisq(chi2, abs(dfs_diff), lower.tail = FALSE) out <- data.frame( - Dev = devs, + Deviance = devs, df = dfs, df_diff = dfs_diff, Chi2 = chi2, From 3cda2dc39b4c605d5058b56bc21ceb384b9d4490 Mon Sep 17 00:00:00 2001 From: jmgirard Date: Wed, 26 Mar 2025 21:42:49 -0500 Subject: [PATCH 3/3] add tests --- R/test_likelihoodratio.R | 2 +- tests/testthat/_snaps/check_dag.md | 275 ++++++++++++++++++ tests/testthat/_snaps/mclogit.md | 22 ++ .../testthat/_snaps/model_performance.rma.md | 7 + tests/testthat/test-test_likelihoodratio.R | 3 + 5 files changed, 308 insertions(+), 1 deletion(-) create mode 100644 tests/testthat/_snaps/check_dag.md create mode 100644 tests/testthat/_snaps/mclogit.md create mode 100644 tests/testthat/_snaps/model_performance.rma.md diff --git a/R/test_likelihoodratio.R b/R/test_likelihoodratio.R index 403c8dea2..06f6b1cc3 100644 --- a/R/test_likelihoodratio.R +++ b/R/test_likelihoodratio.R @@ -126,9 +126,9 @@ test_likelihoodratio.ListNestedRegressions <- function(objects, estimator = "ML" p <- stats::pchisq(chi2, abs(dfs_diff), lower.tail = FALSE) out <- data.frame( - Deviance = devs, df = dfs, df_diff = dfs_diff, + deviance = devs, Chi2 = chi2, p = p, stringsAsFactors = FALSE diff --git a/tests/testthat/_snaps/check_dag.md b/tests/testthat/_snaps/check_dag.md new file mode 100644 index 000000000..779ec5a71 --- /dev/null +++ b/tests/testthat/_snaps/check_dag.md @@ -0,0 +1,275 @@ +# check_dag + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + + Identification of direct and total effects + + Model is correctly specified. + No adjustment needed to estimate the direct and total effect of `x` on `y`. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + - Adjustment: b + + Identification of direct and total effects + + Model is correctly specified. + All minimal sufficient adjustments to estimate the direct and total effect were done. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, at least adjust for `b`. Currently, the model does not adjust for any variables. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + - Adjustment: c + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, at least adjust for `b` and `c`. Currently, the model only adjusts for `c`. You possibly also need to adjust for `b` to block biasing paths. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + - Adjustment: c + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, at least adjust for `b` and `c`. Currently, the model only adjusts for `c`. You possibly also need to adjust for `b` to block biasing paths. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: mpg + - Exposure: wt + - Adjustments: cyl, disp and gear + + Identification of direct and total effects + + Model is correctly specified. + All minimal sufficient adjustments to estimate the direct and total effect were done. + + +# check_dag, multiple adjustment sets + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: exam + - Exposure: podcast + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, at least adjust for one of the following sets: + - alertness, prepared + - alertness, skills_course + - mood, prepared + - mood, skills_course. + Currently, the model does not adjust for any variables. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: exam + - Exposure: podcast + - Adjustments: alertness and prepared + + Identification of direct and total effects + + Model is correctly specified. + All minimal sufficient adjustments to estimate the direct and total effect were done. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: exam + - Exposure: podcast + - Adjustment: alertness + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, at least adjust for one of the following sets: + - alertness, prepared + - alertness, skills_course. + Currently, the model only adjusts for `alertness`. You possibly also need to adjust for `prepared` and `skills_course` to block biasing paths. + + +# check_dag, different adjustements for total and direct + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: outcome + - Exposure: exposure + + Identification of direct effects + + Incorrectly adjusted! + To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model does not adjust for any variables. + + Identification of total effects + + Incorrectly adjusted! + To estimate the total effect, at least adjust for `x1`. Currently, the model does not adjust for any variables. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: outcome + - Exposure: exposure + - Adjustment: x1 + + Identification of direct effects + + Incorrectly adjusted! + To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model only adjusts for `x1`. You possibly also need to adjust for `x2` to block biasing paths. + + Identification of total effects + + Model is correctly specified. + All minimal sufficient adjustments to estimate the total effect were done. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: outcome + - Exposure: exposure + - Adjustment: x2 + + Identification of direct effects + + Incorrectly adjusted! + To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model only adjusts for `x2`. You possibly also need to adjust for `x1` to block biasing paths. + + Identification of total effects + + Incorrectly adjusted! + To estimate the total effect, do not adjust for `x2`. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: outcome + - Exposure: exposure + - Adjustments: x1 and x2 + + Identification of direct effects + + Model is correctly specified. + All minimal sufficient adjustments to estimate the direct effect were done. + + Identification of total effects + + Incorrectly adjusted! + To estimate the total effect, do not adjust for some or all of `x1` and `x2`. + + +# check_dag, collider bias + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: SMD_ICD11 + - Exposure: agegroup + - Adjustments: edgroup3, gender_kid, pss4_kid_sum_2sd and residence + + Identification of direct effects + + Incorrectly adjusted! + To estimate the direct effect, at least adjust for `edgroup3`, `gender_kid`, `pss4_kid_sum_2sd`, `residence` and `sm_h_total_kid`. Currently, the model only adjusts for `edgroup3`, `gender_kid`, `pss4_kid_sum_2sd` and `residence`. You possibly also need to adjust for `sm_h_total_kid` to block biasing paths. + + Identification of total effects + + Model is correctly specified. + All minimal sufficient adjustments to estimate the total effect were done. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: SMD_ICD11 + - Exposure: agegroup + - Adjustments: edgroup3, gender_kid, pss4_kid_sum_2sd, residence and sm_h_total_kid + - Collider: sm_h_total_kid + + Identification of direct effects + + Incorrectly adjusted! + Your model adjusts for a potential collider. To estimate the direct effect, do not adjust for `sm_h_total_kid` to avoid collider-bias. It is recommended to double-check for the collider-bias on the dagitty-website. + + Identification of total effects + + Incorrectly adjusted! + Your model adjusts for a potential collider. To estimate the total effect, do not adjust for `sm_h_total_kid` to avoid collider-bias. It is recommended to double-check for the collider-bias on the dagitty-website. + + diff --git a/tests/testthat/_snaps/mclogit.md b/tests/testthat/_snaps/mclogit.md new file mode 100644 index 000000000..6a69f5db3 --- /dev/null +++ b/tests/testthat/_snaps/mclogit.md @@ -0,0 +1,22 @@ +# model_performance + + Code + model_performance(mod_mb) + Output + # Indices of model performance + + AIC | BIC | Nagelkerke's R2 | RMSE | Sigma + ------------------------------------------------- + 38.823 | 47.618 | 0.836 | 0.298 | 1.016 + +--- + + Code + model_performance(mod_mc) + Output + # Indices of model performance + + AIC | BIC | Nagelkerke's R2 | RMSE | Sigma + ------------------------------------------------- + 13.228 | 24.424 | 0.998 | 0.009 | 0.068 + diff --git a/tests/testthat/_snaps/model_performance.rma.md b/tests/testthat/_snaps/model_performance.rma.md new file mode 100644 index 000000000..f7087b31b --- /dev/null +++ b/tests/testthat/_snaps/model_performance.rma.md @@ -0,0 +1,7 @@ +# check_outliers.rma + + Code + out + Output + 2 outliers detected: studies 4 (Hart & Sutherland) and 8 (TPT Madras). + diff --git a/tests/testthat/test-test_likelihoodratio.R b/tests/testthat/test-test_likelihoodratio.R index 9967ae3d4..3e249f026 100644 --- a/tests/testthat/test-test_likelihoodratio.R +++ b/tests/testthat/test-test_likelihoodratio.R @@ -63,6 +63,7 @@ test_that("test_likelihoodratio - lme4 ML", { t1 <- test_lrt(m1, m2, m3) t2 <- suppressMessages(anova(m1, m2, m3)) expect_equal(attributes(t1)$estimator, "ml") + expect_equal(t1$deviance, c(202.2215, 116.9578, 116.1164), tolerance = 1e-3) expect_equal(t1$Chi2, c(NA, 85.26365, 0.84141), tolerance = 1e-3) expect_equal(t1$p, c(NA, 0, 0.35899), tolerance = 1e-3) # close, but not the same @@ -80,6 +81,7 @@ test_that("test_likelihoodratio - lme4 OLS", { test_that("test_likelihoodratio - lme4 REML", { expect_warning(t3 <- test_lrt(m1, m2, m3, estimator = "REML")) expect_equal(attributes(t3)$estimator, "reml") + expect_equal(t3$deviance, c(210.9834, 121.6540, 124.5104), tolerance = 1e-3) expect_equal(t3$Chi2, c(NA, 89.32933, 2.85635), tolerance = 1e-3) expect_equal(t3$p, c(NA, 0, 0.09101), tolerance = 1e-3) }) @@ -91,6 +93,7 @@ m3 <- glm(am ~ mpg + hp + vs, data = mtcars, family = binomial()) test_that("test_likelihoodratio - glm", { t1 <- anova(m1, m2, m3, test = "LRT") t2 <- test_lrt(m1, m2, m3) + expect_equal(t1$`Resid. Dev`, t2$deviance, tolerance = 1e-3) expect_equal(t1$`Pr(>Chi)`, t2$p, tolerance = 1e-3) expect_equal(t1$Deviance, t2$Chi2, tolerance = 1e-3) expect_equal(attributes(t2)$estimator, "ml")