fix

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.13.1
+Version: 0.13.1.1
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",

NEWS.md

@@ -1,3 +1,9 @@
+# modelbased (devel)
+
+## Bug fixes
+
+* Fixed issue in `estimate_slope()` when `p_adjust = "esarey"`.
+
 # modelbased 0.13.1
 
 ## Changes

R/p_adjust.R

@@ -24,13 +24,22 @@
   # we have to print a different error message.
   emmeans_options <- c("scheffe", "mvt", "dunnettx")
 
-  all_methods <- c(tolower(stats::p.adjust.methods), emmeans_options, "tukey", "sidak", "esarey", "sup-t")
+  all_methods <- c(
+    tolower(stats::p.adjust.methods),
+    emmeans_options,
+    "tukey",
+    "sidak",
+    "esarey",
+    "sup-t"
+  )
   insight::validate_argument(p_adjust, all_methods)
 
   # emmeans methods? Then tell user
   if (p_adjust %in% emmeans_options) {
     insight::format_error(paste0(
-      "`p_adjust = \"", p_adjust, "\"` is only available when `backend = \"emmeans\"."
+      "`p_adjust = \"",
+      p_adjust,
+      "\"` is only available when `backend = \"emmeans\"."
     ))
   }
 
@@ -96,12 +105,16 @@
   # columns.
   coef_column <- intersect(c(.valid_coefficient_names(), "estimate"), colnames(params))[1]
   if (is.na(coef_column)) {
-    insight::format_alert("Could not find coefficient column to apply `sup-t` adjustment.")
+    insight::format_alert(
+      "Could not find coefficient column to apply `sup-t` adjustment."
+    )
     return(params)
   }
   se_column <- intersect(c("SE", "std.error"), colnames(params))[1]
   if (is.na(se_column)) {
-    insight::format_alert("Could not extract standard errors to apply `sup-t` adjustment.")
+    insight::format_alert(
+      "Could not extract standard errors to apply `sup-t` adjustment."
+    )
     return(params)
   }
   p_column <- intersect(c("p", "p.value"), colnames(params))[1]
@@ -112,7 +125,9 @@
   ci_low_column <- intersect(c("CI_low", "conf.low"), colnames(params))[1]
   ci_high_column <- intersect(c("CI_high", "conf.high"), colnames(params))[1]
   if (is.na(ci_low_column) || is.na(ci_high_column)) {
-    insight::format_alert("Could not extract confidence intervals to apply `sup-t` adjustment.")
+    insight::format_alert(
+      "Could not extract confidence intervals to apply `sup-t` adjustment."
+    )
     return(params)
   }
   df_column <- intersect(c("df", "df_error"), colnames(params))[1]
@@ -122,18 +137,30 @@
   }
   # 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
+  crit <- mvtnorm::qmvt(
+    ci_level,
+    df = params[[df_column]][1],
+    tail = "both.tails",
+    corr = vc
+  )$quantile
   # update confidence intervals
   params[[ci_low_column]] <- params[[coef_column]] - crit * params[[se_column]]
   params[[ci_high_column]] <- params[[coef_column]] + crit * params[[se_column]]
   # update p-values
   for (i in 1:nrow(params)) {
-    params[[p_column]][i] <- 1 - mvtnorm::pmvt(
-      lower = rep(-abs(stats::qt(params[[p_column]][i] / 2, df = params[[df_column]][i])), nrow(vc)),
-      upper = rep(abs(stats::qt(params[[p_column]][i] / 2, df = params[[df_column]][i])), nrow(vc)),
-      corr = vc,
-      df = params[[df_column]][i]
-    )
+    params[[p_column]][i] <- 1 -
+      mvtnorm::pmvt(
+        lower = rep(
+          -abs(stats::qt(params[[p_column]][i] / 2, df = params[[df_column]][i])),
+          nrow(vc)
+        ),
+        upper = rep(
+          abs(stats::qt(params[[p_column]][i] / 2, df = params[[df_column]][i])),
+          nrow(vc)
+        ),
+        corr = vc,
+        df = params[[df_column]][i]
+      )
   }
   # clean up - remove temporary column
   params[[".sup_df"]] <- NULL
