diff --git a/DESCRIPTION b/DESCRIPTION index 62b83fa1e..b1dc40506 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -192,4 +192,17 @@ Language: en-US RoxygenNote: 7.2.0 Config/testthat/edition: 3 Roxygen: list(markdown = TRUE) -Remotes: easystats/insight +Remotes: + easystats/insight +RemoteType: github +RemoteHost: api.github.com +RemoteRepo: glmmTMB +RemoteUsername: glmmTMB +RemoteRef: ci_tweaks +RemoteSha: fd8bb78acd8b198147285ce94b1f67043349f570 +RemoteSubdir: glmmTMB +GithubRepo: glmmTMB +GithubUsername: glmmTMB +GithubRef: ci_tweaks +GithubSHA1: fd8bb78acd8b198147285ce94b1f67043349f570 +GithubSubdir: glmmTMB diff --git a/NAMESPACE b/NAMESPACE index 2e3aafac0..9f6bcade2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +S3method(as.data.frame,VarCorr.lme) S3method(as.data.frame,glmmTMB) S3method(as.data.frame,lm) S3method(as.data.frame,merMod) diff --git a/R/1_model_parameters.R b/R/1_model_parameters.R index fc87a4878..2a3d74089 100644 --- a/R/1_model_parameters.R +++ b/R/1_model_parameters.R @@ -442,6 +442,7 @@ model_parameters.default <- function(model, verbose = verbose ) + # extract model parameters table, as data frame out <- tryCatch( { .model_parameters_generic( @@ -470,6 +471,7 @@ model_parameters.default <- function(model, } ) + # tell user if something went wrong... if (length(out) == 1 && isTRUE(is.na(out))) { msg <- insight::format_message( paste0("Sorry, `model_parameters()` failed with the following error (possible class '", class(model)[1], "' not supported):\n"), @@ -485,6 +487,12 @@ model_parameters.default <- function(model, } + + +# helper function for the composition of the parameters table, +# including a bunch of attributes required for further processing +# (like printing etc.) + .model_parameters_generic <- function(model, ci = .95, bootstrap = FALSE, @@ -527,6 +535,7 @@ model_parameters.default <- function(model, ) args <- c(args, dots) params <- do.call("bootstrap_parameters", args) + } else { # set default method for CI if (is.null(ci_method) || missing(ci_method)) { @@ -552,10 +561,12 @@ model_parameters.default <- function(model, params <- do.call(".extract_parameters_generic", args) } + # exponentiate coefficients and SE/CI, if requested if (isTRUE(exponentiate) || identical(exponentiate, "nongaussian")) { params <- .exponentiate_parameters(params, model, exponentiate) } + # add further information as attributes params <- .add_model_parameters_attributes( params, model, diff --git a/R/extract_parameters.R b/R/extract_parameters.R index f172351a1..d427234c7 100644 --- a/R/extract_parameters.R +++ b/R/extract_parameters.R @@ -180,6 +180,7 @@ # ==== degrees of freedom + if (!is.null(ci_method)) { df_error <- degrees_of_freedom(model, method = ci_method, verbose = FALSE) } else { @@ -310,24 +311,12 @@ .add_sigma_residual_df <- function(params, model) { if (is.null(params$Component) || !"sigma" %in% params$Component) { - sig <- tryCatch( - { - suppressWarnings(insight::get_sigma(model, ci = NULL, verbose = FALSE)) - }, - error = function(e) { - NULL - } - ) + sig <- tryCatch(suppressWarnings(insight::get_sigma(model, ci = NULL, verbose = FALSE)), + error = function(e) NULL) attr(params, "sigma") <- as.numeric(sig) - resdf <- tryCatch( - { - suppressWarnings(insight::get_df(model, type = "residual")) - }, - error = function(e) { - NULL - } - ) + resdf <- tryCatch(suppressWarnings(insight::get_df(model, type = "residual")), + error = function(e) NULL) attr(params, "residual_df") <- as.numeric(resdf) } params diff --git a/R/extract_random_variances.R b/R/extract_random_variances.R index 85c9671b7..b01ac5d8f 100644 --- a/R/extract_random_variances.R +++ b/R/extract_random_variances.R @@ -107,6 +107,8 @@ + + # workhorse ------------------------ .extract_random_variances_helper <- function(model, @@ -116,147 +118,69 @@ ci_method = NULL, verbose = FALSE, ...) { - varcorr <- .get_variance_information(model, component) - - ran_intercept <- tryCatch(data.frame(.random_intercept_variance(varcorr)), - error = function(e) NULL) - - ran_slope <- tryCatch(data.frame(.random_slope_variance(model, varcorr)), - error = function(e) NULL) - - ran_corr <- tryCatch(data.frame(.random_slope_intercept_corr(model, varcorr)), - error = function(e) NULL) - - ran_slopes_corr <- tryCatch(data.frame(.random_slopes_corr(model, varcorr)), - error = function(e) NULL) - - - # sigma/dispersion only once, - if (component == "conditional") { - ran_sigma <- data.frame(insight::get_sigma(model, ci = NULL, verbose = FALSE)) - } else { - ran_sigma <- NULL + if (!inherits(model, "lme")) { + class(varcorr) <- "VarCorr.merMod" } + # return varcorr matrix + re_data <- as.data.frame(varcorr, order = "lower.tri") + + # extract parameters from SD and COR separately, for sorting + re_sd_intercept <- re_data$var1 == "(Intercept)" & is.na(re_data$var2) & re_data$grp != "Residual" + re_sd_slope <- re_data$var1 != "(Intercept)" & is.na(re_data$var2) & re_data$grp != "Residual" + re_cor_intercept <- re_data$var1 == "(Intercept)" & !is.na(re_data$var2) & re_data$grp != "Residual" + re_cor_slope <- re_data$var1 != "(Intercept)" & !is.na(re_data$var2) & re_data$grp != "Residual" + re_sigma <- re_data$grp == "Residual" + + # merge to sorted data frame + out <- rbind( + re_data[re_sd_intercept, ], + re_data[re_sd_slope, ], + re_data[re_cor_intercept, ], + re_data[re_cor_slope, ], + re_data[re_sigma, ] + ) + out$Parameter <- NA - # random intercept - tau00 - if (!is.null(ran_intercept) && nrow(ran_intercept) > 0) { - colnames(ran_intercept) <- "Coefficient" - ran_intercept$Group <- rownames(ran_intercept) - ran_intercept$Parameter <- "SD (Intercept)" - } - ran_groups_int <- ran_intercept$Group - - - # random slope - tau11 - if (!is.null(ran_slope) && nrow(ran_slope) > 0) { - colnames(ran_slope) <- "Coefficient" - ran_slope$Group <- rownames(ran_slope) - ran_groups_slp <- gsub("\\..*", "", ran_slope$Group) - ran_slope$Parameter <- NA - for (i in unique(ran_groups_slp)) { - slopes <- which(grepl(paste0("^\\Q", i, "\\E"), ran_slope$Group)) - if (length(slopes)) { - ran_slope$Parameter[slopes] <- paste0( - "SD (", gsub("^\\.", "", gsub(i, "", ran_slope$Group[slopes], fixed = TRUE)), ")" - ) - ran_slope$Group[slopes] <- i - } - } - } else { - ran_groups_slp <- NULL + # rename SD + sds <- !is.na(out$var1) & is.na(out$var2) + if (any(sds)) { + out$Parameter[sds] <- paste0("SD (", out$var1[sds], ")") } - ran_groups <- unique(c(ran_groups_int, ran_groups_slp)) - - - # random slope-intercept correlation - rho01 - if (!is.null(ran_corr) && nrow(ran_corr) > 0) { - if (ncol(ran_corr) > 1 && all(colnames(ran_corr) %in% ran_groups)) { - ran_corr <- datawizard::reshape_longer(ran_corr, colnames_to = "Group", values_to = "Coefficient", rows_to = "Slope") - ran_corr$Parameter <- paste0("Cor (Intercept~", ran_corr$Slope, ": ", ran_intercept$Group, ")") - ran_corr <- datawizard::data_reorder(ran_corr, select = c("Parameter", "Coefficient", "Group")) - ran_corr$Slope <- NULL - } else if (!is.null(ran_intercept$Group) && colnames(ran_corr)[1] == ran_intercept$Group[1]) { - colnames(ran_corr)[1] <- "Coefficient" - ran_corr$Parameter <- paste0("Cor (Intercept~", row.names(ran_corr), ": ", ran_intercept$Group[1], ")") - ran_corr$Group <- ran_intercept$Group[1] - } else if (!is.null(ran_groups) && colnames(ran_corr)[1] == ran_groups[1]) { - # colnames(ran_corr)[1] <- "Coefficient" - # ran_corr$Parameter <- paste0("Cor (Reference~", row.names(ran_corr), ")") - # ran_corr$Group <- ran_groups[1] - # this occurs when model has no random intercept - # will be captures by random slope-slope-correlation - ran_corr <- NULL - } else { - colnames(ran_corr) <- "Coefficient" - ran_corr$Group <- rownames(ran_corr) - ran_corr$Parameter <- NA - for (i in unique(ran_groups)) { - corrs <- which(grepl(paste0("^\\Q", i, "\\E"), ran_corr$Group)) - if (length(corrs)) { - param_name <- i - cor_slopes <- which(grepl(paste0("^\\Q", i, "\\E"), ran_slope$Group)) - if (length(cor_slopes)) { - param_name <- paste0( - gsub("SD \\((.*)\\)", "\\1", ran_slope$Parameter[cor_slopes]), - ": ", - i - ) - } - ran_corr$Parameter[corrs] <- paste0("Cor (Intercept~", param_name, ")") - ran_corr$Group[corrs] <- i - } - } - } + # rename correlations + corrs <- !is.na(out$var2) + if (any(corrs)) { + out$Parameter[corrs] <- paste0("Cor (", out$var1[corrs], "~", out$var2[corrs], ")") } - - # random slope correlation - rho00 - if (!is.null(ran_slopes_corr) && nrow(ran_slopes_corr) > 0) { - term <- rownames(ran_slopes_corr) - colnames(ran_slopes_corr) <- "Coefficient" - rg <- paste0("(", paste0(c(ran_groups, paste0(ran_groups, "\\.\\d+")), collapse = "|"), ")") - grp <- gsub(paste0(rg, "\\.(.*)-(.*)"), "\\1", term) - slope1 <- gsub(paste0(rg, "\\.(.*)-(.*)"), "\\2", term) - slope2 <- gsub(paste0(rg, "\\.(.*)-(.*)"), "\\3", term) - ran_slopes_corr$Parameter <- paste0("Cor (", slope1, "~", slope2, ": ", grp, ")") - ran_slopes_corr$Group <- grp - } - - - # residuals - sigma - if (!is.null(ran_sigma) && nrow(ran_sigma) > 0) { - colnames(ran_sigma) <- "Coefficient" - ran_sigma$Group <- "Residual" - ran_sigma$Parameter <- "SD (Observations)" + # rename sigma + sigma <- out$grp == "Residual" + if (any(sigma)) { + out$Parameter[sigma] <- "SD (Observations)" } - # row bind all random effect variances, if possible - out <- tryCatch( - { - out_list <- insight::compact_list(list(ran_intercept, ran_slope, ran_corr, ran_slopes_corr, ran_sigma)) - do.call(rbind, out_list) - }, - error = function(e) { - NULL - } + # rename columns + out <- datawizard::data_rename( + out, + pattern = c("grp", "sdcor"), + replacement = c("Group", "Coefficient") ) - if (is.null(out)) { - return(NULL) - } + # fix names for uncorrelated slope-intercepts + pattern <- paste0("(", paste0(insight::find_random(model, flatten = TRUE), collapse = "|"), ")\\.\\d+$") + out$Group <- gsub(pattern, "\\1", out$Group) - rownames(out) <- NULL + # remove non-used columns + out$var1 <- NULL + out$var2 <- NULL + out$grp <- NULL + out$vcov <- NULL + out$sdcor <- NULL - # variances to SD (sqrt), except correlations and Sigma - corr_param <- grepl("^Cor (.*)", out$Parameter) - sigma_param <- out$Parameter == "SD (Observations)" - not_cor_and_sigma <- !corr_param & !sigma_param - if (any(not_cor_and_sigma)) { - out$Coefficient[not_cor_and_sigma] <- sqrt(out$Coefficient[not_cor_and_sigma]) - } + # fix intercept names + out$Parameter <- gsub("(Intercept)", "Intercept", out$Parameter, fixed = TRUE) stat_column <- gsub("-statistic", "", insight::find_statistic(model), fixed = TRUE) @@ -282,6 +206,10 @@ } out[ci_cols] <- NA + # variances to SD (sqrt), except correlations and Sigma + corr_param <- grepl("^Cor (.*)", out$Parameter) + sigma_param <- out$Parameter == "SD (Observations)" + # add confidence intervals? if (!is.null(ci) && !all(is.na(ci)) && length(ci) == 1) { out <- .random_sd_ci(model, out, ci_method, ci, corr_param, sigma_param, component, verbose = verbose) @@ -299,6 +227,78 @@ +#' @export +as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ...) { + # retrieve RE SD and COR + stddevs <- sapply(x[, "StdDev"], as.numeric) + if ("Corr" %in% colnames(x)) { + corrs <- suppressWarnings(sapply(x[, "Corr"], as.numeric)) + } else { + corrs <- NULL + } + grps <- grepl(" =$", names(stddevs)) + + # for multiple grouping factors, split at each group + if (any(grps)) { + from <- which(grps) + to <- c(which(grps) - 1, length(grps))[-1] + out_sd <- do.call(rbind, lapply(1:length(from), function(i) { + values <- stddevs[from[i]:to[i]] + data.frame( + grp = gsub("(.*) =$", "\\1", names(values[1])), + var1 = names(values[-1]), + var2 = NA_character_, + sdcor = unname(values[-1]), + stringsAsFactors = FALSE + ) + })) + if (!is.null(corrs)) { + out_cor <- do.call(rbind, lapply(1:length(from), function(i) { + values <- corrs[from[i]:to[i]] + data.frame( + grp = gsub("(.*) =$", "\\1", names(values[1])), + var1 = "(Intercept)", + var2 = names(values[-1]), + sdcor = unname(values[-1]), + stringsAsFactors = FALSE + ) + })) + } else { + out_cor <- NULL + } + } else { + out_sd <- data.frame( + grp = gsub("(.*) =(.*)", "\\1", attributes(x)$title), + var1 = names(stddevs), + var2 = NA_character_, + sdcor = unname(stddevs), + stringsAsFactors = FALSE + ) + if (!is.null(corrs)) { + out_cor <- data.frame( + grp = gsub("(.*) =(.*)", "\\1", attributes(x)$title), + var1 = "(Intercept)", + var2 = names(corrs), + sdcor = unname(corrs), + stringsAsFactors = FALSE + ) + } else { + out_cor <- NULL + } + } + + out_sd$grp[out_sd$var1 == "Residual"] <- "Residual" + out_sd$var1[out_sd$grp == "Residual"] <- NA_character_ + out_sd$var2[out_sd$grp == "Residual"] <- NA_character_ + out_cor <- out_cor[!is.na(out_cor$sdcor), ] + + rbind(out_sd, out_cor) +} + + + + + # extract CI for random SD ------------------------ .random_sd_ci <- function(model, out, ci_method, ci, corr_param, sigma_param, component = NULL, verbose = FALSE) { @@ -309,38 +309,58 @@ } if (inherits(model, c("merMod", "glmerMod", "lmerMod"))) { - if (!is.null(ci_method) && ci_method %in% c("profile", "boot")) { - 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) - rn <- gsub("sd_(.*)(\\|)(.*)", "\\1: \\3", rn) - rn <- gsub("|", ":", rn, fixed = TRUE) - rn <- gsub("[\\(\\)]", "", rn) - rn <- gsub("cor_(.*)\\.(.*)", "cor \\2", rn) + # lme4 - boot and profile - var_ci_corr_param <- grepl("^cor ", rn) - var_ci_sigma_param <- rn == "sigma" + 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))) + colnames(var_ci) <- c("CI_low", "CI_high") + + rn <- row.names(var_ci) + rn <- gsub("sd_(.*)(\\|)(.*)", "\\1: \\3", rn) + rn <- gsub("|", ":", rn, fixed = TRUE) + rn <- gsub("[\\(\\)]", "", rn) + rn <- gsub("cor_(.*)\\.(.*)", "cor \\2", rn) + + var_ci_corr_param <- grepl("^cor ", rn) + var_ci_sigma_param <- rn == "sigma" + + out$CI_low[!corr_param & !sigma_param] <- var_ci$CI_low[!var_ci_corr_param & !var_ci_sigma_param] + out$CI_high[!corr_param & !sigma_param] <- var_ci$CI_high[!var_ci_corr_param & !var_ci_sigma_param] + + if (any(sigma_param) && any(var_ci_sigma_param)) { + out$CI_low[sigma_param] <- var_ci$CI_low[var_ci_sigma_param] + out$CI_high[sigma_param] <- var_ci$CI_high[var_ci_sigma_param] + } - out$CI_low[!corr_param & !sigma_param] <- var_ci$CI_low[!var_ci_corr_param & !var_ci_sigma_param] - out$CI_high[!corr_param & !sigma_param] <- var_ci$CI_high[!var_ci_corr_param & !var_ci_sigma_param] + if (any(corr_param) && any(var_ci_corr_param)) { + out$CI_low[corr_param] <- var_ci$CI_low[var_ci_corr_param] + out$CI_high[corr_param] <- var_ci$CI_high[var_ci_corr_param] + } + out + }, + error = function(e) { + if (isTRUE(verbose)) { + message(insight::format_message( + "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')." + )) + } + out + } + ) + } else if (!is.null(ci_method)) { - if (any(sigma_param) && any(var_ci_sigma_param)) { - out$CI_low[sigma_param] <- var_ci$CI_low[var_ci_sigma_param] - out$CI_high[sigma_param] <- var_ci$CI_high[var_ci_sigma_param] - } + # lme4 - wald / normal CI - if (any(corr_param) && any(var_ci_corr_param)) { - out$CI_low[corr_param] <- var_ci$CI_low[var_ci_corr_param] - out$CI_high[corr_param] <- var_ci$CI_high[var_ci_corr_param] - } - } else if (!is.null(ci_method)) { # Wald based CIs # see https://stat.ethz.ch/pipermail/r-sig-mixed-models/2022q1/029985.html if (all(insight::check_if_installed(c("merDeriv", "lme4"), quietly = TRUE))) { # this may fail, so wrap in try-catch - tryCatch( + out <- tryCatch( { # vcov from full model. the parameters from vcov have a different # order, so we need to restore the "original" order of random effect @@ -381,7 +401,7 @@ var_ci_corr_param <- grepl("(.*)\\.\\(Intercept\\)", var_ci$Parameter) if (any(var_ci_corr_param)) { rnd_slope_terms <- gsub("(.*)\\.\\(Intercept\\)", "\\1", var_ci$Parameter[var_ci_corr_param]) - var_ci$Parameter[var_ci_corr_param] <- paste0("Cor (Intercept~", rnd_slope_terms, ": ", var_ci$Group[var_ci_corr_param], ")") + var_ci$Parameter[var_ci_corr_param] <- paste0("Cor (Intercept~", rnd_slope_terms, ")") } # correlations w/o intercept? usually only for factors @@ -391,7 +411,7 @@ if (any(rnd_slope_corr)) { 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]) + 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] } } @@ -437,6 +457,7 @@ "Some of the standard errors and confidence intervals of the random effects parameters are probably not meaningful!" )) } + out }, error = function(e) { if (isTRUE(verbose)) { @@ -452,6 +473,7 @@ )) } } + out } ) } else if (isTRUE(verbose)) { @@ -459,6 +481,9 @@ } } } else if (inherits(model, "glmmTMB")) { + + # glmmTMB random-effects-CI + ## TODO "profile" seems to be less stable, so only wald? out <- tryCatch( { @@ -466,53 +491,61 @@ as.data.frame(suppressWarnings(stats::confint(model, parm = "theta_", method = "wald", level = ci))), as.data.frame(suppressWarnings(stats::confint(model, parm = "sigma", method = "wald", level = ci))) ) + colnames(var_ci) <- c("CI_low", "CI_high", "not_used") var_ci$Component <- "conditional" - # # regex-pattern to find conditional and ZI components - group_factor <- insight::find_random(model, flatten = TRUE) - group_factor2 <- paste0("(", paste(group_factor, collapse = "|"), ")") - var_ci$Parameter <- row.names(var_ci) - pattern <- paste0("^(zi\\.|", group_factor2, "\\.zi\\.)") - zi_rows <- grepl(pattern, var_ci$Parameter) - if (any(zi_rows)) { - var_ci$Component[zi_rows] <- "zi" - } - # 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"]) + if (utils::packageVersion("glmmTMB") > "1.1.3") { + var_ci$Component[grepl("^zi\\.", var_ci$Parameter)] <- "zi" + # remove cond/zi prefix + var_ci$Parameter <- gsub("^(cond\\.|zi\\.)(.*)", "\\2", var_ci$Parameter) + # copy RE group + var_ci$Group <- gsub("(.*)\\|(.*)$", "\\2", var_ci$Parameter) + var_ci$Parameter <- gsub("(.*)\\|(.*)$", "\\1", var_ci$Parameter) + var_ci$Group[rownames(var_ci) == "sigma"] <- "Residual" } else { - var_ci$Group <- group_factor - # check if sigma was properly identified - if (!"sigma" %in% var_ci$Group && "sigma" %in% rownames(var_ci)) { - var_ci$Group[rownames(var_ci) == "sigma"] <- "Residual" + # regex-pattern to find conditional and ZI components + group_factor <- insight::find_random(model, flatten = TRUE) + group_factor2 <- paste0("(", paste(group_factor, collapse = "|"), ")") + + pattern <- paste0("^(zi\\.|", group_factor2, "\\.zi\\.)") + zi_rows <- grepl(pattern, var_ci$Parameter) + if (any(zi_rows)) { + var_ci$Component[zi_rows] <- "zi" + } + + # 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"]) + } else { + var_ci$Group <- group_factor + # check if sigma was properly identified + if (!"sigma" %in% var_ci$Group && "sigma" %in% rownames(var_ci)) { + var_ci$Group[rownames(var_ci) == "sigma"] <- "Residual" + } } - } - var_ci$Group[var_ci$Group == "sigma"] <- "Residual" - # remove cond/zi prefix - pattern <- paste0("^(cond\\.|zi\\.|", group_factor, "\\.cond\\.|", group_factor, "\\.zi\\.)(.*)") - for (p in pattern) { - var_ci$Parameter <- gsub(p, "\\2", var_ci$Parameter) + # remove cond/zi prefix + pattern <- paste0("^(cond\\.|zi\\.|", group_factor, "\\.cond\\.|", group_factor, "\\.zi\\.)(.*)") + for (p in pattern) { + var_ci$Parameter <- gsub(p, "\\2", var_ci$Parameter) + } } + # fix SD and Cor names var_ci$Parameter <- gsub(".Intercept.", "(Intercept)", var_ci$Parameter, fixed = TRUE) var_ci$Parameter <- gsub("^(Std\\.Dev\\.)(.*)", "SD \\(\\2\\)", var_ci$Parameter) - var_ci$Parameter <- gsub("^Cor\\.(.*)\\.(.*)", "Cor \\(\\2~\\1:", var_ci$Parameter) + var_ci$Parameter <- gsub("^Cor\\.(.*)\\.(.*)", "Cor \\(\\2~\\1\\)", var_ci$Parameter) # minor cleaning var_ci$Parameter <- gsub("((", "(", var_ci$Parameter, fixed = TRUE) var_ci$Parameter <- gsub("))", ")", var_ci$Parameter, fixed = TRUE) var_ci$Parameter <- gsub(")~", "~", var_ci$Parameter, fixed = TRUE) # fix sigma var_ci$Parameter[var_ci$Parameter == "sigma"] <- "SD (Observations)" - # add name of group factor to cor - cor_params <- grepl("^Cor ", var_ci$Parameter) - if (any(cor_params)) { - var_ci$Parameter[cor_params] <- paste0(var_ci$Parameter[cor_params], " ", group_factor, ")") - } + var_ci$Group[var_ci$Group == "sigma"] <- "Residual" # remove unused columns (that are added back after merging) out$CI_low <- NULL @@ -523,83 +556,35 @@ var_ci$not_used <- NULL var_ci$Component <- NULL - merge(out, var_ci, sort = FALSE, all.x = TRUE) - # - # groups <- utils::stack(insight::find_random(model, flatten = FALSE)) - # colnames(groups) <- c("Group", "Component") - # groups$Component <- ifelse(groups$Component == "random", "conditional", "zi") - # - # # regex-pattern to find conditional and ZI components - # group_factor <- insight::find_random(model, flatten = TRUE) - # group_factor2 <- paste0("(", paste(group_factor, collapse = "|"), ")") - # - # thetas <- as.data.frame(suppressWarnings(stats::confint(model, parm = "theta_", method = "wald", level = ci))) - # thetas$Parameter <- row.names(thetas) - # thetas$Component <- "conditional" - # # find zi-prefix, to set correct component value - # pattern <- paste0("^(zi\\.|", group_factor2, "\\.zi\\.)") - # thetas$Component[grepl(pattern, row.names(thetas))] <- "zi" - # - # if (nrow(thetas) == nrow(groups)) { - # thetas <- cbind(thetas, groups) - # } else { - # thetas <- merge(thetas, groups, sort = FALSE) - # } - # - # # reorder columns - # thetas <- datawizard::data_relocate(thetas, cols = "Component", after = "Group") - # thetas <- datawizard::data_relocate(thetas, cols = "Parameter") - # - # sigma <- as.data.frame(suppressWarnings(stats::confint(model, parm = "sigma", method = "wald", level = ci))) - # - # # check for sigma component - # if (nrow(sigma) > 0) { - # sigma$Parameter <- row.names(sigma) - # sigma$Group <- "Residual" - # sigma$Component <- "conditional" - # sigma <- datawizard::data_relocate(sigma, cols = "Parameter") - # var_ci <- rbind(thetas, sigma) - # } else { - # var_ci <- thetas - # } - # - # colnames(var_ci) <- c("Parameter", "CI_low", "CI_high", "not_used", "Group", "Component") - # - # # remove cond/zi prefix - # pattern <- paste0("^(cond\\.|zi\\.|", group_factor2, "\\.cond\\.|", group_factor2, "\\.zi\\.)") - # var_ci$Parameter <- gsub(pattern, "", var_ci$Parameter) - # # fix SD and Cor names - # var_ci$Parameter <- gsub(".Intercept.", "(Intercept)", var_ci$Parameter, fixed = TRUE) - # var_ci$Parameter <- gsub("^(Std\\.Dev\\.)(.*)", "SD \\(\\2\\)", var_ci$Parameter) - # var_ci$Parameter <- gsub("^Cor\\.(.*)\\.(.*)", "Cor \\(\\2~\\1:", var_ci$Parameter) - # # minor cleaning - # var_ci$Parameter <- gsub("((", "(", var_ci$Parameter, fixed = TRUE) - # var_ci$Parameter <- gsub("))", ")", var_ci$Parameter, fixed = TRUE) - # var_ci$Parameter <- gsub(")~", "~", var_ci$Parameter, fixed = TRUE) - # # fix sigma - # var_ci$Parameter[var_ci$Parameter == "sigma"] <- "SD (Observations)" - # # add name of group factor to cor - # cor_params <- grepl("^Cor ", var_ci$Parameter) - # if (any(cor_params)) { - # # this might break if length(group_factor) > 1; I don't have a test case handy - # var_ci$Parameter[cor_params] <- paste0(var_ci$Parameter[cor_params], " ", group_factor, ")") - # } - # - # # remove unused columns (that are added back after merging) - # out$CI_low <- NULL - # out$CI_high <- NULL - # - # # filter component - # var_ci <- var_ci[var_ci$Component == component, ] - # var_ci$not_used <- NULL - # var_ci$Component <- NULL - # - # merge(out, var_ci, sort = FALSE, all.x = TRUE) + # 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)) + + if (singular_fit) { + message(insight::format_message( + "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!" + )) + } else if (missing_ci) { + message(insight::format_message( + "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!" + )) + } + } + + # merge and sort + out$.sort_id <- 1:nrow(out) + out <- merge(out, var_ci, sort = FALSE, all.x = TRUE) + out <- out[order(out$.sort_id), ] + out$.sort_id <- NULL + out }, error = function(e) { if (isTRUE(verbose)) { message(insight::format_message( - "Cannot compute standard errors and confidence intervals for random effects parameters.", + "Cannot compute confidence intervals for random effects parameters.", "Your model may suffer from singularity (see '?lme4::isSingular' and '?performance::check_singularity')." )) } @@ -613,6 +598,8 @@ + + # Extract Variance and Correlation Components ---- # store essential information about variance components... @@ -675,11 +662,13 @@ 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, "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, "useSc") <- FALSE } varcorr <- datawizard::compact_list(list(vc1, vc2)) @@ -695,22 +684,18 @@ varcorr <- list(varcorr) names(varcorr) <- re_names[1] attr(varcorr, "sc") <- model$coef$sigma2[[1]] + attr(varcorr, "useSc") <- TRUE # nlme # --------------------------- } else if (inherits(model, "lme")) { - re_names <- insight::find_random(model, split_nested = TRUE, flatten = TRUE) - if (.is_nested_lme(model)) { - varcorr <- .get_nested_lme_varcorr(model) - } else { - varcorr <- list(nlme::getVarCov(model)) - } - names(varcorr) <- re_names + varcorr <- lme4::VarCorr(model) # ordinal # --------------------------- } else if (inherits(model, "clmm")) { varcorr <- ordinal::VarCorr(model) + attr(varcorr, "useSc") <- FALSE # glmmadmb # --------------------------- @@ -744,6 +729,7 @@ # --------------------------- } else if (inherits(model, "cpglmm")) { varcorr <- cplm::VarCorr(model) + attr(varcorr, "useSc") <- FALSE # lme4 / glmmTMB # --------------------------- @@ -767,64 +753,6 @@ -# Caution! this is somewhat experimental... -# It retrieves the variance-covariance matrix of random effects -# from nested lme-models. -.get_nested_lme_varcorr <- function(model) { - # installed? - insight::check_if_installed("lme4") - - vcor <- lme4::VarCorr(model) - class(vcor) <- "matrix" - - re_index <- (which(rownames(vcor) == "(Intercept)") - 1)[-1] - vc_list <- split(data.frame(vcor, stringsAsFactors = FALSE), findInterval(1:nrow(vcor), re_index)) - vc_rownames <- split(rownames(vcor), findInterval(1:nrow(vcor), re_index)) - re_pars <- unique(unlist(insight::find_parameters(model)["random"])) - re_names <- insight::find_random(model, split_nested = TRUE, flatten = TRUE) - - names(vc_list) <- re_names - - mapply( - function(x, y) { - if ("Corr" %in% colnames(x)) { - g_cor <- suppressWarnings(stats::na.omit(as.numeric(x[, "Corr"]))) - } else { - g_cor <- NULL - } - row.names(x) <- as.vector(y) - vl <- rownames(x) %in% re_pars - x <- suppressWarnings(apply(x[vl, vl, drop = FALSE], MARGIN = c(1, 2), FUN = as.numeric)) - m1 <- matrix(, nrow = nrow(x), ncol = ncol(x)) - m1[1:nrow(m1), 1:ncol(m1)] <- as.vector(x[, 1]) - rownames(m1) <- rownames(x) - colnames(m1) <- rownames(x) - - if (!is.null(g_cor)) { - m1_cov <- sqrt(prod(diag(m1))) * g_cor - for (j in 1:ncol(m1)) { - m1[j, nrow(m1) - j + 1] <- m1_cov[1] - } - } - - attr(m1, "cor_slope_intercept") <- g_cor - m1 - }, - vc_list, - vc_rownames, - SIMPLIFY = FALSE - ) -} - - -.is_nested_lme <- function(model) { - re <- insight::find_random(model) - if (is.null(re)) { - return(FALSE) - } - sapply(re, function(i) any(grepl(":", i, fixed = TRUE))) -} - # glmmTMB returns a list of model information, one for conditional @@ -846,184 +774,3 @@ x } } - - - - -#### helper to extract various random effect variances ----------------------- - - -# random slope-variances (tau 11) ---- -# ---------------------------------------------- -.random_slope_variance <- function(model, varcorr) { - if (inherits(model, "lme")) { - unlist(lapply(varcorr, function(x) diag(x)[-1])) - } else { - # random slopes for correlated slope-intercept - out <- unlist(lapply(varcorr, function(x) diag(x)[-1])) - # check for uncorrelated random slopes-intercept - non_intercepts <- which(sapply(varcorr, function(i) !grepl("^\\(Intercept\\)", dimnames(i)[[1]][1]))) - if (length(non_intercepts)) { - if (length(non_intercepts) == length(varcorr)) { - out <- unlist(lapply(varcorr, function(x) diag(x))) - } else { - dn <- unlist(lapply(varcorr, function(i) dimnames(i)[1])[non_intercepts]) - rndslopes <- unlist(lapply(varcorr, function(i) { - if (is.null(dim(i)) || identical(dim(i), c(1, 1))) { - as.vector(i) - } else { - as.vector(diag(i)) - } - })[non_intercepts]) - # random slopes for uncorrelated slope-intercept - names(rndslopes) <- gsub("(.*)\\.\\d+$", "\\1", names(rndslopes)) - rndslopes <- stats::setNames(rndslopes, paste0(names(rndslopes), ".", dn)) - # anything missing? (i.e. correlated slope-intercept slopes) - missig_rnd_slope <- setdiff(names(out), names(rndslopes)) - if (length(missig_rnd_slope)) { - # sanity check - to_remove <- c() - for (j in 1:length(out)) { - # identical random slopes might have different names, so - # we here check if random slopes from correlated and uncorrelated - # are duplicated (i.e. their difference is 0 - including a tolerance) - # and then remove duplicated elements - the_same <- which(abs(outer(out[j], rndslopes, `-`)) < .0001) - if (length(the_same) && grepl(dn[the_same], names(out[j]), fixed = TRUE)) { - to_remove <- c(to_remove, j) - } - } - if (length(to_remove)) { - out <- out[-to_remove] - } - out <- c(out, rndslopes) - } else { - out <- rndslopes - } - } - } - out - } -} - - - -# random intercept-variances, i.e. -# between-subject-variance (tau 00) ---- -# ---------------------------------------------- -.random_intercept_variance <- function(varcorr) { - vars <- lapply(varcorr, function(i) i[1]) - # check for uncorrelated random slopes-intercept - non_intercepts <- which(sapply(varcorr, function(i) !grepl("^\\(Intercept\\)", dimnames(i)[[1]][1]))) - if (length(non_intercepts)) { - vars <- vars[-non_intercepts] - } - - sapply(vars, function(i) i) -} - - - -# slope-intercept-correlations (rho 01) ---- -# ---------------------------------------------- -.random_slope_intercept_corr <- function(model, varcorr) { - if (inherits(model, "lme")) { - rho01 <- unlist(sapply(varcorr, function(i) attr(i, "cor_slope_intercept"))) - if (is.null(rho01)) { - vc <- lme4::VarCorr(model) - if ("Corr" %in% colnames(vc)) { - re_name <- insight::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 - } - } - } - rho01 - } else { - corrs <- lapply(varcorr, attr, "correlation") - rho01 <- sapply(corrs, function(i) { - if (!is.null(i) && colnames(i)[1] == "(Intercept)") { - i[-1, 1] - } else { - NULL - } - }) - unlist(rho01) - } -} - - - -# slope-slope-correlations (rho 00) ---- -# ---------------------------------------------- -.random_slopes_corr <- function(model, varcorr) { - corrs <- lapply(varcorr, attr, "correlation") - rnd_slopes <- unlist(insight::find_random_slopes(model)) - - # check if any categorical random slopes. we then have - # correlation among factor levels - cat_random_slopes <- tryCatch( - { - d <- insight::get_data(model)[rnd_slopes] - any(sapply(d, is.factor)) - }, - error = function(e) { - NULL - } - ) - - # check if any polynomial / I term in random slopes. - # we then have correlation among levels - rs_names <- unique(unlist(lapply(corrs, colnames))) - pattern <- paste0("(I|poly)(.*)(", paste0(rnd_slopes, collapse = "|"), ")") - poly_random_slopes <- any(grepl(pattern, rs_names)) - - if (length(rnd_slopes) < 2 && !isTRUE(cat_random_slopes) && !isTRUE(poly_random_slopes)) { - return(NULL) - } - - rho00 <- tryCatch( - { - insight::compact_list(lapply(corrs, function(d) { - d[upper.tri(d, diag = TRUE)] <- NA - d <- as.data.frame(d) - - d <- datawizard::reshape_longer(d, colnames_to = "Parameter1", rows_to = "Parameter2", verbose = FALSE) - d <- d[stats::complete.cases(d), ] - d <- d[!d$Parameter1 %in% c("Intercept", "(Intercept)"), ] - - if (nrow(d) == 0) { - return(NULL) - } - - d$Parameter <- paste0(d$Parameter1, "-", d$Parameter2) - d$Parameter1 <- d$Parameter2 <- NULL - stats::setNames(d$Value, d$Parameter) - })) - }, - error = function(e) { - NULL - } - ) - - # rho01 <- tryCatch( - # { - # sapply(corrs, function(i) { - # if (!is.null(i)) { - # slope_pairs <- utils::combn(x = rnd_slopes, m = 2, simplify = FALSE) - # lapply(slope_pairs, function(j) { - # stats::setNames(i[j[1], j[2]], paste0("..", paste0(j, collapse = "-"))) - # }) - # } else { - # NULL - # } - # }) - # }, - # error = function(e) { - # NULL - # } - # ) - - unlist(rho00) -} diff --git a/R/format.R b/R/format.R index b51632203..acacd2142 100644 --- a/R/format.R +++ b/R/format.R @@ -73,6 +73,8 @@ format.parameters_model <- function(x, ran_pars <- which(x$Effects == "random") stddevs <- grepl("^SD \\(", x$Parameter[ran_pars]) x$Parameter[ran_pars[stddevs]] <- paste0(gsub("(.*)\\)", "\\1", x$Parameter[ran_pars[stddevs]]), ": ", x$Group[ran_pars[stddevs]], ")") + corrs <- grepl("^Cor \\(", x$Parameter[ran_pars]) + x$Parameter[ran_pars[corrs]] <- paste0(gsub("(.*)\\)", "\\1", x$Parameter[ran_pars[corrs]]), ": ", x$Group[ran_pars[corrs]], ")") x$Parameter[x$Parameter == "SD (Observations: Residual)"] <- "SD (Residual)" x$Group <- NULL } diff --git a/R/utils_format.R b/R/utils_format.R index a9acd5c61..5b12416f4 100644 --- a/R/utils_format.R +++ b/R/utils_format.R @@ -36,6 +36,9 @@ x[[param_col]] <- insight::trim_ws(paste0(x[[param_col]], linesep, "(", x$SE, ")")) x <- x[c(param_col, "p")] colnames(x) <- paste0(colnames(x), " (", modelname, ")") + } else if (style %in% c("est", "coef")) { + x <- x[1] + colnames(x) <- paste0(colnames(x), " (", modelname, ")") } x[[1]][x[[1]] == "()"] <- "" x diff --git a/R/utils_model_parameters.R b/R/utils_model_parameters.R index 77033abfd..d18ebaea5 100644 --- a/R/utils_model_parameters.R +++ b/R/utils_model_parameters.R @@ -29,9 +29,7 @@ # for simplicity, we just use the model information from the first formula # when we have multivariate response models... - if (!is.null(info) && - insight::is_multivariate(model) && - !"is_zero_inflated" %in% names(info)) { + if (insight::is_multivariate(model) && !"is_zero_inflated" %in% names(info)) { info <- info[[1]] } @@ -81,10 +79,11 @@ } - # Models for which titles should be removed - - # here we add exceptions for objects that should - # not have a table headline - if (inherits(model, c("emmGrid", "emm_list", "lm", "glm", "coxph", "bfsl", "deltaMethod"))) { + # Models for which titles should be removed - here we add exceptions for + # objects that should not have a table headline like "# Fixed Effects", when + # there is nothing else than fixed effects (redundant title) + if (inherits(model, c("mediate", "emmGrid", "emm_list", "lm", "glm", "coxph", "bfsl", "deltaMethod"))) { + attr(params, "no_caption") <- TRUE attr(params, "title") <- "" } diff --git a/tests/testthat/test-GLMMadaptive.R b/tests/testthat/test-GLMMadaptive.R index b0d01580e..bb3aa8df9 100644 --- a/tests/testthat/test-GLMMadaptive.R +++ b/tests/testthat/test-GLMMadaptive.R @@ -167,8 +167,8 @@ if (requiet("testthat") && expect_equal( params$Parameter, c( - "SD (Intercept)", "SD (DOY)", "Cor (Intercept~DOY: site)", "SD (Observations)", - "SD (Intercept)", "SD (DOP)", "Cor (Intercept~DOP: site)" + "SD (Intercept)", "SD (DOY)", "Cor (Intercept~DOY)", "SD (Observations)", + "SD (Intercept)", "SD (DOP)", "Cor (Intercept~DOP)" ) ) expect_equal( @@ -180,7 +180,7 @@ if (requiet("testthat") && ) expect_equal( params$Coefficient, - c(0.56552, 0.29951, 0.06307, 0, 1.02233, 0.38209, -0.17162), + c(0.56552, 0.29951, 0.06307, 1.61936, 1.02233, 0.38209, -0.17162), tolerance = 1e-2 ) }) diff --git a/tests/testthat/test-glmer.R b/tests/testthat/test-glmer.R index 776d18ac6..30a28c6c9 100644 --- a/tests/testthat/test-glmer.R +++ b/tests/testthat/test-glmer.R @@ -54,8 +54,7 @@ if (.runThisTest && "", "Parameter | Coefficient", "----------------------------------", - "SD (Intercept: herd) | 0.64", - "SD (Residual) | 1.00" + "SD (Intercept: herd) | 0.64" ) ) @@ -85,8 +84,7 @@ if (.runThisTest && "", "Parameter | Coefficient | SE | 95% CI", "--------------------------------------------------------", - "SD (Intercept: herd) | 0.64 | 0.18 | [0.37, 1.11]", - "SD (Residual) | 1.00 | | " + "SD (Intercept: herd) | 0.64 | 0.18 | [0.37, 1.11]" ) ) }) diff --git a/tests/testthat/test-glmmTMB.R b/tests/testthat/test-glmmTMB.R index f6a640b4f..3d1c92482 100644 --- a/tests/testthat/test-glmmTMB.R +++ b/tests/testthat/test-glmmTMB.R @@ -152,10 +152,7 @@ if (.runThisTest && ) expect_equal( model_parameters(m1, effects = "all")$Coefficient, - c( - 1.2628, -1.14165, 0.73354, -0.38939, 2.05407, -1.00823, 0.9312, - 1, 1.17399 - ), + c(1.2628, -1.14165, 0.73354, -0.38939, 2.05407, -1.00823, 0.9312, 1.17399), tolerance = 1e-3 ) expect_equal( @@ -176,17 +173,13 @@ if (.runThisTest && model_parameters(m1)$Component, c( "conditional", "conditional", "conditional", "zero_inflated", - "zero_inflated", "zero_inflated", "conditional", "conditional", - "zero_inflated" + "zero_inflated", "zero_inflated", "conditional", "zero_inflated" ) ) expect_null(model_parameters(m2, effects = "fixed")$Component) expect_equal( model_parameters(m2)$Component, - c( - "conditional", "conditional", "conditional", "conditional", - "conditional" - ) + c("conditional", "conditional", "conditional", "conditional") ) expect_equal( model_parameters(m3, effects = "fixed")$Component, @@ -241,29 +234,29 @@ if (.runThisTest && test_that("model_parameters.mixed-ran_pars", { params <- model_parameters(m1, effects = "random") - expect_equal(c(nrow(params), ncol(params)), c(3, 9)) + expect_equal(c(nrow(params), ncol(params)), c(2, 9)) expect_equal( colnames(params), c("Parameter", "Coefficient", "SE", "CI", "CI_low", "CI_high", "Effects", "Group", "Component") ) expect_equal( params$Parameter, - c("SD (Intercept)", "SD (Observations)", "SD (Intercept)") + c("SD (Intercept)", "SD (Intercept)") ) expect_equal( params$Component, - c("conditional", "conditional", "zero_inflated") + c("conditional", "zero_inflated") ) expect_equal( params$Coefficient, - c(0.9312, 1, 1.17399), + c(0.9312, 1.17399), tolerance = 1e-2 ) }) test_that("model_parameters.mixed-all_pars", { params <- model_parameters(m1, effects = "all") - expect_equal(c(nrow(params), ncol(params)), c(9, 12)) + expect_equal(c(nrow(params), ncol(params)), c(8, 12)) expect_equal( colnames(params), c( @@ -275,20 +268,19 @@ if (.runThisTest && params$Parameter, c( "(Intercept)", "child", "camper1", "(Intercept)", "child", - "camper1", "SD (Intercept)", "SD (Observations)", "SD (Intercept)" + "camper1", "SD (Intercept)", "SD (Intercept)" ) ) expect_equal( params$Component, c( "conditional", "conditional", "conditional", "zero_inflated", - "zero_inflated", "zero_inflated", "conditional", "conditional", - "zero_inflated" + "zero_inflated", "zero_inflated", "conditional", "zero_inflated" ) ) expect_equal( params$Coefficient, - c(1.2628, -1.14165, 0.73354, -0.38939, 2.05407, -1.00823, 0.9312, 1, 1.17399), + c(1.2628, -1.14165, 0.73354, -0.38939, 2.05407, -1.00823, 0.9312, 1.17399), tolerance = 1e-2 ) }) @@ -359,7 +351,7 @@ if (.runThisTest && test_that("model_parameters.mixed-ran_pars", { params <- model_parameters(m4, effects = "random") - expect_equal(c(nrow(params), ncol(params)), c(7, 9)) + expect_equal(c(nrow(params), ncol(params)), c(6, 9)) expect_equal( colnames(params), c("Parameter", "Coefficient", "SE", "CI", "CI_low", "CI_high", "Effects", "Group", "Component") @@ -367,20 +359,20 @@ if (.runThisTest && expect_equal( params$Parameter, c( - "SD (Intercept)", "SD (xb)", "Cor (Intercept~xb: persons)", "SD (Observations)", - "SD (Intercept)", "SD (zg)", "Cor (Intercept~zg: persons)" + "SD (Intercept)", "SD (xb)", "Cor (Intercept~xb)", + "SD (Intercept)", "SD (zg)", "Cor (Intercept~zg)" ) ) expect_equal( params$Component, c( - "conditional", "conditional", "conditional", "conditional", + "conditional", "conditional", "conditional", "zero_inflated", "zero_inflated", "zero_inflated" ) ) expect_equal( params$Coefficient, - c(3.40563, 1.21316, -1, 1, 2.73583, 1.56833, 1), + c(3.40563, 1.21316, -1, 2.73583, 1.56833, 1), tolerance = 1e-2 ) }) @@ -403,7 +395,7 @@ if (.runThisTest && test_that("model_parameters, exp, glmmTMB", { mp1 <- model_parameters(m_exp, exponentiate = TRUE) mp2 <- model_parameters(m_exp, exponentiate = FALSE) - expect_equal(mp1$Coefficient, c(0.49271, 6.75824, 1.14541, 5.56294), tolerance = 1e-3) + expect_equal(mp1$Coefficient, c(0.49271, 6.75824, 5.56294, 1.14541), tolerance = 1e-3) expect_equal(mp1$Coefficient[3:4], mp2$Coefficient[3:4], tolerance = 1e-3) }) @@ -444,8 +436,7 @@ if (.runThisTest && "---------------------------------------------------------", "SD (Intercept: persons) | 3.41 | [ 1.67, 6.93]", "SD (xb: persons) | 1.21 | [ 0.60, 2.44]", - "Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00]", - "SD (Residual) | 1.00 | " + "Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00]" ) ) @@ -496,8 +487,7 @@ if (.runThisTest && "---------------------------------------------------------", "SD (Intercept: persons) | 3.41 | [ 1.67, 6.93]", "SD (xb: persons) | 1.21 | [ 0.60, 2.44]", - "Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00]", - "SD (Residual) | 1.00 | " + "Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00]" ) ) @@ -549,8 +539,8 @@ if (.runThisTest && "SD (Intercept: persons) | 3.41 | [ 1.67, 6.93]", "SD (xb: persons) | 1.21 | [ 0.60, 2.44]", "Cor (Intercept~xb: persons) | -1.00 | [-1.00, 1.00]", - "SD (Residual) | 1.00 | ", - "", "# Random Effects (Zero-Inflated Model)", "", "Parameter | Coefficient | 95% CI", + "", + "# Random Effects (Zero-Inflated Model)", "", "Parameter | Coefficient | 95% CI", "---------------------------------------------------------", "SD (Intercept: persons) | 2.74 | [ 1.16, 6.43]", "SD (zg: persons) | 1.57 | [ 0.64, 3.82]", @@ -589,7 +579,6 @@ if (.runThisTest && "SD (Intercept: persons) | 3.4056 | [ 1.67398, 6.92858]", "SD (xb: persons) | 1.2132 | [ 0.60210, 2.44440]", "Cor (Intercept~xb: persons) | -1.0000 | [-1.00000, 1.00000]", - "SD (Residual) | 1.0000 | ", "", "# Random Effects (Zero-Inflated Model)", "", @@ -627,8 +616,8 @@ if (.runThisTest && "SD (Intercept: persons) | 3.4056 | [ 1.67398, 6.92858]", "SD (xb: persons) | 1.2132 | [ 0.60210, 2.44440]", "Cor (Intercept~xb: persons) | -1.0000 | [-1.00000, 1.00000]", - "SD (Residual) | 1.0000 | ", - "", "# Random Effects (Zero-Inflated Model)", "", "Parameter | Coefficient | 95% CI", + "", + "# Random Effects (Zero-Inflated Model)", "", "Parameter | Coefficient | 95% CI", "---------------------------------------------------------------", "SD (Intercept: persons) | 2.7358 | [ 1.16473, 6.42622]", "SD (zg: persons) | 1.5683 | [ 0.64351, 3.82225]", @@ -669,7 +658,6 @@ if (.runThisTest && "----------------------------------------------------------------", "SD (Intercept: Session:Participant) | 0.27 | [0.08, 0.87]", "SD (Intercept: Participant) | 0.38 | [0.16, 0.92]", - "SD (Residual) | 2.06 | [1.30, 3.27]", "", "# Random Effects: zero_inflated", "", @@ -679,6 +667,39 @@ if (.runThisTest && "SD (Intercept: Participant) | 2.39 | [1.25, 4.57]" ) ) + mp <- model_parameters(model_pr, effects = "fixed", component = "all") + out <- utils::capture.output(print(mp)) + expect_equal( + out, + c("# Fixed Effects", + "", + "Parameter | Log-Mean | SE | 95% CI | z | p", + "---------------------------------------------------------------------", + "(Intercept) | 2.12 | 0.30 | [ 1.53, 2.71] | 7.05 | < .001", + "Surface [Lingual] | 0.01 | 0.29 | [-0.56, 0.58] | 0.04 | 0.971 ", + "Surface [Occlusal] | 0.54 | 0.22 | [ 0.10, 0.98] | 2.43 | 0.015 ", + "Side [Anterior] | 0.04 | 0.32 | [-0.58, 0.66] | 0.14 | 0.889 ", + "Side [Left] | -0.04 | 0.20 | [-0.44, 0.37] | -0.17 | 0.862 ", + "Jaw [Maxillar] | -0.10 | 0.21 | [-0.51, 0.30] | -0.51 | 0.612 ", + "", + "# Zero-Inflated", + "", + "Parameter | Log-Odds | SE | 95% CI | z | p", + "----------------------------------------------------------------------", + "(Intercept) | 4.87 | 0.93 | [ 3.04, 6.69] | 5.23 | < .001", + "Surface [Lingual] | 0.93 | 0.34 | [ 0.27, 1.60] | 2.75 | 0.006 ", + "Surface [Occlusal] | -1.01 | 0.29 | [-1.59, -0.44] | -3.45 | < .001", + "Side [Anterior] | -0.20 | 0.37 | [-0.93, 0.52] | -0.55 | 0.583 ", + "Side [Left] | -0.38 | 0.27 | [-0.91, 0.14] | -1.44 | 0.151 ", + "Jaw [Maxillar] | 0.59 | 0.24 | [ 0.11, 1.07] | 2.42 | 0.016 ", + "", + "# Dispersion", + "", + "Parameter | Coefficient | 95% CI", + "----------------------------------------", + "(Intercept) | 2.06 | [1.30, 3.27]" + ) + ) }) } } @@ -687,7 +708,7 @@ if (.runThisTest && test_that("model_parameters.mixed-all", { params <- model_parameters(m4, effects = "all") - expect_equal(c(nrow(params), ncol(params)), c(13, 12)) + expect_equal(c(nrow(params), ncol(params)), c(12, 12)) expect_equal( colnames(params), c( @@ -699,8 +720,8 @@ if (.runThisTest && params$Parameter, c( "(Intercept)", "child", "camper1", "(Intercept)", "child", - "camper1", "SD (Intercept)", "SD (xb)", "Cor (Intercept~xb: persons)", - "SD (Observations)", "SD (Intercept)", "SD (zg)", "Cor (Intercept~zg: persons)" + "camper1", "SD (Intercept)", "SD (xb)", "Cor (Intercept~xb)", + "SD (Intercept)", "SD (zg)", "Cor (Intercept~zg)" ) ) expect_equal( @@ -708,15 +729,14 @@ if (.runThisTest && c( "conditional", "conditional", "conditional", "zero_inflated", "zero_inflated", "zero_inflated", "conditional", "conditional", - "conditional", "conditional", "zero_inflated", "zero_inflated", - "zero_inflated" + "conditional", "zero_inflated", "zero_inflated", "zero_inflated" ) ) expect_equal( params$Coefficient, c( 2.54713, -1.08747, 0.2723, 1.88964, 0.15712, -0.17007, 3.40563, - 1.21316, -1, 1, 2.73583, 1.56833, 1 + 1.21316, -1, 2.73583, 1.56833, 1 ), tolerance = 1e-2 ) @@ -792,7 +812,6 @@ if (.runThisTest && "Parameter | Coefficient | 95% CI", "----------------------------------------------------", "SD (Intercept: persons) | 0.93 | [0.46, 1.89]", - "SD (Residual) | 1.00 | ", "", "# Random Effects (Zero-Inflated Model)", "", diff --git a/tests/testthat/test-model_parameters.blmerMod.R b/tests/testthat/test-model_parameters.blmerMod.R index cd865775c..e489074f6 100644 --- a/tests/testthat/test-model_parameters.blmerMod.R +++ b/tests/testthat/test-model_parameters.blmerMod.R @@ -22,7 +22,7 @@ if (requiet("testthat") && requiet("parameters") && requiet("blme")) { ) expect_equal( params$Parameter, - c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "Cor (Intercept~Days: Subject)", "SD (Observations)") + c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "Cor (Intercept~Days)", "SD (Observations)") ) }) } diff --git a/tests/testthat/test-model_parameters_random_pars.R b/tests/testthat/test-model_parameters_random_pars.R index 6be63db2d..7835b7760 100644 --- a/tests/testthat/test-model_parameters_random_pars.R +++ b/tests/testthat/test-model_parameters_random_pars.R @@ -23,7 +23,7 @@ if (.runThisTest && expect_equal(mp$Coefficient, as.data.frame(lme4::VarCorr(model))$sdcor, tolerance = 1e-3) expect_equal(mp$SE, c(5.83626, 1.24804, 0.31859, 1.50801), tolerance = 1e-3) expect_equal(mp$CI_low, c(15.5817, 3.91828, -0.50907, 22.80044), tolerance = 1e-3) - expect_equal(mp$Parameter, c("SD (Intercept)", "SD (Days)", "Cor (Intercept~Days: Subject)", "SD (Observations)")) + expect_equal(mp$Parameter, c("SD (Intercept)", "SD (Days)", "Cor (Intercept~Days)", "SD (Observations)")) }) model <- lmer(Reaction ~ Days + (1 + Days || Subject), data = sleepstudy) @@ -88,8 +88,8 @@ if (.runThisTest && expect_equal(mp$CI_low, c(16.7131, 21.12065, 24.1964, -0.36662, -0.59868, -0.93174, 24.18608), tolerance = 1e-3) expect_equal( mp$Parameter, - c("SD (Intercept)", "SD (Days2(3,6])", "SD (Days2(6,10])", "Cor (Intercept~Days2(3,6]: Subject)", - "Cor (Intercept~Days2(6,10]: Subject)", "Cor (Days2(3,6]~Days2(6,10]: Subject)", + c("SD (Intercept)", "SD (Days2(3,6])", "SD (Days2(6,10])", "Cor (Intercept~Days2(3,6])", + "Cor (Intercept~Days2(6,10])", "Cor (Days2(3,6]~Days2(6,10])", "SD (Observations)") ) }) @@ -104,8 +104,8 @@ if (.runThisTest && expect_equal( mp$Parameter, c("SD (Days2(-1,3])", "SD (Days2(3,6])", "SD (Days2(6,10])", - "Cor (Days2(-1,3]~Days2(3,6]: Subject)", "Cor (Days2(-1,3]~Days2(6,10]: Subject)", - "Cor (Days2(3,6]~Days2(6,10]: Subject)", "SD (Observations)") + "Cor (Days2(-1,3]~Days2(3,6])", "Cor (Days2(-1,3]~Days2(6,10])", + "Cor (Days2(3,6]~Days2(6,10])", "SD (Observations)") ) }) @@ -113,12 +113,13 @@ if (.runThisTest && mp <- model_parameters(model, effects = "random") test_that("model_parameters-random pars 9", { - expect_equal(mp$Coefficient, as.data.frame(lme4::VarCorr(model))$sdcor[-6:-5], tolerance = 1e-3) + expect_equal(mp$Coefficient, as.data.frame(lme4::VarCorr(model))$sdcor, tolerance = 1e-3) expect_true(all(is.na(mp$SE))) expect_equal( mp$Parameter, c("SD (Intercept)", "SD (Days2(-1,3])", "SD (Days2(3,6])", "SD (Days2(6,10])", - "Cor (Days2(3,6]~Days2(6,10]: Subject.1)", "SD (Observations)") + "Cor (Days2(-1,3]~Days2(3,6])", "Cor (Days2(-1,3]~Days2(6,10])", + "Cor (Days2(3,6]~Days2(6,10])", "SD (Observations)") ) }) @@ -132,8 +133,8 @@ if (.runThisTest && expect_equal( mp$Parameter, c("SD (Days2(-1,3])", "SD (Days2(3,6])", "SD (Days2(6,10])", - "Cor (Days2(-1,3]~Days2(3,6]: Subject)", "Cor (Days2(-1,3]~Days2(6,10]: Subject)", - "Cor (Days2(3,6]~Days2(6,10]: Subject)", "SD (Observations)") + "Cor (Days2(-1,3]~Days2(3,6])", "Cor (Days2(-1,3]~Days2(6,10])", + "Cor (Days2(3,6]~Days2(6,10])", "SD (Observations)") ) }) } diff --git a/tests/testthat/test-random_effects_ci-glmmTMB.R b/tests/testthat/test-random_effects_ci-glmmTMB.R index d4a01fb0c..0ada1aebf 100644 --- a/tests/testthat/test-random_effects_ci-glmmTMB.R +++ b/tests/testthat/test-random_effects_ci-glmmTMB.R @@ -1,10 +1,14 @@ -.runThisTest <- Sys.getenv("RunAllparametersTests") == "yes" +# Test Setup -------------------------------- -osx <- tryCatch( +# only run for dev-version, not on CRAN +.runThisTest <- length(strsplit(packageDescription("parameters")$Version, "\\.")[[1]]) > 3 + +# only run on windows - there are some minor rounding issues on other OS +win_os <- tryCatch( { si <- Sys.info() if (!is.null(si["sysname"])) { - si["sysname"] == "Darwin" || grepl("^darwin", R.version$os) + si["sysname"] == "Windows" || grepl("^mingw", R.version$os) } else { FALSE } @@ -14,26 +18,33 @@ osx <- tryCatch( } ) + + +# tests -------------------------------- + ## TODO also check messages for profiled CI -if (.runThisTest && !osx && +if (.runThisTest && win_os && requiet("testthat") && requiet("parameters") && requiet("glmmTMB") && - requiet("lme4")) { + requiet("lme4") && + packageVersion("glmmTMB") > "1.1.3") { data(sleepstudy) data(cake) set.seed(123) sleepstudy$Months <- sample(1:4, nrow(sleepstudy), TRUE) + set.seed(123) m1 <- suppressWarnings(glmmTMB(angle ~ temperature + (temperature | recipe) + (temperature | replicate), data = cake)) m2 <- glmmTMB(Reaction ~ Days + (Days | Subject), data = sleepstudy) m3 <- suppressWarnings(glmmTMB(angle ~ temperature + (temperature | recipe), data = cake)) m4 <- suppressWarnings(glmmTMB(angle ~ temperature + (temperature | replicate), data = cake)) m5 <- suppressWarnings(glmmTMB(Reaction ~ Days + (Days + Months | Subject), data = sleepstudy)) - expect_message(expect_message(mp1 <- model_parameters(m1), "singularity"), "singularity") + set.seed(123) + expect_message(mp1 <- model_parameters(m1), "singularity") mp2 <- model_parameters(m2) # works expect_message(mp3 <- model_parameters(m3), "singularity") # no SE/CI expect_message(mp4 <- model_parameters(m4), "singularity") # no SE/CI @@ -42,11 +53,13 @@ if (.runThisTest && !osx && test_that("random effects CIs, two slopes, categorical", { expect_equal( mp1$CI_low, - c(28.9123, 5.03115, -1.87304, -2.42081, -3.2708, -2.57695, NA, - NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, - NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, - NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), - tolerance = 1e-3, + c(28.9123, 5.03115, -1.87304, -2.42081, -3.2708, -2.57695, 0.21571, + 4.17466, NaN, 0, 0.26247, 0.34089, 0.02477, 0.65731, 0.3902, + 0.14685, 0.01322, 0.62182, 0.99915, NaN, NaN, NaN, -0.31609, + -0.48806, NaN, -0.8346, NaN, -0.60153, NA, NA, NA, NA, NA, NA, + NA, NA, NA, NaN, NA, NA, NA, NA, NA, NA, NA, NA, NA, NaN, 4.12529 + ), + tolerance = 1e-2, ignore_attr = TRUE ) @@ -57,22 +70,34 @@ if (.runThisTest && !osx && "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", "SD (temperature^5)", "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", - "SD (temperature^5)", "Cor (Intercept~temperature.C: recipe)", - "Cor (Intercept~temperature.C: replicate)", "Cor (Intercept~temperature.L: recipe)", - "Cor (Intercept~temperature.L: replicate)", "Cor (Intercept~temperature.Q: recipe)", - "Cor (Intercept~temperature.Q: replicate)", "Cor (Intercept~temperature^4: recipe)", - "Cor (Intercept~temperature^4: replicate)", "Cor (Intercept~temperature^5: recipe)", - "Cor (Intercept~temperature^5: replicate)", "Cor (temperature.L~temperature.C: recipe)", - "Cor (temperature.Q~temperature.C: recipe)", "Cor (temperature.L~temperature.Q: recipe)", - "Cor (temperature.L~temperature^4: recipe)", "Cor (temperature.Q~temperature^4: recipe)", - "Cor (temperature.C~temperature^4: recipe)", "Cor (temperature.L~temperature^5: recipe)", - "Cor (temperature.Q~temperature^5: recipe)", "Cor (temperature.C~temperature^5: recipe)", - "Cor (temperature^4~temperature^5: recipe)", "Cor (temperature.L~temperature.C: replicate)", - "Cor (temperature.Q~temperature.C: replicate)", "Cor (temperature.L~temperature.Q: replicate)", - "Cor (temperature.L~temperature^4: replicate)", "Cor (temperature.Q~temperature^4: replicate)", - "Cor (temperature.C~temperature^4: replicate)", "Cor (temperature.L~temperature^5: replicate)", - "Cor (temperature.Q~temperature^5: replicate)", "Cor (temperature.C~temperature^5: replicate)", - "Cor (temperature^4~temperature^5: replicate)", "SD (Observations)" + "SD (temperature^5)", "Cor (Intercept~temperature.L)", "Cor (Intercept~temperature.Q)", + "Cor (Intercept~temperature.C)", "Cor (Intercept~temperature^4)", + "Cor (Intercept~temperature^5)", "Cor (Intercept~temperature.L)", + "Cor (Intercept~temperature.Q)", "Cor (Intercept~temperature.C)", + "Cor (Intercept~temperature^4)", "Cor (Intercept~temperature^5)", + "Cor (temperature.L~temperature.Q)", "Cor (temperature.L~temperature.C)", + "Cor (temperature.L~temperature^4)", "Cor (temperature.L~temperature^5)", + "Cor (temperature.Q~temperature.C)", "Cor (temperature.Q~temperature^4)", + "Cor (temperature.Q~temperature^5)", "Cor (temperature.C~temperature^4)", + "Cor (temperature.C~temperature^5)", "Cor (temperature^4~temperature^5)", + "Cor (temperature.L~temperature.Q)", "Cor (temperature.L~temperature.C)", + "Cor (temperature.L~temperature^4)", "Cor (temperature.L~temperature^5)", + "Cor (temperature.Q~temperature.C)", "Cor (temperature.Q~temperature^4)", + "Cor (temperature.Q~temperature^5)", "Cor (temperature.C~temperature^4)", + "Cor (temperature.C~temperature^5)", "Cor (temperature^4~temperature^5)", + "SD (Observations)") + ) + + expect_equal( + mp1$Group, + c("", "", "", "", "", "", "recipe", "replicate", "recipe", "recipe", + "recipe", "recipe", "recipe", "replicate", "replicate", "replicate", + "replicate", "replicate", "recipe", "recipe", "recipe", "recipe", + "recipe", "replicate", "replicate", "replicate", "replicate", + "replicate", "recipe", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "recipe", "recipe", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "Residual" ) ) }) @@ -82,13 +107,13 @@ if (.runThisTest && !osx && expect_equal( mp2$CI_low, c(238.40611, 7.52295, 15.01709, 3.80546, -0.48781, 22.80047), - tolerance = 1e-3, + tolerance = 1e-2, ignore_attr = TRUE ) expect_equal( mp2$Parameter, - c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "Cor (Intercept~Days: Subject)", + c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "Cor (Intercept~Days)", "SD (Observations)") ) }) @@ -98,9 +123,9 @@ if (.runThisTest && !osx && expect_equal( mp3$CI_low, c(31.20278, 4.35879, -2.63767, -2.80041, -3.54983, -3.16627, - NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, - NA, NA, NA, NA, NA, NA), - tolerance = 1e-3, + 0, NaN, NaN, 0, NaN, NaN, NaN, NaN, -0.49203, -0.41167, NaN, + NA, NA, NA, NA, NA, NA, NA, NA, NA, NaN, 7.08478), + tolerance = 1e-2, ignore_attr = TRUE ) @@ -109,15 +134,22 @@ if (.runThisTest && !osx && c("(Intercept)", "temperature.L", "temperature.Q", "temperature.C", "temperature^4", "temperature^5", "SD (Intercept)", "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", - "SD (temperature^5)", "Cor (Intercept~temperature.L: recipe)", - "Cor (Intercept~temperature.Q: recipe)", "Cor (Intercept~temperature.C: recipe)", - "Cor (Intercept~temperature^4: recipe)", "Cor (Intercept~temperature^5: recipe)", - "Cor (temperature.L~temperature.C: recipe)", "Cor (temperature.Q~temperature.C: recipe)", - "Cor (temperature.L~temperature.Q: recipe)", "Cor (temperature.L~temperature^4: recipe)", - "Cor (temperature.Q~temperature^4: recipe)", "Cor (temperature.C~temperature^4: recipe)", - "Cor (temperature.L~temperature^5: recipe)", "Cor (temperature.Q~temperature^5: recipe)", - "Cor (temperature.C~temperature^5: recipe)", "Cor (temperature^4~temperature^5: recipe)", - "SD (Observations)") + "SD (temperature^5)", "Cor (Intercept~temperature.L)", "Cor (Intercept~temperature.Q)", + "Cor (Intercept~temperature.C)", "Cor (Intercept~temperature^4)", + "Cor (Intercept~temperature^5)", "Cor (temperature.L~temperature.Q)", + "Cor (temperature.L~temperature.C)", "Cor (temperature.L~temperature^4)", + "Cor (temperature.L~temperature^5)", "Cor (temperature.Q~temperature.C)", + "Cor (temperature.Q~temperature^4)", "Cor (temperature.Q~temperature^5)", + "Cor (temperature.C~temperature^4)", "Cor (temperature.C~temperature^5)", + "Cor (temperature^4~temperature^5)", "SD (Observations)") + ) + + expect_equal( + mp3$Group, + c("", "", "", "", "", "", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "Residual") ) }) @@ -126,9 +158,10 @@ if (.runThisTest && !osx && expect_equal( mp4$CI_low, c(29.01131, 5.01247, -1.89444, -1.96275, -2.66798, -2.50892, - NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, - NA, NA, NA, NA, NA, NA), - tolerance = 1e-3, + 4.23497, 0.62985, 0.36934, 0.1398, 0.01133, 0.60758, 0.56678, + 0.26866, NaN, NaN, NaN, NA, NA, NA, NA, NA, NA, NA, NA, NA, NaN, + 4.23582), + tolerance = 1e-2, ignore_attr = TRUE ) @@ -137,31 +170,41 @@ if (.runThisTest && !osx && c("(Intercept)", "temperature.L", "temperature.Q", "temperature.C", "temperature^4", "temperature^5", "SD (Intercept)", "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", - "SD (temperature^5)", "Cor (Intercept~temperature.L: replicate)", - "Cor (Intercept~temperature.Q: replicate)", "Cor (Intercept~temperature.C: replicate)", - "Cor (Intercept~temperature^4: replicate)", "Cor (Intercept~temperature^5: replicate)", - "Cor (temperature.L~temperature.C: replicate)", "Cor (temperature.Q~temperature.C: replicate)", - "Cor (temperature.L~temperature.Q: replicate)", "Cor (temperature.L~temperature^4: replicate)", - "Cor (temperature.Q~temperature^4: replicate)", "Cor (temperature.C~temperature^4: replicate)", - "Cor (temperature.L~temperature^5: replicate)", "Cor (temperature.Q~temperature^5: replicate)", - "Cor (temperature.C~temperature^5: replicate)", "Cor (temperature^4~temperature^5: replicate)", - "SD (Observations)") + "SD (temperature^5)", "Cor (Intercept~temperature.L)", "Cor (Intercept~temperature.Q)", + "Cor (Intercept~temperature.C)", "Cor (Intercept~temperature^4)", + "Cor (Intercept~temperature^5)", "Cor (temperature.L~temperature.Q)", + "Cor (temperature.L~temperature.C)", "Cor (temperature.L~temperature^4)", + "Cor (temperature.L~temperature^5)", "Cor (temperature.Q~temperature.C)", + "Cor (temperature.Q~temperature^4)", "Cor (temperature.Q~temperature^5)", + "Cor (temperature.C~temperature^4)", "Cor (temperature.C~temperature^5)", + "Cor (temperature^4~temperature^5)", "SD (Observations)") + ) + + expect_equal( + mp4$Group, + c("", "", "", "", "", "", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "Residual") ) }) + test_that("random effects CIs, double slope", { expect_equal( mp5$CI_low, - c(238.40607, 7.52296, NA, NA, NA, NA, NA, NA, NA), - tolerance = 1e-3, + c(238.40607, 7.52296, 15.01708, 3.80547, NaN, -0.48781, NaN, + NaN, 22.80046), + tolerance = 1e-2, ignore_attr = TRUE ) expect_equal( mp5$Parameter, c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "SD (Months)", - "Cor (Intercept~Days: Subject)", "Cor (Intercept~Months: Subject)", - "Cor (Days~Months: Subject)", "SD (Observations)") + "Cor (Intercept~Days)", "Cor (Intercept~Months)", + "Cor (Days~Months)", "SD (Observations)") ) }) @@ -174,9 +217,11 @@ if (.runThisTest && !osx && set.seed(123) sleepstudy$Months <- sample(1:4, nrow(sleepstudy), TRUE) + set.seed(123) m2 <- glmmTMB(Reaction ~ Days + (0 + Days | Subject), data = sleepstudy) m5 <- suppressWarnings(glmmTMB(Reaction ~ Days + (0 + Days + Months | Subject), data = sleepstudy)) + set.seed(123) mp2 <- model_parameters(m2) expect_message(mp5 <- model_parameters(m5), "singularity") # no SE/CI @@ -184,7 +229,7 @@ if (.runThisTest && !osx && expect_equal( mp2$CI_low, c(243.55046, 6.89554, 4.98429, 25.94359), - tolerance = 1e-3, + tolerance = 1e-2, ignore_attr = TRUE ) @@ -197,14 +242,14 @@ if (.runThisTest && !osx && test_that("random effects CIs, simple slope", { expect_equal( mp5$CI_low, - c(237.03695, 9.04139, NA, NA, NA, NA), - tolerance = 1e-3, + c(237.03695, 9.04139, NaN, 8.95755, NaN, 30.67054), + tolerance = 1e-2, ignore_attr = TRUE ) expect_equal( mp5$Parameter, - c("(Intercept)", "Days", "SD (Days)", "SD (Months)", "Cor (Days~Months: Subject)", + c("(Intercept)", "Days", "SD (Days)", "SD (Months)", "Cor (Days~Months)", "SD (Observations)") ) }) diff --git a/tests/testthat/test-random_effects_ci.R b/tests/testthat/test-random_effects_ci.R index 8df3b91db..7c97d072f 100644 --- a/tests/testthat/test-random_effects_ci.R +++ b/tests/testthat/test-random_effects_ci.R @@ -38,15 +38,18 @@ if (.runThisTest && !osx && expect_message(mp4 <- model_parameters(m4), "meaningful") expect_message(mp5 <- model_parameters(m5), "meaningful") + + # model 1 --------------------- + test_that("random effects CIs, two slopes, categorical", { expect_equal( mp1$CI_low, c(28.75568, 4.97893, -1.95002, -2.69995, -3.62201, -2.69102, 4.28558, 0.21474, 0.40062, 0.10169, 0.04953, 1e-05, 0.55398, - 0, 2e-05, 0.6333, 1.09851, 0.00944, -1, -1, -0.65406, -1, -0.69103, - -1, -0.95271, -1, -0.90617, -1, -0.99802, -0.99836, -1, -1, -1, - -1, -0.75274, -0.96895, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, - -1, -1, 4.07985), + 0, 2e-05, 0.6333, 1.09851, 0.00944, -0.65406, -0.69103, -1, -0.95271, + -0.90617, -1, -1, -1, -1, -1, -1, -0.99802, -1, -0.75274, -0.99836, + -1, -0.96895, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, 4.07985), tolerance = 1e-3, ignore_attr = TRUE ) @@ -58,27 +61,41 @@ if (.runThisTest && !osx && "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", "SD (temperature^5)", "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", - "SD (temperature^5)", "Cor (Intercept~temperature.C: replicate)", - "Cor (Intercept~temperature.C: recipe)", "Cor (Intercept~temperature.L: replicate)", - "Cor (Intercept~temperature.L: recipe)", "Cor (Intercept~temperature.Q: replicate)", - "Cor (Intercept~temperature.Q: recipe)", "Cor (Intercept~temperature^4: replicate)", - "Cor (Intercept~temperature^4: recipe)", "Cor (Intercept~temperature^5: replicate)", - "Cor (Intercept~temperature^5: recipe)", "Cor (temperature.L~temperature.C: replicate)", - "Cor (temperature.Q~temperature.C: replicate)", "Cor (temperature.L~temperature.Q: replicate)", - "Cor (temperature.L~temperature^4: replicate)", "Cor (temperature.Q~temperature^4: replicate)", - "Cor (temperature.C~temperature^4: replicate)", "Cor (temperature.L~temperature^5: replicate)", - "Cor (temperature.Q~temperature^5: replicate)", "Cor (temperature.C~temperature^5: replicate)", - "Cor (temperature^4~temperature^5: replicate)", "Cor (temperature.L~temperature.C: recipe)", - "Cor (temperature.Q~temperature.C: recipe)", "Cor (temperature.L~temperature.Q: recipe)", - "Cor (temperature.L~temperature^4: recipe)", "Cor (temperature.Q~temperature^4: recipe)", - "Cor (temperature.C~temperature^4: recipe)", "Cor (temperature.L~temperature^5: recipe)", - "Cor (temperature.Q~temperature^5: recipe)", "Cor (temperature.C~temperature^5: recipe)", - "Cor (temperature^4~temperature^5: recipe)", "SD (Observations)" - ) + "SD (temperature^5)", "Cor (Intercept~temperature.L)", "Cor (Intercept~temperature.Q)", + "Cor (Intercept~temperature.C)", "Cor (Intercept~temperature^4)", + "Cor (Intercept~temperature^5)", "Cor (Intercept~temperature.L)", + "Cor (Intercept~temperature.Q)", "Cor (Intercept~temperature.C)", + "Cor (Intercept~temperature^4)", "Cor (Intercept~temperature^5)", + "Cor (temperature.L~temperature.Q)", "Cor (temperature.L~temperature.C)", + "Cor (temperature.L~temperature^4)", "Cor (temperature.L~temperature^5)", + "Cor (temperature.Q~temperature.C)", "Cor (temperature.Q~temperature^4)", + "Cor (temperature.Q~temperature^5)", "Cor (temperature.C~temperature^4)", + "Cor (temperature.C~temperature^5)", "Cor (temperature^4~temperature^5)", + "Cor (temperature.L~temperature.Q)", "Cor (temperature.L~temperature.C)", + "Cor (temperature.L~temperature^4)", "Cor (temperature.L~temperature^5)", + "Cor (temperature.Q~temperature.C)", "Cor (temperature.Q~temperature^4)", + "Cor (temperature.Q~temperature^5)", "Cor (temperature.C~temperature^4)", + "Cor (temperature.C~temperature^5)", "Cor (temperature^4~temperature^5)", + "SD (Observations)") + ) + + expect_equal( + mp1$Group, + c("", "", "", "", "", "", "replicate", "recipe", "replicate", + "replicate", "replicate", "replicate", "replicate", "recipe", + "recipe", "recipe", "recipe", "recipe", "replicate", "replicate", + "replicate", "replicate", "replicate", "recipe", "recipe", "recipe", + "recipe", "recipe", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "recipe", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "recipe", "recipe", "Residual") ) }) + + # model 2 --------------------- + test_that("random effects CIs, simple slope", { expect_equal( mp2$CI_low, @@ -89,12 +106,20 @@ if (.runThisTest && !osx && expect_equal( mp2$Parameter, - c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "Cor (Intercept~Days: Subject)", + c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "Cor (Intercept~Days)", "SD (Observations)") ) + + expect_equal( + mp2$Group, + c("", "", "Subject", "Subject", "Subject", "Residual") + ) }) + + # model 3 --------------------- + test_that("random effects CIs, categorical slope-1", { expect_equal( mp3$CI_low, @@ -110,26 +135,36 @@ if (.runThisTest && !osx && c("(Intercept)", "temperature.L", "temperature.Q", "temperature.C", "temperature^4", "temperature^5", "SD (Intercept)", "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", - "SD (temperature^5)", "Cor (Intercept~temperature.L: recipe)", - "Cor (Intercept~temperature.Q: recipe)", "Cor (Intercept~temperature.C: recipe)", - "Cor (Intercept~temperature^4: recipe)", "Cor (Intercept~temperature^5: recipe)", - "Cor (temperature.L~temperature.C: recipe)", "Cor (temperature.Q~temperature.C: recipe)", - "Cor (temperature.L~temperature.Q: recipe)", "Cor (temperature.L~temperature^4: recipe)", - "Cor (temperature.Q~temperature^4: recipe)", "Cor (temperature.C~temperature^4: recipe)", - "Cor (temperature.L~temperature^5: recipe)", "Cor (temperature.Q~temperature^5: recipe)", - "Cor (temperature.C~temperature^5: recipe)", "Cor (temperature^4~temperature^5: recipe)", - "SD (Observations)") + "SD (temperature^5)", "Cor (Intercept~temperature.L)", "Cor (Intercept~temperature.Q)", + "Cor (Intercept~temperature.C)", "Cor (Intercept~temperature^4)", + "Cor (Intercept~temperature^5)", "Cor (temperature.L~temperature.Q)", + "Cor (temperature.L~temperature.C)", "Cor (temperature.L~temperature^4)", + "Cor (temperature.L~temperature^5)", "Cor (temperature.Q~temperature.C)", + "Cor (temperature.Q~temperature^4)", "Cor (temperature.Q~temperature^5)", + "Cor (temperature.C~temperature^4)", "Cor (temperature.C~temperature^5)", + "Cor (temperature^4~temperature^5)", "SD (Observations)") + ) + + expect_equal( + mp3$Group, + c("", "", "", "", "", "", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "Residual") ) }) + + # model 4 --------------------- + test_that("random effects CIs, categorical slope-2", { expect_equal( mp4$CI_low, c(28.88523, 4.96796, -1.93239, -1.98597, -2.68858, -2.5524, 4.27899, 0.35378, 0.08109, 0.03419, 0, 0.49982, -0.68893, -0.71984, -1, - -0.96725, -0.92158, -0.99894, -0.99924, -1, -1, -1, -1, -0.80378, - -0.9778, -1, -1, 4.21143), + -0.96725, -0.92158, -1, -0.99894, -1, -0.80378, -0.99924, -1, + -0.9778, -1, -1, -1, 4.21143), tolerance = 1e-3, ignore_attr = TRUE ) @@ -139,18 +174,30 @@ if (.runThisTest && !osx && c("(Intercept)", "temperature.L", "temperature.Q", "temperature.C", "temperature^4", "temperature^5", "SD (Intercept)", "SD (temperature.L)", "SD (temperature.Q)", "SD (temperature.C)", "SD (temperature^4)", - "SD (temperature^5)", "Cor (Intercept~temperature.L: replicate)", - "Cor (Intercept~temperature.Q: replicate)", "Cor (Intercept~temperature.C: replicate)", - "Cor (Intercept~temperature^4: replicate)", "Cor (Intercept~temperature^5: replicate)", - "Cor (temperature.L~temperature.C: replicate)", "Cor (temperature.Q~temperature.C: replicate)", - "Cor (temperature.L~temperature.Q: replicate)", "Cor (temperature.L~temperature^4: replicate)", - "Cor (temperature.Q~temperature^4: replicate)", "Cor (temperature.C~temperature^4: replicate)", - "Cor (temperature.L~temperature^5: replicate)", "Cor (temperature.Q~temperature^5: replicate)", - "Cor (temperature.C~temperature^5: replicate)", "Cor (temperature^4~temperature^5: replicate)", - "SD (Observations)") + "SD (temperature^5)", "Cor (Intercept~temperature.L)", "Cor (Intercept~temperature.Q)", + "Cor (Intercept~temperature.C)", "Cor (Intercept~temperature^4)", + "Cor (Intercept~temperature^5)", "Cor (temperature.L~temperature.Q)", + "Cor (temperature.L~temperature.C)", "Cor (temperature.L~temperature^4)", + "Cor (temperature.L~temperature^5)", "Cor (temperature.Q~temperature.C)", + "Cor (temperature.Q~temperature^4)", "Cor (temperature.Q~temperature^5)", + "Cor (temperature.C~temperature^4)", "Cor (temperature.C~temperature^5)", + "Cor (temperature^4~temperature^5)", "SD (Observations)") + ) + + expect_equal( + mp4$Group, + c("", "", "", "", "", "", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "replicate", "replicate", + "replicate", "replicate", "replicate", "Residual") ) }) + + + # model 5 --------------------- + test_that("random effects CIs, double slope", { expect_equal( mp5$CI_low, @@ -162,15 +209,21 @@ if (.runThisTest && !osx && expect_equal( mp5$Parameter, c("(Intercept)", "Days", "SD (Intercept)", "SD (Days)", "SD (Months)", - "Cor (Intercept~Days: Subject)", "Cor (Intercept~Months: Subject)", - "Cor (Days~Months: Subject)", "SD (Observations)") + "Cor (Intercept~Days)", "Cor (Intercept~Months)", + "Cor (Days~Months)", "SD (Observations)") ) - }) + expect_equal( + mp5$Group, + c("", "", "Subject", "Subject", "Subject", "Subject", "Subject", + "Subject", "Residual") + ) + }) + # no random intercept -------------------------- data(sleepstudy) set.seed(123) @@ -206,12 +259,16 @@ if (.runThisTest && !osx && expect_equal( mp5$Parameter, - c("(Intercept)", "Days", "SD (Days)", "SD (Months)", "Cor (Days~Months: Subject)", + c("(Intercept)", "Days", "SD (Days)", "SD (Months)", "Cor (Days~Months)", "SD (Observations)") ) }) + + + # poly random slope -------------------------- + data(cake) m <- lmer(angle ~ poly(temp, 2) + (poly(temp, 2) | replicate) + (1 | recipe), data = cake) mp <- model_parameters(m) @@ -229,10 +286,59 @@ if (.runThisTest && !osx && mp$Parameter, c("(Intercept)", "poly(temp, 2)1", "poly(temp, 2)2", "SD (Intercept)", "SD (Intercept)", "SD (poly(temp, 2)1)", "SD (poly(temp, 2)2)", - "Cor (Intercept~poly(temp, 2)1: replicate)", "Cor (Intercept~poly(temp, 2)2: replicate)", - "Cor (poly(temp, 2)1~poly(temp, 2)2: replicate)", "SD (Observations)" + "Cor (Intercept~poly(temp, 2)1)", "Cor (Intercept~poly(temp, 2)2)", + "Cor (poly(temp, 2)1~poly(temp, 2)2)", "SD (Observations)" ) ) }) + + + + # poly and categorical random slope -------------------------- + + m <- lmer(angle ~ poly(temp, 2) + (poly(temp, 2) | replicate) + (temperature | recipe), + data = cake) + mp <- model_parameters(m, effects = "random") + + test_that("random effects CIs, poly categorical slope", { + ## NOTE check back every now and then and see if tests still work + skip("works interactively") + + expect_equal( + mp$CI_low, + c(4.27846, 0.22005, 8.22659, 1.17579, 0, 5e-05, 0.37736, 1.24258, + 0, -0.77207, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, 4.22056), + tolerance = 1e-3, + ignore_attr = TRUE + ) + + expect_equal( + mp$Parameter, + c("SD (Intercept)", "SD (Intercept)", "SD (poly(temp, 2)1)", + "SD (poly(temp, 2)2)", "SD (temperature.L)", "SD (temperature.Q)", + "SD (temperature.C)", "SD (temperature^4)", "SD (temperature^5)", + "Cor (Intercept~poly(temp, 2)1)", "Cor (Intercept~poly(temp, 2)2)", + "Cor (Intercept~temperature.L)", "Cor (Intercept~temperature.Q)", + "Cor (Intercept~temperature.C)", "Cor (Intercept~temperature^4)", + "Cor (Intercept~temperature^5)", "Cor (poly(temp, 2)1~poly(temp, 2)2)", + "Cor (temperature.L~temperature.Q)", "Cor (temperature.L~temperature.C)", + "Cor (temperature.L~temperature^4)", "Cor (temperature.L~temperature^5)", + "Cor (temperature.Q~temperature.C)", "Cor (temperature.Q~temperature^4)", + "Cor (temperature.Q~temperature^5)", "Cor (temperature.C~temperature^4)", + "Cor (temperature.C~temperature^5)", "Cor (temperature^4~temperature^5)", + "SD (Observations)") + ) + + expect_equal( + mp$Group, + c("replicate", "recipe", "replicate", "replicate", "recipe", + "recipe", "recipe", "recipe", "recipe", "replicate", "replicate", + "recipe", "recipe", "recipe", "recipe", "recipe", "replicate", + "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", "recipe", + "recipe", "recipe", "recipe", "Residual") + ) + }) + }