Merge branch 'master' into upkeep

jacob-long · Jul 28, 2024 · 64bd43a · 64bd43a
2 parents de7f328 + 6a9b659
commit 64bd43a
Show file tree

Hide file tree

Showing 77 changed files with 580 additions and 516 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -42,4 +42,4 @@ Enhances:
   rstanarm
 VignetteBuilder: knitr
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.2.9000
diff --git a/NEWS.md b/NEWS.md
@@ -3,12 +3,15 @@
 Enhancements:
 
 * `sim_slopes()` now supports non-continuous variables in the `pred` argument.
+* `sim_slopes()` now has an `at` argument, allowing you to specify an exact,
+perhaps non-centered, level for variables not involved in the interaction.
 * `interact_plot()` now has provisional support for factor predictors (`pred`).
 Users will receive a message because of the possibility for unexpected 
 behavior. `cat_plot()` likewise has support for continuous moderators. (#54)
 * Website and some documentation have been revamped and upgraded.
 * Users can now change the axis labels for `johnson_neyman()` plots via the 
 arguments `y.label` and `modx.label`. (#56)
+* Models produced by the `panelr` package are better supported.
 
 Bug fixes:
 

diff --git a/R/cat_plot.R b/R/cat_plot.R
@@ -97,7 +97,7 @@
 #'
 #'   The function is designed for two and three-way interactions. For
 #'   additional terms, the
-#'   \code{\link[effects]{effects}} package may be better suited to the task.
+#'   \code{effects} package may be better suited to the task.
 #'
 #'   This function supports nonlinear and generalized linear models and by
 #'   default will plot them on

diff --git a/R/int_utils.R b/R/int_utils.R
@@ -354,18 +354,19 @@ auto_mod_vals <-
 ## Centering
 
 center_ss <- function(d, weights, facvars = NULL, fvars, pred, resp, modx,
-                        survey, design = NULL, mod2, wname, offname, centered) {
+                      survey, design = NULL, mod2, wname, offname, centered,
+                      at = NULL) {
 
   # Just need to pick a helper function based on survey vs no survey
   if (survey == TRUE) {
 
     out <- center_ss_survey(d, weights, facvars, fvars, pred, resp, modx,
-                              survey, design, mod2, wname, offname, centered)
+                            survey, design, mod2, wname, offname, centered, at)
 
   } else {
 
     out <- center_ss_non_survey(d, weights, facvars, fvars, pred, resp, modx,
-                                  mod2, wname, offname, centered)
+                                mod2, wname, offname, centered, at)
 
   }
 
@@ -377,7 +378,7 @@ center_ss <- function(d, weights, facvars = NULL, fvars, pred, resp, modx,
 ## If not svydesign, centering is fairly straightforward
 
 center_ss_non_survey <- function(d, weights, facvars = NULL, fvars, pred,
-                                   resp, modx, mod2, wname, offname, centered) {
+                                 resp, modx, mod2, wname, offname, centered, at) {
 
   omitvars <- c(pred, resp, modx, mod2, wname, offname)
 
@@ -389,8 +390,8 @@ center_ss_non_survey <- function(d, weights, facvars = NULL, fvars, pred,
   if (centered[1] != "all" && centered[1] != "none") {
 
     if (any(omitvars %in% centered)) {
-      warning("Moderators, outcome variables, and weights/offsets",
-              " cannot be centered.")
+      warn_wrap("Moderators, outcome variables, and weights/offsets
+                cannot be centered.")
       centered <- centered[centered %nin% omitvars]
     }
     if (length(centered) > 0) {
@@ -422,7 +423,7 @@ center_ss_non_survey <- function(d, weights, facvars = NULL, fvars, pred,
     # Dealing with two-level factors that aren't part
     # of an interaction/focal pred
     for (v in fv2) {
-      if (is.factor(d[[v]]) & length(unique(d[[v]])) == 2) {
+      if (is.factor(d[[v]]) && length(unique(d[[v]])) == 2) {
 
         facvars <- c(facvars, v)
 
@@ -431,8 +432,9 @@ center_ss_non_survey <- function(d, weights, facvars = NULL, fvars, pred,
 
   }
 
-  # Fixes a data type error with predict() later
-  d <- as.data.frame(d)
+  if (!is.null(at)) {
+    d <- set_at(at = at, d = d)
+  }
 
   out <- list(d = d, facvars = facvars, fvars = fvars, design = NULL)
 
@@ -444,9 +446,9 @@ center_ss_non_survey <- function(d, weights, facvars = NULL, fvars, pred,
 
 center_ss_survey <- function(d, weights, facvars = NULL, fvars, pred, resp,
                              modx, survey, design, mod2, wname, offname,
-                             centered) {
+                             centered, at) {
 
-  omitvars <- c(pred, resp, modx, mod2, wname, offname)
+  omitvars <- c(pred, resp, modx, mod2, wname, offname, names(at))
 
   # Dealing with two-level factors that aren't part of an interaction
   # /focal pred
@@ -456,8 +458,8 @@ center_ss_survey <- function(d, weights, facvars = NULL, fvars, pred, resp,
   if (centered[1] != "all" && centered[1] != "none") {
 
     if (any(omitvars %in% centered)) {
-      warning("Moderators, outcome variables, and weights/offsets",
-              " cannot be centered.")
+      warn_wrap("Moderators, outcome variables, and weights/offsets cannot be 
+                 centered.")
       centered <- centered[centered %nin% omitvars]
     }
     design <- gscale(vars = centered, data = design, center.only = TRUE)
@@ -475,33 +477,48 @@ center_ss_survey <- function(d, weights, facvars = NULL, fvars, pred, resp,
     }
 
   } else if (centered == "none") {
-
     # Dealing with two-level factors that aren't part
     # of an interaction/focal pred
     for (v in fv2) {
       if (is.factor(d[[v]]) && length(unique(d[[v]])) == 2) {
-
         facvars <- c(facvars, v)
-
       }
     }
 
   } else if (centered == "all") {
-
     # Center all non-focal
     ndfvars <- fvars[fvars %nin% omitvars]
 
     if (length(ndfvars) > 0) {
       design <- gscale(vars = ndfvars, data = design, center.only = TRUE)
       d <- design$variables
     }
+  }
 
+  if (!is.null(at)) {
+    d <- set_at(at = at, d = d)
   }
 
   out <- list(d = d, design = design, facvars = facvars, fvars = fvars)
 
   return(out)
+}
 
+#### Deal with at variables #################################################
+set_at <- function(at, d) {
+  for (v in names(at)) {
+    if (v %nin% names(d)) stop_wrap("`at` variable ", v, " not found in data.")
+    if (!is.numeric(d[[v]])) {
+      warn_wrap("Inclusion of non-numeric variable ", v, " in `at` argument
+                is not currently supported. As an alternative, treat the 
+                variable as a factor and use the relevel() function to
+                set this value as its reference level before fitting your
+                model.")
+    } else {
+      d[[v]] <- d[[v]] - at[[v]]
+    }
+  }
+  return(d)
 }
 
 #### Send deprecation warnings ##############################################
@@ -923,7 +940,6 @@ drop_factor_levels <- function(d, var, values, labels) {
 
 }
 
-
 # get_contrasts <- function(model) {
 #   form <- as.formula(formula(model))
 #   as.data.frame(t(attr(terms(form), "factors")))

diff --git a/R/interact_plot.R b/R/interact_plot.R
@@ -257,10 +257,13 @@
 #'
 #' @author Jacob Long \email{jacob.long@@sc.edu}
 #'
-#' @seealso \code{\link[rockchalk]{plotSlopes}} from \pkg{rockchalk} performs a
+#' @seealso \code{plotSlopes} from \code{rockchalk} performs a
 #'   similar function, but
 #'   with R's base graphics---this function is meant, in part, to emulate
 #'   its features.
+#' 
+#'   Functions from the `margins` and `sjPlot` packages may also be useful
+#'   if this one isn't working for you.
 #'
 #'   \code{\link{sim_slopes}} performs a simple slopes analysis with a similar
 #'   argument syntax to this function.

diff --git a/R/johnson_neyman.R b/R/johnson_neyman.R
@@ -190,7 +190,7 @@ johnson_neyman <- function(model, pred, modx, vmat = NULL, alpha = 0.05,
   # Construct interaction term
   ## Create helper function to use either fixef() or coef() depending on input
   get_coef <- function(mod) {
-    if (inherits(mod, "merMod") | inherits(mod, "brmsfit")) {
+    if (inherits(mod, "merMod") || inherits(mod, "brmsfit")) {
       coef <- lme4::fixef(model)
       if (inherits(mod, "brmsfit")) {
         coefs <- coef[,1, drop = TRUE]
@@ -209,25 +209,44 @@ johnson_neyman <- function(model, pred, modx, vmat = NULL, alpha = 0.05,
     } else {
       pred_names <- pred
   }
-  ## Hard to predict which order lm() will have the predictors in
+
+  ## Old comment: Hard to predict which order lm() will have the predictors in
+  ## New: Note that this information is actually in the terms object, but 
+  ## I'm going to leave it this way for now since I need to work around 
+  ## panelr compatibility anyway — it manually calculates the interaction 
+  ## term rather than specifying it via the design matrix
   # first possible ordering
-  intterm1 <- paste(pred_names, ":", modx, sep = "")
+  intterm1 <- paste0(pred_names, ":", modx)
   # is it in the coef names?
   intterm1tf <- any(intterm1 %in% names(get_coef(model)))
   # second possible ordering
-  intterm2 <- paste(modx, ":", pred_names, sep = "")
+  intterm2 <- paste0(modx, ":", pred_names)
   # is it in the coef names?
   intterm2tf <- any(intterm2 %in% names(get_coef(model)))
   # Taking care of other business, creating coefs object for later
   coefs <- get_coef(model)
 
-
   ## Now we know which of the two is found in the coefficents
   # Using this to get the index of the TRUE one
   inttermstf <- c(intterm1tf, intterm2tf)
   intterms <- c(intterm1, intterm2) # Both names, want to keep one
   intterm <- intterms[which(inttermstf)] # Keep the index that is TRUE
-
+  # See if we can recover if intterm isn't found, this is motivated by 
+  # desire for compatibility with panelr
+  if (length(intterm) == 0) {
+    intterm1 <- paste0("`", pred_names, ":", modx, "`")
+    intterm2 <- paste0("`", modx, ":", pred_names, "`")
+    intterm1tf <- any(intterm1 %in% names(get_coef(model)))
+    intterm2tf <- any(intterm2 %in% names(get_coef(model)))
+    inttermstf <- c(intterm1tf, intterm2tf)
+    intterms <- c(intterm1, intterm2) # Both names, want to keep one
+    intterm <- intterms[which(inttermstf)] 
+
+    if (length(intterm) == 0) {
+      stop_wrap("Could not find interaction term in the model, so 
+                Johnson-Neyman interval could not be calculated.")
+    }
+  }
   # Getting the range of the moderator
   modrange <- range(model.frame(model)[,un_bt(modx)])
   modrangeo <- range(model.frame(model)[,un_bt(modx)]) # for use later
@@ -288,7 +307,7 @@ johnson_neyman <- function(model, pred, modx, vmat = NULL, alpha = 0.05,
     test <- 0
     i <- 1 + length(marginal_effects)
 
-    while (test == 0 & i > 1) {
+    while (test == 0 && i > 1) {
 
       i <- i - 1
       test <- min(ps[ps_o][1:i] <= multipliers[i] * (alpha * 2))