Skip to content

Revise compare_models() for Bayesian models #780

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
strengejacke wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Revise compare_models() for Bayesian models #780

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6752208
Revise compare_models() for Bayesian models
strengejacke Nov 24, 2024
4be74cf
Merge branch 'main' into strengejacke/issue716
strengejacke Dec 2, 2024
a1c129d
Merge branch 'main' into strengejacke/issue716
strengejacke Dec 29, 2024
f856b5f
Merge branch 'main' into strengejacke/issue716
strengejacke Jan 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 29 additions & 22 deletions R/compare_performance.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,10 @@ compare_performance <- function(..., metrics = "all", rank = FALSE, estimator =
model_objects <- insight::ellipsis_info(..., only_models = TRUE)

# ensure proper object names
model_objects <- .check_objectnames(model_objects, sapply(match.call(expand.dots = FALSE)[["..."]], as.character))
model_objects <- .check_objectnames(
model_objects,
sapply(match.call(expand.dots = FALSE)[["..."]], as.character)
)

# drop unsupport models
supported_models <- sapply(model_objects, function(i) insight::is_model_supported(i) | inherits(i, "lavaan"))
Expand Down Expand Up @@ -119,14 +122,14 @@ compare_performance <- function(..., metrics = "all", rank = FALSE, estimator =
})
dfs <- Reduce(function(x, y) merge(x, y, all = TRUE, sort = FALSE), m)

if (any(c("AIC", "AICc", "BIC", "WAIC") %in% names(dfs))) {
if (any(c("AIC", "AICc", "BIC", "WAIC") %in% colnames(dfs))) {
dfs$AIC_wt <- .ic_weight(dfs[["AIC"]])
dfs$AICc_wt <- .ic_weight(dfs[["AICc"]])
dfs$BIC_wt <- .ic_weight(dfs[["BIC"]])
dfs$WAIC_wt <- .ic_weight(dfs[["WAIC"]])
}

if ("LOOIC" %in% names(dfs)) {
if ("LOOIC" %in% colnames(dfs)) {
lpd_point <- do.call(cbind, lapply(attri, function(x) x$loo$pointwise[, "elpd_loo"]))
dfs$LOOIC_wt <- as.numeric(loo::stacking_weights(lpd_point))
}
Expand All @@ -144,49 +147,53 @@ compare_performance <- function(..., metrics = "all", rank = FALSE, estimator =
}

# Reorder columns
if (all(c("BIC", "BF") %in% names(dfs))) {
idx1 <- grep("^BIC$", names(dfs))
idx2 <- grep("BF", names(dfs), fixed = TRUE)
if (all(c("BIC", "BF") %in% colnames(dfs))) {
idx1 <- grep("^BIC$", colnames(dfs))
idx2 <- grep("BF", colnames(dfs), fixed = TRUE)
last_part <- (idx1 + 1):ncol(dfs)
dfs <- dfs[, c(1:idx1, idx2, last_part[last_part != idx2])]
}
if (all(c("AIC", "AIC_wt") %in% names(dfs))) {
idx1 <- grep("^AIC$", names(dfs))
idx2 <- grep("AIC_wt", names(dfs), fixed = TRUE)
if (all(c("AIC", "AIC_wt") %in% colnames(dfs))) {
idx1 <- grep("^AIC$", colnames(dfs))
idx2 <- grep("AIC_wt", colnames(dfs), fixed = TRUE)
last_part <- (idx1 + 1):ncol(dfs)
dfs <- dfs[, c(1:idx1, idx2, last_part[last_part != idx2])]
}
if (all(c("BIC", "BIC_wt") %in% names(dfs))) {
idx1 <- grep("^BIC$", names(dfs))
idx2 <- grep("BIC_wt", names(dfs), fixed = TRUE)
if (all(c("BIC", "BIC_wt") %in% colnames(dfs))) {
idx1 <- grep("^BIC$", colnames(dfs))
idx2 <- grep("BIC_wt", colnames(dfs), fixed = TRUE)
last_part <- (idx1 + 1):ncol(dfs)
dfs <- dfs[, c(1:idx1, idx2, last_part[last_part != idx2])]
}
if (all(c("AICc", "AICc_wt") %in% names(dfs))) {
idx1 <- grep("^AICc$", names(dfs))
idx2 <- grep("AICc_wt", names(dfs), fixed = TRUE)
if (all(c("AICc", "AICc_wt") %in% colnames(dfs))) {
idx1 <- grep("^AICc$", colnames(dfs))
idx2 <- grep("AICc_wt", colnames(dfs), fixed = TRUE)
last_part <- (idx1 + 1):ncol(dfs)
dfs <- dfs[, c(1:idx1, idx2, last_part[last_part != idx2])]
}
if (all(c("WAIC", "WAIC_wt") %in% names(dfs))) {
idx1 <- grep("^WAIC$", names(dfs))
idx2 <- grep("WAIC_wt", names(dfs), fixed = TRUE)
if (all(c("WAIC", "WAIC_wt") %in% colnames(dfs))) {
idx1 <- grep("^WAIC$", colnames(dfs))
idx2 <- grep("WAIC_wt", colnames(dfs), fixed = TRUE)
last_part <- (idx1 + 1):ncol(dfs)
dfs <- dfs[, c(1:idx1, idx2, last_part[last_part != idx2])]
}
if (all(c("LOOIC", "LOOIC_wt") %in% names(dfs))) {
idx1 <- grep("^LOOIC$", names(dfs))
idx2 <- grep("LOOIC_wt", names(dfs), fixed = TRUE)
if (all(c("LOOIC", "LOOIC_wt") %in% colnames(dfs))) {
idx1 <- grep("^LOOIC$", colnames(dfs))
idx2 <- grep("LOOIC_wt", colnames(dfs), fixed = TRUE)
last_part <- (idx1 + 1):ncol(dfs)
dfs <- dfs[, c(1:idx1, idx2, last_part[last_part != idx2])]
}
if (any(startsWith(colnames(dfs), "R2"))) {
idx <- which(startsWith(colnames(dfs), "R2"))
dfs <- dfs[, c(1, 2, idx, setdiff(seq_len(ncol(dfs)), idx))]
}

# for REML fits, warn user
if (isTRUE(verbose) &&
# only warn for REML fit
identical(estimator, "REML") &&
# only for IC comparison
any(grepl("(AIC|BIC)", names(dfs))) &&
any(grepl("(AIC|BIC)", colnames(dfs))) &&
# only when mixed models are involved, others probably don't have problems with REML fit
any(sapply(model_objects, insight::is_mixed_model)) &&
# only if not all models have same fixed effects (else, REML is ok)
Expand Down