Merge branch 'main' into rc_0_27_0

Author: strengejacke

DESCRIPTION

@@ -175,6 +175,7 @@ Suggests:
     mmrm,
     multcomp,
     MuMIn,
+    mvtnorm,
     NbClust,
     nFactors,
     nestedLogit,

NEWS.md

@@ -26,6 +26,9 @@
 * Improved efficiency in `model_parameters()` for models from packages *brms*
   and *rstanarm*.
 
+* `p_adjust` for `model_parameters()` gets a new options, `"sup-t"`, to calculate
+  simultaneous confidence intervals.
+
 # parameters 0.26.0
 
 ## Changes

R/1_model_parameters.R

@@ -364,6 +364,10 @@
 #'     Indices of Effect Existence and Significance in the Bayesian Framework.
 #'     Frontiers in Psychology, 10, 2767. \doi{10.3389/fpsyg.2019.02767}
 #'
+#'   - Montiel Olea, J. L., and Plagborg-Møller, M. (2019). Simultaneous
+#'     confidence bands: Theory, implementation, and an application to SVARs.
+#'     Journal of Applied Econometrics, 34(1), 1–17. \doi{10.1002/jae.2656}
+#'
 #'   - Neter, J., Wasserman, W., and Kutner, M. H. (1989). Applied linear
 #'     regression models.
 #'
@@ -463,11 +467,12 @@ parameters <- model_parameters
 #'   coefficient. The transformed confidence interval more clearly captures this
 #'   uncertainty. For `compare_parameters()`, `exponentiate = "nongaussian"`
 #'   will only exponentiate coefficients from non-Gaussian families.
-#' @param p_adjust Character vector, if not `NULL`, indicates the method to
-#'   adjust p-values. See [`stats::p.adjust()`] for details. Further
-#'   possible adjustment methods are `"tukey"`, `"scheffe"`,
-#'   `"sidak"` and `"none"` to explicitly disable adjustment for
-#'   `emmGrid` objects (from **emmeans**).
+#' @param p_adjust String value, if not `NULL`, indicates the method to adjust
+#'   p-values. See [`stats::p.adjust()`] for details. Further possible
+#'   adjustment methods are `"tukey"`, `"scheffe"`, `"sidak"`, `"sup-t"`, and
+#'   `"none"` to explicitly disable adjustment for `emmGrid` objects (from
+#'   **emmeans**). `"sup-t"` computes simultaneous confidence bands, also called
+#'   sup-t confidence band (Montiel Olea & Plagborg-Møller, 2019).
 #' @param ci_method Method for computing degrees of freedom for
 #'   confidence intervals (CI) and the related p-values. Allowed are following
 #'   options (which vary depending on the model class): `"residual"`,

R/format_p_adjust.R

@@ -25,36 +25,26 @@ format_p_adjust <- function(method) {
     tukey = "Tukey",
     scheffe = "Scheffe",
     sidak = "Sidak",
+    `sup-t` = "Simultaneous confidence bands",
     method
   )
 }
 
 
