From a7547a991304c28f8d9f506851ef48f3e66981a2 Mon Sep 17 00:00:00 2001 From: Daniel Date: Fri, 5 Apr 2024 15:20:52 +0200 Subject: [PATCH] Fix issue in compare_parameters() for blme (#960) --- DESCRIPTION | 2 +- NEWS.md | 4 + R/compare_parameters.R | 2 +- R/dof.R | 29 ++- R/extract_random_variances.R | 78 ++++--- R/methods_lme4.R | 35 ++- inst/WORDLIST | 6 + man/model_parameters.merMod.Rd | 36 ++- tests/testthat/_snaps/glmmTMB-profile_CI.md | 4 +- tests/testthat/_snaps/glmmTMB.md | 169 ++------------- tests/testthat/test-compare_parameters.R | 16 ++ tests/testthat/test-glmmTMB.R | 205 +++++++++++++++--- .../testthat/test-random_effects_ci-glmmTMB.R | 37 +++- 13 files changed, 390 insertions(+), 233 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index e70ec2421..dc59633cf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: parameters Title: Processing of Model Parameters -Version: 0.21.6.1 +Version: 0.21.6.2 Authors@R: c(person(given = "Daniel", family = "Lüdecke", diff --git a/NEWS.md b/NEWS.md index 161302af7..0ebc98645 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # parameters 0.21.7 +## Bug fixes + +* Fixes issue in `compare_parameters()` for models from package *blme*. + # parameters 0.21.6 ## New supported models diff --git a/R/compare_parameters.R b/R/compare_parameters.R index 2c6fa7943..05305e244 100644 --- a/R/compare_parameters.R +++ b/R/compare_parameters.R @@ -169,7 +169,7 @@ compare_parameters <- function(..., dat <- model } else { # set default-ci_type for Bayesian models - if (.is_bayesian_model(model) && !ci_method %in% c("hdi", "quantile", "ci", "eti", "si", "bci", "bcai")) { + if (.is_bayesian_model(model, exclude = c("bmerMod", "bayesx", "blmerMod", "bglmerMod")) && !ci_method %in% c("hdi", "quantile", "ci", "eti", "si", "bci", "bcai")) { # nolint ci_method_tmp <- "eti" } else { ci_method_tmp <- ci_method diff --git a/R/dof.R b/R/dof.R index 502714d34..f0848eb2c 100644 --- a/R/dof.R +++ b/R/dof.R @@ -188,7 +188,7 @@ dof <- degrees_of_freedom #' @keywords internal .degrees_of_freedom_residual <- function(model, verbose = TRUE) { - if (.is_bayesian_model(model) && !inherits(model, c("bayesx", "blmerMod", "bglmerMod"))) { + if (.is_bayesian_model(model, exclude = c("bmerMod", "bayesx", "blmerMod", "bglmerMod"))) { model <- bayestestR::bayesian_as_frequentist(model) } @@ -277,7 +277,7 @@ dof <- degrees_of_freedom return(TRUE) } else { if (verbose) { - insight::format_alert(sprintf("`%s` must be one of \"wald\", \"residual\" or \"profile\". Using \"wald\" now.", type)) + insight::format_alert(sprintf("`%s` must be one of \"wald\", \"residual\" or \"profile\". Using \"wald\" now.", type)) # nolint } return(FALSE) } @@ -289,7 +289,7 @@ dof <- degrees_of_freedom return(TRUE) } else { if (verbose) { - insight::format_alert(sprintf("`%s` must be one of \"wald\", \"normal\" or \"boot\". Using \"wald\" now.", type)) + insight::format_alert(sprintf("`%s` must be one of \"wald\", \"normal\" or \"boot\". Using \"wald\" now.", type)) # nolint } return(FALSE) } @@ -298,24 +298,24 @@ dof <- degrees_of_freedom info <- insight::model_info(model, verbose = FALSE) if (!is.null(info) && isFALSE(info$is_mixed) && method == "boot") { if (verbose) { - insight::format_alert(sprintf("`%s=boot` only works for mixed models of class `merMod`. To bootstrap this model, use `bootstrap=TRUE, ci_method=\"bcai\"`.", type)) + insight::format_alert(sprintf("`%s=boot` only works for mixed models of class `merMod`. To bootstrap this model, use `bootstrap=TRUE, ci_method=\"bcai\"`.", type)) # nolint } return(TRUE) } if (is.null(info) || !info$is_mixed) { - if (!(method %in% c("analytical", "any", "fit", "betwithin", "nokr", "wald", "ml1", "profile", "boot", "uniroot", "residual", "normal"))) { + if (!(method %in% c("analytical", "any", "fit", "betwithin", "nokr", "wald", "ml1", "profile", "boot", "uniroot", "residual", "normal"))) { # nolint if (verbose) { - insight::format_alert(sprintf("`%s` must be one of \"residual\", \"wald\", \"normal\", \"profile\", \"boot\", \"uniroot\", \"betwithin\" or \"ml1\". Using \"wald\" now.", type)) + insight::format_alert(sprintf("`%s` must be one of \"residual\", \"wald\", \"normal\", \"profile\", \"boot\", \"uniroot\", \"betwithin\" or \"ml1\". Using \"wald\" now.", type)) # nolint } return(FALSE) } return(TRUE) } - if (!(method %in% c("analytical", "any", "fit", "satterthwaite", "betwithin", "kenward", "kr", "nokr", "wald", "ml1", "profile", "boot", "uniroot", "residual", "normal"))) { + if (!(method %in% c("analytical", "any", "fit", "satterthwaite", "betwithin", "kenward", "kr", "nokr", "wald", "ml1", "profile", "boot", "uniroot", "residual", "normal"))) { # nolint if (verbose) { - insight::format_alert(sprintf("`%s` must be one of \"residual\", \"wald\", \"normal\", \"profile\", \"boot\", \"uniroot\", \"kenward\", \"satterthwaite\", \"betwithin\" or \"ml1\". Using \"wald\" now.", type)) + insight::format_alert(sprintf("`%s` must be one of \"residual\", \"wald\", \"normal\", \"profile\", \"boot\", \"uniroot\", \"kenward\", \"satterthwaite\", \"betwithin\" or \"ml1\". Using \"wald\" now.", type)) # nolint } return(FALSE) } @@ -335,12 +335,17 @@ dof <- degrees_of_freedom # helper -.is_bayesian_model <- function(x) { - inherits(x, c( +.is_bayesian_model <- function(x, exclude = NULL) { + bayes_classes <- c( "brmsfit", "stanfit", "MCMCglmm", "stanreg", "stanmvreg", "bmerMod", "BFBayesFactor", "bamlss", "bayesx", "mcmc", "bcplm", "bayesQR", "BGGM", "meta_random", "meta_fixed", "meta_bma", "blavaan", - "blrm" - )) + "blrm", "blmerMod" + ) + # if exclude is not NULL, remove elements in exclude from bayes_class + if (!is.null(exclude)) { + bayes_classes <- bayes_classes[!bayes_classes %in% exclude] + } + inherits(x, bayes_classes) } diff --git a/R/extract_random_variances.R b/R/extract_random_variances.R index df9eef30f..620aabfa0 100644 --- a/R/extract_random_variances.R +++ b/R/extract_random_variances.R @@ -63,14 +63,14 @@ # check for errors if (is.null(out)) { if (isTRUE(verbose)) { - insight::format_warning("Something went wrong when calculating random effects parameters. Only showing model's fixed effects now. You may use `effects=\"fixed\"` to speed up the call to `model_parameters()`.") + insight::format_warning("Something went wrong when calculating random effects parameters. Only showing model's fixed effects now. You may use `effects=\"fixed\"` to speed up the call to `model_parameters()`.") # nolint } return(NULL) } out$Component <- "conditional" - if (insight::model_info(model, verbose = FALSE)$is_zero_inflated && !is.null(insight::find_random(model)$zero_inflated_random)) { + if (insight::model_info(model, verbose = FALSE)$is_zero_inflated && !is.null(insight::find_random(model)$zero_inflated_random)) { # nolint zi_var <- suppressWarnings( .extract_random_variances_helper( model, @@ -254,7 +254,9 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... sdcor = unname(values[-1]) ) })) - if (!is.null(corrs)) { + if (is.null(corrs)) { + out_cor <- NULL + } else { out_cor <- do.call(rbind, lapply(seq_along(from), function(i) { values <- corrs[from[i]:to[i]] .data_frame( @@ -264,8 +266,6 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... sdcor = unname(values[-1]) ) })) - } else { - out_cor <- NULL } } else { out_sd <- .data_frame( @@ -274,15 +274,15 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... var2 = NA_character_, sdcor = unname(stddevs) ) - if (!is.null(corrs)) { + if (is.null(corrs)) { + out_cor <- NULL + } else { out_cor <- .data_frame( grp = gsub("(.*) =(.*)", "\\1", attributes(x)$title), var1 = "(Intercept)", var2 = names(corrs), sdcor = unname(corrs) ) - } else { - out_cor <- NULL } } @@ -300,7 +300,14 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... # extract CI for random SD ------------------------ -.random_sd_ci <- function(model, out, ci_method, ci, ci_random, corr_param, sigma_param, component = NULL, verbose = FALSE) { +.random_sd_ci <- function(model, + out, + ci_method, + ci, ci_random, + corr_param, + sigma_param, + component = NULL, + verbose = FALSE) { ## TODO needs to be removed once MCM > 0.1.5 is on CRAN if (startsWith(insight::safe_deparse(insight::get_call(model)), "mcm_lmer")) { return(out) @@ -340,7 +347,13 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... if (!is.null(ci_method) && ci_method %in% c("profile", "boot")) { out <- tryCatch( { - var_ci <- as.data.frame(suppressWarnings(stats::confint(model, parm = "theta_", oldNames = FALSE, method = ci_method, level = ci))) + var_ci <- as.data.frame(suppressWarnings(stats::confint( + model, + parm = "theta_", + oldNames = FALSE, + method = ci_method, + level = ci + ))) colnames(var_ci) <- c("CI_low", "CI_high") rn <- row.names(var_ci) @@ -370,7 +383,8 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... if (isTRUE(verbose)) { insight::format_alert( "Cannot compute profiled standard errors and confidence intervals for random effects parameters.", - "Your model may suffer from singularity (see '?lme4::isSingular' and '?performance::check_singularity')." + "Your model may suffer from singularity (see '?lme4::isSingular' and '?performance::check_singularity').", + "You may try to impose a prior on the random effects parameters, e.g. using the {.pkg glmmTMB} package." ) } out @@ -445,7 +459,7 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... for (gr in setdiff(unique(out$Group), "Residual")) { rnd_slope_corr_grp <- rnd_slope_corr & out$Group == gr dummy <- gsub("Cor \\((.*)~(.*)\\)", "\\2.\\1", out$Parameter[rnd_slope_corr_grp]) - var_ci$Parameter[var_ci$Group == gr][match(dummy, var_ci$Parameter[var_ci$Group == gr])] <- out$Parameter[rnd_slope_corr_grp] + var_ci$Parameter[var_ci$Group == gr][match(dummy, var_ci$Parameter[var_ci$Group == gr])] <- out$Parameter[rnd_slope_corr_grp] # nolint } } @@ -484,10 +498,11 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... out$CI_high[!var_ci_corr_param] <- exp(log(coefs) + stats::qnorm(0.975) * delta_se) # warn if singular fit - if (isTRUE(verbose) && insight::check_if_installed("performance", quietly = TRUE) && isTRUE(performance::check_singularity(model))) { + if (isTRUE(verbose) && insight::check_if_installed("performance", quietly = TRUE) && isTRUE(performance::check_singularity(model))) { # nolint insight::format_alert( - "Your model may suffer from singularity (see see `?lme4::isSingular` and `?performance::check_singularity`).", - "Some of the standard errors and confidence intervals of the random effects parameters are probably not meaningful!" + "Your model may suffer from singularity (see see `?lme4::isSingular` and `?performance::check_singularity`).", # nolint + "Some of the standard errors and confidence intervals of the random effects parameters are probably not meaningful!", # nolint + "You may try to impose a prior on the random effects parameters, e.g. using the {.pkg glmmTMB} package." # nolint ) } out @@ -495,17 +510,18 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... error = function(e) { if (isTRUE(verbose)) { if (grepl("nAGQ of at least 1 is required", e$message, fixed = TRUE)) { - insight::format_alert("Argument `nAGQ` needs to be larger than 0 to compute confidence intervals for random effect parameters.") + insight::format_alert("Argument `nAGQ` needs to be larger than 0 to compute confidence intervals for random effect parameters.") # nolint } if (grepl("Multiple cluster variables detected.", e$message, fixed = TRUE)) { - insight::format_alert("Confidence intervals for random effect parameters are currently not supported for multiple grouping variables.") + insight::format_alert("Confidence intervals for random effect parameters are currently not supported for multiple grouping variables.") # nolint } if (grepl("exactly singular", e$message, fixed = TRUE) || grepl("computationally singular", e$message, fixed = TRUE) || grepl("Exact singular", e$message, fixed = TRUE)) { insight::format_alert( "Cannot compute standard errors and confidence intervals for random effects parameters.", - "Your model may suffer from singularity (see see `?lme4::isSingular` and `?performance::check_singularity`)." + "Your model may suffer from singularity (see see `?lme4::isSingular` and `?performance::check_singularity`).", # nolint + "You may try to impose a prior on the random effects parameters, e.g. using the {.pkg glmmTMB} package." # nolint ) } } @@ -513,7 +529,7 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... } ) } else if (isTRUE(verbose)) { - insight::format_alert("Package 'merDeriv' needs to be installed to compute confidence intervals for random effect parameters.") + insight::format_alert("Package 'merDeriv' needs to be installed to compute confidence intervals for random effect parameters.") # nolint } } } else if (inherits(model, "glmmTMB")) { @@ -553,8 +569,8 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... # add Group var_ci$Group <- NA if (length(group_factor) > 1) { - var_ci$Group[var_ci$Component == "conditional"] <- gsub(paste0("^", group_factor2, "\\.cond\\.(.*)"), "\\1", var_ci$Parameter[var_ci$Component == "conditional"]) - var_ci$Group[var_ci$Component == "zi"] <- gsub(paste0("^", group_factor2, "\\.zi\\.(.*)"), "\\1", var_ci$Parameter[var_ci$Component == "zi"]) + var_ci$Group[var_ci$Component == "conditional"] <- gsub(paste0("^", group_factor2, "\\.cond\\.(.*)"), "\\1", var_ci$Parameter[var_ci$Component == "conditional"]) # nolint + var_ci$Group[var_ci$Component == "zi"] <- gsub(paste0("^", group_factor2, "\\.zi\\.(.*)"), "\\1", var_ci$Parameter[var_ci$Component == "zi"]) # nolint } else { var_ci$Group <- group_factor # check if sigma was properly identified @@ -594,17 +610,19 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... # check results - warn user if (isTRUE(verbose)) { missing_ci <- any(is.na(var_ci$CI_low) | is.na(var_ci$CI_high)) - singular_fit <- insight::check_if_installed("performance", quietly = TRUE) & isTRUE(performance::check_singularity(model)) + singular_fit <- insight::check_if_installed("performance", quietly = TRUE) & isTRUE(performance::check_singularity(model)) # nolint if (singular_fit) { insight::format_alert( "Your model may suffer from singularity (see `?lme4::isSingular` and `?performance::check_singularity`).", - "Some of the confidence intervals of the random effects parameters are probably not meaningful!" + "Some of the confidence intervals of the random effects parameters are probably not meaningful!", + "You may try to impose a prior on the random effects parameters, e.g. using the {.pkg glmmTMB} package." # nolint ) } else if (missing_ci) { insight::format_alert( "Your model may suffer from singularity (see `?lme4::isSingular` and `?performance::check_singularity`).", - "Some of the confidence intervals of the random effects parameters could not be calculated or are probably not meaningful!" + "Some of the confidence intervals of the random effects parameters could not be calculated or are probably not meaningful!", # nolint + "You may try to impose a prior on the random effects parameters, e.g. using the {.pkg glmmTMB} package." # nolint ) } } @@ -695,13 +713,13 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... vc1 <- list(vc1) names(vc1) <- re_names[[1]] - attr(vc1, "sc") <- sqrt(insight::get_deviance(model, verbose = FALSE) / insight::get_df(model, type = "residual", verbose = FALSE)) + attr(vc1, "sc") <- sqrt(insight::get_deviance(model, verbose = FALSE) / insight::get_df(model, type = "residual", verbose = FALSE)) # nolint attr(vc1, "useSc") <- TRUE if (!is.null(vc2)) { vc2 <- list(vc2) names(vc2) <- re_names[[2]] - attr(vc2, "sc") <- sqrt(insight::get_deviance(model, verbose = FALSE) / insight::get_df(model, type = "residual", verbose = FALSE)) + attr(vc2, "sc") <- sqrt(insight::get_deviance(model, verbose = FALSE) / insight::get_df(model, type = "residual", verbose = FALSE)) # nolint attr(vc2, "useSc") <- FALSE } @@ -742,12 +760,12 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ... varcorr <- 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$cov[, 1, ])) - colnames(out) <- rownames(out) <- gsub("Intercept", "(Intercept)", rownames(element$cov), fixed = TRUE) - } else { + if (is.null(element$cov)) { out <- as.matrix(drop(element$sd[, 1])^2) colnames(out) <- rownames(out) <- gsub("Intercept", "(Intercept)", rownames(element$sd), fixed = TRUE) + } else { + out <- as.matrix(drop(element$cov[, 1, ])) + colnames(out) <- rownames(out) <- gsub("Intercept", "(Intercept)", rownames(element$cov), fixed = TRUE) } attr(out, "sttdev") <- element$sd[, 1] } else { diff --git a/R/methods_lme4.R b/R/methods_lme4.R index 6eb50c384..854edd29a 100644 --- a/R/methods_lme4.R +++ b/R/methods_lme4.R @@ -40,7 +40,7 @@ #' #' @inheritSection model_parameters Confidence intervals and approximation of degrees of freedom #' -#' @section Confidence intervals for random effect variances: +#' @section Confidence intervals for random effects variances: #' For models of class `merMod` and `glmmTMB`, confidence intervals for random #' effect variances can be calculated. #' @@ -61,6 +61,34 @@ #' - For models of class `glmmTMB`, confidence intervals for random effect #' variances always use a Wald t-distribution approximation. #' +#' @section Singular fits (random effects variances near zero): +#' If a model is "singular", this means that some dimensions of the +#' variance-covariance matrix have been estimated as exactly zero. This +#' often occurs for mixed models with complex random effects structures. +#' +#' There is no gold-standard about how to deal with singularity and which +#' random-effects specification to choose. One way is to fully go Bayesian +#' (with informative priors). Other proposals are listed in the documentation +#' of [`performance::check_singularity()`]. However, since version 1.1.9, the +#' **glmmTMB** package allows to use priors in a frequentist framework, too. One +#' recommendation is to use a Gamma prior (_Chung et al. 2013_). The mean may +#' vary from 1 to very large values (like `1e8`), and the shape parameter should +#' be set to a value of 2.5. You can then `update()` your model with the specified +#' prior. In **glmmTMB**, the code would look like this: +#' +#' ``` +#' # "model" is an object of class gmmmTMB +#' prior <- data.frame( +#' prior = "gamma(1, 2.5)", # mean can be 1, but even 1e8 +#' class = "ranef" # for random effects +#' ) +#' model_with_priors <- update(model, priors = prior) +#' ``` +#' +#' Large values for the mean parameter of the Gamma prior have no large impact +#' on the random effects variances in terms of a "bias". Thus, if `1` doesn't +#' fix the singular fit, you can safely try larger values. +#' #' @section Dispersion parameters in *glmmTMB*: #' For some models from package **glmmTMB**, both the dispersion parameter and #' the residual variance from the random effects parameters are shown. Usually, @@ -82,6 +110,11 @@ #' use `effects = "fixed"`. There is also a [`plot()`-method](https://easystats.github.io/see/articles/parameters.html) #' implemented in the [**see**-package](https://easystats.github.io/see/). #' +#' @references +#' Chung Y, Rabe-Hesketh S, Dorie V, Gelman A, and Liu J. 2013. "A Nondegenerate +#' Penalized Likelihood Estimator for Variance Parameters in Multilevel Models." +#' Psychometrika 78 (4): 685–709. \doi{10.1007/s11336-013-9328-2} +#' #' @examplesIf require("lme4") && require("glmmTMB") #' library(parameters) #' data(mtcars) diff --git a/inst/WORDLIST b/inst/WORDLIST index 7637348ab..a7e977782 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -37,6 +37,7 @@ DirichletReg DoF DoFs Dom +Dorie Dupont EFA EGAnet @@ -64,6 +65,7 @@ HDI HEXACO HLM Heisig +Hesketh Heteroskedasticity Higgs Hinkley @@ -86,6 +88,7 @@ LMMs Lakens Laparra Lawley +Liu MADs MCMCglmm MLM @@ -103,6 +106,7 @@ NL Neter Neyman Nieto +Nondegenerate Nonresponse ORCID Olkin @@ -116,6 +120,7 @@ PloS Psychometrika REWB ROPE's +Rabe Rafi Rchoice Revelle @@ -174,6 +179,7 @@ behaviours betareg biserial blavaan +blme bmwiernik brglm brms diff --git a/man/model_parameters.merMod.Rd b/man/model_parameters.merMod.Rd index 71c2874fb..8b1b8e75c 100644 --- a/man/model_parameters.merMod.Rd +++ b/man/model_parameters.merMod.Rd @@ -344,7 +344,7 @@ If the calculation of random effects parameters takes too long, you may use \code{effects = "fixed"}. There is also a \href{https://easystats.github.io/see/articles/parameters.html}{\code{plot()}-method} implemented in the \href{https://easystats.github.io/see/}{\strong{see}-package}. } -\section{Confidence intervals for random effect variances}{ +\section{Confidence intervals for random effects variances}{ For models of class \code{merMod} and \code{glmmTMB}, confidence intervals for random effect variances can be calculated. @@ -366,6 +366,35 @@ variances always use a Wald t-distribution approximation. } } +\section{Singular fits (random effects variances near zero)}{ + +If a model is "singular", this means that some dimensions of the +variance-covariance matrix have been estimated as exactly zero. This +often occurs for mixed models with complex random effects structures. + +There is no gold-standard about how to deal with singularity and which +random-effects specification to choose. One way is to fully go Bayesian +(with informative priors). Other proposals are listed in the documentation +of \code{\link[performance:check_singularity]{performance::check_singularity()}}. However, since version 1.1.9, the +\strong{glmmTMB} package allows to use priors in a frequentist framework, too. One +recommendation is to use a Gamma prior (\emph{Chung et al. 2013}). The mean may +vary from 1 to very large values (like \code{1e8}), and the shape parameter should +be set to a value of 2.5. You can then \code{update()} your model with the specified +prior. In \strong{glmmTMB}, the code would look like this: + +\if{html}{\out{
}}\preformatted{# "model" is an object of class gmmmTMB +prior <- data.frame( + prior = "gamma(1, 2.5)", # mean can be 1, but even 1e8 + class = "ranef" # for random effects +) +model_with_priors <- update(model, priors = prior) +}\if{html}{\out{
}} + +Large values for the mean parameter of the Gamma prior have no large impact +on the random effects variances in terms of a "bias". Thus, if \code{1} doesn't +fix the singular fit, you can safely try larger values. +} + \section{Dispersion parameters in \emph{glmmTMB}}{ For some models from package \strong{glmmTMB}, both the dispersion parameter and @@ -569,6 +598,11 @@ model_parameters(model, bootstrap = TRUE, iterations = 50, verbose = FALSE) } \dontshow{\}) # examplesIf} } +\references{ +Chung Y, Rabe-Hesketh S, Dorie V, Gelman A, and Liu J. 2013. "A Nondegenerate +Penalized Likelihood Estimator for Variance Parameters in Multilevel Models." +Psychometrika 78 (4): 685–709. \doi{10.1007/s11336-013-9328-2} +} \seealso{ \code{\link[insight:standardize_names]{insight::standardize_names()}} to rename columns into a consistent, standardized naming scheme. diff --git a/tests/testthat/_snaps/glmmTMB-profile_CI.md b/tests/testthat/_snaps/glmmTMB-profile_CI.md index 0041f484f..b4499193a 100644 --- a/tests/testthat/_snaps/glmmTMB-profile_CI.md +++ b/tests/testthat/_snaps/glmmTMB-profile_CI.md @@ -7,7 +7,7 @@ Parameter | Coefficient | SE | 95% CI | z | p -------------------------------------------------------------------- - (Intercept) | 251.41 | 6.63 | [237.68, 265.13] | 37.91 | < .001 + (Intercept) | 251.40 | 6.63 | [237.68, 265.13] | 37.91 | < .001 Days | 10.47 | 1.50 | [ 7.36, 13.58] | 6.97 | < .001 # Random Effects @@ -34,7 +34,7 @@ Parameter | Coefficient | SE | 95% CI | z | p -------------------------------------------------------------------- - (Intercept) | 251.41 | 6.63 | [237.68, 265.13] | 37.91 | < .001 + (Intercept) | 251.40 | 6.63 | [237.68, 265.13] | 37.91 | < .001 Days | 10.47 | 1.50 | [ 7.36, 13.58] | 6.97 | < .001 # Random Effects diff --git a/tests/testthat/_snaps/glmmTMB.md b/tests/testthat/_snaps/glmmTMB.md index 3879d85b3..33db59721 100644 --- a/tests/testthat/_snaps/glmmTMB.md +++ b/tests/testthat/_snaps/glmmTMB.md @@ -10,127 +10,6 @@ [5] "child | -1.09 | 0.10 | [-1.28, -0.90] | -11.09 | < .001" [6] "camper [1] | 0.27 | 0.10 | [ 0.07, 0.47] | 2.70 | 0.007 " ---- - - Code - print(mp) - Output - # Random Effects - - Parameter | Coefficient | 95% CI - --------------------------------------------------------- - SD (Intercept: persons) | 3.41 | [ 1.68, 6.92] - SD (xb: persons) | 1.21 | [ 0.60, 2.44] - Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00] - Message - - Uncertainty intervals for random effect variances computed using a - Wald z-distribution approximation. - ---- - - Code - out[-6] - Output - [1] "# Fixed Effects (Zero-Inflation Component)" - [2] "" - [3] "Parameter | Log-Mean | SE | 95% CI | z | p" - [4] "-------------------------------------------------------------" - [5] "(Intercept) | 1.89 | 0.66 | [ 0.59, 3.19] | 2.85 | 0.004" - [6] "camper [1] | -0.17 | 0.39 | [-0.93, 0.59] | -0.44 | 0.660" - ---- - - Code - print(mp) - Output - # Random Effects (Zero-Inflation Component) - - Parameter | Coefficient | 95% CI - -------------------------------------------------------- - SD (Intercept: persons) | 2.74 | [1.17, 6.40] - SD (zg: persons) | 1.57 | [0.65, 3.81] - Cor (Intercept~zg: persons) | 1.00 | - Message - - Uncertainty intervals for random effect variances computed using a - Wald z-distribution approximation. - ---- - - Code - out[-5] - Output - [1] "# Fixed Effects" - [2] "" - [3] "Parameter | Log-Mean | SE | 95% CI | z | p" - [4] "----------------------------------------------------------------" - [5] "child | -1.09 | 0.10 | [-1.28, -0.90] | -11.09 | < .001" - [6] "camper [1] | 0.27 | 0.10 | [ 0.07, 0.47] | 2.70 | 0.007 " - [7] "" - [8] "# Random Effects" - [9] "" - [10] "Parameter | Coefficient | 95% CI" - [11] "---------------------------------------------------------" - [12] "SD (Intercept: persons) | 3.41 | [ 1.68, 6.92]" - [13] "SD (xb: persons) | 1.21 | [ 0.60, 2.44]" - [14] "Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00]" - ---- - - Code - out[-6] - Output - [1] "# Fixed Effects (Zero-Inflation Component)" - [2] "" - [3] "Parameter | Log-Mean | SE | 95% CI | z | p" - [4] "-------------------------------------------------------------" - [5] "(Intercept) | 1.89 | 0.66 | [ 0.59, 3.19] | 2.85 | 0.004" - [6] "camper [1] | -0.17 | 0.39 | [-0.93, 0.59] | -0.44 | 0.660" - [7] "" - [8] "# Random Effects (Zero-Inflation Component)" - [9] "" - [10] "Parameter | Coefficient | 95% CI" - [11] "--------------------------------------------------------" - [12] "SD (Intercept: persons) | 2.74 | [1.17, 6.40]" - [13] "SD (zg: persons) | 1.57 | [0.65, 3.81]" - [14] "Cor (Intercept~zg: persons) | 1.00 | " - ---- - - Code - out[-c(5, 14)] - Output - [1] "# Fixed Effects (Count Model)" - [2] "" - [3] "Parameter | Log-Mean | SE | 95% CI | z | p" - [4] "----------------------------------------------------------------" - [5] "child | -1.09 | 0.10 | [-1.28, -0.90] | -11.09 | < .001" - [6] "camper [1] | 0.27 | 0.10 | [ 0.07, 0.47] | 2.70 | 0.007 " - [7] "" - [8] "# Fixed Effects (Zero-Inflation Component)" - [9] "" - [10] "Parameter | Log-Odds | SE | 95% CI | z | p" - [11] "-------------------------------------------------------------" - [12] "(Intercept) | 1.89 | 0.66 | [ 0.59, 3.19] | 2.85 | 0.004" - [13] "camper [1] | -0.17 | 0.39 | [-0.93, 0.59] | -0.44 | 0.660" - [14] "" - [15] "# Random Effects Variances" - [16] "" - [17] "Parameter | Coefficient | 95% CI" - [18] "---------------------------------------------------------" - [19] "SD (Intercept: persons) | 3.41 | [ 1.68, 6.92]" - [20] "SD (xb: persons) | 1.21 | [ 0.60, 2.44]" - [21] "Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00]" - [22] "" - [23] "# Random Effects (Zero-Inflation Component)" - [24] "" - [25] "Parameter | Coefficient | 95% CI" - [26] "--------------------------------------------------------" - [27] "SD (Intercept: persons) | 2.74 | [1.17, 6.40]" - [28] "SD (zg: persons) | 1.57 | [0.65, 3.81]" - [29] "Cor (Intercept~zg: persons) | 1.00 | " - # print-model_parameters glmmTMB digits Code @@ -140,31 +19,31 @@ [2] "" [3] "Parameter | Log-Mean | SE | 95% CI | z | p" [4] "--------------------------------------------------------------------------" - [5] "child | -1.0875 | 0.0981 | [-1.27966, -0.89529] | -11.0903 | < .001" - [6] "camper [1] | 0.2723 | 0.1009 | [ 0.07462, 0.46998] | 2.6997 | 0.007 " + [5] "child | -1.0875 | 0.0981 | [-1.27967, -0.89528] | -11.0901 | < .001" + [6] "camper [1] | 0.2723 | 0.1009 | [ 0.07461, 0.46999] | 2.6997 | 0.007 " [7] "" [8] "# Fixed Effects (Zero-Inflation Component)" [9] "" [10] "Parameter | Log-Odds | SE | 95% CI | z | p" [11] "-----------------------------------------------------------------------" - [12] "(Intercept) | 1.8896 | 0.6642 | [ 0.58788, 3.19139] | 2.8451 | 0.004" - [13] "camper [1] | -0.1701 | 0.3868 | [-0.92809, 0.58796] | -0.4397 | 0.660" + [12] "(Intercept) | 1.8896 | 0.6642 | [ 0.58780, 3.19147] | 2.8449 | 0.004" + [13] "camper [1] | -0.1701 | 0.3869 | [-0.92836, 0.58822] | -0.4396 | 0.660" [14] "" [15] "# Random Effects Variances" [16] "" [17] "Parameter | Coefficient | 95% CI" [18] "---------------------------------------------------------------" - [19] "SD (Intercept: persons) | 3.4056 | [ 1.67675, 6.91715]" - [20] "SD (xb: persons) | 1.2132 | [ 0.60312, 2.44025]" + [19] "SD (Intercept: persons) | 3.4056 | [ 1.64567, 7.04777]" + [20] "SD (xb: persons) | 1.2132 | [ 0.59190, 2.48650]" [21] "Cor (Intercept~xb: persons) | -1.0000 | [-1.00000, 1.00000]" [22] "" [23] "# Random Effects (Zero-Inflation Component)" [24] "" - [25] "Parameter | Coefficient | 95% CI" - [26] "--------------------------------------------------------------" - [27] "SD (Intercept: persons) | 2.7358 | [1.16901, 6.40268]" - [28] "SD (zg: persons) | 1.5683 | [0.64628, 3.80586]" - [29] "Cor (Intercept~zg: persons) | 1.0000 | " + [25] "Parameter | Coefficient | 95% CI" + [26] "---------------------------------------------------------------" + [27] "SD (Intercept: persons) | 2.7358 | [ 1.16329, 6.43414]" + [28] "SD (zg: persons) | 1.5683 | [ 0.64246, 3.82852]" + [29] "Cor (Intercept~zg: persons) | 1.0000 | [-1.00000, 1.00000]" --- @@ -175,31 +54,31 @@ [2] "" [3] "Parameter | Log-Mean | SE | 95% CI | z | p" [4] "--------------------------------------------------------------------------" - [5] "child | -1.0875 | 0.0981 | [-1.27966, -0.89529] | -11.0903 | < .001" - [6] "camper [1] | 0.2723 | 0.1009 | [ 0.07462, 0.46998] | 2.6997 | 0.007 " + [5] "child | -1.0875 | 0.0981 | [-1.27967, -0.89528] | -11.0901 | < .001" + [6] "camper [1] | 0.2723 | 0.1009 | [ 0.07461, 0.46999] | 2.6997 | 0.007 " [7] "" [8] "# Fixed Effects (Zero-Inflation Component)" [9] "" [10] "Parameter | Log-Odds | SE | 95% CI | z | p" [11] "-----------------------------------------------------------------------" - [12] "(Intercept) | 1.8896 | 0.6642 | [ 0.58788, 3.19139] | 2.8451 | 0.004" - [13] "camper [1] | -0.1701 | 0.3868 | [-0.92809, 0.58796] | -0.4397 | 0.660" + [12] "(Intercept) | 1.8896 | 0.6642 | [ 0.58780, 3.19147] | 2.8449 | 0.004" + [13] "camper [1] | -0.1701 | 0.3869 | [-0.92836, 0.58822] | -0.4396 | 0.660" [14] "" [15] "# Random Effects Variances" [16] "" [17] "Parameter | Coefficient | 95% CI" [18] "---------------------------------------------------------------" - [19] "SD (Intercept: persons) | 3.4056 | [ 1.67675, 6.91715]" - [20] "SD (xb: persons) | 1.2132 | [ 0.60312, 2.44025]" + [19] "SD (Intercept: persons) | 3.4056 | [ 1.64567, 7.04777]" + [20] "SD (xb: persons) | 1.2132 | [ 0.59190, 2.48650]" [21] "Cor (Intercept~xb: persons) | -1.0000 | [-1.00000, 1.00000]" [22] "" [23] "# Random Effects (Zero-Inflation Component)" [24] "" - [25] "Parameter | Coefficient | 95% CI" - [26] "--------------------------------------------------------------" - [27] "SD (Intercept: persons) | 2.7358 | [1.16901, 6.40268]" - [28] "SD (zg: persons) | 1.5683 | [0.64628, 3.80586]" - [29] "Cor (Intercept~zg: persons) | 1.0000 | " + [25] "Parameter | Coefficient | 95% CI" + [26] "---------------------------------------------------------------" + [27] "SD (Intercept: persons) | 2.7358 | [ 1.16329, 6.43414]" + [28] "SD (zg: persons) | 1.5683 | [ 0.64246, 3.82852]" + [29] "Cor (Intercept~zg: persons) | 1.0000 | [-1.00000, 1.00000]" # print-model_parameters glmmTMB CI alignment @@ -219,7 +98,7 @@ ---------------------------------------------------------------- SD (Intercept: Session:Participant) | 0.69 | [0.40, 1.19] SD (Intercept: Participant) | 2.39 | [1.25, 4.57] - Message + Message Uncertainty intervals for random effect variances computed using a Wald z-distribution approximation. @@ -256,7 +135,7 @@ Parameter | Coefficient | 95% CI ---------------------------------------- (Intercept) | 2.06 | [1.30, 3.27] - Message + Message Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald z-distribution approximation. diff --git a/tests/testthat/test-compare_parameters.R b/tests/testthat/test-compare_parameters.R index fe261d8d9..ce299dc56 100644 --- a/tests/testthat/test-compare_parameters.R +++ b/tests/testthat/test-compare_parameters.R @@ -209,3 +209,19 @@ withr::with_options( }) } ) + +skip_on_cran() +skip_if_not_installed("blme") +skip_if_not_installed("glmmTMB") +skip_if_not_installed("lme4") + +test_that("compare_parameters, works with blmer and glmmTMB", { + data(sleepstudy, package = "lme4") + control <- lme4::lmerControl(check.conv.grad = "ignore") + fm1 <- blme::blmer(Reaction ~ Days + (0 + Days | Subject), sleepstudy, + control = control, + cov.prior = gamma + ) + fm2 <- glmmTMB::glmmTMB(Reaction ~ Days + (1 + Days | Subject), sleepstudy) + expect_silent(compare_parameters(fm1, fm2)) +}) diff --git a/tests/testthat/test-glmmTMB.R b/tests/testthat/test-glmmTMB.R index fcfb533db..a8d80751f 100644 --- a/tests/testthat/test-glmmTMB.R +++ b/tests/testthat/test-glmmTMB.R @@ -200,9 +200,9 @@ withr::with_options( test_that("model_parameters.mixed-random", { params <- model_parameters(m1, effects = "random", group_level = TRUE) - expect_equal(c(nrow(params), ncol(params)), c(8, 10)) - expect_identical( - colnames(params), + expect_identical(c(nrow(params), ncol(params)), c(8L, 10L)) + expect_named( + params, c( "Parameter", "Level", "Coefficient", "SE", "CI", "CI_low", "CI_high", "Component", "Effects", "Group" @@ -231,9 +231,9 @@ withr::with_options( test_that("model_parameters.mixed-ran_pars", { params <- model_parameters(m1, effects = "random") - expect_equal(c(nrow(params), ncol(params)), c(2, 9)) - expect_identical( - colnames(params), + expect_identical(c(nrow(params), ncol(params)), c(2L, 9L)) + expect_named( + params, c("Parameter", "Coefficient", "SE", "CI", "CI_low", "CI_high", "Effects", "Group", "Component") ) expect_identical( @@ -253,9 +253,9 @@ withr::with_options( test_that("model_parameters.mixed-all_pars", { params <- model_parameters(m1, effects = "all") - expect_equal(c(nrow(params), ncol(params)), c(8, 12)) - expect_identical( - colnames(params), + expect_identical(c(nrow(params), ncol(params)), c(8L, 12L)) + expect_named( + params, c( "Parameter", "Coefficient", "SE", "CI", "CI_low", "CI_high", "z", "df_error", "p", "Effects", "Group", "Component" @@ -285,8 +285,8 @@ withr::with_options( test_that("model_parameters.mixed-all", { params <- model_parameters(m1, effects = "all", group_level = TRUE) expect_identical(c(nrow(params), ncol(params)), c(14L, 13L)) - expect_identical( - colnames(params), + expect_named( + params, c( "Parameter", "Level", "Coefficient", "SE", "CI", "CI_low", "CI_high", "z", "df_error", "p", "Component", "Effects", @@ -347,10 +347,15 @@ withr::with_options( )) test_that("model_parameters.mixed-ran_pars", { - params <- model_parameters(m4, effects = "random") + expect_message( + { + params <- model_parameters(m4, effects = "random") + }, + regex = "Your model may" + ) expect_identical(c(nrow(params), ncol(params)), c(6L, 9L)) - expect_identical( - colnames(params), + expect_named( + params, c("Parameter", "Coefficient", "SE", "CI", "CI_low", "CI_high", "Effects", "Group", "Component") ) expect_identical( @@ -408,33 +413,173 @@ withr::with_options( test_that("print-model_parameters glmmTMB", { skip_on_os(c("mac", "linux", "solaris")) - skip_if_not(getRversion() < "4.3.0") + skip_if_not(getRversion() >= "4.3.3") mp <- model_parameters(m4, effects = "fixed", component = "conditional") out <- utils::capture.output(print(mp)) expect_snapshot(out[-5]) - mp <- model_parameters(m4, ci_random = TRUE, effects = "random", component = "conditional") - expect_snapshot(print(mp)) + + mp <- model_parameters(m4, ci_random = TRUE, effects = "random", component = "conditional", verbose = FALSE) + out <- utils::capture.output(print(mp)) + expect_identical( + attributes(mp)$pretty_labels, + c( + `SD (Intercept)` = "SD (Intercept)", `SD (xb)` = "SD (xb)", + `Cor (Intercept~xb)` = "Cor (Intercept~xb)" + ) + ) + expect_identical( + substr(out, 1, 30), + c( + "# Random Effects", + "", + "Parameter | ", + "------------------------------", + "SD (Intercept: persons) | ", + "SD (xb: persons) | ", + "Cor (Intercept~xb: persons) | " + ) + ) + expect_equal(mp$Coefficient, c(3.40563, 1.21316, -1), tolerance = 1e-3) + expect_equal(mp$CI_low, c(1.64567, 0.5919, -1), tolerance = 1e-3) + mp <- model_parameters(m4, ci_random = TRUE, effects = "fixed", component = "zero_inflated") out <- utils::capture.output(print(mp)) - expect_snapshot(out[-6]) + expect_identical( + attributes(mp)$pretty_labels, + c(`(Intercept)` = "(Intercept)", child = "child", camper1 = "camper [1]") + ) + expect_identical( + substr(out, 1, 12), + c( + "# Fixed Effe", "", "Parameter ", "------------", "(Intercept) ", + "child ", "camper [1] " + ) + ) + expect_equal(mp$Coefficient, c(1.88964, 0.15712, -0.17007), tolerance = 1e-3) + expect_equal(mp$CI_low, c(0.5878, -0.78781, -0.92836), tolerance = 1e-3) - mp <- model_parameters(m4, ci_random = TRUE, effects = "random", component = "zero_inflated") - expect_snapshot(print(mp)) - mp <- model_parameters(m4, ci_random = TRUE, effects = "all", component = "conditional") + mp <- model_parameters(m4, ci_random = TRUE, effects = "random", component = "zero_inflated", verbose = FALSE) out <- utils::capture.output(print(mp)) - expect_snapshot(out[-5]) + expect_identical( + attributes(mp)$pretty_labels, + c( + `SD (Intercept)` = "SD (Intercept)", `SD (zg)` = "SD (zg)", + `Cor (Intercept~zg)` = "Cor (Intercept~zg)" + ) + ) + expect_identical( + substr(out, 1, 30), + c( + "# Random Effects (Zero-Inflati", "", "Parameter | ", + "------------------------------", "SD (Intercept: persons) | ", + "SD (zg: persons) | ", "Cor (Intercept~zg: persons) | " + ) + ) + expect_equal(mp$Coefficient, c(2.73583, 1.56833, 1), tolerance = 1e-3) + expect_equal(mp$CI_low, c(1.16329, 0.64246, -1), tolerance = 1e-3) - mp <- model_parameters(m4, effects = "all", ci_random = TRUE, component = "zero_inflated") + + mp <- model_parameters(m4, ci_random = TRUE, effects = "all", component = "conditional", verbose = FALSE) out <- utils::capture.output(print(mp)) - expect_snapshot(out[-6]) + expect_identical( + attributes(mp)$pretty_labels, + c( + `(Intercept)` = "(Intercept)", child = "child", camper1 = "camper [1]", + `SD (Intercept)` = "SD (Intercept)", `SD (xb)` = "SD (xb)", + `Cor (Intercept~xb)` = "Cor (Intercept~xb)" + ) + ) + expect_identical( + substr(out, 1, 30), + c( + "# Fixed Effects", "", "Parameter | Log-Mean | SE ", + "------------------------------", + "(Intercept) | 2.55 | 0.25 ", "child | -1.09 | 0.10 ", + "camper [1] | 0.27 | 0.10 ", "", "# Random Effects", "", + "Parameter | ", "------------------------------", + "SD (Intercept: persons) | ", "SD (xb: persons) | ", + "Cor (Intercept~xb: persons) | " + ) + ) + expect_equal(mp$Coefficient, c(2.54713, -1.08747, 0.2723, 3.40563, 1.21316, -1), tolerance = 1e-3) + expect_equal(mp$CI_low, c(2.06032, -1.27967, 0.07461, 1.64567, 0.5919, -1), tolerance = 1e-3) - mp <- model_parameters(m4, effects = "all", component = "all", ci_random = TRUE) + + mp <- model_parameters(m4, effects = "all", ci_random = TRUE, component = "zero_inflated", verbose = FALSE) out <- utils::capture.output(print(mp)) - expect_snapshot(out[-c(5, 14)]) + expect_identical( + attributes(mp)$pretty_labels, + c( + `(Intercept)` = "(Intercept)", child = "child", camper1 = "camper [1]", + `SD (Intercept)` = "SD (Intercept)", `SD (zg)` = "SD (zg)", + `Cor (Intercept~zg)` = "Cor (Intercept~zg)" + ) + ) + expect_identical( + substr(out, 1, 30), + c( + "# Fixed Effects (Zero-Inflatio", "", "Parameter | Log-Mean | SE ", + "------------------------------", "(Intercept) | 1.89 | 0.66 ", + "child | 0.16 | 0.48 ", "camper [1] | -0.17 | 0.39 ", + "", "# Random Effects (Zero-Inflati", "", "Parameter | ", + "------------------------------", "SD (Intercept: persons) | ", + "SD (zg: persons) | ", "Cor (Intercept~zg: persons) | " + ) + ) + expect_equal(mp$Coefficient, c(1.88964, 0.15712, -0.17007, 2.73583, 1.56833, 1), tolerance = 1e-3) + expect_equal(mp$CI_low, c(0.5878, -0.78781, -0.92836, 1.16329, 0.64246, -1), tolerance = 1e-3) + + + mp <- model_parameters(m4, effects = "all", component = "all", ci_random = TRUE, verbose = FALSE) + out <- utils::capture.output(print(mp)) + expect_identical( + attributes(mp)$pretty_labels, + c( + `(Intercept)` = "(Intercept)", child = "child", camper1 = "camper [1]", + `(Intercept)` = "(Intercept)", child = "child", camper1 = "camper1", # nolint + `SD (Intercept)` = "SD (Intercept)", `SD (xb)` = "SD (xb)", + `Cor (Intercept~xb)` = "Cor (Intercept~xb)", + `SD (Intercept)` = "SD (Intercept)", `SD (zg)` = "SD (zg)", # nolint + `Cor (Intercept~zg)` = "Cor (Intercept~zg)" + ) + ) + expect_identical( + substr(out, 1, 30), + c( + "# Fixed Effects (Count Model)", "", "Parameter | Log-Mean | SE ", + "------------------------------", "(Intercept) | 2.55 | 0.25 ", + "child | -1.09 | 0.10 ", "camper [1] | 0.27 | 0.10 ", + "", "# Fixed Effects (Zero-Inflatio", "", "Parameter | Log-Odds | SE ", + "------------------------------", "(Intercept) | 1.89 | 0.66 ", + "child | 0.16 | 0.48 ", "camper [1] | -0.17 | 0.39 ", + "", "# Random Effects Variances", "", "Parameter | ", + "------------------------------", "SD (Intercept: persons) | ", + "SD (xb: persons) | ", "Cor (Intercept~xb: persons) | ", + "", "# Random Effects (Zero-Inflati", "", "Parameter | ", + "------------------------------", "SD (Intercept: persons) | ", + "SD (zg: persons) | ", "Cor (Intercept~zg: persons) | " + ) + ) + expect_equal( + mp$Coefficient, + c( + 2.54713, -1.08747, 0.2723, 1.88964, 0.15712, -0.17007, 3.40563, + 1.21316, -1, 2.73583, 1.56833, 1 + ), + tolerance = 1e-3 + ) + expect_equal( + mp$CI_low, + c( + 2.06032, -1.27967, 0.07461, 0.5878, -0.78781, -0.92836, 1.64567, + 0.5919, -1, 1.16329, 0.64246, -1 + ), + tolerance = 1e-3 + ) }) @@ -442,7 +587,7 @@ withr::with_options( test_that("print-model_parameters glmmTMB digits", { skip_on_os(c("mac", "linux", "solaris")) - skip_if_not(getRversion() < "4.3.0") + skip_if_not(getRversion() >= "4.3.3") mp <- model_parameters(m4, ci_random = TRUE, effects = "all", component = "all") out <- utils::capture.output(print(mp, digits = 4, ci_digits = 5)) @@ -459,7 +604,7 @@ withr::with_options( skip_if_not_installed("curl") skip_if_offline() skip_on_os(c("mac", "linux", "solaris")) - skip_if_not(getRversion() < "4.3.0") + skip_if_not(getRversion() >= "4.3.3") model_pr <- tryCatch( { @@ -488,7 +633,7 @@ withr::with_options( test_that("model_parameters.mixed-all", { skip_on_os(c("mac", "linux", "solaris")) - skip_if_not(getRversion() < "4.3.0") + skip_if_not(getRversion() >= "4.3.3") params <- model_parameters(m4, effects = "all") expect_identical(c(nrow(params), ncol(params)), c(12L, 12L)) @@ -527,7 +672,7 @@ withr::with_options( test_that("print-model_parameters", { skip_on_os(c("mac", "linux", "solaris")) - skip_if_not(getRversion() < "4.3.0") + skip_if_not(getRversion() >= "4.3.3") mp <- model_parameters(m1, effects = "fixed", verbose = FALSE) expect_snapshot(mp) diff --git a/tests/testthat/test-random_effects_ci-glmmTMB.R b/tests/testthat/test-random_effects_ci-glmmTMB.R index c5c1d3eed..7f3279969 100644 --- a/tests/testthat/test-random_effects_ci-glmmTMB.R +++ b/tests/testthat/test-random_effects_ci-glmmTMB.R @@ -13,7 +13,7 @@ skip_on_cran() data(sleepstudy, package = "lme4") data(cake, package = "lme4") set.seed(123) -sleepstudy$Months <- sample(1:4, nrow(sleepstudy), TRUE) +sleepstudy$Months <- sample.int(4, nrow(sleepstudy), TRUE) set.seed(123) m1 <- suppressWarnings(glmmTMB::glmmTMB( @@ -26,11 +26,26 @@ m4 <- suppressWarnings(glmmTMB::glmmTMB(angle ~ temperature + (temperature | rep m5 <- suppressWarnings(glmmTMB::glmmTMB(Reaction ~ Days + (Days + Months | Subject), data = sleepstudy)) set.seed(123) -expect_message(mp1 <- model_parameters(m1, ci_random = TRUE), "singularity") +expect_message( + { + mp1 <- model_parameters(m1, ci_random = TRUE) + }, + "singularity" +) mp2 <- model_parameters(m2, ci_random = TRUE) # works -expect_message(mp3 <- model_parameters(m3, ci_random = TRUE), "singularity") # no SE/CI +expect_message( + { + mp3 <- model_parameters(m3, ci_random = TRUE) + }, + "singularity" +) # no SE/CI mp4 <- model_parameters(m4, ci_random = TRUE) -expect_message(mp5 <- model_parameters(m5, ci_random = TRUE), "singularity") # no SE/CI +expect_message( + { + mp5 <- model_parameters(m5, ci_random = TRUE) + }, + "singularity" +) # no SE/CI test_that("random effects CIs, two slopes, categorical", { ## FIXME: Results differ across R versions, no idea why... @@ -198,10 +213,7 @@ test_that("random effects CIs, categorical slope-2", { test_that("random effects CIs, double slope", { expect_equal( mp5$CI_low, - c( - 238.40607, 7.52296, 15.01708, 3.80547, NaN, -0.48781, NaN, - NaN, 22.80046 - ), + c(238.40606, 7.52296, 15.0171, 3.80547, 0, -0.48781, NaN, NaN, 22.80045), tolerance = 1e-2, ignore_attr = TRUE ) @@ -220,7 +232,7 @@ test_that("random effects CIs, double slope", { test_that("random effects CIs, simple slope", { data(sleepstudy, package = "lme4") set.seed(123) - sleepstudy$Months <- sample(1:4, nrow(sleepstudy), TRUE) + sleepstudy$Months <- sample.int(4, nrow(sleepstudy), TRUE) set.seed(123) m2 <- glmmTMB::glmmTMB(Reaction ~ Days + (0 + Days | Subject), data = sleepstudy) @@ -228,7 +240,12 @@ test_that("random effects CIs, simple slope", { set.seed(123) mp2 <- model_parameters(m2, ci_random = TRUE) - expect_message(mp5 <- model_parameters(m5, ci_random = TRUE), "singularity") # no SE/CI + expect_message( + { + mp5 <- model_parameters(m5, ci_random = TRUE) + }, + "singularity" + ) # no SE/CI expect_equal( mp2$CI_low, c(243.55046, 6.89554, 4.98429, 25.94359),