|
11 | 11 | focal <- attributes(params)$contrast |
12 | 12 | # Use .safe to handle cases where no statistic is extracted |
13 | 13 | statistic <- .safe(insight::get_statistic(model)$Statistic) |
14 | | - dof <- insight::get_df(model, type = "wald", verbose = FALSE) |
| 14 | + # extract degrees of freedom |
| 15 | + dof <- .safe(params$df[1]) |
| 16 | + if (is.null(dof)) { |
| 17 | + dof <- insight::get_df(model, type = "wald", verbose = FALSE) |
| 18 | + } |
15 | 19 |
|
16 | 20 | # harmonize argument |
17 | 21 | p_adjust <- tolower(p_adjust) |
|
20 | 24 | # we have to print a different error message. |
21 | 25 | emmeans_options <- c("scheffe", "mvt", "dunnettx") |
22 | 26 |
|
23 | | - all_methods <- c(tolower(stats::p.adjust.methods), emmeans_options, "tukey", "sidak", "esarey") |
| 27 | + all_methods <- c(tolower(stats::p.adjust.methods), emmeans_options, "tukey", "sidak", "esarey", "sup-t") |
24 | 28 | insight::validate_argument(p_adjust, all_methods) |
25 | 29 |
|
26 | 30 | # emmeans methods? Then tell user |
|
35 | 39 | return(.p_adjust_esarey(params)) |
36 | 40 | } |
37 | 41 |
|
| 42 | + # sup-t is a longer subroutine, so we handle it separately |
| 43 | + if (p_adjust == "sup-t") { |
| 44 | + return(.p_adjust_supt(params, model)) |
| 45 | + } |
| 46 | + |
38 | 47 | # needed for rank adjustment |
39 | 48 | focal_terms <- datagrid[focal] |
40 | 49 | rank_adjust <- prod(vapply(focal_terms, insight::n_unique, numeric(1))) |
|
68 | 77 | } |
69 | 78 |
|
70 | 79 |
|
| 80 | +.p_adjust_supt <- function(params, model) { |
| 81 | + insight::check_if_installed("mvtnorm") |
| 82 | + # get correlation matrix, based on the covariance matrix |
| 83 | + vc <- .safe(stats::cov2cor(attributes(params)$vcov)) |
| 84 | + if (is.null(vc)) { |
| 85 | + insight::format_alert("Could not calculate covariance matrix for `sup-t` adjustment.") |
| 86 | + return(params) |
| 87 | + } |
| 88 | + # get confidence interval level, or set default |
| 89 | + ci_level <- attributes(params)$ci |
| 90 | + if (is.null(ci_level)) { |
| 91 | + ci_level <- 0.95 |
| 92 | + } |
| 93 | + # several sanity checks - we can either have a marginaleffects object, when |
| 94 | + # `estimate_slopes()` was called, or a modelbased object, when processing / |
| 95 | + # formatting was already done. So we check for both, and extract the required |
| 96 | + # columns. |
| 97 | + df_column <- intersect(c("df", "df_error"), colnames(params))[1] |
| 98 | + if (is.na(df_column)) { |
| 99 | + df_column <- ".sup_df" |
| 100 | + params[[df_column]] <- Inf |
| 101 | + } |
| 102 | + coef_column <- intersect(c(.valid_coefficient_names(), "estimate"), colnames(params))[1] |
| 103 | + if (is.na(coef_column)) { |
| 104 | + insight::format_alert("Could not find coefficient column to apply `sup-t` adjustment.") |
| 105 | + return(params) |
| 106 | + } |
| 107 | + se_column <- intersect(c("SE", "std.error"), colnames(params))[1] |
| 108 | + if (is.na(se_column)) { |
| 109 | + insight::format_alert("Could not extract standard errors to apply `sup-t` adjustment.") |
| 110 | + return(params) |
| 111 | + } |
| 112 | + p_column <- intersect(c("p", "p.value"), colnames(params))[1] |
| 113 | + if (is.na(p_column)) { |
| 114 | + insight::format_alert("Could not extract p-values to apply `sup-t` adjustment.") |
| 115 | + return(params) |
| 116 | + } |
| 117 | + ci_low_column <- intersect(c("CI_low", "conf.low"), colnames(params))[1] |
| 118 | + ci_high_column <- intersect(c("CI_high", "conf.high"), colnames(params))[1] |
| 119 | + if (is.na(ci_low_column) || is.na(ci_high_column)) { |
| 120 | + insight::format_alert("Could not extract confidence intervals to apply `sup-t` adjustment.") |
| 121 | + return(params) |
| 122 | + } |
| 123 | + # calculate updated confidence interval level, based on simultaenous |
| 124 | + # confidence intervals (https://onlinelibrary.wiley.com/doi/10.1002/jae.2656) |
| 125 | + crit <- mvtnorm::qmvt(ci_level, df = params[[df_column]][1], tail = "both.tails", corr = vc)$quantile |
| 126 | + # update confidence intervals |
| 127 | + params[[ci_low_column]] <- params[[coef_column]] - crit * params[[se_column]] |
| 128 | + params[[ci_high_column]] <- params[[coef_column]] + crit * params[[se_column]] |
| 129 | + # update p-values |
| 130 | + for (i in 1:nrow(params)) { |
| 131 | + params[[p_column]][i] <- 1 - mvtnorm::pmvt( |
| 132 | + lower = rep(-abs(stats::qt(params[[p_column]][i] / 2, df = params[[df_column]][i])), nrow(vc)), |
| 133 | + upper = rep(abs(stats::qt(params[[p_column]][i] / 2, df = params[[df_column]][i])), nrow(vc)), |
| 134 | + corr = vc, |
| 135 | + df = params[[df_column]][i] |
| 136 | + ) |
| 137 | + } |
| 138 | + # clean up - remove temporary column |
| 139 | + params[[".sup_df"]] <- NULL |
| 140 | + |
| 141 | + params |
| 142 | +} |
| 143 | + |
| 144 | + |
71 | 145 | .p_adjust_esarey <- function(x) { |
72 | 146 | # only for slopes |
73 | 147 | if (!inherits(x, c("estimate_slopes", "marginaleffects_slopes"))) { |
|
0 commit comments