+# p-value adjustment -----
+
 .p_adjust <- function(params, p_adjust, model = NULL, verbose = TRUE) {
   # check if we have any adjustment at all, and a p-column
   if (!is.null(p_adjust) && "p" %in% colnames(params) && p_adjust != "none") {
     ## TODO add "mvt" method from emmeans
 
     # prepare arguments
-    all_methods <- c(stats::p.adjust.methods, "tukey", "scheffe", "sidak")
+    all_methods <- c(stats::p.adjust.methods, "tukey", "scheffe", "sidak", "sup-t")
 
     # for interaction terms, e.g. for "by" argument in emmeans
     # pairwise comparison, we have to adjust the rank resp. the
     # number of estimates in a comparison family
-    rank_adjust <- tryCatch(
-      {
-        correction <- 1
-        by_vars <- model@misc$by.vars
-        if (!is.null(by_vars) && by_vars %in% colnames(params)) {
-          correction <- insight::n_unique(params[[by_vars]])
-        }
-        correction
-      },
-      error = function(e) {
-        1
-      }
-    )
-
+    rank_adjust <- .p_adjust_rank(model, params)
 
     # only proceed if valid argument-value
     if (tolower(p_adjust) %in% tolower(all_methods)) {
@@ -69,44 +59,18 @@ format_p_adjust <- function(method) {
         params$p <- stats::p.adjust(params$p, method = p_adjust)
       } else if (tolower(p_adjust) == "tukey") {
         # tukey adjustment
-        if ("df" %in% colnames(params) && length(stat_column) > 0) {
-          params$p <- suppressWarnings(stats::ptukey(
-            sqrt(2) * abs(params[[stat_column]]),
-            nrow(params) / rank_adjust,
-            params$df,
-            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(params[[stat_column]]), df = params$df, lower.tail = FALSE)
-            verbose <- FALSE
-          }
-        }
+        result <- .p_adjust_tukey(params, stat_column, rank_adjust, verbose)
+        params <- result$params
+        verbose <- result$verbose
       } else if (tolower(p_adjust) == "scheffe" && !is.null(model)) {
         # scheffe adjustment
-        if ("df" %in% colnames(params) && length(stat_column) > 0) {
-          # 1st try
-          scheffe_ranks <- try(qr(model@linfct)$rank, silent = TRUE)
-
-          # 2nd try
-          if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
-            scheffe_ranks <- try(model$qr$rank, silent = TRUE)
-          }
-
-          if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
-            scheffe_ranks <- nrow(params)
-          }
-          scheffe_ranks <- scheffe_ranks / rank_adjust
-          params$p <- stats::pf(params[[stat_column]]^2 / scheffe_ranks,
-            df1 = scheffe_ranks,
-            df2 = params$df,
-            lower.tail = FALSE
-          )
-        }
+        params <- .p_adjust_scheffe(model, params, stat_column, rank_adjust)
       } else if (tolower(p_adjust) == "sidak") {
         # sidak adjustment
         params$p <- 1 - (1 - params$p)^(nrow(params) / rank_adjust)
+      }  else if (tolower(p_adjust) == "sup-t") {
+        # sup-t adjustment
+        params <- .p_adjust_supt(model, params)
       }
 
       if (isTRUE(all(old_p_vals == params$p)) && !identical(p_adjust, "none") && verbose) {
@@ -118,3 +82,133 @@ format_p_adjust <- function(method) {
   }
   params
 }
+
+
+# calculate rank adjustment -----
+
+.p_adjust_rank <- function(model, params) {
+  tryCatch(
+    {
+      correction <- 1
+      by_vars <- model@misc$by.vars
+      if (!is.null(by_vars) && by_vars %in% colnames(params)) {
+        correction <- insight::n_unique(params[[by_vars]])
+      }
+      correction
+    },
+    error = function(e) {
+      1
+    }
+  )
+}
+
+
+# tukey adjustment -----
+
+.p_adjust_tukey <- function(params, stat_column, rank_adjust = 1, verbose = TRUE) {
+  df_column <- colnames(params)[stats::na.omit(match(c("df", "df_error"), colnames(params)))][1]
+  if (length(df_column) && length(stat_column)) {
+    params$p <- suppressWarnings(stats::ptukey(
+      sqrt(2) * abs(params[[stat_column]]),
+      nmeans = nrow(params) / rank_adjust,
+      df = params[[df_column]],
+      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(params[[stat_column]]),
+        df = params[[df_column]],
+        lower.tail = FALSE
+      )
+      verbose <- FALSE
+    }
+  }
+  list(params = params, verbose = verbose)
+}
+
+
+# scheffe adjustment -----
+
+.p_adjust_scheffe <- function(model, params, stat_column, rank_adjust = 1) {
+  df_column <- colnames(params)[stats::na.omit(match(c("df", "df_error"), colnames(params)))][1]
+  if (length(df_column) && length(stat_column)) {
+    # 1st try
+    scheffe_ranks <- try(qr(model@linfct)$rank, silent = TRUE)
+
+    # 2nd try
+    if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
+      scheffe_ranks <- try(model$qr$rank, silent = TRUE)
+    }
+
+    if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
+      scheffe_ranks <- nrow(params)
+    }
+    scheffe_ranks <- scheffe_ranks / rank_adjust
+    params$p <- stats::pf(params[[stat_column]]^2 / scheffe_ranks,
+      df1 = scheffe_ranks,
+      df2 = params[[df_column]],
+      lower.tail = FALSE
+    )
+  }
+  params
+}
+
+
+# sup-t adjustment -----
+
+.p_adjust_supt <- function(model, params) {
+  if ("Component" %in% colnames(params) && insight::n_unique(params$Component) > 1) {
+    # if we have multiple components, we adjust each component separately
+    for (component in unique(params$Component)) {
+      if (!is.na(component)) {
+        params[which(params$Component == component), ] <- .supt_adjust(
+          params[which(params$Component == component), ],
+          model,
+          component = component
+        )
+      }
+    }
+    params
+  } else {
+    .supt_adjust(params, model)
+  }
+}
+
+
+.supt_adjust <- function(params, model, component = NULL) {
+  insight::check_if_installed("mvtnorm")
+  # get correlation matrix, based on the covariance matrix
+  vc <- .safe(stats::cov2cor(insight::get_varcov(model, component = component)))
+  if (is.null(vc)) {
+    insight::format_warning("Could not calculate covariance matrix for `sup-t` adjustment.")
+    return(params)
+  }
+  # get confidence interval level, or set default
+  ci_level <- .safe(params$CI[1])
+  if (is.null(ci_level)) {
+    ci_level <- 0.95
+  }
+  # find degrees of freedom column, if available
+  df_column <- colnames(params)[stats::na.omit(match(c("df", "df_error"), colnames(params)))][1]
+  if (length(df_column) == 0) {
+    return(params)
+  }
+  # calculate updated confidence interval level, based on simultaenous
+  # confidence intervals (https://onlinelibrary.wiley.com/doi/10.1002/jae.2656)
+  crit <- mvtnorm::qmvt(ci_level, df = params[[df_column]][1], tail = "both.tails", corr = vc)$quantile
+  # update confidence intervals
+  params$CI_low <- params$Coefficient - crit * params$SE
+  params$CI_high <- params$Coefficient + crit * params$SE
+  # udpate p-values
+  for (i in 1:nrow(params)) {
+    params$p[i] <- 1 - mvtnorm::pmvt(
+      lower = rep(-abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])), nrow(vc)),
+      upper = rep(abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])), nrow(vc)),
+      corr = vc,
+      df = params[[df_column]][i]
+    )
+  }
+  params
+}

inst/WORDLIST

@@ -117,6 +117,8 @@ McNemar
 Merkle
 Metaclustering
 Monti
+Møller
+Montiel
 Mundlak
 NHST
 NL
@@ -127,6 +129,7 @@ Nievergelt
 Nondegenerate
 Nonresponse
 ORCID
+Olea
 Olkin
 PCoA
 PHQ
@@ -135,6 +138,7 @@ PMCMRplus
 PeerJ
 Pernet
 Pettersson
+Plagborg
 PloS
 Psychometrika
 REWB
@@ -174,6 +178,7 @@ Stata
 Stigum
 Struyf
 Synthese
+SVARs
 TOST
 Tabachnick
 Thiyagarajan