From 0ec79fe9a02d15a8948a0c5f991b4995afedae72 Mon Sep 17 00:00:00 2001 From: Daniel Date: Fri, 14 Jun 2024 18:48:10 +0200 Subject: [PATCH 01/21] Test get_variance for zero-inflation models --- R/compute_variances.R | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/R/compute_variances.R b/R/compute_variances.R index 95b01f78f..273dad798 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -831,7 +831,28 @@ cvsquared <- tryCatch( { if (faminfo$link_function == "tweedie") { + # Tweedie models ------------------------------------------------------ + # --------------------------------------------------------------------- dispersion_param <- .variance_family_tweedie(x, mu, sig) + + } else if (faminfo$is_zero_inflated) { + # Zero Inflated models ------------------------------------------------ + # --------------------------------------------------------------------- + dispersion_param <- switch(faminfo$family, + # (zero-inflated) poisson ---- + # ---------------------------- + poisson = .variance_family_poisson(x, mu, faminfo), + + # hurdle-poisson ---- + # ------------------- + `hurdle poisson` = , + truncated_poisson = stats::family(x)$variance(sig), + + # others ---- + # ----------- + sig + ) + } else { dispersion_param <- switch(faminfo$family, From 53b005285381c529494d1d8ff74c313ffd7eddf5 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 11:25:15 +0200 Subject: [PATCH 02/21] docs --- R/compute_variances.R | 76 ++++++++++++++++++++++++++++++++----------- R/get_variances.R | 4 ++- man/get_variance.Rd | 4 ++- 3 files changed, 63 insertions(+), 21 deletions(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index 273dad798..d064e532a 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -4,7 +4,7 @@ name_full = NULL, verbose = TRUE, tolerance = 1e-8, - model_component = "conditional", + model_component = "full", model_null = NULL, approximation = "lognormal") { ## Original code taken from GitGub-Repo of package glmmTMB @@ -22,6 +22,7 @@ # check argument approx_method <- match.arg(approximation, c("lognormal", "delta", "trigamma", "observation_level")) + # sanity checks - distribution supported? if (any(faminfo$family == "truncated_nbinom1")) { if (verbose) { format_warning(sprintf( @@ -32,6 +33,19 @@ return(NA) } + # check whether R2 should be calculated for the full model, or the + # conditional model only + if (is.null(model_component) || model_component %in% c("zi", "zero_inflated")) { + model_component <- "full" + } + + # zero-inflated model, but not conditioning on full model? + if (!identical(model_component, "full") && (faminfo$is_zero_inflated || faminfo$is_hurdle) && verbose) { + format_alert( + "Zero-inflation part of the model is not considered for variance decomposition. Use `model_component = \"full\"` to take both the conditional and the zero-inflation model into account.", # nolint + ) + } + # rename lme4 neg-binom family if (startsWith(faminfo$family, "Negative Binomial")) { faminfo$family <- "negative binomial" @@ -43,8 +57,7 @@ x, faminfo = faminfo, name_fun = name_fun, - verbose = verbose, - model_component = model_component + verbose = verbose ) # we also need necessary model information, like fixed and random effects, @@ -56,8 +69,7 @@ model_null, faminfo = faminfo, name_fun = name_fun, - verbose = verbose, - model_component = model_component + verbose = verbose ) # Test for non-zero random effects ((near) singularity) @@ -133,6 +145,7 @@ revar_null = var.random_null, approx_method = approximation, name = name_full, + model_component = model_component, verbose = verbose ) } @@ -557,6 +570,7 @@ revar_null = NULL, name, approx_method = "lognormal", + model_component = NULL, verbose = TRUE) { # get overdispersion parameter / sigma sig <- .safe(get_sigma(x)) @@ -589,6 +603,7 @@ revar_null = revar_null, approx_method = approx_method, name = name, + model_component = model_component, verbose = verbose ), .badlink(faminfo$link_function, faminfo$family, verbose = verbose) @@ -653,6 +668,7 @@ revar_null = revar_null, approx_method = approx_method, name = name, + model_component = model_component, verbose = verbose ), sqrt = 0.25 * sig, @@ -685,6 +701,7 @@ revar_null = revar_null, name = name, approx_method = approx_method, + model_component = model_component, verbose = verbose ), .badlink(faminfo$link_function, faminfo$family, verbose = verbose) @@ -702,6 +719,7 @@ revar_null = revar_null, approx_method = approx_method, name = name, + model_component = model_component, verbose = verbose ), .badlink(faminfo$link_function, faminfo$family, verbose = verbose) @@ -767,6 +785,7 @@ revar_null = NULL, name, approx_method = "lognormal", + model_component = NULL, verbose = TRUE) { check_if_installed("lme4", "to compute variances for mixed models") @@ -828,6 +847,14 @@ ) } + # ------------------------------------------------------------------ + # we have to make exception and check for tweedie-link_function here. + # Some models from tweedie family have saved this information in + # "faminfo$family", but some also in "faminfo$link_function" (with a different + # family). So we have to check for both. + # Same applies to zero-inflated models, that's why we better double-check here. + # ------------------------------------------------------------------ + cvsquared <- tryCatch( { if (faminfo$link_function == "tweedie") { @@ -835,42 +862,54 @@ # --------------------------------------------------------------------- dispersion_param <- .variance_family_tweedie(x, mu, sig) - } else if (faminfo$is_zero_inflated) { + } else if (identical(model_component, "full") && (faminfo$is_zero_inflated || faminfo$is_hurdle)) { # Zero Inflated models ------------------------------------------------ # --------------------------------------------------------------------- dispersion_param <- switch(faminfo$family, + # (zero-inflated) poisson ---- # ---------------------------- - poisson = .variance_family_poisson(x, mu, faminfo), + poisson = , + `zero-inflated poisson` = .variance_family_poisson(x, mu, faminfo), # hurdle-poisson ---- # ------------------- `hurdle poisson` = , truncated_poisson = stats::family(x)$variance(sig), + # (zero-inflated) negative binomial ---- + # -------------------------------------- + nbinom = , + nbinom1 = , + nbinom2 = , + negbinomial = , + `negative binomial` = , + `zero-inflated negative binomial` = .variance_family_nbinom(model, mu, sig, faminfo), + + # hurdle negative binomial ---- + # ----------------------------- + truncated_nbinom2 = stats::family(x)$variance(mu, sig), + # others ---- # ----------- sig ) } else { + # All other models ------------------------------------------------ + # ----------------------------------------------------------------- dispersion_param <- switch(faminfo$family, - # (zero-inflated) poisson ---- + # (generalized, compoised) poisson ---- # ---------------------------- - `zero-inflated poisson` = .variance_family_poisson(x, mu, faminfo), poisson = 1, - - # hurdle-poisson ---- - # ------------------- - `hurdle poisson` = , - truncated_poisson = stats::family(x)$variance(sig), + genpois = .variance_family_nbinom(x, mu, sig, faminfo), # Gamma, exponential ---- # ----------------------- Gamma = stats::family(x)$variance(sig), - # (zero-inflated) negative binomial ---- + # negative binomial ---- # -------------------------------------- nbinom = , nbinom1 = , @@ -878,17 +917,16 @@ quasipoisson = , negbinomial = , `negative binomial` = sig, - `zero-inflated negative binomial` = , - genpois = .variance_family_nbinom(x, mu, sig, faminfo), - truncated_nbinom2 = stats::family(x)$variance(mu, sig), # beta-alike ---- # --------------- beta = .variance_family_beta(x, mu, sig), ordbeta = .variance_family_orderedbeta(x, mu), + betabinomial = .variance_family_betabinom(x, mu, sig), + + ## TODO: check alternatives, but probably less accurate # betabinomial = .variance_family_beta(x, mu, sig), # betabinomial = stats::family(x)$variance(mu, sig), - betabinomial = .variance_family_betabinom(x, mu, sig), # other distributions ---- # ------------------------ diff --git a/R/get_variances.R b/R/get_variances.R index 93e6493b2..a5518a3b2 100644 --- a/R/get_variances.R +++ b/R/get_variances.R @@ -28,7 +28,9 @@ #' specific variance, however, it can also be `"observation_level"`. See #' _Nakagawa et al. 2017_, in particular supplement 2, for details. #' @param model_component For models that can have a zero-inflation component, -#' specify for which component variances should be returned. +#' specify for which component variances should be returned. If `NULL` or `"full"` +#' (the default), both the conditional and the zero-inflation component are taken +#' into account. If `"conditional"`, only the conditional component is considered. #' @param ... Currently not used. #' #' @return A list with following elements: diff --git a/man/get_variance.Rd b/man/get_variance.Rd index afa4c6d65..e23e00484 100644 --- a/man/get_variance.Rd +++ b/man/get_variance.Rd @@ -87,7 +87,9 @@ specific variance, however, it can also be \code{"observation_level"}. See \item{verbose}{Toggle off warnings.} \item{model_component}{For models that can have a zero-inflation component, -specify for which component variances should be returned.} +specify for which component variances should be returned. If \code{NULL} or \code{"full"} +(the default), both the conditional and the zero-inflation component are taken +into account. If \code{"conditional"}, only the conditional component is considered.} } \value{ A list with following elements: From a08e75f62708d6665087c63046fdb63d015a1fc8 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 11:41:42 +0200 Subject: [PATCH 03/21] add warning --- R/compute_variances.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index d064e532a..a33dd3c87 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -42,7 +42,7 @@ # zero-inflated model, but not conditioning on full model? if (!identical(model_component, "full") && (faminfo$is_zero_inflated || faminfo$is_hurdle) && verbose) { format_alert( - "Zero-inflation part of the model is not considered for variance decomposition. Use `model_component = \"full\"` to take both the conditional and the zero-inflation model into account.", # nolint + "Zero-inflation part of the model is not considered for variance decomposition. Use `model_component = \"full\"` to take both the conditional and the zero-inflation model into account." # nolint ) } From d572671ba7f783508ff90d00a60d55af12a7df1d Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 13:52:19 +0200 Subject: [PATCH 04/21] news --- DESCRIPTION | 2 +- NEWS.md | 14 ++++++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6b8f80dce..3f7b55c4f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: insight Title: Easy Access to Model Information for Various Model Objects -Version: 0.20.1.5 +Version: 0.20.1.6 Authors@R: c(person(given = "Daniel", family = "Lüdecke", diff --git a/NEWS.md b/NEWS.md index afdee4658..b1788dac8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,20 @@ * Improved accuracy of singularity-checks in `get_variance()`. +* `get_variance()` gets a few new arguments: + + * `null_model`, to provide a null-model to be used for the calculation of + random effect variances. If `NULL`, the null-model is computed internally. + This argument is optional, but may be useful to save time, or when the + null-model cannot be calculated internally. + + * `approximation`, indicating the approximation method for the distribution-specific + (observation level, or residual) variance. + + * `model_component`, for models that can have a zero-inflation component, + specify for which component variances should be returned. By default, both + the conditional and the zero-inflation component are taken into account. + ## Bug fixes * `null_model()` now correctly handles zero-inflated models from package From 531882cbd2a3ab6fd8ea396ca00575912efa5ab8 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 14:31:12 +0200 Subject: [PATCH 05/21] rename variables --- R/compute_variances.R | 373 +++++++++++++++++++++--------------------- R/get_variances.R | 6 +- 2 files changed, 191 insertions(+), 188 deletions(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index a33dd3c87..75779850a 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -1,4 +1,4 @@ -.compute_variances <- function(x, +.compute_variances <- function(model, component, name_fun = NULL, name_full = NULL, @@ -17,7 +17,7 @@ # needed for singularity check check_if_installed("performance", reason = "to check for singularity") - faminfo <- model_info(x, verbose = FALSE) + faminfo <- model_info(model, verbose = FALSE) # check argument approx_method <- match.arg(approximation, c("lognormal", "delta", "trigamma", "observation_level")) @@ -54,7 +54,7 @@ # get necessary model information, like fixed and random effects, # variance-covariance matrix etc. vals <- .get_variance_information( - x, + model, faminfo = faminfo, name_fun = name_fun, verbose = verbose @@ -63,7 +63,7 @@ # we also need necessary model information, like fixed and random effects, # variance-covariance matrix etc. for the null model if (is.null(model_null)) { - model_null <- .safe(null_model(x, verbose = FALSE)) + model_null <- .safe(null_model(model, verbose = FALSE)) } vals_null <- .get_variance_information( model_null, @@ -74,7 +74,7 @@ # Test for non-zero random effects ((near) singularity) no_random_variance <- FALSE - singular_fit <- isTRUE(.safe(performance::check_singularity(x, tolerance = tolerance))) + singular_fit <- isTRUE(.safe(performance::check_singularity(model, tolerance = tolerance))) if (singular_fit && !(component %in% c("slope", "intercept"))) { if (verbose) { @@ -104,21 +104,24 @@ } # Are random slopes present as fixed effects? Warn. - if (!.random_slopes_in_fixed(x) && verbose) { + if (!.random_slopes_in_fixed(model) && verbose) { format_warning( sprintf("Random slopes not present as fixed effects. This artificially inflates the conditional %s.", name_full), "Solution: Respecify fixed structure!" ) } - # Separate observation variance from variance of random effects + # Separate observation variance from variance of random effects, e.g. + # if we have a random effects term that has one observation per group + # (i.e. number of random effects "groups" is the same as number of + # observations in the model) nr <- vapply(vals$re, nrow, numeric(1)) - not.obs.terms <- names(nr[nr != n_obs(x)]) - obs.terms <- names(nr[nr == n_obs(x)]) + not.obs.terms <- names(nr[nr != n_obs(model)]) + obs.terms <- names(nr[nr == n_obs(model)]) # Variance of random effects if (component %in% c("random", "all") && isFALSE(no_random_variance)) { - var.random <- .compute_variance_random(not.obs.terms, x = x, vals = vals) + var.random <- .compute_variance_random(model, not.obs.terms, vals = vals) } # Variance of random effects for NULL model @@ -127,8 +130,8 @@ nr <- vapply(vals_null$re, nrow, numeric(1)) not.obs.terms_null <- names(nr[nr != n_obs(model_null)]) var.random_null <- .compute_variance_random( + model = model_null, not.obs.terms_null, - x = model_null, vals = vals_null ) } @@ -138,7 +141,7 @@ if (component %in% c("residual", "distribution", "all")) { var.distribution <- .compute_variance_distribution( - x = x, + model, var.cor = vals$vc, faminfo, model_null = model_null, @@ -152,7 +155,7 @@ if (component %in% c("residual", "dispersion", "all")) { var.dispersion <- .compute_variance_dispersion( - x = x, + model, vals = vals, faminfo = faminfo, obs.terms = obs.terms @@ -163,21 +166,21 @@ var.residual <- var.distribution + var.dispersion } - if (isTRUE(faminfo$is_mixed) || inherits(x, c("wblm", "wbgee"))) { + if (isTRUE(faminfo$is_mixed) || inherits(model, c("wblm", "wbgee"))) { if (component %in% c("intercept", "all")) { - var.intercept <- .between_subject_variance(vals, x) + var.intercept <- .between_subject_variance(vals) } if (component %in% c("slope", "all")) { - var.slope <- .random_slope_variance(vals, x) + var.slope <- .random_slope_variance(model, vals) } if (component %in% c("rho01", "all")) { - cor.slope_intercept <- .random_slope_intercept_corr(vals, x) + cor.slope_intercept <- .random_slope_intercept_corr(model, vals) } if (component %in% c("rho00", "all")) { - cor.slopes <- .random_slopes_corr(vals, x) + cor.slopes <- .random_slopes_corr(model, vals) } } else { var.intercept <- NULL @@ -217,119 +220,119 @@ # # basically, this function should return a list that has the same # structure for any mixed models like this code for lme4: -# beta = lme4::fixef(x), -# X = lme4::getME(x, "X"), -# vc = lme4::VarCorr(x), -# re = lme4::ranef(x) +# beta = lme4::fixef(model), +# X = lme4::getME(model, "X"), +# vc = lme4::VarCorr(model), +# re = lme4::ranef(model) # -.get_variance_information <- function(x, +.get_variance_information <- function(model, faminfo, name_fun = "get_variances", verbose = TRUE, model_component = "conditional") { # sanity check - if (is.null(x)) { + if (is.null(model)) { return(NULL) } # installed? check_if_installed("lme4", reason = "to compute variances for mixed models") - if (inherits(x, "lme")) { + if (inherits(model, "lme")) { check_if_installed("nlme", reason = "to compute variances for mixed models") } - if (inherits(x, "rstanarm")) { + if (inherits(model, "rstanarm")) { check_if_installed("rstanarm", reason = "to compute variances for mixed models") } # stanreg # --------------------------- - if (inherits(x, "stanreg")) { + if (inherits(model, "stanreg")) { vals <- list( - beta = lme4::fixef(x), - X = rstanarm::get_x(x), - vc = lme4::VarCorr(x), - re = lme4::ranef(x) + beta = lme4::fixef(model), + X = rstanarm::get_x(model), + vc = lme4::VarCorr(model), + re = lme4::ranef(model) ) # GLMMapdative # --------------------------- - } else if (inherits(x, "MixMod")) { + } else if (inherits(model, "MixMod")) { vc1 <- vc2 <- NULL - re_names <- find_random(x) + re_names <- find_random(model) - vc_cond <- !startsWith(colnames(x$D), "zi_") + vc_cond <- !startsWith(colnames(model$D), "zi_") if (any(vc_cond)) { - vc1 <- x$D[vc_cond, vc_cond, drop = FALSE] + vc1 <- model$D[vc_cond, vc_cond, drop = FALSE] attr(vc1, "stddev") <- sqrt(diag(vc1)) - attr(vc1, "correlation") <- stats::cov2cor(x$D[vc_cond, vc_cond, drop = FALSE]) + attr(vc1, "correlation") <- stats::cov2cor(model$D[vc_cond, vc_cond, drop = FALSE]) } - vc_zi <- startsWith(colnames(x$D), "zi_") + vc_zi <- startsWith(colnames(model$D), "zi_") if (any(vc_zi)) { - colnames(x$D) <- gsub("^zi_(.*)", "\\1", colnames(x$D)) - rownames(x$D) <- colnames(x$D) - vc2 <- x$D[vc_zi, vc_zi, drop = FALSE] + colnames(model$D) <- gsub("^zi_(.*)", "\\1", colnames(model$D)) + rownames(model$D) <- colnames(model$D) + vc2 <- model$D[vc_zi, vc_zi, drop = FALSE] attr(vc2, "stddev") <- sqrt(diag(vc2)) - attr(vc2, "correlation") <- stats::cov2cor(x$D[vc_zi, vc_zi, drop = FALSE]) + attr(vc2, "correlation") <- stats::cov2cor(model$D[vc_zi, vc_zi, drop = FALSE]) } vc1 <- list(vc1) names(vc1) <- re_names[[1]] - attr(vc1, "sc") <- sqrt(get_deviance(x, verbose = FALSE) / get_df(x, type = "residual", verbose = FALSE)) + attr(vc1, "sc") <- sqrt(get_deviance(model, verbose = FALSE) / get_df(model, type = "residual", verbose = FALSE)) if (!is.null(vc2)) { vc2 <- list(vc2) names(vc2) <- re_names[[2]] - attr(vc2, "sc") <- sqrt(get_deviance(x, verbose = FALSE) / get_df(x, type = "residual", verbose = FALSE)) + attr(vc2, "sc") <- sqrt(get_deviance(model, verbose = FALSE) / get_df(model, type = "residual", verbose = FALSE)) } vcorr <- compact_list(list(vc1, vc2)) names(vcorr) <- c("cond", "zi")[seq_along(vcorr)] vals <- list( - beta = lme4::fixef(x), - X = get_modelmatrix(x), + beta = lme4::fixef(model), + X = get_modelmatrix(model), vc = vcorr, - re = list(lme4::ranef(x)) + re = list(lme4::ranef(model)) ) - names(vals$re) <- x$id_name + names(vals$re) <- model$id_name # joineRML # --------------------------- - } else if (inherits(x, "mjoint")) { - re_names <- find_random(x, flatten = TRUE) - vcorr <- summary(x)$D + } else if (inherits(model, "mjoint")) { + re_names <- find_random(model, flatten = TRUE) + vcorr <- summary(model)$D attr(vcorr, "stddev") <- sqrt(diag(vcorr)) attr(vcorr, "correlation") <- stats::cov2cor(vcorr) vcorr <- list(vcorr) names(vcorr) <- re_names[1] - attr(vcorr, "sc") <- x$coef$sigma2[[1]] + attr(vcorr, "sc") <- model$coef$sigma2[[1]] vals <- list( - beta = lme4::fixef(x), - X = matrix(1, nrow = n_obs(x), dimnames = list(NULL, "(Intercept)_1")), + beta = lme4::fixef(model), + X = matrix(1, nrow = n_obs(model), dimnames = list(NULL, "(Intercept)_1")), vc = vcorr, - re = list(lme4::ranef(x)) + re = list(lme4::ranef(model)) ) names(vals$re) <- re_names[seq_along(vals$re)] # nlme / glmmPQL # --------------------------- - } else if (inherits(x, "lme")) { - re_names <- find_random(x, split_nested = TRUE, flatten = TRUE) - comp_x <- get_modelmatrix(x) + } else if (inherits(model, "lme")) { + re_names <- find_random(model, split_nested = TRUE, flatten = TRUE) + comp_x <- get_modelmatrix(model) rownames(comp_x) <- seq_len(nrow(comp_x)) - if (.is_nested_lme(x)) { - vals_vc <- .get_nested_lme_varcorr(x) - vals_re <- lme4::ranef(x) + if (.is_nested_lme(model)) { + vals_vc <- .get_nested_lme_varcorr(model) + vals_re <- lme4::ranef(model) } else { - vals_vc <- list(nlme::getVarCov(x)) - vals_re <- list(lme4::ranef(x)) + vals_vc <- list(nlme::getVarCov(model)) + vals_re <- list(lme4::ranef(model)) } vals <- list( - beta = lme4::fixef(x), + beta = lme4::fixef(model), X = comp_x, vc = vals_vc, re = vals_re @@ -339,34 +342,34 @@ # ordinal # --------------------------- - } else if (inherits(x, "clmm")) { + } else if (inherits(model, "clmm")) { if (requireNamespace("ordinal", quietly = TRUE)) { - mm <- get_modelmatrix(x) + mm <- get_modelmatrix(model) vals <- list( - beta = c("(Intercept)" = 1, stats::coef(x)[intersect(names(stats::coef(x)), colnames(mm))]), + beta = c("(Intercept)" = 1, stats::coef(model)[intersect(names(stats::coef(model)), colnames(mm))]), X = mm, - vc = ordinal::VarCorr(x), - re = ordinal::ranef(x) + vc = ordinal::VarCorr(model), + re = ordinal::ranef(model) ) } # glmmadmb # --------------------------- - } else if (inherits(x, "glmmadmb")) { + } else if (inherits(model, "glmmadmb")) { vals <- list( - beta = lme4::fixef(x), - X = get_modelmatrix(x), - vc = lme4::VarCorr(x), - re = lme4::ranef(x) + beta = lme4::fixef(model), + X = get_modelmatrix(model), + vc = lme4::VarCorr(model), + re = lme4::ranef(model) ) # brms # --------------------------- - } else if (inherits(x, "brmsfit")) { - comp_x <- get_modelmatrix(x) + } else if (inherits(model, "brmsfit")) { + comp_x <- get_modelmatrix(model) rownames(comp_x) <- seq_len(nrow(comp_x)) - vc <- lapply(names(lme4::VarCorr(x)), function(i) { - element <- lme4::VarCorr(x)[[i]] + vc <- lapply(names(lme4::VarCorr(model)), function(i) { + element <- lme4::VarCorr(model)[[i]] if (i != "residual__") { if (is.null(element$cov)) { out <- as.matrix(drop(element$sd[, 1])^2) @@ -382,13 +385,13 @@ out }) vc <- compact_list(vc) - names(vc) <- setdiff(names(lme4::VarCorr(x)), "residual__") - attr(vc, "sc") <- lme4::VarCorr(x)$residual__$sd[1, 1] + names(vc) <- setdiff(names(lme4::VarCorr(model)), "residual__") + attr(vc, "sc") <- lme4::VarCorr(model)$residual__$sd[1, 1] vals <- list( - beta = lme4::fixef(x)[, 1], + beta = lme4::fixef(model)[, 1], X = comp_x, vc = vc, - re = lapply(lme4::ranef(x), function(re) { + re = lapply(lme4::ranef(model), function(re) { reval <- as.data.frame(drop(re[, 1, ])) colnames(reval) <- gsub("Intercept", "(Intercept)", dimnames(re)[[3]], fixed = TRUE) reval @@ -398,32 +401,32 @@ # cpglmm # --------------------------- - } else if (inherits(x, "cpglmm")) { + } else if (inherits(model, "cpglmm")) { # installed? check_if_installed("cplm") vals <- list( - beta = cplm::fixef(x), - X = cplm::model.matrix(x), - vc = cplm::VarCorr(x), - re = cplm::ranef(x) + beta = cplm::fixef(model), + X = cplm::model.matrix(model), + vc = cplm::VarCorr(model), + re = cplm::ranef(model) ) # lme4 / glmmTMB # --------------------------- } else { vals <- list( - beta = lme4::fixef(x), - X = lme4::getME(x, "X"), - vc = lme4::VarCorr(x), - re = lme4::ranef(x) + beta = lme4::fixef(model), + X = lme4::getME(model, "X"), + vc = lme4::VarCorr(model), + re = lme4::ranef(model) ) } # for glmmTMB, tell user that dispersion model is ignored - if (inherits(x, c("glmmTMB", "MixMod"))) { + if (inherits(model, c("glmmTMB", "MixMod"))) { if (is.null(model_component) || model_component == "conditional") { vals <- lapply(vals, .collapse_cond) } else { @@ -431,7 +434,7 @@ } } - if (!is.null(find_formula(x)[["dispersion"]]) && verbose) { + if (!is.null(find_formula(model)[["dispersion"]]) && verbose) { format_warning(sprintf("%s ignores effects of dispersion model.", name_fun)) } @@ -461,21 +464,21 @@ # glmmTMB returns a list of model information, one for conditional # and one for zero-inflated part, so here we "unlist" it, returning # only the conditional part. -.collapse_cond <- function(x) { - if (is.list(x) && "cond" %in% names(x)) { - x[["cond"]] +.collapse_cond <- function(model) { + if (is.list(model) && "cond" %in% names(model)) { + model[["cond"]] } else { - x + model } } # same as above, but for the zero-inflation component -.collapse_zi <- function(x) { - if (is.list(x) && "zi" %in% names(x)) { - x[["zi"]] +.collapse_zi <- function(model) { + if (is.list(model) && "zi" %in% names(model)) { + model[["zi"]] } else { - x + model } } @@ -497,13 +500,13 @@ # dispersion-specific variance ---- # --------------------------------- -.compute_variance_dispersion <- function(x, vals, faminfo, obs.terms) { +.compute_variance_dispersion <- function(model, vals, faminfo, obs.terms) { if (faminfo$is_linear) { 0 } else if (length(obs.terms) == 0) { 0 } else { - .compute_variance_random(obs.terms, x = x, vals = vals) + .compute_variance_random(model = model, obs.terms, vals = vals) } } @@ -514,7 +517,7 @@ # This is in line with Nakagawa et al. 2017, Suppl. 2, and package MuMIm # see https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2017.0213&file=rsif20170213supp2.pdf # ---------------------------------------------------------------------------- -.compute_variance_random <- function(terms, x, vals) { +.compute_variance_random <- function(model, terms, vals) { if (is.null(terms)) { return(NULL) } @@ -536,7 +539,7 @@ Z <- vals$X[, rn, drop = FALSE] Z.m <- Z %*% Sigma - sum(diag(crossprod(Z.m, Z))) / n_obs(x) + sum(diag(crossprod(Z.m, Z))) / n_obs(model) } # if (inherits(x, "MixMod")) { @@ -563,7 +566,7 @@ # despite being correctly re-formulated in "null_model()", returns slightly # different values for the log/delta/trigamma approximation. # ----------------------------------------------------------------------------- -.compute_variance_distribution <- function(x, +.compute_variance_distribution <- function(model, var.cor, faminfo, model_null = NULL, @@ -573,7 +576,7 @@ model_component = NULL, verbose = TRUE) { # get overdispersion parameter / sigma - sig <- .safe(get_sigma(x)) + sig <- .safe(get_sigma(model)) if (is.null(sig)) { sig <- 1 @@ -596,7 +599,7 @@ probit = , cloglog = , clogloglink = .variance_distributional( - x, + model, faminfo = faminfo, sig = sig, model_null = model_null, @@ -613,7 +616,7 @@ # -------------------------- # we need this to adjust for "cbind()" outcomes - y_factor <- .binomial_response_weight(x) + y_factor <- .binomial_response_weight(model) # for observation level approximation, when we don't want the "fixed" # residual variance, pi^2/3, but the variance based on the distribution @@ -646,7 +649,7 @@ # ------------- resid.variance <- .variance_distributional( - x, + model, faminfo = faminfo, sig = sig, model_null = model_null, @@ -661,7 +664,7 @@ resid.variance <- switch(faminfo$link_function, log = .variance_distributional( - x, + model, faminfo = faminfo, sig = sig, model_null = model_null, @@ -680,7 +683,7 @@ resid.variance <- switch(faminfo$link_function, inverse = , - identity = stats::family(x)$variance(sig), + identity = stats::family(model)$variance(sig), log = switch(approx_method, delta = 1 / sig^-2, trigamma = trigamma(sig^-2), @@ -694,7 +697,7 @@ resid.variance <- switch(faminfo$link_function, logit = .variance_distributional( - x, + model, faminfo = faminfo, sig = sig, model_null = model_null, @@ -712,7 +715,7 @@ resid.variance <- switch(faminfo$link_function, logit = .variance_distributional( - x, + model, faminfo = faminfo, sig = sig, model_null = model_null, @@ -736,12 +739,12 @@ # # calculate weight if we have a "cbind()" response for binomial models # needed to weight the residual variance -.binomial_response_weight <- function(x) { +.binomial_response_weight <- function(model) { # we need this to adjust for "cbind()" outcomes - resp_value <- .safe(stats::model.frame(x)[, find_response(x)]) + resp_value <- .safe(stats::model.frame(model)[, find_response(model)]) # sanity check if (is.null(resp_value)) { - resp_value <- .safe(get_data(x, source = "frame")[, find_response(x)]) + resp_value <- .safe(get_data(model, source = "frame")[, find_response(model)]) } if (!is.null(resp_value) && !is.null(ncol(resp_value)) && ncol(resp_value) > 1) { mean(rowSums(resp_value, na.rm = TRUE)) @@ -778,7 +781,7 @@ # VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation # VarOtF <- trigamma((1 / lambda + 1 / thetaF)^(-1)) # trigamma function # ----------------------------------------------------------------------------- -.variance_distributional <- function(x, +.variance_distributional <- function(model, faminfo, sig, model_null = NULL, @@ -826,7 +829,7 @@ # transform expected mean of the null model if (is.null(faminfo$family)) { - mu <- link_inverse(x)(mu) + mu <- link_inverse(model)(mu) } else { # transform mu mu <- switch(faminfo$family, @@ -842,7 +845,7 @@ negbinomial = , tweedie = , `negative binomial` = exp(mu + 0.5 * as.vector(revar_null)), - link_inverse(x)(mu) ## TODO: check if this is better than "exp(mu)" + link_inverse(model)(mu) ## TODO: check if this is better than "exp(mu)" # exp(mu) ) } @@ -860,7 +863,7 @@ if (faminfo$link_function == "tweedie") { # Tweedie models ------------------------------------------------------ # --------------------------------------------------------------------- - dispersion_param <- .variance_family_tweedie(x, mu, sig) + dispersion_param <- .variance_family_tweedie(model, mu, sig) } else if (identical(model_component, "full") && (faminfo$is_zero_inflated || faminfo$is_hurdle)) { # Zero Inflated models ------------------------------------------------ @@ -870,12 +873,12 @@ # (zero-inflated) poisson ---- # ---------------------------- poisson = , - `zero-inflated poisson` = .variance_family_poisson(x, mu, faminfo), + `zero-inflated poisson` = .variance_family_poisson(model, mu, faminfo), # hurdle-poisson ---- # ------------------- `hurdle poisson` = , - truncated_poisson = stats::family(x)$variance(sig), + truncated_poisson = stats::family(model)$variance(sig), # (zero-inflated) negative binomial ---- # -------------------------------------- @@ -888,7 +891,7 @@ # hurdle negative binomial ---- # ----------------------------- - truncated_nbinom2 = stats::family(x)$variance(mu, sig), + truncated_nbinom2 = stats::family(model)$variance(mu, sig), # others ---- # ----------- @@ -903,11 +906,11 @@ # (generalized, compoised) poisson ---- # ---------------------------- poisson = 1, - genpois = .variance_family_nbinom(x, mu, sig, faminfo), + genpois = .variance_family_nbinom(model, mu, sig, faminfo), # Gamma, exponential ---- # ----------------------- - Gamma = stats::family(x)$variance(sig), + Gamma = stats::family(model)$variance(sig), # negative binomial ---- # -------------------------------------- @@ -920,21 +923,21 @@ # beta-alike ---- # --------------- - beta = .variance_family_beta(x, mu, sig), - ordbeta = .variance_family_orderedbeta(x, mu), - betabinomial = .variance_family_betabinom(x, mu, sig), + beta = .variance_family_beta(model, mu, sig), + ordbeta = .variance_family_orderedbeta(model, mu), + betabinomial = .variance_family_betabinom(model, mu, sig), ## TODO: check alternatives, but probably less accurate - # betabinomial = .variance_family_beta(x, mu, sig), - # betabinomial = stats::family(x)$variance(mu, sig), + # betabinomial = .variance_family_beta(model, mu, sig), + # betabinomial = stats::family(model)$variance(mu, sig), # other distributions ---- # ------------------------ - tweedie = .variance_family_tweedie(x, mu, sig), + tweedie = .variance_family_tweedie(model, mu, sig), # default variance for non-captured distributions ---- # ---------------------------------------------------- - .variance_family_default(x, mu, verbose) + .variance_family_default(model, mu, verbose) ) } @@ -974,11 +977,11 @@ ) } }, - error = function(x) { + error = function(e) { if (verbose) { format_warning( "Can't calculate model's distribution-specific variance. Results are not reliable.", - paste("The following error occured: ", x$message) + paste("The following error occured: ", e$message) ) } 0 @@ -996,15 +999,15 @@ # Get distributional variance for poisson-family # ---------------------------------------------- -.variance_family_poisson <- function(x, mu, faminfo) { +.variance_family_poisson <- function(model, mu, faminfo) { if (faminfo$is_zero_inflated) { - .variance_zip(x, faminfo, family_var = mu) - } else if (inherits(x, "MixMod")) { + .variance_zip(model, faminfo, family_var = mu) + } else if (inherits(model, "MixMod")) { mu - } else if (inherits(x, "cpglmm")) { - .get_cplm_family(x)$variance(mu) + } else if (inherits(model, "cpglmm")) { + .get_cplm_family(model)$variance(mu) } else { - stats::family(x)$variance(mu) + stats::family(model)$variance(mu) } } @@ -1012,10 +1015,10 @@ # Get distributional variance for beta-family # ---------------------------------------------- -.variance_family_beta <- function(x, mu, phi) { - stats::family(x)$variance(mu) - # if (inherits(x, "MixMod")) { - # stats::family(x)$variance(mu) +.variance_family_beta <- function(model, mu, phi) { + stats::family(model)$variance(mu) + # if (inherits(model, "MixMod")) { + # stats::family(model)$variance(mu) # } else { # # was: # # mu * (1 - mu) / (1 + phi) @@ -1028,9 +1031,9 @@ # Get distributional variance for ordered beta-family # ---------------------------------------------- -.variance_family_orderedbeta <- function(x, mu) { - if (inherits(x, "MixMod")) { - stats::family(x)$variance(mu) +.variance_family_orderedbeta <- function(model, mu) { + if (inherits(model, "MixMod")) { + stats::family(model)$variance(mu) } else { mu * (1 - mu) } @@ -1040,9 +1043,9 @@ # Get distributional variance for beta-family # ---------------------------------------------- -.variance_family_betabinom <- function(x, mu, phi) { - if (inherits(x, "MixMod")) { - stats::family(x)$variance(mu) +.variance_family_betabinom <- function(model, mu, phi) { + if (inherits(model, "MixMod")) { + stats::family(model)$variance(mu) } else { mu * (1 - mu) / (1 + phi) # n <- .binomial_response_weight(x) @@ -1054,15 +1057,15 @@ # Get distributional variance for tweedie-family # ---------------------------------------------- -.variance_family_tweedie <- function(x, mu, phi) { - if (inherits(x, "cpglmm")) { - phi <- x@phi - p <- x@p - 2 +.variance_family_tweedie <- function(model, mu, phi) { + if (inherits(model, "cpglmm")) { + phi <- model@phi + p <- model@p - 2 } else { - if ("psi" %in% names(x$fit$par)) { - psi <- x$fit$par["psi"] # glmmmTMB >= 1.1.5 + if ("psi" %in% names(model$fit$par)) { + psi <- model$fit$par["psi"] # glmmmTMB >= 1.1.5 } else { - psi <- x$fit$par["thetaf"] + psi <- model$fit$par["thetaf"] } p <- unname(stats::plogis(psi) + 1) } @@ -1073,17 +1076,17 @@ # Get distributional variance for nbinom-family # ---------------------------------------------- -.variance_family_nbinom <- function(x, mu, sig, faminfo) { +.variance_family_nbinom <- function(model, mu, sig, faminfo) { if (faminfo$is_zero_inflated) { if (missing(sig)) sig <- 0 - .variance_zinb(x, sig, faminfo, family_var = mu * (1 + sig)) - } else if (inherits(x, "MixMod")) { + .variance_zinb(model, sig, faminfo, family_var = mu * (1 + sig)) + } else if (inherits(model, "MixMod")) { if (missing(sig)) { return(rep(1e-16, length(mu))) } mu * (1 + sig) } else { - stats::family(x)$variance(mu, sig) + stats::family(model)$variance(mu, sig) } } @@ -1151,26 +1154,26 @@ # Get distribution-specific variance for general and # undefined families / link-functions # ---------------------------------------------- -.variance_family_default <- function(x, mu, verbose) { +.variance_family_default <- function(model, mu, verbose) { # installed? check_if_installed("lme4") tryCatch( - if (inherits(x, "merMod")) { - mu * (1 + mu / lme4::getME(x, "glmer.nb.theta")) - } else if (inherits(x, "MixMod")) { - stats::family(x)$variance(mu) - } else if (inherits(x, "glmmTMB")) { - if (is.null(x$theta)) { - theta <- lme4::getME(x, "theta") + if (inherits(model, "merMod")) { + mu * (1 + mu / lme4::getME(xmodel, "glmer.nb.theta")) + } else if (inherits(model, "MixMod")) { + stats::family(model)$variance(mu) + } else if (inherits(model, "glmmTMB")) { + if (is.null(model$theta)) { + theta <- lme4::getME(model, "theta") } else { - theta <- x$theta + theta <- model$theta } mu * (1 + mu / theta) } else { - mu * (1 + mu / x$theta) + mu * (1 + mu / model$theta) }, - error = function(x) { + error = function(e) { if (verbose) { format_warning("Can't calculate model's distribution-specific variance. Results are not reliable.") } @@ -1216,7 +1219,7 @@ # random intercept-variances, i.e. # between-subject-variance (tau 00) ---- # ---------------------------------------------- -.between_subject_variance <- function(vals, x) { +.between_subject_variance <- function(vals) { vars <- lapply(vals$vc, function(i) i[1]) # check for uncorrelated random slopes-intercept non_intercepts <- which(sapply(vals$vc, function(i) !startsWith(dimnames(i)[[1]][1], "(Intercept)"))) @@ -1233,12 +1236,12 @@ # random slope-variances (tau 11) ---- # ---------------------------------------------- -.random_slope_variance <- function(vals, x) { - if (inherits(x, "lme")) { - unlist(lapply(vals$vc, function(x) diag(x)[-1])) +.random_slope_variance <- function(model, vals) { + if (inherits(model, "lme")) { + unlist(lapply(vals$vc, function(i) diag(i)[-1])) } else { # random slopes for correlated slope-intercept - out <- unlist(lapply(vals$vc, function(x) diag(x)[-1])) + out <- unlist(lapply(vals$vc, function(i) diag(i)[-1])) # check for uncorrelated random slopes-intercept non_intercepts <- which(sapply(vals$vc, function(i) !startsWith(dimnames(i)[[1]][1], "(Intercept)"))) if (length(non_intercepts) == length(vals$vc)) { @@ -1288,13 +1291,13 @@ # slope-intercept-correlations (rho 01) ---- # ---------------------------------------------- -.random_slope_intercept_corr <- function(vals, x) { - if (inherits(x, "lme")) { +.random_slope_intercept_corr <- function(model, vals) { + if (inherits(model, "lme")) { rho01 <- unlist(sapply(vals$vc, attr, which = "cor_slope_intercept")) if (is.null(rho01)) { - vc <- lme4::VarCorr(x) + vc <- lme4::VarCorr(model) if ("Corr" %in% colnames(vc)) { - re_name <- find_random(x, split_nested = FALSE, flatten = TRUE) + re_name <- find_random(model, split_nested = FALSE, flatten = TRUE) rho01 <- as.vector(suppressWarnings(stats::na.omit(as.numeric(vc[, "Corr"])))) if (length(re_name) == length(rho01)) { names(rho01) <- re_name @@ -1321,15 +1324,15 @@ # slope-slope-correlations (rho 01) ---- # ---------------------------------------------- -.random_slopes_corr <- function(vals, x) { +.random_slopes_corr <- function(model, vals) { corrs <- lapply(vals$vc, attr, "correlation") - rnd_slopes <- unlist(find_random_slopes(x)) + rnd_slopes <- unlist(find_random_slopes(model)) # check if any categorical random slopes. we then have # correlation among factor levels cat_random_slopes <- tryCatch( { - d <- get_data(x, verbose = FALSE)[rnd_slopes] + d <- get_data(model, verbose = FALSE)[rnd_slopes] any(sapply(d, is.factor)) }, error = function(e) { @@ -1373,7 +1376,7 @@ # { # sapply(corrs, function(i) { # if (!is.null(i)) { - # slope_pairs <- utils::combn(x = rnd_slopes, m = 2, simplify = FALSE) + # slope_pairs <- utils::combn(model = rnd_slopes, m = 2, simplify = FALSE) # lapply(slope_pairs, function(j) { # stats::setNames(i[j[1], j[2]], paste0("..", paste0(j, collapse = "-"))) # }) diff --git a/R/get_variances.R b/R/get_variances.R index a5518a3b2..7ef656a1e 100644 --- a/R/get_variances.R +++ b/R/get_variances.R @@ -173,7 +173,7 @@ get_variance.merMod <- function(x, ...) { component <- match.arg(component) .safe(.compute_variances( - x, + model = x, component = component, name_fun = "get_variance", name_full = "random effect variances", @@ -232,7 +232,7 @@ get_variance.glmmTMB <- function(x, ...) { component <- match.arg(component) .safe(.compute_variances( - x, + model = x, component = component, name_fun = "get_variance", name_full = "random effect variances", @@ -262,7 +262,7 @@ get_variance.mixed <- function(x, ...) { component <- match.arg(component) .safe(.compute_variances( - x$full_model, + model = x$full_model, component = component, name_fun = "get_variance", name_full = "random effect variances", From 9e9c65d0615a90cdb84dce82df609bd8b37842a3 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 14:41:27 +0200 Subject: [PATCH 06/21] rename variablkes --- R/compute_variances.R | 130 ++++++++++++++++++++++-------------------- 1 file changed, 67 insertions(+), 63 deletions(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index 75779850a..f49eb2e2e 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -53,7 +53,7 @@ # get necessary model information, like fixed and random effects, # variance-covariance matrix etc. - vals <- .get_variance_information( + mixed_effects_info <- .get_variance_information( model, faminfo = faminfo, name_fun = name_fun, @@ -65,7 +65,7 @@ if (is.null(model_null)) { model_null <- .safe(null_model(model, verbose = FALSE)) } - vals_null <- .get_variance_information( + me_info_null <- .get_variance_information( model_null, faminfo = faminfo, name_fun = name_fun, @@ -100,7 +100,7 @@ # Get variance of fixed effects: multiply coefs by design matrix if (component %in% c("fixed", "all")) { - var.fixed <- .compute_variance_fixed(vals) + var.fixed <- .compute_variance_fixed(mixed_effects_info) } # Are random slopes present as fixed effects? Warn. @@ -115,24 +115,24 @@ # if we have a random effects term that has one observation per group # (i.e. number of random effects "groups" is the same as number of # observations in the model) - nr <- vapply(vals$re, nrow, numeric(1)) + nr <- vapply(mixed_effects_info$re, nrow, numeric(1)) not.obs.terms <- names(nr[nr != n_obs(model)]) obs.terms <- names(nr[nr == n_obs(model)]) # Variance of random effects if (component %in% c("random", "all") && isFALSE(no_random_variance)) { - var.random <- .compute_variance_random(model, not.obs.terms, vals = vals) + var.random <- .compute_variance_random(model, not.obs.terms, mixed_effects_info) } # Variance of random effects for NULL model - if (!singular_fit && !is.null(vals_null)) { + if (!singular_fit && !is.null(me_info_null)) { # Separate observation variance from variance of random effects - nr <- vapply(vals_null$re, nrow, numeric(1)) + nr <- vapply(me_info_null$re, nrow, numeric(1)) not.obs.terms_null <- names(nr[nr != n_obs(model_null)]) var.random_null <- .compute_variance_random( model = model_null, not.obs.terms_null, - vals = vals_null + mixed_effects_info = me_info_null ) } @@ -142,7 +142,7 @@ if (component %in% c("residual", "distribution", "all")) { var.distribution <- .compute_variance_distribution( model, - var.cor = vals$vc, + var.cor = mixed_effects_info$vc, faminfo, model_null = model_null, revar_null = var.random_null, @@ -156,7 +156,7 @@ if (component %in% c("residual", "dispersion", "all")) { var.dispersion <- .compute_variance_dispersion( model, - vals = vals, + mixed_effects_info = mixed_effects_info, faminfo = faminfo, obs.terms = obs.terms ) @@ -168,19 +168,19 @@ if (isTRUE(faminfo$is_mixed) || inherits(model, c("wblm", "wbgee"))) { if (component %in% c("intercept", "all")) { - var.intercept <- .between_subject_variance(vals) + var.intercept <- .between_subject_variance(mixed_effects_info) } if (component %in% c("slope", "all")) { - var.slope <- .random_slope_variance(model, vals) + var.slope <- .random_slope_variance(model, mixed_effects_info) } if (component %in% c("rho01", "all")) { - cor.slope_intercept <- .random_slope_intercept_corr(model, vals) + cor.slope_intercept <- .random_slope_intercept_corr(model, mixed_effects_info) } if (component %in% c("rho00", "all")) { - cor.slopes <- .random_slopes_corr(model, vals) + cor.slopes <- .random_slopes_corr(model, mixed_effects_info) } } else { var.intercept <- NULL @@ -249,7 +249,7 @@ # stanreg # --------------------------- if (inherits(model, "stanreg")) { - vals <- list( + mixed_effects_info <- list( beta = lme4::fixef(model), X = rstanarm::get_x(model), vc = lme4::VarCorr(model), @@ -291,13 +291,13 @@ vcorr <- compact_list(list(vc1, vc2)) names(vcorr) <- c("cond", "zi")[seq_along(vcorr)] - vals <- list( + mixed_effects_info <- list( beta = lme4::fixef(model), X = get_modelmatrix(model), vc = vcorr, re = list(lme4::ranef(model)) ) - names(vals$re) <- model$id_name + names(mixed_effects_info$re) <- model$id_name # joineRML # --------------------------- @@ -310,13 +310,13 @@ names(vcorr) <- re_names[1] attr(vcorr, "sc") <- model$coef$sigma2[[1]] - vals <- list( + mixed_effects_info <- list( beta = lme4::fixef(model), X = matrix(1, nrow = n_obs(model), dimnames = list(NULL, "(Intercept)_1")), vc = vcorr, re = list(lme4::ranef(model)) ) - names(vals$re) <- re_names[seq_along(vals$re)] + names(mixed_effects_info$re) <- re_names[seq_along(mixed_effects_info$re)] # nlme / glmmPQL # --------------------------- @@ -331,21 +331,21 @@ vals_vc <- list(nlme::getVarCov(model)) vals_re <- list(lme4::ranef(model)) } - vals <- list( + mixed_effects_info <- list( beta = lme4::fixef(model), X = comp_x, vc = vals_vc, re = vals_re ) - names(vals$re) <- re_names - names(vals$vc) <- re_names + names(mixed_effects_info$re) <- re_names + names(mixed_effects_info$vc) <- re_names # ordinal # --------------------------- } else if (inherits(model, "clmm")) { if (requireNamespace("ordinal", quietly = TRUE)) { mm <- get_modelmatrix(model) - vals <- list( + mixed_effects_info <- list( beta = c("(Intercept)" = 1, stats::coef(model)[intersect(names(stats::coef(model)), colnames(mm))]), X = mm, vc = ordinal::VarCorr(model), @@ -356,7 +356,7 @@ # glmmadmb # --------------------------- } else if (inherits(model, "glmmadmb")) { - vals <- list( + mixed_effects_info <- list( beta = lme4::fixef(model), X = get_modelmatrix(model), vc = lme4::VarCorr(model), @@ -387,7 +387,7 @@ vc <- compact_list(vc) names(vc) <- setdiff(names(lme4::VarCorr(model)), "residual__") attr(vc, "sc") <- lme4::VarCorr(model)$residual__$sd[1, 1] - vals <- list( + mixed_effects_info <- list( beta = lme4::fixef(model)[, 1], X = comp_x, vc = vc, @@ -397,7 +397,7 @@ reval }) ) - names(vals$beta) <- gsub("Intercept", "(Intercept)", names(vals$beta), fixed = TRUE) + names(mixed_effects_info$beta) <- gsub("Intercept", "(Intercept)", names(mixed_effects_info$beta), fixed = TRUE) # cpglmm # --------------------------- @@ -405,7 +405,7 @@ # installed? check_if_installed("cplm") - vals <- list( + mixed_effects_info <- list( beta = cplm::fixef(model), X = cplm::model.matrix(model), vc = cplm::VarCorr(model), @@ -415,7 +415,7 @@ # lme4 / glmmTMB # --------------------------- } else { - vals <- list( + mixed_effects_info <- list( beta = lme4::fixef(model), X = lme4::getME(model, "X"), vc = lme4::VarCorr(model), @@ -428,9 +428,9 @@ if (inherits(model, c("glmmTMB", "MixMod"))) { if (is.null(model_component) || model_component == "conditional") { - vals <- lapply(vals, .collapse_cond) + mixed_effects_info <- lapply(mixed_effects_info, .collapse_cond) } else { - vals <- lapply(vals, .collapse_zi) + mixed_effects_info <- lapply(mixed_effects_info, .collapse_zi) } } @@ -439,13 +439,13 @@ } # fix rank deficiency - rankdef <- is.na(vals$beta) + rankdef <- is.na(mixed_effects_info$beta) if (any(rankdef)) { - rankdef_names <- names(vals$beta)[rankdef] - vals$beta <- vals$beta[setdiff(names(vals$beta), rankdef_names)] + rankdef_names <- names(mixed_effects_info$beta)[rankdef] + mixed_effects_info$beta <- mixed_effects_info$beta[setdiff(names(mixed_effects_info$beta), rankdef_names)] } - vals + mixed_effects_info } @@ -492,21 +492,21 @@ # However, package MuMIn differs and uses "fitted()" instead, leading to minor # deviations # ----------------------------------------------------------------------------- -.compute_variance_fixed <- function(vals) { - with(vals, stats::var(as.vector(beta %*% t(X)))) +.compute_variance_fixed <- function(mixed_effects_info) { + with(mixed_effects_info, stats::var(as.vector(beta %*% t(X)))) } # dispersion-specific variance ---- # --------------------------------- -.compute_variance_dispersion <- function(model, vals, faminfo, obs.terms) { +.compute_variance_dispersion <- function(model, mixed_effects_info, faminfo, obs.terms) { if (faminfo$is_linear) { 0 } else if (length(obs.terms) == 0) { 0 } else { - .compute_variance_random(model = model, obs.terms, vals = vals) + .compute_variance_random(model, obs.terms, mixed_effects_info) } } @@ -517,7 +517,7 @@ # This is in line with Nakagawa et al. 2017, Suppl. 2, and package MuMIm # see https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2017.0213&file=rsif20170213supp2.pdf # ---------------------------------------------------------------------------- -.compute_variance_random <- function(model, terms, vals) { +.compute_variance_random <- function(model, terms, mixed_effects_info) { if (is.null(terms)) { return(NULL) } @@ -525,29 +525,29 @@ rn <- rownames(Sigma) # fix for models w/o intercept - if (!any(startsWith(colnames(vals$X), "(Intercept)"))) { - vals$X <- cbind("(Intercept)" = 1, vals$X) + if (!any(startsWith(colnames(mixed_effects_info$X), "(Intercept)"))) { + mixed_effects_info$X <- cbind("(Intercept)" = 1, mixed_effects_info$X) } if (!is.null(rn)) { - valid <- rownames(Sigma) %in% colnames(vals$X) + valid <- rownames(Sigma) %in% colnames(mixed_effects_info$X) if (!all(valid)) { rn <- rn[valid] Sigma <- Sigma[valid, valid, drop = FALSE] } } - Z <- vals$X[, rn, drop = FALSE] + Z <- mixed_effects_info$X[, rn, drop = FALSE] Z.m <- Z %*% Sigma sum(diag(crossprod(Z.m, Z))) / n_obs(model) } # if (inherits(x, "MixMod")) { - # .sigma_sum(vals$vc) + # .sigma_sum(mixed_effects_info$vc) # } else { - # sum(sapply(vals$vc[terms], .sigma_sum)) + # sum(sapply(mixed_effects_info$vc[terms], .sigma_sum)) # } - sum(sapply(vals$vc[terms], .sigma_sum)) + sum(sapply(mixed_effects_info$vc[terms], .sigma_sum)) } @@ -1219,10 +1219,12 @@ # random intercept-variances, i.e. # between-subject-variance (tau 00) ---- # ---------------------------------------------- -.between_subject_variance <- function(vals) { - vars <- lapply(vals$vc, function(i) i[1]) +.between_subject_variance <- function(mixed_effects_info) { + vars <- lapply(mixed_effects_info$vc, function(i) i[1]) # check for uncorrelated random slopes-intercept - non_intercepts <- which(sapply(vals$vc, function(i) !startsWith(dimnames(i)[[1]][1], "(Intercept)"))) + non_intercepts <- which(sapply(mixed_effects_info$vc, function(i) { + !startsWith(dimnames(i)[[1]][1], "(Intercept)") + })) if (length(non_intercepts)) { vars <- vars[-non_intercepts] } @@ -1236,19 +1238,21 @@ # random slope-variances (tau 11) ---- # ---------------------------------------------- -.random_slope_variance <- function(model, vals) { +.random_slope_variance <- function(model, mixed_effects_info) { if (inherits(model, "lme")) { - unlist(lapply(vals$vc, function(i) diag(i)[-1])) + unlist(lapply(mixed_effects_info$vc, function(i) diag(i)[-1])) } else { # random slopes for correlated slope-intercept - out <- unlist(lapply(vals$vc, function(i) diag(i)[-1])) + out <- unlist(lapply(mixed_effects_info$vc, function(i) diag(i)[-1])) # check for uncorrelated random slopes-intercept - non_intercepts <- which(sapply(vals$vc, function(i) !startsWith(dimnames(i)[[1]][1], "(Intercept)"))) - if (length(non_intercepts) == length(vals$vc)) { - out <- unlist(lapply(vals$vc, diag)) + non_intercepts <- which(sapply(mixed_effects_info$vc, function(i) { + !startsWith(dimnames(i)[[1]][1], "(Intercept)") + })) + if (length(non_intercepts) == length(mixed_effects_info$vc)) { + out <- unlist(lapply(mixed_effects_info$vc, diag)) } else if (length(non_intercepts)) { - dn <- unlist(lapply(vals$vc, function(i) dimnames(i)[1])[non_intercepts]) - rndslopes <- unlist(lapply(vals$vc, function(i) { + dn <- unlist(lapply(mixed_effects_info$vc, function(i) dimnames(i)[1])[non_intercepts]) + rndslopes <- unlist(lapply(mixed_effects_info$vc, function(i) { if (is.null(dim(i)) || identical(dim(i), c(1, 1))) { as.vector(i) } else { @@ -1291,9 +1295,9 @@ # slope-intercept-correlations (rho 01) ---- # ---------------------------------------------- -.random_slope_intercept_corr <- function(model, vals) { +.random_slope_intercept_corr <- function(model, mixed_effects_info) { if (inherits(model, "lme")) { - rho01 <- unlist(sapply(vals$vc, attr, which = "cor_slope_intercept")) + rho01 <- unlist(sapply(mixed_effects_info$vc, attr, which = "cor_slope_intercept")) if (is.null(rho01)) { vc <- lme4::VarCorr(model) if ("Corr" %in% colnames(vc)) { @@ -1306,7 +1310,7 @@ } rho01 } else { - corrs <- lapply(vals$vc, attr, "correlation") + corrs <- lapply(mixed_effects_info$vc, attr, "correlation") rho01 <- sapply(corrs, function(i) { if (!is.null(i) && colnames(i)[1] == "(Intercept)") { i[-1, 1] @@ -1324,8 +1328,8 @@ # slope-slope-correlations (rho 01) ---- # ---------------------------------------------- -.random_slopes_corr <- function(model, vals) { - corrs <- lapply(vals$vc, attr, "correlation") +.random_slopes_corr <- function(model, mixed_effects_info) { + corrs <- lapply(mixed_effects_info$vc, attr, "correlation") rnd_slopes <- unlist(find_random_slopes(model)) # check if any categorical random slopes. we then have @@ -1333,7 +1337,7 @@ cat_random_slopes <- tryCatch( { d <- get_data(model, verbose = FALSE)[rnd_slopes] - any(sapply(d, is.factor)) + any(vapply(d, is.factor, logical(1))) }, error = function(e) { NULL From db490b3ffe43b464339a946279f151b9310bd374 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 14:51:23 +0200 Subject: [PATCH 07/21] rename --- R/compute_variances.R | 57 +++++++++++++++++++++---------------------- 1 file changed, 28 insertions(+), 29 deletions(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index f49eb2e2e..aabfc80b8 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -11,7 +11,7 @@ ## Author: Ben Bolker, who used an cleaned-up/adapted ## version of Jon Lefcheck's code from SEMfit - ## Major revisions and adaption to more complex models and other packages + ## Revisions and adaption to more complex models and other packages ## by Daniel Lüdecke # needed for singularity check @@ -116,22 +116,22 @@ # (i.e. number of random effects "groups" is the same as number of # observations in the model) nr <- vapply(mixed_effects_info$re, nrow, numeric(1)) - not.obs.terms <- names(nr[nr != n_obs(model)]) - obs.terms <- names(nr[nr == n_obs(model)]) + not_obs_terms <- names(nr[nr != n_obs(model)]) + obs_terms <- names(nr[nr == n_obs(model)]) # Variance of random effects if (component %in% c("random", "all") && isFALSE(no_random_variance)) { - var.random <- .compute_variance_random(model, not.obs.terms, mixed_effects_info) + var.random <- .compute_variance_random(model, not_obs_terms, mixed_effects_info) } # Variance of random effects for NULL model if (!singular_fit && !is.null(me_info_null)) { # Separate observation variance from variance of random effects nr <- vapply(me_info_null$re, nrow, numeric(1)) - not.obs.terms_null <- names(nr[nr != n_obs(model_null)]) + not_obs_terms_null <- names(nr[nr != n_obs(model_null)]) var.random_null <- .compute_variance_random( model = model_null, - not.obs.terms_null, + not_obs_terms_null, mixed_effects_info = me_info_null ) } @@ -142,7 +142,7 @@ if (component %in% c("residual", "distribution", "all")) { var.distribution <- .compute_variance_distribution( model, - var.cor = mixed_effects_info$vc, + var_cor = mixed_effects_info$vc, faminfo, model_null = model_null, revar_null = var.random_null, @@ -158,7 +158,7 @@ model, mixed_effects_info = mixed_effects_info, faminfo = faminfo, - obs.terms = obs.terms + obs_terms = obs_terms ) } @@ -228,8 +228,7 @@ .get_variance_information <- function(model, faminfo, name_fun = "get_variances", - verbose = TRUE, - model_component = "conditional") { + verbose = TRUE) { # sanity check if (is.null(model)) { return(NULL) @@ -424,16 +423,16 @@ } - # for glmmTMB, tell user that dispersion model is ignored - + # for models with zero-inflation, we only want the conditional part if (inherits(model, c("glmmTMB", "MixMod"))) { - if (is.null(model_component) || model_component == "conditional") { - mixed_effects_info <- lapply(mixed_effects_info, .collapse_cond) - } else { - mixed_effects_info <- lapply(mixed_effects_info, .collapse_zi) - } + mixed_effects_info <- lapply(mixed_effects_info, .collapse_cond) } + # currently, we don't support calculating all variance components + # for the zero-inflated part of the model only. This is not fully implemented + # mixed_effects_info <- lapply(mixed_effects_info, .collapse_zi) + + # for glmmTMB, tell user that dispersion model is ignored if (!is.null(find_formula(model)[["dispersion"]]) && verbose) { format_warning(sprintf("%s ignores effects of dispersion model.", name_fun)) } @@ -464,21 +463,21 @@ # glmmTMB returns a list of model information, one for conditional # and one for zero-inflated part, so here we "unlist" it, returning # only the conditional part. -.collapse_cond <- function(model) { - if (is.list(model) && "cond" %in% names(model)) { - model[["cond"]] +.collapse_cond <- function(x) { + if (is.list(x) && "cond" %in% names(x)) { + x[["cond"]] } else { - model + x } } # same as above, but for the zero-inflation component -.collapse_zi <- function(model) { - if (is.list(model) && "zi" %in% names(model)) { - model[["zi"]] +.collapse_zi <- function(x) { + if (is.list(x) && "zi" %in% names(x)) { + x[["zi"]] } else { - model + x } } @@ -500,13 +499,13 @@ # dispersion-specific variance ---- # --------------------------------- -.compute_variance_dispersion <- function(model, mixed_effects_info, faminfo, obs.terms) { +.compute_variance_dispersion <- function(model, mixed_effects_info, faminfo, obs_terms) { if (faminfo$is_linear) { 0 - } else if (length(obs.terms) == 0) { + } else if (length(obs_terms) == 0) { 0 } else { - .compute_variance_random(model, obs.terms, mixed_effects_info) + .compute_variance_random(model, obs_terms, mixed_effects_info) } } @@ -567,7 +566,7 @@ # different values for the log/delta/trigamma approximation. # ----------------------------------------------------------------------------- .compute_variance_distribution <- function(model, - var.cor, + var_cor, faminfo, model_null = NULL, revar_null = NULL, From b2c8ef7e12bbeef0932e80d27a04083981ed9ca7 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 14:58:38 +0200 Subject: [PATCH 08/21] docs --- R/get_variances.R | 9 ++++++++- man/get_variance.Rd | 9 ++++++++- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/R/get_variances.R b/R/get_variances.R index 7ef656a1e..e51a02429 100644 --- a/R/get_variances.R +++ b/R/get_variances.R @@ -111,7 +111,14 @@ #' from **blme**), `clmm`, `cpglmm`, `glmmadmb`, `glmmTMB`, `MixMod`, `lme`, #' `mixed`, `rlmerMod`, `stanreg`, `brmsfit` or `wbm`. Support for objects of #' class `MixMod` (**GLMMadaptive**), `lme` (**nlme**) or `brmsfit` (**brms**) is -#' experimental and may not work for all models. +#' not fully implemented or tested, and therefore may not work for all models +#' of the aforementioned classes. +#' +#' Extracting variance components for models with zero-inflation part is not +#' straightforward, because it is not definitly clear how the distribution-specific +#' variance should be calculated. Therefore, it is recommended to carefully +#' inspect the results, and probably validate against other models, e.g. Bayesian +#' models (although results may be only roughly comparable). #' #' @references #' - Johnson, P. C. D. (2014). Extension of Nakagawa & Schielzeth’s R2 GLMM to diff --git a/man/get_variance.Rd b/man/get_variance.Rd index e23e00484..e21334c06 100644 --- a/man/get_variance.Rd +++ b/man/get_variance.Rd @@ -120,7 +120,14 @@ This function supports models of class \code{merMod} (including models from \strong{blme}), \code{clmm}, \code{cpglmm}, \code{glmmadmb}, \code{glmmTMB}, \code{MixMod}, \code{lme}, \code{mixed}, \code{rlmerMod}, \code{stanreg}, \code{brmsfit} or \code{wbm}. Support for objects of class \code{MixMod} (\strong{GLMMadaptive}), \code{lme} (\strong{nlme}) or \code{brmsfit} (\strong{brms}) is -experimental and may not work for all models. +not fully implemented or tested, and therefore may not work for all models +of the aforementioned classes. + +Extracting variance components for models with zero-inflation part is not +straightforward, because it is not definitly clear how the distribution-specific +variance should be calculated. Therefore, it is recommended to carefully +inspect the results, and probably validate against other models, e.g. Bayesian +models (although results may be only roughly comparable). } \section{Fixed effects variance}{ From 5dff93f2530115e1ebf458867c708596eff6f074 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 15:09:16 +0200 Subject: [PATCH 09/21] fix --- R/compute_variances.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index aabfc80b8..cbfdd8bfb 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -1159,7 +1159,7 @@ tryCatch( if (inherits(model, "merMod")) { - mu * (1 + mu / lme4::getME(xmodel, "glmer.nb.theta")) + mu * (1 + mu / lme4::getME(model, "glmer.nb.theta")) } else if (inherits(model, "MixMod")) { stats::family(model)$variance(mu) } else if (inherits(model, "glmmTMB")) { From 3d50ddc2d300a385eff8a568cbd6e94743cd8141 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 15:41:56 +0200 Subject: [PATCH 10/21] fix test --- tests/testthat/test-r2_nakagawa_negbin_zi.R | 27 ++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-r2_nakagawa_negbin_zi.R b/tests/testthat/test-r2_nakagawa_negbin_zi.R index 5bc50ada6..582adde60 100644 --- a/tests/testthat/test-r2_nakagawa_negbin_zi.R +++ b/tests/testthat/test-r2_nakagawa_negbin_zi.R @@ -85,24 +85,45 @@ test_that("glmmTMB, Nbinom1 zero-inflated", { VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation VarOtF <- trigamma((1 / lambda + 1 / thetaF)^-1) # trigamma function + # test message + expect_message( + performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, model_component = "conditional"), + regex = "Zero-inflation part of" + ) + # lognormal R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF) R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF) - out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr) + out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, model_component = "conditional", verbose = FALSE) expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE) expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE) + # full model + out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr) + expect_equal(out$R2_conditional, 0.7543869, tolerance = 1e-4, ignore_attr = TRUE) + expect_equal(out$R2_marginal, 0.6448776, tolerance = 1e-4, ignore_attr = TRUE) + # delta R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF) R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF) - out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "delta") + out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "delta", verbose = FALSE, model_component = "conditional") expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE) expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE) + # full model + out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "delta") + expect_equal(out$R2_conditional, 0.6942946, tolerance = 1e-4, ignore_attr = TRUE) + expect_equal(out$R2_marginal, 0.5935086, tolerance = 1e-4, ignore_attr = TRUE) + # trigamma R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF) R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF) - out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "trigamma") + out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "trigamma", verbose = FALSE, model_component = "conditional") expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE) expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE) + + # full model + out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "trigamma") + expect_equal(out$R2_conditional, 0.6051817, tolerance = 1e-4, ignore_attr = TRUE) + expect_equal(out$R2_marginal, 0.5173316, tolerance = 1e-4, ignore_attr = TRUE) }) From 894bd35901985d14a6fe0e06ec9ff4f14ebdc5a2 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 15:51:57 +0200 Subject: [PATCH 11/21] fix test --- tests/testthat/test-r2_nakagawa_poisson_zi.R | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-r2_nakagawa_poisson_zi.R b/tests/testthat/test-r2_nakagawa_poisson_zi.R index b0bb27a3b..d565d81b3 100644 --- a/tests/testthat/test-r2_nakagawa_poisson_zi.R +++ b/tests/testthat/test-r2_nakagawa_poisson_zi.R @@ -19,13 +19,18 @@ test_that("glmmTMB, Poisson zero-inflated", { family = poisson(), data = Salamanders ) out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m)) - out2 <- performance::r2_nakagawa(m) + out2 <- performance::r2_nakagawa(m, model_component = "conditional", verbose = FALSE) # matches theoretical values expect_equal(out2$R2_marginal, 0.4636197, ignore_attr = TRUE, tolerance = 1e-4) expect_equal(out2$R2_conditional, 0.5751936, ignore_attr = TRUE, tolerance = 1e-4) expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-1) expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-1) + # full model + out <- performance::r2_nakagawa(m) + expect_equal(out2$R2_marginal, 0.4636197, ignore_attr = TRUE, tolerance = 1e-4) + expect_equal(out2$R2_conditional, 0.5751936, ignore_attr = TRUE, tolerance = 1e-4) + # glmmTMB, sqrt, no random slope ------------------------------------------------- m <- glmmTMB::glmmTMB(count ~ mined + (1 | site), ziformula = ~mined, @@ -43,7 +48,7 @@ test_that("glmmTMB, Poisson zero-inflated", { family = poisson(), data = Salamanders )) out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m)) - out2 <- suppressWarnings(performance::r2_nakagawa(m, tolerance = 1e-8)) + out2 <- suppressWarnings(performance::r2_nakagawa(m, tolerance = 1e-8, model_component = "conditional", verbose = FALSE)) # we have slight differences here: MuMIn uses "var(fitted())" to exctract fixed # effects variances, while insight uses "var(beta %*% t(mm))". The latter gives # different values when random slopes are involved From 9fd98e22b539506c7de4babdc5cc292d9060682f Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 15:54:16 +0200 Subject: [PATCH 12/21] fix test --- tests/testthat/test-r2_nakagawa_poisson_zi.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-r2_nakagawa_poisson_zi.R b/tests/testthat/test-r2_nakagawa_poisson_zi.R index d565d81b3..199dd51d7 100644 --- a/tests/testthat/test-r2_nakagawa_poisson_zi.R +++ b/tests/testthat/test-r2_nakagawa_poisson_zi.R @@ -28,8 +28,8 @@ test_that("glmmTMB, Poisson zero-inflated", { # full model out <- performance::r2_nakagawa(m) - expect_equal(out2$R2_marginal, 0.4636197, ignore_attr = TRUE, tolerance = 1e-4) - expect_equal(out2$R2_conditional, 0.5751936, ignore_attr = TRUE, tolerance = 1e-4) + expect_equal(out$R2_marginal, 0.2923215, ignore_attr = TRUE, tolerance = 1e-4) + expect_equal(out$R2_conditional, 0.362671, ignore_attr = TRUE, tolerance = 1e-4) # glmmTMB, sqrt, no random slope ------------------------------------------------- m <- glmmTMB::glmmTMB(count ~ mined + (1 | site), @@ -57,6 +57,12 @@ test_that("glmmTMB, Poisson zero-inflated", { expect_equal(out2$R2_marginal, 0.524714, ignore_attr = TRUE, tolerance = 1e-1) expect_equal(out2$R2_conditional, 0.6498465, ignore_attr = TRUE, tolerance = 1e-1) + # full model + out <- suppressWarnings(performance::r2_nakagawa(m)) + expect_equal(out$R2_marginal, 0.3344432, ignore_attr = TRUE, tolerance = 1e-4) + expect_equal(out$R2_conditional, 0.4142003, ignore_attr = TRUE, tolerance = 1e-4) + + # glmmTMB, sqrt, random slope ------------------------------------------------- m <- suppressWarnings(glmmTMB::glmmTMB(count ~ mined + cover + (1 + cover | site), ziformula = ~mined, From 58e0d640ef62fbcba4c3693a6407be14e9d39519 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 16:05:53 +0200 Subject: [PATCH 13/21] typo --- R/get_variances.R | 2 +- man/get_variance.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/get_variances.R b/R/get_variances.R index e51a02429..1ed61919d 100644 --- a/R/get_variances.R +++ b/R/get_variances.R @@ -115,7 +115,7 @@ #' of the aforementioned classes. #' #' Extracting variance components for models with zero-inflation part is not -#' straightforward, because it is not definitly clear how the distribution-specific +#' straightforward, because it is not definitely clear how the distribution-specific #' variance should be calculated. Therefore, it is recommended to carefully #' inspect the results, and probably validate against other models, e.g. Bayesian #' models (although results may be only roughly comparable). diff --git a/man/get_variance.Rd b/man/get_variance.Rd index e21334c06..61193f50f 100644 --- a/man/get_variance.Rd +++ b/man/get_variance.Rd @@ -124,7 +124,7 @@ not fully implemented or tested, and therefore may not work for all models of the aforementioned classes. Extracting variance components for models with zero-inflation part is not -straightforward, because it is not definitly clear how the distribution-specific +straightforward, because it is not definitely clear how the distribution-specific variance should be calculated. Therefore, it is recommended to carefully inspect the results, and probably validate against other models, e.g. Bayesian models (although results may be only roughly comparable). From d9cce15890aa8cd34dd7b53f3ce5185549043dd2 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 16:06:56 +0200 Subject: [PATCH 14/21] styler --- R/compute_variances.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index cbfdd8bfb..cabd7e636 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -863,7 +863,6 @@ # Tweedie models ------------------------------------------------------ # --------------------------------------------------------------------- dispersion_param <- .variance_family_tweedie(model, mu, sig) - } else if (identical(model_component, "full") && (faminfo$is_zero_inflated || faminfo$is_hurdle)) { # Zero Inflated models ------------------------------------------------ # --------------------------------------------------------------------- @@ -896,7 +895,6 @@ # ----------- sig ) - } else { # All other models ------------------------------------------------ # ----------------------------------------------------------------- From 868f0217f0cfd3fe51dd437256e2e66a0a125189 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 16:35:53 +0200 Subject: [PATCH 15/21] add test --- R/compute_variances.R | 26 +++++++++++---- .../test-r2_nakagawa_truncated_poisson.R | 32 +++++++++++++++++++ 2 files changed, 51 insertions(+), 7 deletions(-) create mode 100644 tests/testthat/test-r2_nakagawa_truncated_poisson.R diff --git a/R/compute_variances.R b/R/compute_variances.R index cabd7e636..244501481 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -22,6 +22,12 @@ # check argument approx_method <- match.arg(approximation, c("lognormal", "delta", "trigamma", "observation_level")) + # check whether R2 should be calculated for the full model, or the + # conditional model only + if (is.null(model_component) || model_component %in% c("zi", "zero_inflated")) { + model_component <- "full" + } + # sanity checks - distribution supported? if (any(faminfo$family == "truncated_nbinom1")) { if (verbose) { @@ -33,10 +39,12 @@ return(NA) } - # check whether R2 should be calculated for the full model, or the - # conditional model only - if (is.null(model_component) || model_component %in% c("zi", "zero_inflated")) { - model_component <- "full" + # sanity checks - distribution supported with requested component? + if ((any(startsWith(faminfo$family, "truncated_")) || any(startsWith(faminfo$family, "hurdle"))) && model_component != "full") { # nolint + if (verbose) { + format_warning("Truncated or hurdle families are only supported for `model_component = \"full\"`.") + } + return(NA) } # zero-inflated model, but not conditioning on full model? @@ -855,6 +863,8 @@ # "faminfo$family", but some also in "faminfo$link_function" (with a different # family). So we have to check for both. # Same applies to zero-inflated models, that's why we better double-check here. + # If variances for conditional model only are requested, we check the same + # families below again. # ------------------------------------------------------------------ cvsquared <- tryCatch( @@ -876,7 +886,7 @@ # hurdle-poisson ---- # ------------------- `hurdle poisson` = , - truncated_poisson = stats::family(model)$variance(sig), + truncated_poisson = stats::family(model)$variance(mu), # (zero-inflated) negative binomial ---- # -------------------------------------- @@ -896,13 +906,14 @@ sig ) } else { - # All other models ------------------------------------------------ + # All other families (or conditional model only) ------------------ # ----------------------------------------------------------------- dispersion_param <- switch(faminfo$family, # (generalized, compoised) poisson ---- # ---------------------------- poisson = 1, + `zero-inflated poisson` = 1, genpois = .variance_family_nbinom(model, mu, sig, faminfo), # Gamma, exponential ---- @@ -916,7 +927,8 @@ nbinom2 = , quasipoisson = , negbinomial = , - `negative binomial` = sig, + `negative binomial` = , + `zero-inflated negative binomial` = sig, # beta-alike ---- # --------------- diff --git a/tests/testthat/test-r2_nakagawa_truncated_poisson.R b/tests/testthat/test-r2_nakagawa_truncated_poisson.R new file mode 100644 index 000000000..f6b3f29d3 --- /dev/null +++ b/tests/testthat/test-r2_nakagawa_truncated_poisson.R @@ -0,0 +1,32 @@ +skip_on_cran() + +skip_if_not_installed("glmmTMB") +skip_if_not_installed("MuMIn") +skip_if_not_installed("performance") + + +# ============================================================================== +# neg-binomial1 zero-inflated mixed models, glmmTMB +# ============================================================================== + +test_that("glmmTMB, truncated_poisson", { + data(Salamanders, package = "glmmTMB") + + m <- glmmTMB::glmmTMB(count ~ spp + mined + (1 | site), + ziformula = ~ spp + mined, + family = glmmTMB::truncated_poisson(), data = glmmTMB::Salamanders + ) + + # truncated only works for full model + expect_warning(performance::r2_nakagawa(m, model_component = "conditional")) + + # full model + out <- performance::r2_nakagawa(m) + expect_equal(out$R2_conditional, 0.6116589, tolerance = 1e-4, ignore_attr = TRUE) + expect_equal(out$R2_marginal, 0.5411228, tolerance = 1e-4, ignore_attr = TRUE) + + out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m)) + # matches delta values + expect_equal(out1[1, "R2m"], out$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4) + expect_equal(out1[1, "R2c"], out$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4) +}) From d2d1a8c101ccf636287bac715c16014cb8e804e8 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 16:44:28 +0200 Subject: [PATCH 16/21] minor --- NEWS.md | 3 +++ R/format_message.R | 21 ++++++++++++--------- man/format_message.Rd | 11 ++++++++--- 3 files changed, 23 insertions(+), 12 deletions(-) diff --git a/NEWS.md b/NEWS.md index b1788dac8..7d441feb1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -26,6 +26,9 @@ specify for which component variances should be returned. By default, both the conditional and the zero-inflation component are taken into account. +* `format_alert()` and `format_warning()` get an `immediate` argument, to + ouput warnings immediately. + ## Bug fixes * `null_model()` now correctly handles zero-inflated models from package diff --git a/R/format_message.R b/R/format_message.R index 41cdeb146..a89015f8c 100644 --- a/R/format_message.R +++ b/R/format_message.R @@ -98,10 +98,13 @@ format_message <- function(string, #' @param type Type of exception alert to raise. #' Can be `"message"` for `message()`, `"warning"` for `warning()`, #' or `"error"` for `stop()`. -#' @param call. Logical. Indicating if the call should be included in the the +#' @param call Logical. Indicating if the call should be included in the the #' error message. This is usually confusing for users when the function #' producing the warning or error is deep within another function, so the #' default is `FALSE`. +#' @param immediate Logical. Indicating if the *warning* should be printed +#' immediately. Only applies to `format_warning()` or `format_alert()` with +#' `type = "warning"`. The default is `FALSE`. #' #' @examplesIf identical(Sys.getenv("NOT_CRAN"), "true") #' # message @@ -121,7 +124,8 @@ format_alert <- function(string, line_length = 0.9 * getOption("width", 80), indent = " ", type = "message", - call. = FALSE) { + call = FALSE, + immediate = FALSE) { type <- match.arg(type, choices = c("message", "warning", "error")) if (type == "message") { message(format_message( @@ -132,20 +136,20 @@ format_alert <- function(string, warning(format_message( string = string, ..., line_length = line_length, indent = indent - ), call. = call.) + ), call. = call, immediate. = immediate) } else { stop(format_message( string = string, ..., line_length = line_length, indent = indent - ), call. = call.) + ), call. = call) } } #' @name format_warning #' @rdname format_message #' @export -format_warning <- function(...) { - format_alert(..., type = "warning") +format_warning <- function(..., immediate = FALSE) { + format_alert(..., type = "warning", immediate. = immediate) } #' @name format_error @@ -252,9 +256,8 @@ format_error <- function(...) { .find_tokens <- function(string) { tokens <- c("{.b ", "{.i ", "{.url ", "{.pkg ") matches <- vapply(tokens, grepl, TRUE, string, fixed = TRUE) - if (any(matches)) { - matches - } else { + if (!any(matches)) { return(NULL) } + matches } diff --git a/man/format_message.Rd b/man/format_message.Rd index 79b5b9d8d..9bd51a991 100644 --- a/man/format_message.Rd +++ b/man/format_message.Rd @@ -20,10 +20,11 @@ format_alert( line_length = 0.9 * getOption("width", 80), indent = " ", type = "message", - call. = FALSE + call = FALSE, + immediate = FALSE ) -format_warning(...) +format_warning(..., immediate = FALSE) format_error(...) } @@ -44,10 +45,14 @@ line, two white space characters are inserted.} Can be \code{"message"} for \code{message()}, \code{"warning"} for \code{warning()}, or \code{"error"} for \code{stop()}.} -\item{call.}{Logical. Indicating if the call should be included in the the +\item{call}{Logical. Indicating if the call should be included in the the error message. This is usually confusing for users when the function producing the warning or error is deep within another function, so the default is \code{FALSE}.} + +\item{immediate}{Logical. Indicating if the \emph{warning} should be printed +immediately. Only applies to \code{format_warning()} or \code{format_alert()} with +\code{type = "warning"}. The default is \code{FALSE}.} } \value{ For \code{format_message()}, a formatted string. From 859fe25fb88c2b60479531a594086245ff5d4614 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 17:08:59 +0200 Subject: [PATCH 17/21] typo --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 7d441feb1..8532e6871 100644 --- a/NEWS.md +++ b/NEWS.md @@ -27,7 +27,7 @@ the conditional and the zero-inflation component are taken into account. * `format_alert()` and `format_warning()` get an `immediate` argument, to - ouput warnings immediately. + output warnings immediately. ## Bug fixes From 0b9cf95e7c8515dcdfb820ef646074e46b1f8d7f Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 18 Jun 2024 23:45:28 +0200 Subject: [PATCH 18/21] lintr --- R/colour_tools.R | 3 +- R/find_interactions.R | 8 ++--- R/find_offset.R | 26 +++++++-------- R/find_parameters_bayesian.R | 40 +++++++++++------------ R/find_parameters_emmeans.R | 8 ++--- R/format_rope.R | 6 ++-- R/format_value.R | 62 +++++++++++++++++------------------- R/get_varcov_sandwich.R | 6 ++-- R/is_nested_models.R | 10 +++--- 9 files changed, 83 insertions(+), 86 deletions(-) diff --git a/R/colour_tools.R b/R/colour_tools.R index f04e775dd..2d46f2165 100644 --- a/R/colour_tools.R +++ b/R/colour_tools.R @@ -3,7 +3,8 @@ return(FALSE) } - if ((cols <- Sys.getenv("RSTUDIO_CONSOLE_COLOR", "")) != "" && !is.na(as.numeric(cols))) { + cols <- Sys.getenv("RSTUDIO_CONSOLE_COLOR", "") + if (cols != "" && !is.na(as.numeric(cols))) { return(TRUE) } diff --git a/R/find_interactions.R b/R/find_interactions.R index a1fa1b71e..582ccfc40 100644 --- a/R/find_interactions.R +++ b/R/find_interactions.R @@ -71,13 +71,13 @@ find_interactions <- function(x, if (is.null(f)) { return(NULL) } - terms <- labels(stats::terms(f)) + model_terms <- labels(stats::terms(f)) if (main_effects) { - terms + model_terms } else { - interaction_terms <- grepl(":", terms, fixed = TRUE) + interaction_terms <- grepl(":", model_terms, fixed = TRUE) if (any(interaction_terms)) { - terms[interaction_terms] + model_terms[interaction_terms] } else { NULL } diff --git a/R/find_offset.R b/R/find_offset.R index 40b289cc0..8c24b0029 100644 --- a/R/find_offset.R +++ b/R/find_offset.R @@ -25,38 +25,38 @@ #' find_offset(m2) #' @export find_offset <- function(x) { - terms <- .safe( + model_terms <- .safe( as.character(attributes(stats::terms(find_formula(x)[[1]]))$variables), find_terms(x) ) - offset <- NULL + model_offset <- NULL - offcol <- grep("offset(", terms, fixed = TRUE) + offcol <- grep("offset(", model_terms, fixed = TRUE) if (length(offcol)) { - offset <- terms[offcol] + model_offset <- model_terms[offcol] } model_call <- get_call(x) - if (is.null(offset) && object_has_names(model_call, "offset")) { - offset <- safe_deparse(model_call$offset) + if (is.null(model_offset) && object_has_names(model_call, "offset")) { + model_offset <- safe_deparse(model_call$offset) } # fixest sometimes returns a weird macro syntax instead of the real offset # if we have to implement too many model-specific workarounds, it may eventually be worth it to do S3 # VAB: no test because I can only replicate in a weird {etwfe} example if (inherits(x, "fixest")) { - if (is.null(offset) || startsWith(offset, "..")) { - offset <- x[["model_info"]][["offset"]] + if (is.null(model_offset) || startsWith(model_offset, "..")) { + model_offset <- x[["model_info"]][["offset"]] } - offset <- sub("^~", "", offset) + model_offset <- sub("^~", "", model_offset) } - offset <- clean_names(offset) + model_offset <- clean_names(model_offset) # sometimes we get an empty list (e.g., fixest with iris dataset) - if (length(offset)) { - return(offset) + if (length(model_offset)) { + model_offset } else { - return(NULL) + NULL } } diff --git a/R/find_parameters_bayesian.R b/R/find_parameters_bayesian.R index 684bad68a..cd770d152 100644 --- a/R/find_parameters_bayesian.R +++ b/R/find_parameters_bayesian.R @@ -187,7 +187,7 @@ find_parameters.bamlss <- function(x, ignore <- grepl("(\\.alpha|logLik|\\.accepted|\\.edf)$", cn) cond <- cn[grepl("^(mu\\.p\\.|pi\\.p\\.)", cn) & !ignore] - sigma <- cn[startsWith(cn, "sigma.p.") & !ignore] + aux <- cn[startsWith(cn, "sigma.p.") & !ignore] smooth_terms <- cn[grepl("^mu\\.s\\.(.*)(\\.tau\\d+|\\.edf)$", cn)] alpha <- cn[endsWith(cn, ".alpha")] @@ -195,7 +195,7 @@ find_parameters.bamlss <- function(x, l <- compact_list(list( conditional = cond, smooth_terms = smooth_terms, - sigma = sigma, + sigma = aux, alpha = alpha )[elements]) @@ -244,10 +244,10 @@ find_parameters.brmsfit <- function(x, car_struc <- fe[fe %in% c("car", "sdcar")] smooth_terms <- fe[startsWith(fe, "sds_")] priors <- fe[startsWith(fe, "prior_")] - sigma <- fe[startsWith(fe, "sigma_") | grepl("sigma", fe, fixed = TRUE)] + sigma_param <- fe[startsWith(fe, "sigma_") | grepl("sigma", fe, fixed = TRUE)] randsigma <- fe[grepl("^r_(.*__sigma)", fe, perl = TRUE)] - beta <- fe[grepl("beta", fe, fixed = TRUE)] - randbeta <- fe[grepl("^r_(.*__beta)", fe, perl = TRUE)] + fixed_beta <- fe[grepl("beta", fe, fixed = TRUE)] + rand_beta <- fe[grepl("^r_(.*__beta)", fe, perl = TRUE)] mix <- fe[grepl("mix", fe, fixed = TRUE)] shiftprop <- fe[grepl("shiftprop", fe, fixed = TRUE)] dispersion <- fe[grepl("dispersion", fe, fixed = TRUE)] @@ -256,8 +256,8 @@ find_parameters.brmsfit <- function(x, # if auxiliary is modelled directly, we need to remove duplicates here # e.g. "b_sigma..." is in "cond" and in "sigma" now, we just need it in "cond". - sigma <- setdiff(sigma, c(cond, rand, rand_sd, rand_cor, randsigma, car_struc, "prior_sigma")) - beta <- setdiff(beta, c(cond, rand, rand_sd, randbeta, rand_cor, car_struc)) + sigma_param <- setdiff(sigma_param, c(cond, rand, rand_sd, rand_cor, randsigma, car_struc, "prior_sigma")) + fixed_beta <- setdiff(fixed_beta, c(cond, rand, rand_sd, rand_beta, rand_cor, car_struc)) auxiliary <- setdiff(auxiliary, c(cond, rand, rand_sd, rand_cor, car_struc)) l <- compact_list(list( @@ -267,10 +267,10 @@ find_parameters.brmsfit <- function(x, zero_inflated_random = c(randzi, randzi_sd, randzi_cor), simplex = simo, smooth_terms = smooth_terms, - sigma = sigma, + sigma = sigma_param, sigma_random = randsigma, - beta = beta, - beta_random = randbeta, + beta = fixed_beta, + beta_random = rand_beta, dispersion = dispersion, mix = mix, shiftprop = shiftprop, @@ -320,15 +320,15 @@ find_parameters.brmsfit <- function(x, } if (object_has_names(l, "sigma")) { - sigma <- l$sigma[grepl(sprintf("^sigma_\\Q%s\\E$", i), l$sigma)] + sigma_param <- l$sigma[grepl(sprintf("^sigma_\\Q%s\\E$", i), l$sigma)] } else { - sigma <- NULL + sigma_param <- NULL } if (object_has_names(l, "beta")) { - beta <- l$beta[grepl(sprintf("^beta_\\Q%s\\E$", i), l$sigma)] + fixed_beta <- l$beta[grepl(sprintf("^beta_\\Q%s\\E$", i), l$beta)] } else { - beta <- NULL + fixed_beta <- NULL } if (object_has_names(l, "dispersion")) { @@ -368,8 +368,8 @@ find_parameters.brmsfit <- function(x, zero_inflated_random = zero_inflated_random, simplex = simplex, smooth_terms = smooth_terms, - sigma = sigma, - beta = beta, + sigma = sigma_param, + beta = fixed_beta, dispersion = dispersion, mix = mix, priors = priors, @@ -441,7 +441,7 @@ find_parameters.stanreg <- function(x, rand <- fe[startsWith(fe, "b[")] rand_sd <- fe[startsWith(fe, "Sigma[")] smooth_terms <- fe[startsWith(fe, "smooth_sd")] - sigma <- fe[grepl("sigma", fe, fixed = TRUE)] + sigma_param <- fe[grepl("sigma", fe, fixed = TRUE)] auxiliary <- fe[grepl("(shape|phi|precision)", fe)] # remove auxiliary from conditional @@ -451,7 +451,7 @@ find_parameters.stanreg <- function(x, conditional = cond, random = c(rand, rand_sd), smooth_terms = smooth_terms, - sigma = sigma, + sigma = sigma_param, auxiliary = auxiliary )) @@ -499,7 +499,7 @@ find_parameters.stanmvreg <- function(x, rand <- fe[startsWith(fe, "b[")] rand_sd <- fe[startsWith(fe, "Sigma[")] smooth_terms <- fe[startsWith(fe, "smooth_sd")] - sigma <- fe[endsWith(fe, "|sigma") & .grep_non_smoothers(fe)] + sigma_param <- fe[endsWith(fe, "|sigma") & .grep_non_smoothers(fe)] auxiliary <- fe[grepl("(shape|phi|precision)", fe)] # remove auxiliary from conditional @@ -509,7 +509,7 @@ find_parameters.stanmvreg <- function(x, conditional = cond, random = c(rand, rand_sd), smooth_terms = smooth_terms, - sigma = sigma, + sigma = sigma_param, auxiliary = auxiliary )) diff --git a/R/find_parameters_emmeans.R b/R/find_parameters_emmeans.R index fde61d8e1..480e4c6c6 100644 --- a/R/find_parameters_emmeans.R +++ b/R/find_parameters_emmeans.R @@ -30,12 +30,10 @@ find_parameters.emmGrid <- function(x, flatten = FALSE, merge_parameters = FALSE } else { out <- stats::setNames(list(col_names), unique(.classify_emmeans(x))) } + } else if ("Component" %in% colnames(params)) { + out <- lapply(split(params, params$Component), function(i) i[[1]]) } else { - if ("Component" %in% colnames(params)) { - out <- lapply(split(params, params$Component), function(i) i[[1]]) - } else { - out <- stats::setNames(list(params[[1]]), unique(.classify_emmeans(x))) - } + out <- stats::setNames(list(params[[1]]), unique(.classify_emmeans(x))) } if (flatten) { diff --git a/R/format_rope.R b/R/format_rope.R index 799759eae..ae63ce4d1 100644 --- a/R/format_rope.R +++ b/R/format_rope.R @@ -11,15 +11,15 @@ #' format_rope(c(0.02, 0.12, 0.357, 0), name = NULL) #' @export format_rope <- function(rope_percentage, name = "in ROPE", digits = 2) { - text <- ifelse(rope_percentage == 0, "0%", + out <- ifelse(rope_percentage == 0, "0%", ifelse(rope_percentage == 1, "100%", # nolint format_value(rope_percentage, digits = digits, as_percent = TRUE) ) ) if (!is.null(name)) { - text <- paste(text, name) + out <- paste(out, name) } - text + out } diff --git a/R/format_value.R b/R/format_value.R index 06abb9300..2f9a48481 100644 --- a/R/format_value.R +++ b/R/format_value.R @@ -246,41 +246,39 @@ format_percent <- function(x, ...) { ) ) } + } else if (is.character(digits) && grepl("scientific", digits, fixed = TRUE)) { + digits <- tryCatch( + expr = { + as.numeric(gsub("scientific", "", digits, fixed = TRUE)) + }, + error = function(e) { + 5 + } + ) + if (is.na(digits)) digits <- 5 + x <- sprintf("%.*e", digits, x) + } else if (is.character(digits) && grepl("signif", digits, fixed = TRUE)) { + digits <- tryCatch( + expr = { + as.numeric(gsub("signif", "", digits, fixed = TRUE)) + }, + error = function(e) { + NA + } + ) + if (is.na(digits)) digits <- 3 + x <- as.character(signif(x, digits)) } else { - if (is.character(digits) && grepl("scientific", digits, fixed = TRUE)) { - digits <- tryCatch( - expr = { - as.numeric(gsub("scientific", "", digits, fixed = TRUE)) - }, - error = function(e) { - 5 - } - ) - if (is.na(digits)) digits <- 5 - x <- sprintf("%.*e", digits, x) - } else if (is.character(digits) && grepl("signif", digits, fixed = TRUE)) { - digits <- tryCatch( - expr = { - as.numeric(gsub("signif", "", digits, fixed = TRUE)) - }, - error = function(e) { - NA - } - ) - if (is.na(digits)) digits <- 3 - x <- as.character(signif(x, digits)) + need_sci <- (abs(x) >= 1e+5 | (log10(abs(x)) < -digits)) & x != 0 + if (.zap_small) { + x <- ifelse(is.na(x), .missing, sprintf("%.*f", digits, x)) } else { - need_sci <- (abs(x) >= 1e+5 | (log10(abs(x)) < -digits)) & x != 0 - if (.zap_small) { - x <- ifelse(is.na(x), .missing, sprintf("%.*f", digits, x)) - } else { - x <- ifelse(is.na(x), .missing, - ifelse(need_sci, # nolint - sprintf("%.*e", digits, x), - sprintf("%.*f", digits, x) - ) + x <- ifelse(is.na(x), .missing, + ifelse(need_sci, # nolint + sprintf("%.*e", digits, x), + sprintf("%.*f", digits, x) ) - } + ) } } } else if (anyNA(x)) { diff --git a/R/get_varcov_sandwich.R b/R/get_varcov_sandwich.R index 9c06d7404..53480df59 100644 --- a/R/get_varcov_sandwich.R +++ b/R/get_varcov_sandwich.R @@ -45,11 +45,11 @@ # vcov_fun is a function if (is.function(vcov_fun)) { if (is.null(vcov_args) || !is.list(vcov_args)) { - args <- list(x) + my_args <- list(x) } else { - args <- c(list(x), vcov_args) + my_args <- c(list(x), vcov_args) } - .vcov <- do.call("vcov_fun", args) + .vcov <- do.call("vcov_fun", my_args) return(.vcov) } diff --git a/R/is_nested_models.R b/R/is_nested_models.R index 03e0eaa49..f5bd76c2c 100644 --- a/R/is_nested_models.R +++ b/R/is_nested_models.R @@ -26,14 +26,14 @@ #' is_nested_models(m1, m2, m3) #' @export is_nested_models <- function(...) { - objects <- list(...) - object_names <- match.call(expand.dots = FALSE)$`...` + model_objects <- list(...) + object_names <- match.call(expand.dots = FALSE)[["..."]] - if (!all(vapply(objects, is_regression_model, TRUE))) { + if (!all(vapply(model_objects, is_regression_model, TRUE))) { format_error("All models must be valid regression model objects.") } - names(objects) <- object_names - info <- ellipsis_info.ListRegressions(objects) + names(model_objects) <- object_names + info <- ellipsis_info.ListRegressions(model_objects) out <- isTRUE(attributes(info)$is_nested) attr(out, "is_nested_increasing") <- attributes(info)$is_nested_increasing From bf1e5dd621263da68f2c91f14d68ac4958fea92f Mon Sep 17 00:00:00 2001 From: Daniel Date: Wed, 19 Jun 2024 00:10:25 +0200 Subject: [PATCH 19/21] refactor --- R/compute_variances.R | 31 ++++++++----------------------- 1 file changed, 8 insertions(+), 23 deletions(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index 244501481..a8df89fe3 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -698,9 +698,9 @@ ), .badlink(faminfo$link_function, faminfo$family, verbose = verbose) ) - } else if (faminfo$family == "beta") { - # Beta ---- - # ---------- + } else if (faminfo$family == "beta" || faminfo$is_orderedbeta) { + # (Ordered) Beta ---- + # -------------------- resid.variance <- switch(faminfo$link_function, logit = .variance_distributional( @@ -716,24 +716,6 @@ ), .badlink(faminfo$link_function, faminfo$family, verbose = verbose) ) - } else if (faminfo$is_orderedbeta) { - # Ordered Beta ---- - # ------------------ - - resid.variance <- switch(faminfo$link_function, - logit = .variance_distributional( - model, - faminfo = faminfo, - sig = sig, - model_null = model_null, - revar_null = revar_null, - approx_method = approx_method, - name = name, - model_component = model_component, - verbose = verbose - ), - .badlink(faminfo$link_function, faminfo$family, verbose = verbose) - ) } else { resid.variance <- sig } @@ -883,8 +865,10 @@ poisson = , `zero-inflated poisson` = .variance_family_poisson(model, mu, faminfo), - # hurdle-poisson ---- - # ------------------- + # hurdle-poisson/Gamma ---- + # ------------------------- + gamma = , + Gamma = , `hurdle poisson` = , truncated_poisson = stats::family(model)$variance(mu), @@ -918,6 +902,7 @@ # Gamma, exponential ---- # ----------------------- + gamma = , Gamma = stats::family(model)$variance(sig), # negative binomial ---- From bb9a6057249c11e0e0e9283922792abe483dd032 Mon Sep 17 00:00:00 2001 From: Daniel Date: Wed, 19 Jun 2024 08:50:25 +0200 Subject: [PATCH 20/21] model_info(lognormal) --- R/utils_model_info.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/utils_model_info.R b/R/utils_model_info.R index 7bf427308..dd384bcaa 100644 --- a/R/utils_model_info.R +++ b/R/utils_model_info.R @@ -345,7 +345,7 @@ dirichlet_fam || is.ordinal || zero.inf || is.censored || is.survival || is_binomtest || is.categorical || hurdle || is.multinomial || is_chi2test || is_proptest || is_xtab) { linear_model <- FALSE - } else if (!(fitfam %in% c("Student's-t", "t Family", "gaussian", "Gaussian")) && !grepl("(\\st)$", fitfam)) { + } else if (!(fitfam %in% c("Student's-t", "t Family", "gaussian", "Gaussian", "lognormal")) && !grepl("(\\st)$", fitfam)) { linear_model <- FALSE } if (!linear_model && is.survival && fitfam == "gaussian") { From 9dd7e06735d2c7a44139869d044f6d0ef0ace16e Mon Sep 17 00:00:00 2001 From: Daniel Date: Wed, 19 Jun 2024 09:22:28 +0200 Subject: [PATCH 21/21] cleanup --- R/compute_variances.R | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/R/compute_variances.R b/R/compute_variances.R index a8df89fe3..3a49a6051 100644 --- a/R/compute_variances.R +++ b/R/compute_variances.R @@ -822,20 +822,22 @@ } else { # transform mu mu <- switch(faminfo$family, + # beta-alike beta = , betabinomial = , ordbeta = stats::plogis(mu), # for count models, Nakagawa et al. suggest this transformation poisson = , quasipoisson = , + genpois = , nbinom = , nbinom1 = , nbinom2 = , negbinomial = , tweedie = , `negative binomial` = exp(mu + 0.5 * as.vector(revar_null)), - link_inverse(model)(mu) ## TODO: check if this is better than "exp(mu)" - # exp(mu) + # all other models + link_inverse(model)(mu) ) } @@ -865,15 +867,9 @@ poisson = , `zero-inflated poisson` = .variance_family_poisson(model, mu, faminfo), - # hurdle-poisson/Gamma ---- - # ------------------------- - gamma = , - Gamma = , - `hurdle poisson` = , - truncated_poisson = stats::family(model)$variance(mu), - # (zero-inflated) negative binomial ---- # -------------------------------------- + genpois = , nbinom = , nbinom1 = , nbinom2 = , @@ -881,8 +877,12 @@ `negative binomial` = , `zero-inflated negative binomial` = .variance_family_nbinom(model, mu, sig, faminfo), - # hurdle negative binomial ---- - # ----------------------------- + # hurdle and Gamma ---- + # ------------------------- + gamma = , + Gamma = , + `hurdle poisson` = , + truncated_poisson = stats::family(model)$variance(mu), truncated_nbinom2 = stats::family(model)$variance(mu, sig), # others ---- @@ -894,11 +894,10 @@ # ----------------------------------------------------------------- dispersion_param <- switch(faminfo$family, - # (generalized, compoised) poisson ---- - # ---------------------------- + # (compoised) poisson ---- + # ------------------------ poisson = 1, `zero-inflated poisson` = 1, - genpois = .variance_family_nbinom(model, mu, sig, faminfo), # Gamma, exponential ---- # ----------------------- @@ -913,6 +912,7 @@ quasipoisson = , negbinomial = , `negative binomial` = , + genpois = , `zero-inflated negative binomial` = sig, # beta-alike ---- @@ -925,7 +925,7 @@ # betabinomial = .variance_family_beta(model, mu, sig), # betabinomial = stats::family(model)$variance(mu, sig), - # other distributions ---- + # tweed distributions ---- # ------------------------ tweedie = .variance_family_tweedie(model, mu, sig), @@ -953,6 +953,7 @@ nbinom1 = , nbinom2 = , negbinomial = , + genpois = , `negative binomial` = ((1 / mu) + (1 / dispersion_param))^-1, poisson = , quasipoisson = mu / dispersion_param, @@ -964,6 +965,7 @@ nbinom1 = , nbinom2 = , negbinomial = , + genpois = , `negative binomial` = (1 / mu) + (1 / dispersion_param), poisson = , quasipoisson = dispersion_param / mu,