@@ -145,15 +172,19 @@
 .p_adjust_esarey <- function(x) {
   # only for slopes
   if (!inherits(x, c("estimate_slopes", "marginaleffects_slopes"))) {
-    insight::format_error("The `esarey` p-value adjustment is only available for Johnson-Neyman intervals, i.e. when calling `estimate_slopes()` with an interaction term of two numeric predictors.") # nolint
+    insight::format_error(
+      "The `esarey` p-value adjustment is only available for Johnson-Neyman intervals, i.e. when calling `estimate_slopes()` with an interaction term of two numeric predictors."
+    ) # nolint
   }
   # get names of interaction terms
   pred <- attributes(x)$trend
   mod <- attributes(x)$by
 
   # check for valid values - all must be numeric
   if (!all(vapply(attributes(x)$datagrid[c(pred, mod)], is.numeric, logical(1)))) {
-    insight::format_error("The `esarey` p-value adjustment is only available for Johnson-Neyman intervals, i.e. when calling `estimate_slopes()` with an interaction term of two numeric predictors.") # nolint
+    insight::format_error(
+      "The `esarey` p-value adjustment is only available for Johnson-Neyman intervals, i.e. when calling `estimate_slopes()` with an interaction term of two numeric predictors."
+    ) # nolint
   }
 
   int <- paste0(pred, ":", mod)
@@ -186,14 +217,17 @@
   # produces a sequence of marginal effects
   marginal_effects <- beta_pred + beta_int * range_sequence
   # SEs of those marginal effects
-  me_ses <- sqrt(vcov_pred + (range_sequence^2) * vcov_int + 2 * range_sequence * vcov_pred_int)
+  me_ses <- sqrt(
+    vcov_pred + (range_sequence^2) * vcov_int + 2 * range_sequence * vcov_pred_int
+  )
 
   # t-values across range of marginal effects
   statistic <- marginal_effects / me_ses
   # degrees of freedom
   dof <- insight::get_df(model, type = "wald")
   # Get the minimum p values used in the adjustment
-  pvalues <- 2 * pmin(stats::pt(statistic, df = dof), (1 - stats::pt(statistic, df = dof)))
+  pvalues <- 2 *
+    pmin(stats::pt(statistic, df = dof), (1 - stats::pt(statistic, df = dof)))
   # Multipliers
   multipliers <- seq_along(marginal_effects) / length(marginal_effects)
   # Order the pvals

R/summary.R

@@ -4,7 +4,9 @@ summary.estimate_slopes <- function(object, verbose = TRUE, ...) {
   by <- attributes(object)$by
 
   if (verbose && nrow(out) < 50) {
-    insight::format_alert("There might be too few data to accurately determine intervals. Consider setting `length = 100` (or larger) in your call to `estimate_slopes()`.") # nolint
+    insight::format_alert(
+      "There might be too few data to accurately determine intervals. Consider setting `length = 100` (or larger) in your call to `estimate_slopes()`."
+    ) # nolint
   }
 
   # Add "Confidence" col based on the sig index present in the data
@@ -22,6 +24,10 @@ summary.estimate_slopes <- function(object, verbose = TRUE, ...) {
     out <- .estimate_slope_parts(out, by)
   }
 
+  if (anyNA(out$Confidence) && verbose) {
+    insight::format_warning("Significance could not be determined for some intervals.")
+  }
+
   attributes(out) <- utils::modifyList(attributes(object), attributes(out))
   class(out) <- c("summary_estimate_slopes", "data.frame")
   attr(out, "table_title") <- c("Johnson-Neymann Intervals", "blue")
@@ -40,15 +46,23 @@ summary.reshape_grouplevel <- function(object, ...) {
 
 # Utilities ===============================================================
 
-
 .estimate_slope_parts <- function(out, by) {
   # mark all "changes" from negative to positive and vice versa
   index <- 1
   out$switch <- index
   index <- index + 1
 
   for (i in 2:nrow(out)) {
-    if (out$Direction[i] != out$Direction[i - 1] || out$Confidence[i] != out$Confidence[i - 1]) {
+    # we only proceed if we have non-missing values
+    do_switch <- (!is.na(out$Direction[i - 1]) & !is.na(out$Direction[i])) &&
+      out$Direction[i] != out$Direction[i - 1]
+
+    if (!do_switch) {
+      do_switch <- (!is.na(out$Confidence[i - 1]) & !is.na(out$Confidence[i])) &&
+        out$Confidence[i] != out$Confidence[i - 1]
+    }
+
+    if (do_switch) {
       out$switch[i:nrow(out)] <- index # styler: off
       index <- index + 1
     }
@@ -57,14 +71,17 @@ summary.reshape_grouplevel <- function(object, ...) {
   # split into "switches"
   parts <- split(out, out$switch)
 
-  do.call(rbind, lapply(parts, function(i) {
-    data.frame(
-      Start = i[[by]][1],
-      End = i[[by]][nrow(i)],
-      Direction = i$Direction[1],
-      Confidence = i$Confidence[1]
-    )
-  }))
+  do.call(
+    rbind,
+    lapply(parts, function(i) {
+      data.frame(
+        Start = i[[by]][1],
+        End = i[[by]][nrow(i)],
+        Direction = i$Direction[1],
+        Confidence = i$Confidence[1]
+      )
+    })
+  )
 }
 
 
@@ -83,18 +100,27 @@ summary.reshape_grouplevel <- function(object, ...) {
 
   if (confidence == "auto") {
     # TODO: make sure all of these work
-    if ("BF" %in% names(x)) confidence <- "BF"
-    if ("p" %in% names(x)) confidence <- "p"
+    if ("BF" %in% names(x)) {
+      confidence <- "BF"
+    }
+    if ("p" %in% names(x)) {
+      confidence <- "p"
+    }
     if ("pd" %in% names(x)) confidence <- "pd"
   }
 
-  switch(confidence,
+  switch(
+    confidence,
     p = tools::toTitleCase(effectsize::interpret_p(x$p, ...)),
     BF = tools::toTitleCase(effectsize::interpret_bf(x$BF, ...)),
     pd = tools::toTitleCase(effectsize::interpret_pd(x$pd, ...)),
     {
       # Based on CI
-      out <- ifelse((x$CI_high < 0 & x$CI_low < 0) | (x$CI_high > 0 & x$CI_low > 0), "Significant", "Uncertain")
+      out <- ifelse(
+        (x$CI_high < 0 & x$CI_low < 0) | (x$CI_high > 0 & x$CI_low > 0),
+        "Significant",
+        "Uncertain"
+      )
       factor(out, levels = c("Uncertain", "Significant"))
     }
   )

vignettes/introduction_comparisons_3.Rmd

@@ -270,12 +270,10 @@ plot(pr) + ggplot2::facet_wrap(~Illiteracy)
 For Johnson-Neyman intervals, which arise from interactions between two numeric predictors, we are essentially conducting multiple comparisons across the values of the moderator. To account for this, p-values can be adjusted to avoid an inflation of type-I errors. For `estimate_slopes()`, you can use the `p_adjust` argument for this purpose. While common methods like `"holm"` or `"bonferroni"` are available, for the specific case of Johnson-Neyman intervals, the `"esarey"` or `"sup-t"` adjustments are particularly recommended. The `"sup-t"` method, for instance, computes simultaneous confidence bands (Montiel Olea & Plagborg-Møller, 2019), providing a more rigorous test for the interval.
 
 ```{r}
-# we will force to calculate slopes at 200 values for "Illiteracy" using `length`
 slopes <- estimate_slopes(
   m_mod,
   "Murder",
   by = "Illiteracy",
-  length = 200,
   p_adjust = "esarey"
 )
 summary(slopes)