diff --git a/NAMESPACE b/NAMESPACE index 622cb423..3e1109af 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,8 @@ S3method(describe_nonlinear,data.frame) S3method(describe_nonlinear,estimate_predicted) S3method(describe_nonlinear,numeric) +S3method(estimate_contrasts,default) +S3method(estimate_contrasts,estimate_predicted) S3method(format,estimate_contrasts) S3method(format,estimate_grouplevel) S3method(format,estimate_means) diff --git a/R/estimate_contrast_methods.R b/R/estimate_contrast_methods.R new file mode 100644 index 00000000..ea543307 --- /dev/null +++ b/R/estimate_contrast_methods.R @@ -0,0 +1,350 @@ +#' @export +estimate_contrasts.estimate_predicted <- function(model, + contrast = NULL, + by = NULL, + predict = "response", + ci = 0.95, + p_adjust = "none", + comparison = "pairwise", + verbose = TRUE, + ...) { + # sanity check + if (inherits(comparison, "formula")) { + comparison <- all.vars(comparison)[1] + } + comparison <- insight::validate_argument(comparison, c("pairwise", "interaction")) + + # sanity check + if (is.null(contrast)) { + insight::format_error("Argument `contrast` must be specified and cannot be `NULL`.") + } + + # the "model" object is an object of class "estimate_predicted", we want + # to copy that into a separate object, for clearer names + predictions <- object <- model + model <- attributes(object)$model + datagrid <- attributes(object)$datagrid + + # vcov matrix, for adjusting se + vcov_matrix <- .safe(stats::vcov(model, verbose = FALSE, ...)) + + minfo <- insight::model_info(model) + + # model df + dof <- insight::get_df(model, type = "wald", verbose = FALSE) + crit_factor <- (1 + ci) / 2 + + ## TODO: For Bayesian models, we always use the returned standard errors + # need to check whether scale is always correct + + # for non-Gaussian models, we need to adjust the standard errors + if (!minfo$is_linear && !minfo$is_bayesian) { + se_from_predictions <- tryCatch( + { + # arguments for predict(), to get SE on response scale for non-Gaussian models + my_args <- list( + model, + newdata = datagrid, + type = predict, + se.fit = TRUE + ) + # for mixed models, need to set re.form to NULL or NA + if (insight::is_mixed_model(model)) { + my_args$re.form <- NULL + } + do.call(stats::predict, my_args) + }, + error = function(e) { + e + } + ) + # check if everything worked as expected + if (inherits(se_from_predictions, "error")) { + insight::format_error( + "This model (family) is probably not supported. The error that occured was:", + se_from_predictions$message + ) + } + # check if we have standard errors + if (is.null(se_from_predictions$se.fit)) { + insight::format_error("Could not extract standard errors from predictions.") + } + preds_with_se <- merge( + predictions, + cbind(datagrid, se_prob = se_from_predictions$se.fit), + sort = FALSE, + all = TRUE + ) + predictions$SE <- preds_with_se$se_prob + } else { + # for linear models, we don't need adjustment of standard errors + vcov_matrix <- NULL + } + + # compute contrasts or comparisons + out <- switch(comparison, + pairwise = .compute_comparisons(predictions, dof, vcov_matrix, datagrid, contrast, by, crit_factor), + interaction = .compute_interactions(predictions, dof, vcov_matrix, datagrid, contrast, by, crit_factor) + ) + + # restore attributes, for formatting + info <- attributes(object) + attributes(out) <- utils::modifyList(attributes(out), info[.info_elements()]) + + # overwrite some of the attributes + attr(out, "contrast") <- contrast + attr(out, "focal_terms") <- c(contrast, by) + attr(out, "by") <- by + + # format output + out <- format.marginaleffects_contrasts(out, model, p_adjust, comparison, ...) + + # p-value adjustment? + if (!is.null(p_adjust)) { + out <- .p_adjust(model, out, p_adjust, verbose, ...) + } + + # Table formatting + attr(out, "table_title") <- c("Model-based Contrasts Analysis", "blue") + attr(out, "table_footer") <- .table_footer( + out, + by = contrast, + type = "contrasts", + model = model, + info = info + ) + + # Add attributes + attr(out, "model") <- model + attr(out, "response") <- insight::find_response(model) + attr(out, "ci") <- ci + attr(out, "p_adjust") <- p_adjust + + # add attributes from workhorse function + attributes(out) <- utils::modifyList(attributes(out), info[.info_elements()]) + + # Output + class(out) <- unique(c("estimate_contrasts", "see_estimate_contrasts", class(out))) + out +} + + +# pairwise comparisons ---------------------------------------------------- +.compute_comparisons <- function(predictions, dof, vcov_matrix, datagrid, contrast, by, crit_factor) { + # we need the focal terms and all unique values from the datagrid + focal_terms <- c(contrast, by) + at_list <- lapply(datagrid[focal_terms], unique) + + # pairwise comparisons are a bit more complicated, as we need to create + # pairwise combinations of the levels of the focal terms. + + # since we split at "." later, we need to replace "." in all levels + # with a unique character combination + at_list <- lapply(at_list, function(i) { + gsub(".", "#_#", as.character(i), fixed = TRUE) + }) + # create pairwise combinations + level_pairs <- interaction(expand.grid(at_list)) + # using the matrix and then removing the lower triangle, we get all + # pairwise combinations, except the ones that are the same + M <- matrix( + 1, + nrow = length(level_pairs), + ncol = length(level_pairs), + dimnames = list(level_pairs, level_pairs) + ) + M[!upper.tri(M)] <- NA + # table() works fine to create variables of this pairwise combinations + pairs_data <- stats::na.omit(as.data.frame(as.table(M))) + pairs_data$Freq <- NULL + pairs_data <- lapply(pairs_data, as.character) + # the levels are combined by ".", we need to split them and then create + # a list of data frames, where each data frames contains the levels of + # the focal terms as variables + pairs_data <- lapply(pairs_data, function(i) { + # split at ".", which is the separator char for levels + pair <- strsplit(i, ".", fixed = TRUE) + # since we replaced "." with "#_#" in original levels, + # we need to replace it back here + pair <- lapply(pair, gsub, pattern = "#_#", replacement = ".", fixed = TRUE) + datawizard::data_rotate(as.data.frame(pair)) + }) + # now we iterate over all pairs and try to find the corresponding predictions + out <- do.call(rbind, lapply(seq_len(nrow(pairs_data[[1]])), function(i) { + pos1 <- predictions[[focal_terms[1]]] == pairs_data[[1]][i, 1] + pos2 <- predictions[[focal_terms[1]]] == pairs_data[[2]][i, 1] + + if (length(focal_terms) > 1) { + pos1 <- pos1 & predictions[[focal_terms[2]]] == pairs_data[[1]][i, 2] + pos2 <- pos2 & predictions[[focal_terms[2]]] == pairs_data[[2]][i, 2] + } + if (length(focal_terms) > 2) { + pos1 <- pos1 & predictions[[focal_terms[3]]] == pairs_data[[1]][i, 3] + pos2 <- pos2 & predictions[[focal_terms[3]]] == pairs_data[[2]][i, 3] + } + # once we have found the correct rows for the pairs, we can calculate + # the contrast. We need the predicted values first + predicted1 <- predictions$Predicted[pos1] + predicted2 <- predictions$Predicted[pos2] + + # we then create labels for the pairs. "result" is a data frame with + # the labels (of the pairwise contrasts) as columns. + result <- data.frame( + Parameter = paste( + paste0("(", paste(pairs_data[[1]][i, ], collapse = " "), ")"), + paste0("(", paste(pairs_data[[2]][i, ], collapse = " "), ")"), + sep = "-" + ), + stringsAsFactors = FALSE + ) + # we then add the contrast and the standard error. for linear models, the + # SE is sqrt(se1^2 + se2^2). + result$Difference <- predicted1 - predicted2 + # sum of squared standard errors + sum_se_squared <- predictions$SE[pos1]^2 + predictions$SE[pos2]^2 + # for non-Gaussian models, we subtract the covariance of the two predictions + # but only if the vcov_matrix is not NULL and has the correct dimensions + correct_row_dims <- nrow(vcov_matrix) > 0 && all(nrow(vcov_matrix) >= which(pos1)) + correct_col_dims <- ncol(vcov_matrix) > 0 && all(ncol(vcov_matrix) >= which(pos2)) + if (is.null(vcov_matrix) || !correct_row_dims || !correct_col_dims) { + vcov_sub <- 0 + } else { + vcov_sub <- vcov_matrix[which(pos1), which(pos2)]^2 + } + # Avoid negative values in sqrt() + if (vcov_sub >= sum_se_squared) { + result$SE <- sqrt(sum_se_squared) + } else { + result$SE <- sqrt(sum_se_squared - vcov_sub) + } + result + })) + # add CI and p-values + out$CI_low <- out$Difference - stats::qt(crit_factor, df = dof) * out$SE + out$CI_high <- out$Difference + stats::qt(crit_factor, df = dof) * out$SE + out$Statistic <- out$Difference / out$SE + out$p <- 2 * stats::pt(abs(out$Statistic), df = dof, lower.tail = FALSE) + + # filter by "by" + if (!is.null(by)) { + idx <- rep_len(TRUE, nrow(out)) + for (filter_by in by) { + # create index with "by" variables for each comparison pair + filter_data <- do.call(cbind, lapply(pairs_data, function(i) { + colnames(i) <- focal_terms + i[filter_by] + })) + # check which pairs have identical values - these are the rows we want to keep + idx <- idx & unname(apply(filter_data, 1, function(r) r[1] == r[2])) + } + # prepare filtered dataset + filter_column <- pairs_data[[1]] + colnames(filter_column) <- focal_terms + # bind the filtered data to the output + out <- cbind(filter_column[idx, by, drop = FALSE], out[idx, , drop = FALSE]) + } + + rownames(out) <- NULL + out +} + + +# interaction contrasts ---------------------------------------------------- +.compute_interactions <- function(predictions, dof, vcov_matrix, datagrid, contrast, by, crit_factor) { + # we need the focal terms and all unique values from the datagrid + focal_terms <- c(contrast, by) + at_list <- lapply(datagrid[focal_terms], unique) + + ## TODO: interaction contrasts currently only work for two focal terms + if (length(focal_terms) != 2) { + insight::format_error("Interaction contrasts currently only work for two focal terms.") + } + + # create pairwise combinations of first focal term + level_pairs <- at_list[[1]] + M <- matrix( + 1, + nrow = length(level_pairs), + ncol = length(level_pairs), + dimnames = list(level_pairs, level_pairs) + ) + M[!upper.tri(M)] <- NA + # table() works fine to create variables of this pairwise combinations + pairs_focal1 <- stats::na.omit(as.data.frame(as.table(M))) + pairs_focal1$Freq <- NULL + + # create pairwise combinations of second focal term + level_pairs <- at_list[[2]] + M <- matrix( + 1, + nrow = length(level_pairs), + ncol = length(level_pairs), + dimnames = list(level_pairs, level_pairs) + ) + M[!upper.tri(M)] <- NA + # table() works fine to create variables of this pairwise combinations + pairs_focal2 <- stats::na.omit(as.data.frame(as.table(M))) + pairs_focal2$Freq <- NULL + + # now we iterate over all pairs and try to find the corresponding predictions + out <- do.call(rbind, lapply(seq_len(nrow(pairs_focal1)), function(i) { + # differences between levels of first focal term + pos1 <- predictions[[focal_terms[1]]] == pairs_focal1[i, 1] + pos2 <- predictions[[focal_terms[1]]] == pairs_focal1[i, 2] + + do.call(rbind, lapply(seq_len(nrow(pairs_focal2)), function(j) { + # difference between levels of first focal term, *within* first + # level of second focal term + pos_1a <- pos1 & predictions[[focal_terms[2]]] == pairs_focal2[j, 1] + pos_1b <- pos2 & predictions[[focal_terms[2]]] == pairs_focal2[j, 1] + # difference between levels of first focal term, *within* second + # level of second focal term + pos_2a <- pos1 & predictions[[focal_terms[2]]] == pairs_focal2[j, 2] + pos_2b <- pos2 & predictions[[focal_terms[2]]] == pairs_focal2[j, 2] + # once we have found the correct rows for the pairs, we can calculate + # the contrast. We need the predicted values first + predicted1 <- predictions$Predicted[pos_1a] - predictions$Predicted[pos_1b] + predicted2 <- predictions$Predicted[pos_2a] - predictions$Predicted[pos_2b] + # we then create labels for the pairs. "result" is a data frame with + # the labels (of the pairwise contrasts) as columns. + result <- data.frame( + a = paste(pairs_focal1[i, 1], pairs_focal1[i, 2], sep = "-"), + b = paste(pairs_focal2[j, 1], pairs_focal2[j, 2], sep = " and "), + stringsAsFactors = FALSE + ) + colnames(result) <- focal_terms + # we then add the contrast and the standard error. for linear models, the + # SE is sqrt(se1^2 + se2^2) + result$Difference <- predicted1 - predicted2 + sum_se_squared <- sum( + predictions$SE[pos_1a]^2, predictions$SE[pos_1b]^2, + predictions$SE[pos_2a]^2, predictions$SE[pos_2b]^2 + ) + # for non-Gaussian models, we subtract the covariance of the two predictions + # but only if the vcov_matrix is not NULL and has the correct dimensions + correct_row_dims <- nrow(vcov_matrix) > 0 && all(nrow(vcov_matrix) >= which(pos_1a)) && all(nrow(vcov_matrix) >= which(pos_2a)) # nolint + correct_col_dims <- ncol(vcov_matrix) > 0 && all(ncol(vcov_matrix) >= which(pos_1b)) && all(ncol(vcov_matrix) >= which(pos_2b)) # nolint + if (is.null(vcov_matrix) || !correct_row_dims || !correct_col_dims) { + vcov_sub <- 0 + } else { + vcov_sub <- sum( + vcov_matrix[which(pos_1a), which(pos_1b)]^2, + vcov_matrix[which(pos_2a), which(pos_2b)]^2 + ) + } + # Avoid negative values in sqrt() + if (vcov_sub >= sum_se_squared) { + result$SE <- sqrt(sum_se_squared) + } else { + result$SE <- sqrt(sum_se_squared - vcov_sub) + } + result + })) + })) + # add CI and p-values + out$CI_low <- out$Difference - stats::qt(crit_factor, df = dof) * out$SE + out$CI_high <- out$Difference + stats::qt(crit_factor, df = dof) * out$SE + out$Statistic <- out$Difference / out$SE + out$p <- 2 * stats::pt(abs(out$Statistic), df = dof, lower.tail = FALSE) + out +} diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R index 9f46f2ac..a54e70fd 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -92,18 +92,25 @@ #' #' @return A data frame of estimated contrasts. #' @export -estimate_contrasts <- function(model, - contrast = NULL, - by = NULL, - predict = NULL, - ci = 0.95, - p_adjust = "none", - comparison = "pairwise", - marginalize = "average", - backend = getOption("modelbased_backend", "marginaleffects"), - transform = NULL, - verbose = TRUE, - ...) { +estimate_contrasts <- function(model, ...) { + UseMethod("estimate_contrasts") +} + + +#' @rdname estimate_contrasts +#' @export +estimate_contrasts.default <- function(model, + contrast = NULL, + by = NULL, + predict = NULL, + ci = 0.95, + p_adjust = "none", + comparison = "pairwise", + marginalize = "average", + backend = getOption("modelbased_backend", "marginaleffects"), + transform = NULL, + verbose = TRUE, + ...) { ## TODO: remove deprecation warning later if (!is.null(transform)) { insight::format_warning("Argument `transform` is deprecated. Please use `predict` instead.") diff --git a/R/estimate_predicted.R b/R/estimate_predicted.R index 053c8652..dd077afb 100644 --- a/R/estimate_predicted.R +++ b/R/estimate_predicted.R @@ -450,6 +450,7 @@ estimate_relation <- function(model, attr(out, "keep_iterations") <- keep_iterations attr(out, "response") <- model_response attr(out, "model") <- model + attr(out, "datagrid") <- data attr(out, "focal_terms") <- grid_specs$at_specs$varname attr(out, "preserve_range") <- grid_specs$preserve_range attr(out, "table_title") <- c("Model-based Predictions", "blue") diff --git a/R/get_marginalcontrasts.R b/R/get_marginalcontrasts.R index 5549fe89..f5c84ec5 100644 --- a/R/get_marginalcontrasts.R +++ b/R/get_marginalcontrasts.R @@ -283,40 +283,40 @@ get_marginalcontrasts <- function(model, return(params) } + # harmonize argument + p_adjust <- tolower(p_adjust) + all_methods <- c(tolower(stats::p.adjust.methods), "tukey", "sidak") + insight::validate_argument(p_adjust, all_methods) # needed for rank adjustment focal_terms <- datagrid[focal] rank_adjust <- prod(vapply(focal_terms, insight::n_unique, numeric(1))) - # only proceed if valid argument-value - if (tolower(p_adjust) %in% all_methods) { - if (tolower(p_adjust) %in% tolower(stats::p.adjust.methods)) { - # base R adjustments - params[["p"]] <- stats::p.adjust(params[["p"]], method = p_adjust) - } else if (tolower(p_adjust) == "tukey") { - if (!is.null(statistic)) { - # tukey adjustment - params[["p"]] <- suppressWarnings(stats::ptukey( - sqrt(2) * abs(statistic), - rank_adjust, - dof, - lower.tail = FALSE - )) - # for specific contrasts, ptukey might fail, and the tukey-adjustement - # could just be simple p-value calculation - if (all(is.na(params[["p"]]))) { - params[["p"]] <- 2 * stats::pt(abs(statistic), df = dof, lower.tail = FALSE) - } - } else if (verbose) { - insight::format_alert("No test-statistic found. P-values were not adjusted.") + if (p_adjust %in% tolower(stats::p.adjust.methods)) { + # base R adjustments + params[["p"]] <- stats::p.adjust(params[["p"]], method = p_adjust) + } else if (p_adjust == "tukey") { + if (!is.null(statistic)) { + # tukey adjustment + params[["p"]] <- suppressWarnings(stats::ptukey( + sqrt(2) * abs(statistic), + rank_adjust, + dof, + lower.tail = FALSE + )) + # for specific contrasts, ptukey might fail, and the tukey-adjustement + # could just be simple p-value calculation + if (all(is.na(params[["p"]]))) { + params[["p"]] <- 2 * stats::pt(abs(statistic), df = dof, lower.tail = FALSE) } - } else if (tolower(p_adjust) == "sidak") { - # sidak adjustment - params[["p"]] <- 1 - (1 - params[["p"]])^rank_adjust + } else if (verbose) { + insight::format_alert("No test-statistic found. P-values were not adjusted.") } - } else { - insight::format_error(paste0("`p_adjust` must be one of ", toString(all_methods))) + } else if (p_adjust == "sidak") { + # sidak adjustment + params[["p"]] <- 1 - (1 - params[["p"]])^rank_adjust } + params } diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd index 52de36fb..105db8b0 100644 --- a/man/estimate_contrasts.Rd +++ b/man/estimate_contrasts.Rd @@ -2,9 +2,12 @@ % Please edit documentation in R/estimate_contrasts.R \name{estimate_contrasts} \alias{estimate_contrasts} +\alias{estimate_contrasts.default} \title{Estimate Marginal Contrasts} \usage{ -estimate_contrasts( +estimate_contrasts(model, ...) + +\method{estimate_contrasts}{default}( model, contrast = NULL, by = NULL, @@ -22,6 +25,23 @@ estimate_contrasts( \arguments{ \item{model}{A statistical model.} +\item{...}{Other arguments passed, for instance, to \code{\link[insight:get_datagrid]{insight::get_datagrid()}}, +to functions from the \strong{emmeans} or \strong{marginaleffects} package, or to process +Bayesian models via \code{\link[bayestestR:describe_posterior]{bayestestR::describe_posterior()}}. Examples: +\itemize{ +\item \code{insight::get_datagrid()}: Argument such as \code{length} or \code{range} can be used +to control the (number of) representative values. +\item \strong{marginaleffects}: Internally used functions are \code{avg_predictions()} for +means and contrasts, and \code{avg_slope()} for slopes. Therefore, arguments +for instance like \code{vcov}, \code{transform}, \code{equivalence} or \code{slope} can be +passed to those functions. +\item \strong{emmeans}: Internally used functions are \code{emmeans()} and \code{emtrends()}. +Additional arguments can be passed to these functions. +\item Bayesian models: For Bayesian models, parameters are cleaned using +\code{describe_posterior()}, thus, arguments like, for example, \code{centrality}, +\code{rope_range}, or \code{test} are passed to that function. +}} + \item{contrast}{A character vector indicating the name of the variable(s) for which to compute the contrasts.} @@ -151,23 +171,6 @@ as default backend.} \item{transform}{Deprecated, please use \code{predict} instead.} \item{verbose}{Use \code{FALSE} to silence messages and warnings.} - -\item{...}{Other arguments passed, for instance, to \code{\link[insight:get_datagrid]{insight::get_datagrid()}}, -to functions from the \strong{emmeans} or \strong{marginaleffects} package, or to process -Bayesian models via \code{\link[bayestestR:describe_posterior]{bayestestR::describe_posterior()}}. Examples: -\itemize{ -\item \code{insight::get_datagrid()}: Argument such as \code{length} or \code{range} can be used -to control the (number of) representative values. -\item \strong{marginaleffects}: Internally used functions are \code{avg_predictions()} for -means and contrasts, and \code{avg_slope()} for slopes. Therefore, arguments -for instance like \code{vcov}, \code{transform}, \code{equivalence} or \code{slope} can be -passed to those functions. -\item \strong{emmeans}: Internally used functions are \code{emmeans()} and \code{emtrends()}. -Additional arguments can be passed to these functions. -\item Bayesian models: For Bayesian models, parameters are cleaned using -\code{describe_posterior()}, thus, arguments like, for example, \code{centrality}, -\code{rope_range}, or \code{test} are passed to that function. -}} } \value{ A data frame of estimated contrasts. diff --git a/tests/testthat/_snaps/estimate_contrasts_methods.md b/tests/testthat/_snaps/estimate_contrasts_methods.md new file mode 100644 index 00000000..3c85cec4 --- /dev/null +++ b/tests/testthat/_snaps/estimate_contrasts_methods.md @@ -0,0 +1,149 @@ +# estimate_contrasts - Random Effects Levels, pairwise + + Code + print(estimate_contrasts(estim, contrast = c("gender", "employed", "age")), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | Statistic | p + ------------------------------------------------------------------------------------------------- + Male, no, -40 | Female, no, -40 | 0.07 | 1.11 | [-2.11, 2.25] | 0.06 | 0.951 + Male, no, -40 | Male, yes, -40 | -0.88 | 1.18 | [-3.18, 1.43] | -0.75 | 0.455 + Female, no, -40 | Male, yes, -40 | -0.95 | 1.04 | [-2.99, 1.10] | -0.91 | 0.364 + Male, no, -40 | Female, yes, -40 | -0.79 | 1.07 | [-2.89, 1.30] | -0.74 | 0.457 + Female, no, -40 | Female, yes, -40 | -0.86 | 0.92 | [-2.67, 0.94] | -0.94 | 0.348 + Male, yes, -40 | Female, yes, -40 | 0.08 | 0.99 | [-1.86, 2.03] | 0.08 | 0.932 + Male, no, -40 | Male, no, 41-64 | 0.30 | 1.10 | [-1.84, 2.45] | 0.28 | 0.781 + Female, no, -40 | Male, no, 41-64 | 0.24 | 0.95 | [-1.63, 2.10] | 0.25 | 0.805 + Male, yes, -40 | Male, no, 41-64 | 1.18 | 1.02 | [-0.83, 3.19] | 1.15 | 0.248 + Female, yes, -40 | Male, no, 41-64 | 1.10 | 0.90 | [-0.66, 2.86] | 1.22 | 0.221 + Male, no, -40 | Female, no, 41-64 | 1.49 | 0.94 | [-0.35, 3.32] | 1.59 | 0.112 + Female, no, -40 | Female, no, 41-64 | 1.42 | 0.76 | [-0.08, 2.91] | 1.86 | 0.063 + Male, yes, -40 | Female, no, 41-64 | 2.37 | 0.85 | [ 0.70, 4.04] | 2.78 | 0.005 + Female, yes, -40 | Female, no, 41-64 | 2.28 | 0.69 | [ 0.92, 3.64] | 3.29 | 0.001 + Male, no, 41-64 | Female, no, 41-64 | 1.18 | 0.74 | [-0.26, 2.63] | 1.60 | 0.109 + Male, no, -40 | Male, yes, 41-64 | -0.23 | 1.04 | [-2.26, 1.81] | -0.22 | 0.827 + Female, no, -40 | Male, yes, 41-64 | -0.30 | 0.88 | [-2.03, 1.44] | -0.33 | 0.738 + Male, yes, -40 | Male, yes, 41-64 | 0.65 | 0.96 | [-1.23, 2.54] | 0.68 | 0.498 + Female, yes, -40 | Male, yes, 41-64 | 0.57 | 0.83 | [-1.05, 2.19] | 0.69 | 0.492 + Male, no, 41-64 | Male, yes, 41-64 | -0.53 | 0.86 | [-2.22, 1.16] | -0.62 | 0.538 + Female, no, 41-64 | Male, yes, 41-64 | -1.71 | 0.65 | [-2.98, -0.44] | -2.65 | 0.008 + Male, no, -40 | Female, yes, 41-64 | 1.15 | 0.94 | [-0.70, 2.99] | 1.22 | 0.224 + Female, no, -40 | Female, yes, 41-64 | 1.08 | 0.77 | [-0.43, 2.59] | 1.40 | 0.162 + Male, yes, -40 | Female, yes, 41-64 | 2.02 | 0.86 | [ 0.34, 3.71] | 2.36 | 0.018 + Female, yes, -40 | Female, yes, 41-64 | 1.94 | 0.70 | [ 0.56, 3.32] | 2.76 | 0.006 + Male, no, 41-64 | Female, yes, 41-64 | 0.84 | 0.74 | [-0.62, 2.30] | 1.13 | 0.258 + Female, no, 41-64 | Female, yes, 41-64 | -0.34 | 0.48 | [-1.28, 0.60] | -0.71 | 0.476 + Male, yes, 41-64 | Female, yes, 41-64 | 1.37 | 0.66 | [ 0.09, 2.66] | 2.09 | 0.036 + Male, no, -40 | Male, no, 65++ | 0.83 | 1.08 | [-1.29, 2.94] | 0.77 | 0.443 + Female, no, -40 | Male, no, 65++ | 0.76 | 0.93 | [-1.07, 2.58] | 0.82 | 0.415 + Male, yes, -40 | Male, no, 65++ | 1.71 | 1.01 | [-0.26, 3.68] | 1.70 | 0.090 + Female, yes, -40 | Male, no, 65++ | 1.62 | 0.88 | [-0.10, 3.34] | 1.85 | 0.064 + Male, no, 41-64 | Male, no, 65++ | 0.52 | 0.91 | [-1.26, 2.31] | 0.57 | 0.565 + Female, no, 41-64 | Male, no, 65++ | -0.66 | 0.71 | [-2.05, 0.73] | -0.93 | 0.354 + Male, yes, 41-64 | Male, no, 65++ | 1.06 | 0.84 | [-0.59, 2.70] | 1.26 | 0.209 + Female, yes, 41-64 | Male, no, 65++ | -0.32 | 0.72 | [-1.72, 1.09] | -0.44 | 0.658 + Male, no, -40 | Female, no, 65++ | 1.71 | 0.98 | [-0.21, 3.63] | 1.75 | 0.081 + Female, no, -40 | Female, no, 65++ | 1.64 | 0.81 | [ 0.04, 3.23] | 2.01 | 0.044 + Male, yes, -40 | Female, no, 65++ | 2.59 | 0.90 | [ 0.83, 4.35] | 2.88 | 0.004 + Female, yes, -40 | Female, no, 65++ | 2.50 | 0.75 | [ 1.03, 3.97] | 3.34 | < .001 + Male, no, 41-64 | Female, no, 65++ | 1.40 | 0.79 | [-0.14, 2.95] | 1.78 | 0.076 + Female, no, 41-64 | Female, no, 65++ | 0.22 | 0.55 | [-0.85, 1.29] | 0.40 | 0.686 + Male, yes, 41-64 | Female, no, 65++ | 1.94 | 0.71 | [ 0.55, 3.32] | 2.74 | 0.006 + Female, yes, 41-64 | Female, no, 65++ | 0.56 | 0.56 | [-0.53, 1.65] | 1.01 | 0.312 + Male, no, 65++ | Female, no, 65++ | 0.88 | 0.77 | [-0.62, 2.38] | 1.15 | 0.250 + Male, no, -40 | Male, yes, 65++ | -0.10 | 1.39 | [-2.83, 2.63] | -0.07 | 0.942 + Female, no, -40 | Male, yes, 65++ | -0.17 | 1.28 | [-2.69, 2.35] | -0.13 | 0.895 + Male, yes, -40 | Male, yes, 65++ | 0.78 | 1.34 | [-1.85, 3.40] | 0.58 | 0.561 + Female, yes, -40 | Male, yes, 65++ | 0.69 | 1.25 | [-1.75, 3.13] | 0.56 | 0.577 + Male, no, 41-64 | Male, yes, 65++ | -0.41 | 1.27 | [-2.89, 2.08] | -0.32 | 0.750 + Female, no, 41-64 | Male, yes, 65++ | -1.59 | 1.13 | [-3.81, 0.64] | -1.40 | 0.162 + Male, yes, 41-64 | Male, yes, 65++ | 0.13 | 1.22 | [-2.26, 2.52] | 0.10 | 0.917 + Female, yes, 41-64 | Male, yes, 65++ | -1.25 | 1.14 | [-3.48, 0.99] | -1.09 | 0.274 + Male, no, 65++ | Male, yes, 65++ | -0.93 | 1.25 | [-3.39, 1.53] | -0.74 | 0.459 + Female, no, 65++ | Male, yes, 65++ | -1.81 | 1.17 | [-4.10, 0.48] | -1.55 | 0.122 + Male, no, -40 | Female, yes, 65++ | 0.39 | 1.36 | [-2.27, 3.05] | 0.29 | 0.775 + Female, no, -40 | Female, yes, 65++ | 0.32 | 1.25 | [-2.12, 2.76] | 0.26 | 0.798 + Male, yes, -40 | Female, yes, 65++ | 1.27 | 1.30 | [-1.28, 3.82] | 0.97 | 0.330 + Female, yes, -40 | Female, yes, 65++ | 1.18 | 1.20 | [-1.18, 3.54] | 0.98 | 0.326 + Male, no, 41-64 | Female, yes, 65++ | 0.08 | 1.23 | [-2.33, 2.49] | 0.07 | 0.946 + Female, no, 41-64 | Female, yes, 65++ | -1.10 | 1.09 | [-3.23, 1.04] | -1.01 | 0.313 + Male, yes, 41-64 | Female, yes, 65++ | 0.62 | 1.18 | [-1.69, 2.92] | 0.52 | 0.601 + Female, yes, 41-64 | Female, yes, 65++ | -0.76 | 1.09 | [-2.90, 1.39] | -0.69 | 0.489 + Male, no, 65++ | Female, yes, 65++ | -0.44 | 1.21 | [-2.82, 1.94] | -0.36 | 0.717 + Female, no, 65++ | Female, yes, 65++ | -1.32 | 1.13 | [-3.53, 0.89] | -1.17 | 0.241 + Male, yes, 65++ | Female, yes, 65++ | 0.49 | 1.50 | [-2.45, 3.43] | 0.33 | 0.745 + + Variable predicted: qol + Predictors contrasted: gender, employed, age + +--- + + Code + print(estimate_contrasts(estim, contrast = c("gender", "employed"), by = "age"), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + Level1 | Level2 | age | Difference | SE | 95% CI | Statistic | p + ----------------------------------------------------------------------------------------- + Male, no | Female, no | -40 | 0.07 | 1.11 | [-2.11, 2.25] | 0.06 | 0.951 + Male, no | Male, yes | -40 | -0.88 | 1.18 | [-3.18, 1.43] | -0.75 | 0.455 + Female, no | Male, yes | -40 | -0.95 | 1.04 | [-2.99, 1.10] | -0.91 | 0.364 + Male, no | Female, yes | -40 | -0.79 | 1.07 | [-2.89, 1.30] | -0.74 | 0.457 + Female, no | Female, yes | -40 | -0.86 | 0.92 | [-2.67, 0.94] | -0.94 | 0.348 + Male, yes | Female, yes | -40 | 0.08 | 0.99 | [-1.86, 2.03] | 0.08 | 0.932 + Male, no | Female, no | 41-64 | 1.18 | 0.74 | [-0.26, 2.63] | 1.60 | 0.109 + Male, no | Male, yes | 41-64 | -0.53 | 0.86 | [-2.22, 1.16] | -0.62 | 0.538 + Female, no | Male, yes | 41-64 | -1.71 | 0.65 | [-2.98, -0.44] | -2.65 | 0.008 + Male, no | Female, yes | 41-64 | 0.84 | 0.74 | [-0.62, 2.30] | 1.13 | 0.258 + Female, no | Female, yes | 41-64 | -0.34 | 0.48 | [-1.28, 0.60] | -0.71 | 0.476 + Male, yes | Female, yes | 41-64 | 1.37 | 0.66 | [ 0.09, 2.66] | 2.09 | 0.036 + Male, no | Female, no | 65+ | 0.88 | 0.77 | [-0.62, 2.38] | 1.15 | 0.250 + Male, no | Male, yes | 65+ | -0.93 | 1.25 | [-3.39, 1.53] | -0.74 | 0.459 + Female, no | Male, yes | 65+ | -1.81 | 1.17 | [-4.10, 0.48] | -1.55 | 0.122 + Male, no | Female, yes | 65+ | -0.44 | 1.21 | [-2.82, 1.94] | -0.36 | 0.717 + Female, no | Female, yes | 65+ | -1.32 | 1.13 | [-3.53, 0.89] | -1.17 | 0.241 + Male, yes | Female, yes | 65+ | 0.49 | 1.50 | [-2.45, 3.43] | 0.33 | 0.745 + + Variable predicted: qol + Predictors contrasted: gender, employed + +--- + + Code + print(estimate_contrasts(estim, contrast = "employed", by = c("age", "gender")), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + Level1 | Level2 | age | gender | Difference | SE | 95% CI | Statistic | p + ---------------------------------------------------------------------------------------- + no | yes | -40 | Female | -0.86 | 0.92 | [-2.67, 0.94] | -0.94 | 0.348 + no | yes | -40 | Male | -0.88 | 1.18 | [-3.18, 1.43] | -0.75 | 0.455 + no | yes | 41-64 | Female | -0.34 | 0.48 | [-1.28, 0.60] | -0.71 | 0.476 + no | yes | 41-64 | Male | -0.53 | 0.86 | [-2.22, 1.16] | -0.62 | 0.538 + no | yes | 65+ | Female | -1.32 | 1.13 | [-3.53, 0.89] | -1.17 | 0.241 + no | yes | 65+ | Male | -0.93 | 1.25 | [-3.39, 1.53] | -0.74 | 0.459 + + Variable predicted: qol + Predictors contrasted: employed + +# estimate_contrasts - Random Effects Levels, interaction + + Code + print(estimate_contrasts(estim, contrast = c("age", "employed"), comparison = "interaction"), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + age | employed | Difference | SE | 95% CI | Statistic | p + ------------------------------------------------------------------------------ + -40-41-64 | no and yes | -0.49 | 0.94 | [-2.34, 1.36] | -0.52 | 0.600 + -40-65+ | no and yes | 0.44 | 1.41 | [-2.31, 3.20] | 0.31 | 0.753 + 41-64-65+ | no and yes | 0.94 | 1.21 | [-1.43, 3.30] | 0.78 | 0.438 + + Variable predicted: qol + Predictors contrasted: age, employed + diff --git a/tests/testthat/test-attributes_estimatefun.R b/tests/testthat/test-attributes_estimatefun.R index 195cf204..f3e7c1a3 100644 --- a/tests/testthat/test-attributes_estimatefun.R +++ b/tests/testthat/test-attributes_estimatefun.R @@ -87,8 +87,9 @@ test_that("attributes_means", { attributes(estim), c( "names", "row.names", "class", "ci", "keep_iterations", "response", - "model", "focal_terms", "preserve_range", "table_title", "table_footer", - "adjusted_for", "at_specs", "at", "by", "reference", "data" + "model", "datagrid", "focal_terms", "preserve_range", "table_title", + "table_footer", "adjusted_for", "at_specs", "at", "by", "reference", + "data" ) ) }) diff --git a/tests/testthat/test-estimate_contrasts.R b/tests/testthat/test-estimate_contrasts.R index 2a783796..785e4025 100644 --- a/tests/testthat/test-estimate_contrasts.R +++ b/tests/testthat/test-estimate_contrasts.R @@ -120,7 +120,8 @@ test_that("estimate_contrasts - Frequentist, Three factors 1", { 7.03333, 14.01667, 0.05, 7.03333, 6.98333, 6.98333, 11.275, 18.25833, 4.29167, 11.275, 6.98333 ), - tolerance = 1e-4) + tolerance = 1e-4 + ) expect_snapshot(print(estimate_contrasts(model, contrast = c("vs", "am"), by = "gear", backend = "marginaleffects"), zap_small = TRUE, table_width = Inf)) # nolint }) @@ -142,12 +143,13 @@ test_that("estimate_contrasts - Frequentist, Three factors 2", { by = c("c161sex", "c172code", "e16sex"), hypothesis = ~ pairwise | e16sex ) - expect_equal(out$estimate[out$e16sex == "male"], estim$Difference, tolerance = 1e-4) + expect_equal(out$estimate[out$e16sex == "male"], estim$Difference, tolerance = 1e-4) expect_identical(dim(estim), c(15L, 10L)) expect_equal( estim$Difference, - c(2.70334, 2.23511, 1.55845, 2.48322, 2.48543, -0.46823, -1.14489, + c( + 2.70334, 2.23511, 1.55845, 2.48322, 2.48543, -0.46823, -1.14489, -0.22012, -0.21791, -0.67666, 0.24811, 0.25032, 0.92477, 0.92698, 0.00221 ), @@ -156,7 +158,8 @@ test_that("estimate_contrasts - Frequentist, Three factors 2", { expect_true(all(estim$e16sex == "male")) expect_identical( as.character(estim$Level1), - c("Male, mid", "Male, high", "Female, low", "Female, mid", "Female, high", + c( + "Male, mid", "Male, high", "Female, low", "Female, mid", "Female, high", "Male, high", "Female, low", "Female, mid", "Female, high", "Female, low", "Female, mid", "Female, high", "Female, mid", "Female, high", "Female, high" @@ -164,7 +167,8 @@ test_that("estimate_contrasts - Frequentist, Three factors 2", { ) expect_identical( as.character(estim$Level2), - c("Male, low", "Male, low", "Male, low", "Male, low", "Male, low", + c( + "Male, low", "Male, low", "Male, low", "Male, low", "Male, low", "Male, mid", "Male, mid", "Male, mid", "Male, mid", "Male, high", "Male, high", "Male, high", "Female, low", "Female, low", "Female, mid" ) diff --git a/tests/testthat/test-estimate_contrasts_methods.R b/tests/testthat/test-estimate_contrasts_methods.R new file mode 100644 index 00000000..99aaa910 --- /dev/null +++ b/tests/testthat/test-estimate_contrasts_methods.R @@ -0,0 +1,70 @@ +skip_on_cran() +skip_on_os("mac") +skip_if_not_installed("glmmTMB") +skip_if_not(interactive()) + + +test_that("estimate_contrasts - Random Effects Levels, pairwise", { + # sample data set + data(efc, package = "modelbased") + + # numeric to factors, set labels as levels + d <- datawizard::to_factor(efc, select = c("c161sex", "c172code", "c175empl")) + # recode age into three groups + d <- datawizard::recode_values( + d, + select = "c160age", + recode = list(`1` = "min:40", `2` = 41:64, `3` = "65:max") + ) + # rename variables + d <- datawizard::data_rename( + d, + select = c("c161sex", "c160age", "quol_5", "c175empl"), + replacement = c("gender", "age", "qol", "employed") + ) + # age into factor, set levels, and change labels for education + d <- datawizard::data_modify(d, age = factor(age, labels = c("-40", "41-64", "65+"))) + dat <<- d + + # Quality of Life score ranges from 0 to 25 + m_null <- glmmTMB::glmmTMB(qol ~ 1 + (1 | gender:employed:age), data = dat) + estim <- estimate_relation(m_null, by = c("gender", "employed", "age")) + + # test errors + expect_error(estimate_contrasts(estim), regex = "must be specified") + expect_error(estimate_contrasts(estim, "employed", comparison = ~reference), regex = "Invalid option for argument") + + # test output + expect_snapshot(print(estimate_contrasts(estim, contrast = c("gender", "employed", "age")), zap_small = TRUE, table_width = Inf)) + expect_snapshot(print(estimate_contrasts(estim, contrast = c("gender", "employed"), by = "age"), zap_small = TRUE, table_width = Inf)) + expect_snapshot(print(estimate_contrasts(estim, contrast = "employed", by = c("age", "gender")), zap_small = TRUE, table_width = Inf)) +}) + + +test_that("estimate_contrasts - Random Effects Levels, interaction", { + # sample data set + data(efc, package = "modelbased") + + # numeric to factors, set labels as levels + d <- datawizard::to_factor(efc, select = c("c161sex", "c172code", "c175empl")) + # recode age into three groups + d <- datawizard::recode_values( + d, + select = "c160age", + recode = list(`1` = "min:40", `2` = 41:64, `3` = "65:max") + ) + # rename variables + d <- datawizard::data_rename( + d, + select = c("c161sex", "c160age", "quol_5", "c175empl"), + replacement = c("gender", "age", "qol", "employed") + ) + # age into factor, set levels, and change labels for education + d <- datawizard::data_modify(d, age = factor(age, labels = c("-40", "41-64", "65+"))) + dat <<- d + + # Quality of Life score ranges from 0 to 25 + m_null <- glmmTMB::glmmTMB(qol ~ 1 + (1 | gender:employed:age), data = dat) + estim <- estimate_relation(m_null, by = c("age", "employed")) + expect_snapshot(print(estimate_contrasts(estim, contrast = c("age", "employed"), comparison = "interaction"), zap_small = TRUE, table_width = Inf)) +})