# R/pseudo_jn.r In manymome: Mediation, Moderation and Moderated-Mediation After Model Fitting

#### Documented in print.pseudo_johnson_neymanpseudo_johnson_neyman

```#' @title Pseudo Johnson-Neyman Probing
#'
#' @description Use the pseudo
#' Johnson-Neyman approach (Hayes, 2022)
#' to find the range of values of a
#' moderator in which the conditional
#' effect is not significant.
#'
#' @details This function uses the
#' pseudo Johnson-Neyman approach
#' proposed by Hayes (2022) to find the
#' values of a moderator at which a
#' conditional effect is
#' "nearly just significant" based on
#' confidence interval. If an effect is
#' moderated, there will be two such
#' points (though one can be very large
#' or small) forming a range.
#' The conditional effect
#' is not significant within this range,
#' and significant outside this range,
#' based on the confidence interval.
#'
#' This function receives the output
#' of [cond_indirect_effects()]
#' and search for, within
#' a specific range, the two values of
#' the moderator at which
#' the conditional effect is "nearly just significant",
#' that is, the confidence interval
#' "nearly touches" zero.
#'
#' Note that numerical method is used
#' to find the points. Therefore,
#' strictly speaking, the effects at
#' the end points are still either
#' significant or not significant, even
#' if the confidence limit is very close
#' to zero.
#'
#' ## Supported Methods
#'
#' This function supports models fitted
#' by [lm()], [lavaan::sem()],
#' and [semTools::sem.mi()]. This function
#' also supports both bootstrapping
#' and Monte Carlo confidence intervals.
#' It also supports conditional
#' direct paths (no mediator) and
#' conditional indirect paths (with one
#' or more mediator), with `x` and/or
#' `y` standardized.
#'
#' ## Requirements
#'
#' To be eligible for using this function,
#' one form of confidence intervals
#' (e.g, bootstrapping or Monte Carlo)
#' must has been requested (e.g.,
#' setting `boot_ci = TRUE` or
#' `mc_ci = TRUE`) when calling
#' [cond_indirect_effects()].
#'
#' The confidence level of the confidence
#' [cond_indirect_effects()] will be used
#' by this function.
#'
#' ## Possible failures
#'
#' Even if a path has only one moderator,
#' it is possible that no solution, or
#' more than one solution, is/are found
#' if the relation between this moderator
#' and the conditional effect is not linear.
#'
#' the conditional effect is significant
#' over a wide range of value of the
#' moderator.
#'
#' It is advised to use [plot_effect_vs_w()]
#' to examine the relation between the
#' effect and the moderator first before
#' calling this function.
#'
#' ## Speed
#'
#' Note that, for conditional indirect
#' effects, the search can be slow
#' because the confidence interval needs
#' to be recomputed for each new value
#' of the moderator.
#'
#' ## Limitations
#'
#' - This function currently only supports
#' a path with only one moderator,
#'
#' - This function does not yet support
#' multigroup models.
#'
#' @return
#' A list of the class `pseudo_johnson_neyman`
#' (with a print method, [print.pseudo_johnson_neyman()]).
#' It has these major elements:
#'
#' - `cond_effects`: An output of
#' [cond_indirect_effects()] for the
#' two levels of the moderator found.
#'
#' - `w_min_valid`: Logical. If `TRUE`,
#'  the conditional effect is just
#'  significant at the lower level of
#'  the moderator found,
#'  and so is significant below this point.
#'  If `FALSE`, then the lower level of
#'  the moderator found is just the
#'  lower bound of the range searched,
#'  that is, `w_lower`.
#'
#' - `w_max_valid`: Logical. If `TRUE`,
#'  the conditional effect is just
#'  significant at the higher level of
#'  the moderator found,
#'  and so is significant above this point.
#'  If `FALSE`, then the higher level of
#'  the moderator found is just the
#'  upper bound of the range searched,
#'  that is, `w_upper`.
#'
#' @param object A
#' `cond_indirect_effects`-class object,
#' which is the output of [cond_indirect_effects()].
#'
#' @param w_lower The smallest value of
#' the moderator when doing the search.
#' If set to `NULL,` the default, it
#' will be 10 standard deviations
#' below mean, which should be small
#' enough.
#'
#' @param w_upper The largest value of
#' the moderator when doing the search.
#' If set to `NULL,` the default, it
#' will be 10 standard deviations
#' above mean, which should be large
#' enough.
#'
#' @param optimize_method The optimization
#' method to be used. Either
#' `"uniroot"` (the default) or
#' `"optimize"`, corresponding to
#' [stats::uniroot()] and
#' [stats::optimize()], respectively.
#'
#' @param extendInt Used by
#' [stats::uniroot()]. If `"no"`, then
#' search will be conducted strictly
#' within `c(w_lower, w_upper)`. Otherwise,
#' the range is extended based on this
#' for details.
#'
#' @param tol The tolerance level used
#' by both [stats::uniroot()] and
#' [stats::optimize()].
#'
#' @param x The output of
#' [pseudo_johnson_neyman()].
#'
#' @param digits Number of digits to
#' display. Default is 3.
#'
#' @param ... Other arguments. Not used.
#'
#' @references
#' Hayes, A. F. (2022). *Introduction to mediation, moderation, and conditional process analysis: A regression-based approach* (Third edition). The Guilford Press.
#'
#'
#' @seealso [cond_indirect_effects()]
#'
#' @examples
#'
#' library(lavaan)
#'
#' dat <- data_med_mod_a
#' dat\$wx <- dat\$x * dat\$w
#' mod <-
#' "
#' m ~ x + w + wx
#' y ~ m + x
#' "
#' fit <- sem(mod, dat)
#'
#' # In real research, R should be 2000 or even 5000
#' # In real research, no need to set parallel and progress to FALSE
#' # Parallel processing is enabled by default and
#' # progress is displayed by default.
#' boot_out <- do_boot(fit,
#'                     R = 50,
#'                     seed = 4314,
#'                     parallel = FALSE,
#'                     progress = FALSE)
#' out <- cond_indirect_effects(x = "x", y = "y", m = "m",
#'                              wlevels = "w",
#'                              fit = fit,
#'                              boot_ci = TRUE,
#'                              boot_out = boot_out)
#'
#' # Visualize the relation first
#' plot_effect_vs_w(out)
#'
#' out_jn <- pseudo_johnson_neyman(out)
#' out_jn
#'
#' # Plot the range
#' plot_effect_vs_w(out_jn\$cond_effects)
#'
#' @export
#'

pseudo_johnson_neyman <- function(object = NULL,
w_lower = NULL,
w_upper = NULL,
optimize_method = c("uniroot", "optimize"),
extendInt = c("no", "yes", "downX", "upX"),
tol = .Machine\$double.eps^0.25) {
optimize_method <- match.arg(optimize_method)
if (inherits(object, "cond_indirect_effects")) {
if (cond_indirect_effects_has_groups(object)) {
stop("Multigroup models are not supported.")
}
# Do not use update() because it is not reliable in this context
x <- attr(object, "x", exact = TRUE)
y <- attr(object, "y", exact = TRUE)
m <- attr(object, "m", exact = TRUE)
wlevels0 <- attr(object, "wlevels", exact = TRUE)
w <- colnames(wlevels0)[1]
fit <- attr(object, "fit", exact = TRUE)
full_output_i <- attr(object, "full_output", exact = TRUE)[[1]]
out_call <- attr(object, "call", exact = TRUE)
boot_out <- attr(object, "boot_out", exact = TRUE)
mc_out <- attr(object, "mc_out", exact = TRUE)
has_boot_out <- !is.null(boot_out)
has_mc_out <- !is.null(mc_out)
ci_type <- ifelse(has_boot_out, "boot",
ifelse(has_mc_out, "mc",
NA))
boot_ci <- has_boot_out
mc_ci <- has_mc_out
standardized_x <- full_output_i\$standardized_x
standardized_y <- full_output_i\$standardized_y
if (ncol(wlevels0) != 1) {
stop("Support only one moderator.")
}
if (!is.numeric(wlevels0[, 1])) {
stop("Support only numeric moderator.")
}
if (!has_boot_out && !has_mc_out) {
stop("Confidence intervals not in 'object'.")
}
} else {
stop("'object' needs to be an output of cond_indirect_effects().")
}
w_tmp <- mod_levels(w = w,
fit = fit,
w_method = "sd",
sd_from_mean = c(-10, 10))
w_tmp <- range(w_tmp[, w])
if (is.null(w_lower)) {
w_lower <- w_tmp[1]
}
if (is.null(w_upper)) {
w_upper <- w_tmp[2]
}
wlevel_interval <- c(w_lower, w_upper)
w_lower_ci <- pseudo_johnson_neyman_one_bound(w0 = w_lower,
object = object,
type = "ci")
w_upper_ci <- pseudo_johnson_neyman_one_bound(w0 = w_upper,
object = object,
type = "ci")
w_increasing <- (w_upper_ci[1, 2] > w_lower_ci[1, 2])
if (w_increasing) {
w_lb_0 <- list(w = w_lower, objective = w_lower_ci[1, 2])
w_ub_0 <- list(w = w_upper, objective = w_upper_ci[1, 1])
} else {
w_lb_0 <- list(w = w_upper, objective = w_upper_ci[1, 2])
w_ub_0 <- list(w = w_lower, objective = w_lower_ci[1, 1])
}
if (optimize_method == "uniroot") {
w_lb <- tryCatch(pseudo_johnson_neyman_uniroot(object = object,
wlevel_interval = wlevel_interval,
which = "lower",
extendInt = extendInt,
tol = tol),
error = function(e) e)
if (inherits(w_lb, "error")) {
w_lb <- w_lb_0
names(w_lb) <- c("root", "f.root")
}
w_ub <- tryCatch(pseudo_johnson_neyman_uniroot(object = object,
wlevel_interval = wlevel_interval,
which = "upper",
extendInt = extendInt,
tol = tol),
error = function(e) e)
if (inherits(w_ub, "error")) {
w_ub <- w_ub_0
names(w_ub) <- c("root", "f.root")
}
w_df <- rbind(data.frame(w_lb[c("root", "f.root")]),
data.frame(w_ub[c("root", "f.root")]))
}
if (optimize_method == "optimize") {
w_lb <- tryCatch(pseudo_johnson_neyman_optimize(object = object,
wlevel_interval = wlevel_interval,
which = "lower",
tol = tol),
error = function(e) e)
w_ub <- tryCatch(pseudo_johnson_neyman_optimize(object = object,
wlevel_interval = wlevel_interval,
which = "upper",
tol = tol),
error = function(e) e)
if (inherits(w_lb, "error") || inherits(w_ub, "error")) {
stop("The search failed. Try setting optimize_method to 'uniroot'.")
}
w_df <- rbind(data.frame(w_lb),
data.frame(w_ub))
}
colnames(w_df) <- c("w", "fx")
w_df <- w_df[order(w_df\$w, decreasing = TRUE), ]
w_min <- min(w_df\$w)
w_max <- max(w_df\$w)
w_min_valid <- isTRUE(all.equal(w_df[2, "fx"], 0, tolerance = 1e-5))
w_max_valid <- isTRUE(all.equal(w_df[1, "fx"], 0, tolerance = 1e-5))
w_tmp <- mod_levels(w = w,
fit = fit,
values = c(Low = w_min, High = w_max))
w_lower <- min(w_min, w_lower)
w_upper <- max(w_max, w_upper)
out_cond <- cond_indirect_effects(wlevels = w_tmp,
x = x,
y = y,
m = m,
fit = fit,
boot_ci = boot_ci,
boot_out = boot_out,
mc_ci = mc_ci,
mc_out = mc_out,
standardized_x = standardized_x,
standardized_y = standardized_y)
out <- list(cond_effects = out_cond,
w_min_valid = w_min_valid,
w_max_valid = w_max_valid,
w_range_lb = w_min,
w_range_ub = w_max,
w_lower = w_lower,
w_upper = w_upper)
class(out) <- c("pseudo_johnson_neyman", class(out))
out
}

#' @noRd

pseudo_johnson_neyman_uniroot <- function(object,
wlevel_interval,
which = c("lower", "upper"),
extendInt = c("no", "yes", "downX", "upX"),
tol = .Machine\$double.eps^0.25,
maxiter = 1000) {
w_b <- stats::uniroot(pseudo_johnson_neyman_one_bound,
interval = wlevel_interval,
object = object,
which = which,
type = "limit",
extendInt = extendInt,
tol = tol,
maxiter = maxiter)
tmp <- pseudo_johnson_neyman_one_bound(w0 = w_b\$root,
object = object,
which = which,
type = "ci")
chk <- switch(which,
lower = (tmp[2] < 0) && isTRUE(all.equal(tmp[2], 0)),
upper = (tmp[1] > 0) && isTRUE(all.equal(tmp[1], 0)))
if (chk) {
lower = -1 * abs(tmp[2]),
upper = abs(tmp[1]))
w_b <- stats::uniroot(pseudo_johnson_neyman_one_bound,
interval = c(w_b\$root * (1 - 1e-2),
w_b\$root * (1 + 1e-2)),
object = object,
which = which,
type = "limit",
extendInt = "no",
tol = tol,
maxiter = maxiter)
}
w_b
}

#' @noRd

pseudo_johnson_neyman_optimize <- function(object,
wlevel_interval,
which = c("lower", "upper"),
tol = .Machine\$double.eps^0.25,
maximum = FALSE) {
w_b <- stats::optimize(pseudo_johnson_neyman_one_bound,
interval = wlevel_interval,
object = object,
which = which,
type = "distance",
maximum = maximum,
tol = tol)
tmp <- pseudo_johnson_neyman_one_bound(w0 = w_b\$minimum,
object = object,
which = which,
type = "ci")
chk <- switch(which,
lower = (tmp[2] < 0) && isTRUE(all.equal(tmp[2], 0)),
upper = (tmp[1] > 0) && isTRUE(all.equal(tmp[1], 0)))
if (chk) {
lower = sign(tmp[2]) * sqrt(abs(tmp[2])),
upper = sign(tmp[1]) * sqrt(abs(tmp[1])))
w_b <- stats::optimize(pseudo_johnson_neyman_one_bound,
interval = c(w_b\$minimum * (1 - 1e-2),
w_b\$minimum * (1 + 1e-2)),
object = object,
which = which,
type = "distance",
maximum = maximum,
tol = tol)
}
w_b
}

#' @describeIn pseudo_johnson_neyman Print
#' method for output of [pseudo_johnson_neyman()].
#'
#' @export

print.pseudo_johnson_neyman <- function(x, digits = 3, ...) {
out_cond <- x\$cond_effects
full_output_i <- attr(out_cond, "full_output", exact = TRUE)[[1]]
w <- names(full_output_i\$wvalues)[1]
ci_level <- full_output_i\$level
sig_level <- 1 - ci_level
w_range <- c(x\$w_lower, x\$w_upper)
w_range_lb <- min(w_range)
w_range_ub <- max(w_range)
w_min_valid <- x\$w_min_valid
w_max_valid <- x\$w_max_valid
w_found_lb <- x\$w_range_lb
w_found_ub <- x\$w_range_ub
w_lb_str <- formatC(w_found_lb, digits = digits, format = "f")
w_ub_str <- formatC(w_found_ub, digits = digits, format = "f")
w_range_lb_str <- formatC(w_range_lb, digits = digits, format = "f")
w_range_ub_str <- formatC(w_range_ub, digits = digits, format = "f")
cat("\n")
cat("== Pseudo Johnson-Neyman Probing ==\n")
cat("\n")
str_tmp <- character(0)
if (w_min_valid && w_max_valid) {
str_tmp <- c(str_tmp,
paste0("The conditional effect is ",
"not significant when ",
w, " is greater than ",
w_lb_str, " and less than ", w_ub_str,
", and is significant when ",
w, " is outside this range, ",
"at ", sig_level, " level of significance."))
}
if (!w_min_valid && w_max_valid) {
str_tmp <- c(str_tmp,
paste0("The conditional effect is ",
"not significant when ",
w, " is greater than ",
w_lb_str, " and less than ", w_ub_str,
", and is significant when ",
w, " is greater than ", w_ub_str, ", ",
"at ", sig_level, " level of significance."))
}
if (w_min_valid && !w_max_valid) {
str_tmp <- c(str_tmp,
paste0("The conditional effect is ",
"not significant when ",
w, " is greater than ",
w_lb_str, " and less than ", w_ub_str,
", and is significant when ",
w, " is less than ", w_lb_str, ", ",
"at ", sig_level, " level of significance."))
}
if (!w_min_valid && !w_max_valid) {
str_tmp <- c(str_tmp,
paste0("The conditional effect is ",
"not significant when ",
w, " is greater than ",
w_lb_str, " and less than ", w_ub_str,
", ",
"at ", sig_level, " level of significance."))
}
if (!w_min_valid || !w_max_valid) {
str_tmp <- c(str_tmp, "", "-- Note --")
}
if (!w_min_valid) {
str_tmp <- c(str_tmp,
paste0("- The lower bound of the range of ",
"nonsignificance is below ",
"the range being searched (",
w_range_lb_str, " to ", w_range_ub_str, "). ",
"Set a lower value for 'w_lower' if necessary."))
}
if (!w_max_valid) {
str_tmp <- c(str_tmp,
paste0("- The upper bound of the range of ",
"nonsignificance is above ",
"the range being searched (",
w_range_lb_str, " to ", w_range_ub_str, "). ",
"Set a higher value for 'w_upper' if necessary."))
}
str_tmp_final <- strwrap(str_tmp,
exdent = 0)
cat(str_tmp_final, sep = "\n")
print(out_cond)
invisible(x)
}

#' @noRd
# Search for one bound of the range
pseudo_johnson_neyman_one_bound <- function(w0,
object = NULL,
x = NULL,
y = NULL,
m = NULL,
w = NULL,
fit = NULL,
boot_ci = FALSE,
boot_out = NULL,
mc_ci = FALSE,
mc_out = NULL,
standardized_x = FALSE,
standardized_y = FALSE,
which = c("lower", "upper"),
type = c("distance", "limit", "ci", "est", "full_out"),
which <- match.arg(which)
type <- match.arg(type)
if (inherits(object, "cond_indirect_effects")) {
# Do not use update() because it is not reliable in this context
x <- attr(object, "x", exact = TRUE)
y <- attr(object, "y", exact = TRUE)
m <- attr(object, "m", exact = TRUE)
wlevels0 <- attr(object, "wlevels", exact = TRUE)
w <- colnames(wlevels0)[1]
fit <- attr(object, "fit", exact = TRUE)
full_output_i <- attr(object, "full_output", exact = TRUE)[[1]]
out_call <- attr(object, "call", exact = TRUE)
boot_out <- attr(object, "boot_out", exact = TRUE)
mc_out <- attr(object, "mc_out", exact = TRUE)
has_boot_out <- !is.null(boot_out)
has_mc_out <- !is.null(mc_out)
ci_type <- ifelse(has_boot_out, "boot",
ifelse(has_mc_out, "mc",
NA))
boot_ci <- has_boot_out
mc_ci <- has_mc_out
standardized_x <- full_output_i\$standardized_x
standardized_y <- full_output_i\$standardized_y
}
w_i <- w0
names(w_i) <- w
out <- cond_indirect(wvalues = w_i,
x = x,
y = y,
m = m,
fit = fit,
standardized_x = standardized_x,
standardized_y = standardized_y,
boot_ci = boot_ci,
boot_out = boot_out,
mc_ci = mc_ci,
mc_out = mc_out)
out1 <- switch(which, lower = stats::confint(out)[1, 2] + adj,
upper = stats::confint(out)[1, 1] - adj)
return(switch(type,
distance = out1^2,
limit = out1,
ci = stats::confint(out),
est = stats::coef(out),
full_out = out))
}
```

## Try the manymome package in your browser

Any scripts or data that you put into this service are public.

manymome documentation built on June 22, 2024, 9:34 a.m.