diff --git a/.Rbuildignore b/.Rbuildignore index 6516e16a..f769702b 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,4 +5,5 @@ ^\.github$ ^data-raw$ ^rebuild_vignettes.R$ -^\.lintr$ \ No newline at end of file +^\.lintr$ +^\.vscode$ \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 00000000..5270aa4f --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "cSpell.words": [ + "Multigroup" + ] +} \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 964db3fb..306b9b63 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: manymome Title: Mediation, Moderation and Moderated-Mediation After Model Fitting -Version: 0.1.14.1 +Version: 0.1.14.5 Authors@R: c(person(given = "Shu Fai", family = "Cheung", diff --git a/NEWS.md b/NEWS.md index d8ad7a8d..d3ce3816 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,47 @@ -# manymome 0.1.14.1 +# manymome 0.1.14.5 + +## New Features + +- Many functions have been updated to + work for multigroup models fitted + by `lavaan`. Most common tasks are + supported. There likely are functions + that may not yet work on + multigroup models. Checks will be + added to them to alert users. + For now, only some + functions (e.g., + `cond_indirect_effect()`) supports + multigroup models which have + one or more moderators within each + group, but these models are rare. + Functions that do not yet support + multigroup models (e.g, + `mod_levels()`) will raise an error + if used on a multigroup model. + (0.1.14.2 to 0.1.14.5) + +- Relaxed the requirement that only + different paths can be used in `+` + and `-`. They can now be used in + these operations, as they may be + paths in different groups in + multigroup models. (0.1.14.2) + +- The `plot`-method of + `cond_indirect_effects`-class objects + will be forced to be a tumble graph + if the lines for different groups + are to be plotted. In these cases, + the data within each group will be used, + including standardization. This + approach, though leading to results + different from those in single-group + model using the group as a moderator, + makes more sense for multigroup + models, in which the distribution of + variables are allowed to be different + between groups. (0.1.14.2) ## Miscellaneous diff --git a/R/all_indirect_paths.R b/R/all_indirect_paths.R index 49e2a060..eb7e6170 100644 --- a/R/all_indirect_paths.R +++ b/R/all_indirect_paths.R @@ -7,6 +7,17 @@ #' @details It makes use of [igraph::all_simple_paths()] #' to identify paths in a model. #' +#' ## Multigroup Models +#' +#' Since Version 0.1.14.2, support for +#' multigroup models has been added for models +#' fitted by `lavaan`. If a model has more +#' than one group and `group` is not +#' specified, than paths in all groups +#' will be returned. If `group` is +#' specified, than only paths in the +#' selected group will be returned. +#' #' @return #' [all_indirect_paths()] returns #' a list of the class `all_paths`. Each argument is a @@ -49,6 +60,19 @@ #' @param all_paths An `all_paths`-class object. For example, #' the output of [all_indirect_paths()]. #' +#' @param group Either the group number +#' as appeared in the [summary()] +#' or [lavaan::parameterEstimates()] +#' output of a [lavaan::lavaan-class] object, +#' or the group label as used in +#' the [lavaan::lavaan-class] object. +#' Used only when the number of +#' groups is greater than one. Default +#' is `NULL`. If not specified by the model +#' has more than one group, than paths +#' that appears in at least one group +#' will be included in the output. +#' #' @author Shu Fai Cheung #' #' @seealso [indirect_effect()], [lm2list()]. @@ -85,6 +109,32 @@ #' out3 #' names(out3) #' +#' # Multigroup models +#' +#' data(data_med_complicated_mg) +#' mod <- +#' " +#' m11 ~ x1 + x2 + c1 + c2 +#' m12 ~ m11 + c1 + c2 +#' m2 ~ x1 + x2 + c1 + c2 +#' y1 ~ m11 + m12 + x1 + x2 + c1 + c2 +#' y2 ~ m2 + x1 + x2 + c1 + c2 +#' " +#' fit <- sem(mod, data_med_complicated_mg, group = "group") +#' summary(fit) +#' +#' all_indirect_paths(fit, +#' x = "x1", +#' y = "y1") +#' all_indirect_paths(fit, +#' x = "x1", +#' y = "y1", +#' group = 1) +#' all_indirect_paths(fit, +#' x = "x1", +#' y = "y1", +#' group = "Group B") +#' #' @describeIn all_indirect_paths Enumerate all indirect paths. #' #' @order 1 @@ -94,26 +144,115 @@ all_indirect_paths <- function(fit = NULL, exclude = NULL, x = NULL, - y = NULL) { + y = NULL, + group = NULL) { fit_type <- cond_indirect_check_fit(fit) if (is.na(fit_type)) { stop("'fit' is not of a supported type.") } + ngroups <- 1 + group_number <- NULL + group_label <- NULL # Create an adjancey matrix if (identical(fit_type, "lavaan")) { - beta <- lavaan::lavInspect(fit)$beta + + ngroups <- lavaan::lavTech(fit, "ngroups") + if ((ngroups > 1) && !is.null(group)) { + group_labels_all <- lavaan::lavTech(fit, + "group.label") + if (is.numeric(group)) { + group_label <- group_labels_all[group] + group_number <- group + } else { + group_number <- match(group, group_labels_all) + group_label <- group + } + } + tmp <- lavaan::lavInspect(fit, + drop.list.single.group = FALSE) + tmp <- lapply(tmp, function(x) x$beta) + beta <- tmp } if (identical(fit_type, "lavaan.mi")) { - beta <- lavaan::lavInspect(fit)$beta + # TODO: + # Add support for multiple group models. + beta <- list(lavaan::lavInspect(fit)$beta) } if (identical(fit_type, "lm")) { - beta <- beta_from_lm(fit) + beta <- list(beta_from_lm(fit)) + } + if ((ngroups > 1) && + (identical(fit_type, "lavaan"))) { + group_labels_all <- lavaan::lavTech(fit, + "group.label") + if (is.null(group)) { + groups <- group_labels_all + group_numbers <- seq_len(ngroups) + } else { + beta <- beta[group_number] + groups <- group + group_numbers <- group_number + group_labels_all <- group_labels_all[group_number] + } + tmpfct <- function(adj_i, + group_i, + group_label_i, + group_number_i, + exclude = exclude, + x = x, + y = y, + fit = fit, + fit_type = fit_type) { + out <- all_indirect_paths_i(adj = adj_i, + exclude = exclude, + x = x, + y = y, + fit = fit, + fit_type = fit_type) + for (i in seq_along(out)) { + out[[i]]$group_label <- group_label_i + out[[i]]$group_number <- group_number_i + } + out + } + out3 <- mapply(tmpfct, + adj_i = beta, + group_i = groups, + group_label_i = group_labels_all, + group_number_i = group_numbers, + MoreArgs = list(exclude = exclude, + x = x, + y = y, + fit = fit, + fit_type = fit_type), + SIMPLIFY = FALSE) + out3 <- unlist(out3, + recursive = FALSE) + } else { + out3 <- all_indirect_paths_i(adj = beta[[1]], + exclude = exclude, + x = x, + y = y, + fit = fit, + fit_type = fit_type) } - adj <- beta + + class(out3) <- c("all_paths", class(out3)) + attr(out3, "call") <- match.call() + out3 + } + +#' @noRd + +all_indirect_paths_i <- function(adj, + exclude = NULL, + x = NULL, + y = NULL, + fit = NULL, + fit_type = NULL) { adj[adj > 0] <- 1 adj <- t(adj) - # Remove excluded variables if (is.character(exclude)) { adj <- adj[!(rownames(adj) %in% exclude), @@ -170,8 +309,6 @@ all_indirect_paths <- function(fit = NULL, # Format the output out3 <- lapply(out2, to_x_y_m) names(out3) <- sapply(out3, path_name) - class(out3) <- c("all_paths", class(out3)) - attr(out3, "call") <- match.call() out3 } @@ -188,9 +325,15 @@ all_paths_to_df <- function(all_paths) { all_y <- sapply(all_paths, function(x) x$y) all_m <- sapply(all_paths, function(x) x$m, simplify = FALSE) + all_group_label <- sapply(all_paths, function(x) x$group_label) + all_group_number <- sapply(all_paths, function(x) x$group_number) out <- data.frame(x = all_x, y = all_y) out$m <- all_m + if (!any(sapply(all_group_label, is.null))) { + out$group_label <- all_group_label + out$group_number <- all_group_number + } out } diff --git a/R/boot2est_lavaan.R b/R/boot2est_lavaan.R index ce986da3..28134d65 100644 --- a/R/boot2est_lavaan.R +++ b/R/boot2est_lavaan.R @@ -181,12 +181,20 @@ fit2boot_out_do_boot <- function(fit, environment(gen_boot_i_lavaan) <- parent.frame() boot_i <- gen_boot_i_lavaan(fit) } - dat_org <- lav_data_used(fit) - n <- nrow(dat_org) - boot_test <- suppressWarnings(boot_i(dat_org)) + dat_org <- lav_data_used(fit, + drop_list_single_group = TRUE) + ngp <- lavaan::lavTech(fit, "ngroups") + if (ngp == 1) { + n <- nrow(dat_org) + } else { + n <- sapply(dat_org, nrow) + } + boot_test <- suppressWarnings(boot_i(dat_org, + start = lavaan::parameterTable(fit)$start)) + # Increase the tolerance for mutliple group model if (!isTRUE(all.equal(unclass(lavaan::coef(fit)), lavaan::coef(boot_test)[names(lavaan::coef(fit))], - tolerance = sqrt(.Machine$double.eps * 1e04)))) { + tolerance = sqrt(.Machine$double.eps * 1e08)))) { stop(paste("Something is wrong.", "This function cannot reproduce the results.", "Please fit the model with se = 'boot'")) @@ -194,7 +202,12 @@ fit2boot_out_do_boot <- function(fit, ft <- lavaan::lavInspect(boot_test, "timing")$total requireNamespace("parallel", quietly = TRUE) if (!is.null(seed)) set.seed(seed) - ids <- replicate(R, sample.int(n, replace = TRUE), simplify = FALSE) + if (ngp == 1) { + ids <- replicate(R, sample.int(n, replace = TRUE), simplify = FALSE) + } else { + ids <- replicate(R, sapply(n, sample.int, replace = TRUE, simplify = TRUE), + simplify = FALSE) + } if (parallel) { if (is.numeric(ncores)) { ncores0 <- parallel::detectCores() @@ -361,8 +374,13 @@ set_est_i_lavaan <- function(est0, fit, p_free, est_df = NULL) { ptable <- as.data.frame(fit@ParTable) if (!is.null(est_df)) { est_df$est <- NULL - est0 <- merge(est_df, ptable[, c("lhs", "op", "rhs", "est")], - sort = FALSE) + if ("group" %in% colnames(est_df)) { + est0 <- merge(est_df, ptable[, c("lhs", "op", "rhs", "block", "group", "est")], + sort = FALSE) + } else { + est0 <- merge(est_df, ptable[, c("lhs", "op", "rhs", "est")], + sort = FALSE) + } class(est0) <- class(est_df) return(est0) } else { @@ -445,22 +463,49 @@ get_implied_i_lavaan <- function(est0, fit, fit_tmp = NULL) { delta = TRUE) } out <- lav_implied_all(fit) + ngroups <- lavaan::lavTech(fit, "ngroups") out_names <- names(out) implied_names <- names(implied) out1 <- out + # For multigroup models, use the format of lavaan::lav_model_implied() + # - Estimates than groups. for (x in out_names) { if (x %in% implied_names) { if (!is.null(implied[[x]][[1]])) { - out1[[x]][] <- implied[[x]][[1]] + if (ngroups > 1) { + for (j in seq_len(ngroups)) { + out1[[x]][[j]][] <- implied[[x]][[j]] + } + } else { + out1[[x]][] <- implied[[x]][[1]] + } } else { - out1[[x]][] <- NA + if (ngroups > 1) { + for (j in seq_len(ngroups)) { + out1[[x]][[j]][] <- NA + } + } else { + out1[[x]][] <- NA + } } } else { - out1[[x]][] <- NA + if (ngroups > 1) { + for (j in seq_len(ngroups)) { + out1[[x]][[j]][] <- NA + } + } else { + out1[[x]][] <- NA + } } } if (has_lv) { - out1$mean_lv <- implied$mean_lv[[1]] + if (ngroups > 1) { + for (j in seq_len(ngroups)) { + out1[["mean_lv"]][[j]][] <- NA + } + } else { + out1$mean_lv <- implied$mean_lv[[1]] + } } out1 } @@ -566,24 +611,34 @@ gen_boot_i_lavaan <- function(fit) { X_old <- fit_data@X - function(d, i = NULL) { + function(d, i = NULL, start = NULL) { force(fit_data) force(fit_model) force(fit_sampstats) force(fit_opts) force(fit_pt) force(fit) + fit_pt1 <- fit_pt + if (!is.null(start)) { + fit_pt1$start <- start + } if (is.null(i)) { return(lavaan::lavaan(slotData = fit_data, slotModel = fit_model, slotSampleStats = fit_sampstats, slotOptions = fit_opts, - slotParTable = fit_pt)) + slotParTable = fit_pt1)) } else { - # Support for multigroup models will be added later. - b_i <- list(i) + # 2024-03-29: Added support for multigroup models + if (!is.list(i)) { + b_i <- list(i) + } else { + b_i <- i + } X_new <- X_old - X_new[[1]] <- X_new[[1]][i, , drop = FALSE] + for (j in seq_along(X_new)) { + X_new[[j]] <- X_new[[j]][b_i[[j]], , drop = FALSE] + } fit_data_new <- lavaan::lav_data_update( lavdata = fit_data, newX = X_new, @@ -615,7 +670,7 @@ gen_boot_i_lavaan <- function(fit) { slotModel = fit_model_i, slotSampleStats = fit_sampstats_new, slotOptions = fit_opts, - slotParTable = fit_pt), + slotParTable = fit_pt1), error = function(e) e, warning = function(e) e) if (inherits(out, "error") || inherits(out, "warning")) { @@ -631,6 +686,8 @@ gen_boot_i_lavaan <- function(fit) { implied_stats = NA, ok = FALSE) } + # If ngroups > 1, + # cov, mean, and mean_lv are lists. implied <- list(cov = lavaan::lavInspect(out, "cov.all"), mean = lavaan::lavInspect(out, "mean.ov"), mean_lv = lavaan::lavInspect(out, "mean.lv")) diff --git a/R/coef_cond_indirect_effects.R b/R/coef_cond_indirect_effects.R index a7a70aac..89cf7778 100644 --- a/R/coef_cond_indirect_effects.R +++ b/R/coef_cond_indirect_effects.R @@ -62,14 +62,24 @@ #' @export coef.cond_indirect_effects <- function(object, ...) { + has_wlevels <- cond_indirect_effects_has_wlevels(object) + has_groups <- cond_indirect_effects_has_groups(object) wlevels <- attr(object, "wlevels") if ("std" %in% colnames(object)) { out <- object$std } else { out <- object$ind } - if (!is.null(wlevels)) { + if (has_wlevels && !has_groups) { names(out) <- rownames(wlevels) } + if (!has_wlevels && has_groups) { + tmp <- paste0(object$Group, " [", object$Group_ID, "]") + names(out) <- tmp + } + if (has_wlevels && has_groups) { + # TODO: + # - Support objects with both wlevels and groups + } out } diff --git a/R/cond_indirect.R b/R/cond_indirect.R index 2170ee6c..460cae10 100644 --- a/R/cond_indirect.R +++ b/R/cond_indirect.R @@ -104,6 +104,37 @@ #' as `R`, `seed`, and `parallel` will #' be ignored. #' +#' ## Multigroup Models +#' +#' Since Version 0.1.14.2, support for +#' multigroup models has been added for models +#' fitted by `lavaan`. Both bootstrapping +#' and Monte Carlo confidence intervals +#' are supported. When used on +#' a multigroup model: +#' +#' - For [cond_indirect()] and +#' [indirect_effect()], users need to +#' specify the `group` argument +#' (by number or label). When using +#' [cond_indirect_effects()], if +#' `group` is not set, all groups wil +#' be used and the indirect effect +#' in each group will be computed, +#' kind of treating group as a moderator. +#' +#' - For [many_indirect_effects()], +#' the paths can be generated from a +#' multigroup models. +#' +#' - Currently, [cond_indirect_effects()] +#' does not support a multigroup model +#' with moderators on the path selected. +#' The function [cond_indirect()] does +#' not have this limitation but users +#' need to manually specify the desired +#' value of the moderator(s). +#' #' @return [indirect_effect()] and #' [cond_indirect()] return an #' `indirect`-class object. @@ -113,8 +144,7 @@ #' #' These two classes of objects have #' their own print methods for printing -#' the results (see [print.indirect()] -#' and [print.cond_indirect_effects()]). +#' the results (see [print.indirect()] and [print.cond_indirect_effects()]). #' They also have a `coef` method for #' extracting the estimates #' ([coef.indirect()] and @@ -212,7 +242,7 @@ #' `TRUE`, `boot_out` is `NULL`, and #' bootstrap standard errors not #' requested if `fit` is a -#' [lavaan-class] object, this function +#' [lavaan::lavaan-class] object, this function #' will do bootstrapping on `fit`. `R` #' is the number of bootstrap samples. #' Default is 100. For Monte Carlo @@ -367,6 +397,26 @@ #' supplied, will override `boot_ci` #' and `mc_ci`. #' +#' @param group Either the group number +#' as appeared in the [summary()] +#' or [lavaan::parameterEstimates()] +#' output of a [lavaan::lavaan-class] object, +#' or the group label as used in +#' the [lavaan::lavaan-class] object. +#' Used only when the number of +#' groups is greater than one. Default +#' is `NULL`. +#' +#' @param groups Either a vector of +#' group numbers +#' as appeared in the [summary()] +#' or [lavaan::parameterEstimates()] +#' output of a [lavaan::lavaan-class] object, +#' or a vector of group labels as used in +#' the [lavaan::lavaan-class] object. +#' Used only when the number of +#' groups is greater than one. Default +#' is `NULL`. #' #' @seealso [mod_levels()] and #' [merge_mod_levels()] for generating @@ -441,7 +491,8 @@ cond_indirect <- function(x, ci_out = NULL, save_ci_full = FALSE, save_ci_out = TRUE, - ci_type = NULL) { + ci_type = NULL, + group = NULL) { fit_type <- cond_indirect_check_fit(fit) chkpath <- check_path(x = x, y = y, m = m, fit = fit, est = est) if (!chkpath) { @@ -569,7 +620,8 @@ cond_indirect <- function(x, standardized_y = standardized_y, get_prods_only = TRUE, data = fit_data, - expand = TRUE) + expand = TRUE, + group = group) } if (get_prods_only) return(prods) out0 <- indirect_i(x = x, @@ -581,7 +633,8 @@ cond_indirect <- function(x, wvalues = wvalues, standardized_x = standardized_x, standardized_y = standardized_y, - prods = prods) + prods = prods, + group = group) if (mc_ci) { out_mc <- mapply(indirect_i, est = lapply(mc_out, function(x) x$est), @@ -594,7 +647,8 @@ cond_indirect <- function(x, standardized_x = standardized_x, standardized_y = standardized_y, warn = FALSE, - prods = prods), + prods = prods, + group = group), SIMPLIFY = FALSE) if (save_mc_full) { out0$mc_full <- out_mc @@ -632,7 +686,8 @@ cond_indirect <- function(x, standardized_x = standardized_x, standardized_y = standardized_y, warn = FALSE, - prods = prods), + prods = prods, + group = group), SIMPLIFY = FALSE) if (save_boot_full) { out0$boot_full <- out_boot @@ -665,6 +720,34 @@ cond_indirect <- function(x, #' @export #' +#' @examples +#' +#' # Multigroup model for indirect_effect() +#' +#' dat <- data_med_mg +#' mod <- +#' " +#' m ~ x + c1 + c2 +#' y ~ m + x + c1 + c2 +#' " +#' fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, +#' group = "group") +#' +#' # If a model has more than one group, +#' # the argument 'group' must be set. +#' ind1 <- indirect_effect(x = "x", +#' y = "y", +#' m = "m", +#' fit = fit, +#' group = "Group A") +#' ind1 +#' ind2 <- indirect_effect(x = "x", +#' y = "y", +#' m = "m", +#' fit = fit, +#' group = 2) +#' ind2 +#' #' @describeIn cond_indirect Compute the #' indirect effect. A wrapper of #' [cond_indirect()]. Can be used when @@ -690,6 +773,7 @@ indirect_effect <- function(x, make_cluster_args = list(), progress = TRUE, save_boot_full = FALSE, + save_boot_out = TRUE, mc_ci = FALSE, mc_out = NULL, save_mc_full = FALSE, @@ -697,7 +781,8 @@ indirect_effect <- function(x, ci_out = NULL, save_ci_full = FALSE, save_ci_out = TRUE, - ci_type = NULL) { + ci_type = NULL, + group = NULL) { cond_indirect(x = x, y = y, m = m, @@ -716,6 +801,7 @@ indirect_effect <- function(x, make_cluster_args = make_cluster_args, progress = progress, save_boot_full = save_boot_full, + save_boot_out = save_boot_out, mc_ci = mc_ci, mc_out = mc_out, save_mc_full = save_mc_full, @@ -723,7 +809,8 @@ indirect_effect <- function(x, ci_out = ci_out, save_ci_full = save_ci_full, save_ci_out = save_ci_out, - ci_type = ci_type) + ci_type = ci_type, + group = group) } #' @param w_type Character. Whether the @@ -788,6 +875,22 @@ indirect_effect <- function(x, #' cond_indirect_effects(x = "x", y = "y", m = "m1", #' wlevels = w1levels, fit = fit) #' +#' # Multigroup models for cond_indirect_effects() +#' +#' dat <- data_med_mg +#' mod <- +#' " +#' m ~ x + c1 + c2 +#' y ~ m + x + c1 + c2 +#' " +#' fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, +#' group = "group") +#' +#' # If a model has more than one group, +#' # it will be used as a 'moderator'. +#' cond_indirect_effects(x = "x", y = "y", m = "m", +#' fit = fit) +#' #' @export #' #' @describeIn cond_indirect Compute the @@ -823,8 +926,31 @@ cond_indirect_effects <- function(wlevels, mc_out = NULL, ci_out = NULL, ci_type = NULL, + groups = NULL, ...) { + # Check the number of groups and handle multiple-group models + has_group <- FALSE + ngroups <- 1 + group_numbers <- NULL + group_labels <- NULL + if (inherits(fit, "lavaan")) { + ngroups <- lavaan::lavTech(fit, "ngroups") + if (ngroups > 1) { + has_group <- TRUE + tmp <- group_labels_and_numbers(groups = groups, + fit = fit) + group_numbers <- tmp$number + group_labels <- tmp$label + } else { + if (!is.null(groups)) { + stop("The model has only one group but 'groups' is set.") + } + } + } + # Check and process the levels of moderators + has_wlevels <- FALSE if (!missing(wlevels)) { + has_wlevels <- TRUE wlevels_check <- check_wlevels(wlevels) if (!is.null(wlevels_check)) { wlevels <- wlevels_check @@ -846,12 +972,21 @@ cond_indirect_effects <- function(wlevels, } } else { - stop("wlevels is required.") + wlevels <- NULL + if (!has_group) { + stop("wlevels is required for single-group models.") + } + } + if (has_group && has_wlevels) { + stop("Multiple group models with moderators not yet supported.", + "Will be supported soon.") + } + if (has_wlevels) { + k <- nrow(wlevels) + wlevels1 <- split(wlevels, seq_len(k)) + wlevels2 <- lapply(wlevels1, unlist) + names(wlevels2) <- rownames(wlevels) } - k <- nrow(wlevels) - wlevels1 <- split(wlevels, seq_len(k)) - wlevels2 <- lapply(wlevels1, unlist) - names(wlevels2) <- rownames(wlevels) fit_type <- cond_indirect_check_fit(fit) if ((fit_type == "lm") && !inherits(fit, "lm_list") && is.list(fit)) { @@ -930,78 +1065,247 @@ cond_indirect_effects <- function(wlevels, progress = progress) } } - prods <- cond_indirect(wvalues = wlevels2[[1]], - x = x, - y = y, - m = m, - fit = fit, - est = est, - implied_stats = implied_stats, - get_prods_only = TRUE, - ...) - out <- lapply(wlevels2, - function(wv, - x, - y, - m, - fit, - est, - implied_stats, - boot_ci, - boot_out, - R, - seed, - prods, - save_boot_out, - mc_ci, - mc_out, - save_mc_out, - ci_type, - ci_out, - save_ci_out, - ...) { - cond_indirect(wvalues = wv, - x = x, - y = y, - m = m, - fit = fit, - est = est, - implied_stats = implied_stats, - boot_ci = boot_ci, - boot_out = boot_out, - R = R, - seed = seed, - prods = prods, - save_boot_out = FALSE, - mc_ci = mc_ci, - mc_out = mc_out, - save_mc_out = FALSE, - ci_type = ci_type, - ci_out = ci_out, - save_ci_out = FALSE, - ...) - }, - x = x, - y = y, - m = m, - fit = fit, - est = est, - implied_stats = implied_stats, - boot_ci = boot_ci, - boot_out = boot_out, - R = R, - seed = seed, - prods = prods, - save_boot_out = FALSE, - mc_ci = mc_ci, - mc_out = mc_out, - save_mc_out = FALSE, - ci_type = ci_type, - ci_out = ci_out, - save_ci_out = FALSE, - ...) + # TODO: + # Revise cond_indirect and friends such that + # no need to have three very similar blocks. + if (has_wlevels && !has_group) { + prods <- cond_indirect(wvalues = wlevels2[[1]], + x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + get_prods_only = TRUE, + ...) + } + if (!has_wlevels && has_group) { + prods <- cond_indirect(x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + get_prods_only = TRUE, + group = 1, + ...) + } + if (has_wlevels && has_group) { + prods <- cond_indirect(wvalues = wlevels2[[1]], + x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + get_prods_only = TRUE, + group = 1, + ...) + } + # TODO: + # Revise cond_indirect and friends such that + # no need to have three very similar blocks. + if (has_wlevels && !has_group) { + out <- lapply(wlevels2, + function(wv, + x, + y, + m, + fit, + est, + implied_stats, + boot_ci, + boot_out, + R, + seed, + prods, + save_boot_out, + mc_ci, + mc_out, + save_mc_out, + ci_type, + ci_out, + save_ci_out, + ...) { + cond_indirect(wvalues = wv, + x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + boot_ci = boot_ci, + boot_out = boot_out, + R = R, + seed = seed, + prods = prods, + save_boot_out = FALSE, + mc_ci = mc_ci, + mc_out = mc_out, + save_mc_out = FALSE, + ci_type = ci_type, + ci_out = ci_out, + save_ci_out = FALSE, + ...) + }, + x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + boot_ci = boot_ci, + boot_out = boot_out, + R = R, + seed = seed, + prods = prods, + save_boot_out = FALSE, + mc_ci = mc_ci, + mc_out = mc_out, + save_mc_out = FALSE, + ci_type = ci_type, + ci_out = ci_out, + save_ci_out = FALSE, + ...) + } + if (!has_wlevels && has_group) { + out <- lapply(group_numbers, + function(gn, + x, + y, + m, + fit, + est, + implied_stats, + boot_ci, + boot_out, + R, + seed, + prods, + save_boot_out, + mc_ci, + mc_out, + save_mc_out, + ci_type, + ci_out, + save_ci_out, + ...) { + indirect_effect(x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + boot_ci = boot_ci, + boot_out = boot_out, + R = R, + seed = seed, + save_boot_out = FALSE, + mc_ci = mc_ci, + mc_out = mc_out, + save_mc_out = FALSE, + ci_type = ci_type, + ci_out = ci_out, + save_ci_out = FALSE, + group = gn, + ...) + }, + x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + boot_ci = boot_ci, + boot_out = boot_out, + R = R, + seed = seed, + save_boot_out = FALSE, + mc_ci = mc_ci, + mc_out = mc_out, + save_mc_out = FALSE, + ci_type = ci_type, + ci_out = ci_out, + save_ci_out = FALSE, + ...) + } + if (has_wlevels && has_group) { + # TODO + # - Not yet supported. + # - Need to use expand.grid to create + # all combinations of group and wlevels + group_numbers_long <- NULL + wlevels2_long <- NULL + out <- mapply(function(gn, + wv, + x, + y, + m, + fit, + est, + implied_stats, + boot_ci, + boot_out, + R, + seed, + prods, + save_boot_out, + mc_ci, + mc_out, + save_mc_out, + ci_type, + ci_out, + save_ci_out, + ...) { + cond_indirect(wvalues = wv, + x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + boot_ci = boot_ci, + boot_out = boot_out, + R = R, + seed = seed, + prods = prods, + save_boot_out = FALSE, + mc_ci = mc_ci, + mc_out = mc_out, + save_mc_out = FALSE, + ci_type = ci_type, + ci_out = ci_out, + save_ci_out = FALSE, + group = gn, + ...) + }, + gn = group_numbers_long, + wv = wlevels2_long, + x = x, + y = y, + m = m, + fit = fit, + est = est, + implied_stats = implied_stats, + boot_ci = boot_ci, + boot_out = boot_out, + R = R, + seed = seed, + prods = prods, + save_boot_out = FALSE, + mc_ci = mc_ci, + mc_out = mc_out, + save_mc_out = FALSE, + ci_type = ci_type, + ci_out = ci_out, + save_ci_out = FALSE, + ...) + } if (output_type == "data.frame") { - out1 <- cond_indirect_effects_to_df(out, wlevels = wlevels) + out1 <- cond_indirect_effects_to_df(out, + wlevels = wlevels, + group_numbers = group_numbers, + group_labels = group_labels) class(out1) <- c("cond_indirect_effects", class(out1)) attr(out1, "call") <- match.call() attr(out1, "full_output") <- out @@ -1015,6 +1319,8 @@ cond_indirect_effects <- function(wlevels, attr(out1, "x") <- x attr(out1, "y") <- y attr(out1, "m") <- m + # TODO: + # - Store the expanded combination of group and wlevels return(out1) } else { return(out) @@ -1042,18 +1348,38 @@ cond_indirect_effects <- function(wlevels, #' " #' fit <- sem(mod, data_serial_parallel, #' fixed.x = FALSE) -#' #' # All indirect paths from x to y #' paths <- all_indirect_paths(fit, #' x = "x", #' y = "y") #' paths -#' #' # Indirect effect estimates #' out <- many_indirect_effects(paths, #' fit = fit) #' out #' +#' # Multigroup models for many_indirect_effects() +#' +#' data(data_med_complicated_mg) +#' mod <- +#' " +#' m11 ~ x1 + x2 + c1 + c2 +#' m12 ~ m11 + c1 + c2 +#' m2 ~ x1 + x2 + c1 + c2 +#' y1 ~ m11 + m12 + x1 + x2 + c1 + c2 +#' y2 ~ m2 + x1 + x2 + c1 + c2 +#' " +#' fit <- sem(mod, data_med_complicated_mg, group = "group") +#' summary(fit) +#' +#' paths <- all_indirect_paths(fit, +#' x = "x1", +#' y = "y1") +#' paths +#' # Indirect effect estimates for all paths in all groups +#' out <- many_indirect_effects(paths, +#' fit = fit) +#' out #' #' @export #' @@ -1067,12 +1393,22 @@ cond_indirect_effects <- function(wlevels, many_indirect_effects <- function(paths, ...) { path_names <- names(paths) xym <- all_paths_to_df(paths) - out <- mapply(indirect_effect, - x = xym$x, - y = xym$y, - m = xym$m, - MoreArgs = list(...), - SIMPLIFY = FALSE) + if ("group_label" %in% colnames(xym)) { + out <- mapply(indirect_effect, + x = xym$x, + y = xym$y, + m = xym$m, + group = xym$group_number, + MoreArgs = list(...), + SIMPLIFY = FALSE) + } else { + out <- mapply(indirect_effect, + x = xym$x, + y = xym$y, + m = xym$m, + MoreArgs = list(...), + SIMPLIFY = FALSE) + } names(out) <- path_names class(out) <- c("indirect_list", class(out)) attr(out, "paths") <- paths @@ -1110,12 +1446,25 @@ cond_indirect_check_fit <- function(fit) { # information on the levels of moderators. #' @noRd -cond_indirect_effects_to_df <- function(x, wlevels) { - k <- nrow(wlevels) - wlevels_label <- attr(wlevels, "wlevels") - colnames(wlevels_label) <- paste0("[", colnames(wlevels_label), "]") - wlevels2 <- wlevels - colnames(wlevels2) <- paste0("(", colnames(wlevels2), ")") +cond_indirect_effects_to_df <- function(x, + wlevels = NULL, + group_numbers = NULL, + group_labels = NULL) { + has_wlevels <- !is.null(wlevels) + has_group <- !is.null(group_numbers) || !is.null(group_labels) + if (has_wlevels) { + # TOFIX + wlevels_label <- attr(wlevels, "wlevels") + colnames(wlevels_label) <- paste0("[", colnames(wlevels_label), "]") + wlevels2 <- wlevels + colnames(wlevels2) <- paste0("(", colnames(wlevels2), ")") + } + if (has_group) { + group_numbers <- sapply(x, function(x) x$group_number) + group_labels <- sapply(x, function(x) x$group_label) + gp_df <- data.frame(Group = group_labels, + Group_ID = group_numbers) + } standardized_x <- x[[1]]$standardized_x standardized_y <- x[[1]]$standardized_y if (standardized_x || standardized_y) { @@ -1167,8 +1516,13 @@ cond_indirect_effects_to_df <- function(x, wlevels) { out <- data.frame(std = indirect_std, cc, ustd = indirect, check.names = FALSE) } } - out1 <- cbind(wlevels_label, wlevels2, out) - out1 + if (has_wlevels) { + out <- cbind(wlevels_label, wlevels2, out) + } + if (has_group) { + out <- cbind(gp_df, out) + } + out } # Check the argument `wlevels` and convert it to a valid data diff --git a/R/cond_indirect_diff.R b/R/cond_indirect_diff.R index 123fb18d..1741c25f 100644 --- a/R/cond_indirect_diff.R +++ b/R/cond_indirect_diff.R @@ -158,6 +158,11 @@ cond_indirect_diff <- function(output, from = NULL, to = NULL, level = .95) { + has_wlevels <- cond_indirect_effects_has_wlevels(output) + has_groups <- cond_indirect_effects_has_groups(output) + if (has_wlevels && has_groups) { + stop("Objects with both wlevels and groups not yet supported") + } if (missing(output)) { stop("'output' is missing.") } @@ -229,9 +234,29 @@ cond_indirect_diff <- function(output, boot_diff_p <- NA boot_diff_se <- NA } - wlevels <- attr(output, "wlevels") - wlevels_from <- wlevels[from, , drop = FALSE] - wlevels_to <- wlevels[to, , drop = FALSE] + if (has_wlevels && !has_groups) { + wlevels <- attr(output, "wlevels") + wlevels_from <- wlevels[from, , drop = FALSE] + wlevels_to <- wlevels[to, , drop = FALSE] + final_from <- wlevels_from + final_to <- wlevels_to + } + if (!has_wlevels && has_groups) { + group_labels <- output$Group + group_numbers <- output$Group_ID + group_from <- paste0(group_labels[from], " [", + group_numbers[from], "]") + group_to <- paste0(group_labels[to], " [", + group_numbers[to], "]") + names(group_from) <- "Group" + names(group_to) <- "Group" + final_from <- group_from + final_to <- group_to + } + if (has_wlevels && has_groups) { + # TODO: + # - Add support for objects with both wlevels and groups. + } out_diff_ci <- c(NA, NA) out_diff_se <- NA if (has_mc) out_diff_ci <- mc_diff_ci @@ -243,11 +268,13 @@ cond_indirect_diff <- function(output, pvalue = boot_diff_p, se = out_diff_se, level = level, - from = wlevels_from, - to = wlevels_to, + from = final_from, + to = final_to, output = output[c(to, from), ], boot_diff = boot_diff, - mc_diff = mc_diff) + mc_diff = mc_diff, + has_groups = has_groups, + has_wlevels = has_wlevels) class(out) <- c("cond_indirect_diff", class(out)) out } @@ -328,6 +355,8 @@ print.cond_indirect_diff <- function(x, se = FALSE, ...) { full_output_attr <- attr(x$output, "full_output")[[1]] + has_groups <- is.numeric(full_output_attr$group_number) + has_wlevels <- !is.null(full_output_attr$wvalues) print(x$output, digits = digits, annotation = FALSE, ...) x_type <- x$type if (!is.null(x_type)) { @@ -368,6 +397,15 @@ print.cond_indirect_diff <- function(x, } else { cat("\nLevels: \n") print(tofrom, quote = FALSE) + if (!has_wlevels && has_groups) { + tmp <- paste0("Levels compared: ", + xto0, + " - ", + xfrom0) + } else { + tmp <- "Levels compared: Row 1 - Row 2" + } + cat("\n", tmp, "\n", sep = "") } index_df <- data.frame(x = full_output_attr$x, y = full_output_attr$y, diff --git a/R/cond_indirect_effects_math.R b/R/cond_indirect_effects_math.R index 49714760..bfdde11e 100644 --- a/R/cond_indirect_effects_math.R +++ b/R/cond_indirect_effects_math.R @@ -10,8 +10,7 @@ #' and `-` operator are supported. These #' operators can be used to estimate and #' test a function of effects between -#' the same pair of variables but along -#' different paths. +#' the same pair of variables. #' #' For example, they can be used to #' compute and test the total effects @@ -28,8 +27,7 @@ #' the same variable, #' #' 2. the two paths do not end at the -#' same variable, (c) a path appears in -#' both objects, +#' same variable, #' #' 3. moderators are involved but they #' are not set to the same values in @@ -42,6 +40,28 @@ #' estimates stored in #' `mc_out`, if any, are not identical. #' +#' ## Multigroup Models +#' +#' Since Version 0.1.14.2, support for +#' multigroup models has been added for models +#' fitted by `lavaan`. Both bootstrapping +#' and Monte Carlo confidence intervals +#' are supported. These operators can +#' be used to compute and test the +#' difference of an indirect effect +#' between two groups. This can also +#' be used to compute and test the +#' difference between a function of +#' effects between groups, for example, +#' the total indirect effects between +#' two groups. +#' +#' The operators are flexible and allow +#' users to do many possible computations. +#' Therefore, users need to make sure +#' that the function of effects is +#' meaningful. +#' #' @return An 'indirect'-class object #' with a list of effects stored. See #' [indirect_effect()] on details for @@ -91,6 +111,35 @@ NULL #' out123 #' coef(out1) + coef(out2) + coef(out3) #' +#' # Multigroup model with indirect effects +#' +#' dat <- data_med_mg +#' mod <- +#' " +#' m ~ x + c1 + c2 +#' y ~ m + x + c1 + c2 +#' " +#' fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, +#' group = "group") +#' +#' # If a model has more than one group, +#' # the argument 'group' must be set. +#' ind1 <- indirect_effect(x = "x", +#' y = "y", +#' m = "m", +#' fit = fit, +#' group = "Group A") +#' ind1 +#' ind2 <- indirect_effect(x = "x", +#' y = "y", +#' m = "m", +#' fit = fit, +#' group = 2) +#' ind2 +#' +#' # Compute the difference in indirect effects between groups +#' ind2 - ind1 +#' #' @export `+.indirect` <- function(e1, e2) { plusminus(e1, e2, op = "+") @@ -107,24 +156,43 @@ NULL plusminus <- function(e1, e2, op = c("+", "-")) { op <- match.arg(op, c("+", "-")) check_xy(e1, e2) + # group_number and group_label can be vectors + group_number_1 <- e1$group_number + group_number_2 <- e2$group_number + group_label_1 <- e1$group_label + group_label_2 <- e2$group_label + if (is.numeric(group_number_1) && is.numeric(group_number_2)) { + has_group <- TRUE + group_labels <- c(group_label_1, group_label_2) + } else { + has_group <- FALSE + } cp1 <- if (is.list(e1$components)) e1$components else list(e1$components) cp2 <- if (is.list(e2$components)) e2$components else list(e2$components) cp0 <- c(cp1, cp2) + if (has_group) names(cp0) <- group_labels cpc1 <- if (is.list(e1$components_conditional)) e1$components_conditional else list(e1$components_conditional) cpc2 <- if (is.list(e2$components_conditional)) e2$components_conditional else list(e2$components_conditional) cpc0 <- c(cpc1, cpc2) + if (has_group) names(cpc0) <- group_labels m1 <- if (is.list(e1$m)) e1$m else list(e1$m) m2 <- if (is.list(e2$m)) e2$m else list(e2$m) m0 <- c(m1, m2) + if (has_group) names(m0) <- group_labels cv1 <- if (is.list(e1$computation_values)) e1$computation_values else list(e1$computation_values) cv2 <- if (is.list(e2$computation_values)) e2$computation_values else list(e2$computation_values) cv0 <- c(cv1, cv2) + if (has_group) names(cv0) <- group_labels cs1 <- if (is.list(e1$computation_symbols)) e1$computation_symbols else list(e1$computation_symbols) cs2 <- if (is.list(e2$computation_symbols)) e2$computation_symbols else list(e2$computation_symbols) cs0 <- c(cs1, cs2) + if (has_group) names(cs0) <- group_labels ca1 <- if (is.list(e1$call)) e1$call else list(e1$call) ca2 <- if (is.list(e2$call)) e2$call else list(e2$call) ca0 <- c(ca1, ca2) + if (has_group) names(ca0) <- group_labels + gnumber0 <- c(group_number_1, group_number_2) + glabel0 <- c(group_label_1, group_label_2) est0 <- switch(op, "+" = e1$indirect + e2$indirect, "-" = e1$indirect - e2$indirect) @@ -182,6 +250,10 @@ plusminus <- function(e1, e2, op = c("+", "-")) { paste(eval(e1$m), collapse = "->"), "->", e1$y) } + if (has_group && (length(group_label_1) == 1)) { + op1 <- paste0(group_label_1, "[", + group_number_1, "]: ", op1) + } } if (is.null(op2)) { if (is.null(e2$m)) { @@ -191,6 +263,10 @@ plusminus <- function(e1, e2, op = c("+", "-")) { paste(eval(e2$m), collapse = "->"), "->", e2$y) } + if (has_group && (length(group_label_2) == 1)) { + op2 <- paste0(group_label_2, "[", + group_number_2, "]: ", op2) + } } op0 <- paste0("(", op1, ")", "\n", op, "(", op2, ")") @@ -243,7 +319,9 @@ plusminus <- function(e1, e2, op = c("+", "-")) { mc_scale_y = e1$mc_scale_y, level = level0, boot_out = e1$boot_out, - mc_out = e1$mc_out + mc_out = e1$mc_out, + group_number = gnumber0, + group_label = glabel0 ) class(out) <- c("indirect", class(out)) out @@ -269,6 +347,17 @@ check_xy <- function(e1, e2) { scy1 <- e1$scale_y scx2 <- e2$scale_x scy2 <- e2$scale_y + group1 <- e1$group_number + group2 <- e2$group_number + if (is.numeric(group1) && is.numeric(group2)) { + has_group <- TRUE + } else { + has_group <- FALSE + } + if ((is.null(group1) && is.numeric(group2)) || + (is.null(group2) && is.numeric(group1))) { + stop("The objects do not agree in the number of groups.") + } if (!identical(x1, x2)) { stop("The objects to be added do not have the same 'x'.") } @@ -282,9 +371,12 @@ check_xy <- function(e1, e2) { m2 <- list(m2) } m1m2_chk <- intersect(m1, m2) - if (length(m1m2_chk) != 0) { - stop("The objects have one or more paths in common.") - } + # Disable this test. + # - The two effects can be two conditional effects. + # - The two effects can be from two different groups. + # if (length(m1m2_chk) != 0) { + # stop("The objects have one or more paths in common.") + # } if (!identical(stdx1, stdx2)) { stop("x is standardized in one object but not in the other.") } @@ -303,11 +395,13 @@ check_xy <- function(e1, e2) { } } } - if (!identical(scx1, scx2)) { - stop("x is not scaled by the same factor (SD) in the two objects.") - } - if (!identical(scy1, scy2)) { - stop("y is not scaled by the same factor (SD) in the two objects.") + if (!has_group) { + if (!identical(scx1, scx2)) { + stop("x is not scaled by the same factor (SD) in the two objects.") + } + if (!identical(scy1, scy2)) { + stop("y is not scaled by the same factor (SD) in the two objects.") + } } if (!identical(e1$level, e2$level)) { stop("The two objects do not have the same level for confidence interval.") diff --git a/R/cond_indirect_effects_subset.R b/R/cond_indirect_effects_subset.R index 3b107e8f..2cac78f9 100644 --- a/R/cond_indirect_effects_subset.R +++ b/R/cond_indirect_effects_subset.R @@ -71,6 +71,10 @@ NULL `[.cond_indirect_effects` <- function(x, i, j, drop = if (missing(i)) TRUE else length(j) == 1) { + # TODO: + # - Support object with both wlevels and groups + has_wlevels <- cond_indirect_effects_has_wlevels(x) + has_groups <- cond_indirect_effects_has_groups(x) out <- NextMethod() if (is.null(dim(out))) { if (!missing(j)) { @@ -87,9 +91,11 @@ NULL } else { fo <- attr(x, "full_output")[i] attr(out, "full_output") <- fo - wl <- attr(x, "wlevels")[i, , drop = FALSE] - attr(wl, "wlevels") <- attr(wl, "wlevels")[i, , drop = FALSE] - attr(out, "wlevels") <- wl + if (has_wlevels) { + wl <- attr(x, "wlevels")[i, , drop = FALSE] + attr(wl, "wlevels") <- attr(wl, "wlevels")[i, , drop = FALSE] + attr(out, "wlevels") <- wl + } return(out) } } \ No newline at end of file diff --git a/R/confint_cond_indirect_effects.R b/R/confint_cond_indirect_effects.R index def9677a..0e2b3485 100644 --- a/R/confint_cond_indirect_effects.R +++ b/R/confint_cond_indirect_effects.R @@ -93,14 +93,16 @@ confint.cond_indirect_effects <- function(object, parm, level = .95, ...) { # warning("Bootstrapping interval not in the object.") # out0 <- c(NA, NA) # } + has_wlevels <- cond_indirect_effects_has_wlevels(object) + has_groups <- cond_indirect_effects_has_groups(object) out0 <- as.data.frame(object) full_output <- attr(object, "full_output") has_ci <- FALSE - if (is.null(full_output[[1]]$boot_ci)) { + if (!is.null(full_output[[1]]$boot_ci)) { has_ci <- TRUE ci_type <- "boot" } - if (is.null(full_output[[1]]$mc_ci)) { + if (!is.null(full_output[[1]]$mc_ci)) { has_ci <- TRUE ci_type <- "mc" } @@ -117,8 +119,18 @@ confint.cond_indirect_effects <- function(object, parm, level = .95, ...) { trim = TRUE, scientific = FALSE, digits = 2), "%") - wlevels <- attr(object, "wlevels") colnames(out) <- cnames - rownames(out) <- rownames(wlevels) + if (has_wlevels && !has_groups) { + wlevels <- attr(object, "wlevels") + rownames(out) <- rownames(wlevels) + } + if (!has_wlevels && has_groups) { + tmp <- paste0(object$Group, " [", object$Group_ID, "]") + rownames(out) <- tmp + } + if (has_wlevels && has_groups) { + # TODO: + # - Support for having both wlevels and groups + } out } diff --git a/R/dat_2_med_mg.R b/R/dat_2_med_mg.R new file mode 100644 index 00000000..df121ad1 --- /dev/null +++ b/R/dat_2_med_mg.R @@ -0,0 +1,32 @@ +#' @title Sample Dataset: Simple +#' Mediation With Two Groups +#' +#' @description A simple mediation +#' model with two groups. +#' +#' @format A data frame with 100 rows +#' and 5 variables: +#' \describe{ +#' \item{x}{Predictor. Numeric.} +#' \item{m}{Mediator. Numeric.} +#' \item{y}{Outcome variable. Numeric.} +#' \item{c1}{Control variable. Numeric.} +#' \item{c2}{Control variable. Numeric.} +#' \item{group}{Group variable. Character. "Group A" or "Group B"} +#' } +#' +#' @examples +#' library(lavaan) +#' data(data_med_mg) +#' mod <- +#' " +#' m ~ c(a1, a2) * x + c1 + c2 +#' y ~ c(b1, b2) * m + x + c1 + c2 +#' a1b1 := a1 * b1 +#' a2b2 := a2 * b2 +#' abdiff := a2b2 - a1b1 +#' " +#' fit <- sem(mod, data_med_mg, fixed.x = FALSE, +#' group = "group") +#' parameterEstimates(fit) +"data_med_mg" diff --git a/R/dat_5_mg.R b/R/dat_5_mg.R new file mode 100644 index 00000000..6057029b --- /dev/null +++ b/R/dat_5_mg.R @@ -0,0 +1,37 @@ +#' @title Sample Dataset: A Complicated +#' Mediation Model With Two Groups +#' +#' @description A mediation model with +#' two predictors, two pathways, and +#' two groups. +#' +#' @format A data frame with 300 rows +#' and 5 variables: +#' \describe{ +#' \item{x1}{Predictor 1. Numeric.} +#' \item{x2}{Predictor 2. Numeric.} +#' \item{m11}{Mediator 1 in Path 1. Numeric.} +#' \item{m12}{Mediator 2 in Path 1. Numeric.} +#' \item{m2}{Mediator in Path 2. Numeric.} +#' \item{y1}{Outcome variable 1. Numeric.} +#' \item{y2}{Outcome variable 2. Numeric.} +#' \item{c1}{Control variable. Numeric.} +#' \item{c2}{Control variable. Numeric.} +#' \item{group}{Group variable. Character. 'Group A' or 'Group B'} +#' } +#' +#' @examples +#' library(lavaan) +#' data(data_med_complicated_mg) +#' dat <- data_med_complicated_mg +#' mod <- +#' " +#' m11 ~ x1 + x2 + c1 + c2 +#' m12 ~ m11 + c1 + c2 +#' m2 ~ x1 + x2 + c1 + c2 +#' y1 ~ m11 + m12 + x1 + x2 + c1 + c2 +#' y2 ~ m2 + x1 + x2 + c1 + c2 +#' " +#' fit <- sem(mod, dat, group = "group") +#' summary(fit) +"data_med_complicated_mg" diff --git a/R/do_boot.R b/R/do_boot.R index 193e2466..02370b1c 100644 --- a/R/do_boot.R +++ b/R/do_boot.R @@ -34,6 +34,15 @@ #' [lm2boot_out()], [fit2boot_out()], or #' [fit2boot_out_do_boot()]. #' +#' ## Multigroup Models +#' +#' Since Version 0.1.14.2, support for +#' multigroup models has been added for models +#' fitted by `lavaan`. The implementation +#' of bootstrapping is identical to +#' that used by `lavaan`, with resampling +#' done within each group. +#' #' @return A `boot_out`-class object #' that can be used for the `boot_out` #' argument of diff --git a/R/do_mc.R b/R/do_mc.R index 136fabe6..655603c8 100644 --- a/R/do_mc.R +++ b/R/do_mc.R @@ -38,6 +38,11 @@ #' estimates is used in all subsequent #' analysis. #' +#' ## Multigroup Models +#' +#' Since Version 0.1.14.2, support for +#' multigroup models has been added for models +#' fitted by `lavaan`. #' #' @return A `mc_out`-class object #' that can be used for the `mc_out` diff --git a/R/find_product.R b/R/find_product.R index 5c823513..eee59a5b 100644 --- a/R/find_product.R +++ b/R/find_product.R @@ -46,6 +46,13 @@ #'@noRd find_product <- function(data, target) { + if (is.list(data) && !is.data.frame(data)) { + ngroups <- length(data) + # Aasume all groups have the same variables + data <- do.call(rbind, data) + } else { + ngroups <- 1 + } a_col <- data[, target] out <- c(NA, NA) q <- 0 @@ -69,6 +76,12 @@ find_product <- function(data, target) { #'@noRd find_all_products <- function(data, expand = TRUE) { + if (is.list(data) && !is.data.frame(data)) { + ngroups <- length(data) + data <- do.call(rbind, data) + } else { + ngroups <- 1 + } out <- sapply(colnames(data), find_product, data = data, USE.NAMES = TRUE, diff --git a/R/get_b.R b/R/get_b.R index e3b0620b..234b2914 100644 --- a/R/get_b.R +++ b/R/get_b.R @@ -30,15 +30,28 @@ get_b <- function(x, y, fit, - est = NULL) { + est = NULL, + group_number = NULL) { if (is.null(est)) { est <- lav_est(fit, se = FALSE, ci = FALSE) } + if (!is.null(group_number)) { + est <- est[est$group == group_number, ] + } i <- (est$lhs == y) & (est$op == "~") & (est$rhs == x) if (isTRUE(any(i))) { - return(est[i, "est"]) + out <- est[i, "est"] + if (length(out) > 1) { + # Multigroup model + if (!is.null(est$group)) { + names(out) <- est$group[i] + } else { + names(out) <- seq_along(out) + } + } + return(out) } else { return(NA) } diff --git a/R/get_prod.R b/R/get_prod.R index d6fd9be9..c2ab3d5a 100644 --- a/R/get_prod.R +++ b/R/get_prod.R @@ -131,10 +131,19 @@ get_prod <- function(x, expand = expand) all_prods_names <- names(all_prods) } + if (!is.null(est$group) || + isTRUE(suppressWarnings(max(est$group) > 1))) { + ngroups <- max(est$group) + } else { + ngroups <- 1 + } + if (ngroups > 1) { + + } i_rhs <- (est$lhs == y) & (est$op == "~") if (isTRUE(any(i_rhs))) { - y_rhs <- est[i_rhs, "rhs"] + y_rhs <- unique(est[i_rhs, "rhs"]) } else { return(NA) } @@ -185,7 +194,15 @@ get_prod <- function(x, } b_prod <- sapply(prod_x, function(x) get_b(x = x, y = y, - est = est)) + est = est), + simplify = FALSE, + USE.NAMES = TRUE) + # Handle multigroup models + if (ngroups == 1) { + b_prod <- unlist(b_prod) + } else { + b_prod <- b_prod + } out <- list(prod = prod_x, b = b_prod, w = w, diff --git a/R/indirect.R b/R/indirect.R index 55bf60b3..b828bdf2 100644 --- a/R/indirect.R +++ b/R/indirect.R @@ -124,6 +124,16 @@ #' latent variables and observed #' variables. Default is `TRUE`. #' +#' @param group Either the group number +#' as appeared in the [summary()] +#' or [lavaan::parameterEstimates()] +#' output of an `lavaan`-class object, +#' or the group label as used in +#' the `lavaan`-class object. +#' Used only when the number of +#' groups is greater than one. Default +#' is NULL. +#' #' @seealso [indirect_effect()], #' [cond_indirect_effects()], and #' [cond_indirect()], the high level @@ -183,10 +193,37 @@ indirect_i <- function(x, data = NULL, expand = TRUE, warn = TRUE, - allow_mixing_lav_and_obs = TRUE) { + allow_mixing_lav_and_obs = TRUE, + group = NULL) { if (is.null(est)) { est <- lav_est(fit) } + ngroups <- 1 + if (!is.null(est$group)) { + if (max(est$group) > 1) { + ngroups <- max(est$group) + } + } + if ((ngroups > 1) && + !is.numeric(group) && + !is.character(group)) { + stop("The model has more than one group but group is not set.") + } + if (ngroups > 1) { + group_labels_all <- lavaan::lavTech(fit, + "group.label") + if (is.numeric(group)) { + group_label <- group_labels_all[group] + group_number <- group + } else { + group_number <- match(group, group_labels_all) + group_label <- group + } + } else { + group_labels_all <- NULL + group_number <- NULL + group_label <- NULL + } chkpath <- check_path(x = x, y = y, m = m, fit = fit, est = est) if (!chkpath) { msg <- paste0("No path from ", sQuote(x), " to ", sQuote(y), @@ -206,7 +243,8 @@ indirect_i <- function(x, bs <- mapply(get_b, x = xs, y = ys, - MoreArgs = list(est = est)) + MoreArgs = list(est = est, + group_number = group_number)) bs_org <- bs names(bs_org) <- bs_names chk_lv <- unique(c(xs, ys)) %in% check_lv_in_est(est) @@ -343,7 +381,9 @@ indirect_i <- function(x, } sum(b_i * wvalues_i) } - b_cond <- sapply(prods, tmpfct) + prods_tmp <- prods_group_i(prods, + group_number = group_number) + b_cond <- sapply(prods_tmp, tmpfct) bs <- bs + b_cond } else { b_cond <- rep(NA, length(bs)) @@ -353,7 +393,8 @@ indirect_i <- function(x, MoreArgs = list(digits = computation_digits, y = y, wvalues = wvalues, - warn = warn), + warn = warn, + group_number = group_number), USE.NAMES = TRUE, SIMPLIFY = FALSE) b_all_str0 <- paste0("(", b_cond_str, ")", collapse = "*") @@ -366,6 +407,10 @@ indirect_i <- function(x, if (is.null(implied_stats)) { implied_stats <- lav_implied_all(fit) } + if (!is.null(group_number)) { + implied_stats <- implied_stats_group_i(implied_stats, + group_number = group_number) + } if (standardized_x) { scale_x <- sqrt(diag(implied_stats$cov)[x]) b_all_str0 <- paste0(b_all_str0, "*(", @@ -398,7 +443,9 @@ indirect_i <- function(x, y = y, m = m, computation_values = b_all_str0, - computation_symbols = b_all_str1) + computation_symbols = b_all_str1, + group_number = group_number, + group_label = group_label) class(out) <- "indirect" return(out) } @@ -406,7 +453,8 @@ indirect_i <- function(x, #' @noRd gen_computation <- function(xi, yi, yiname, digits = 3, y, wvalues = NULL, - warn = TRUE) { + warn = TRUE, + group_number = NULL) { yiname_old <- yiname yiname <- paste0("b.", yiname) if (all(is.na(xi)) || is.null(xi$prod)) { @@ -414,6 +462,11 @@ gen_computation <- function(xi, yi, yiname, digits = 3, y, wvalues = NULL, names(out) <- yiname return(out) } + if (is.numeric(group_number)) { + tmp <- sapply(xi$b, function(xx) xx[group_number]) + names(tmp) <- names(xi$b) + xi$b <- tmp + } b_i <- xi$b b_i0 <- paste0("b.", names(b_i)) w_i <- xi$w @@ -471,15 +524,15 @@ gen_computation <- function(xi, yi, yiname, digits = 3, y, wvalues = NULL, out1 <- paste0(y0, " + ", paste0("(", b_i0, ")*(", w_i2, ")", collapse = " + ")) - - out2 <- paste0("(", formatC(yi, digits = digits, format = "f"), + out2 <- tryCatch(paste0("(", formatC(yi, digits = digits, format = "f"), ") + ", paste0("(", formatC(b_i, digits = digits, format = "f"), ")*(", wvalues_i0, ")", - collapse = " + ")) + collapse = " + ")), error = function(e) e) + if (inherits(out2, "error")) browser() names(out2) <- out1 out2 } @@ -499,14 +552,41 @@ update_prods <- function(prods, est) { if (all(is.na(prods_i))) { return(pout_i) } else { - est_i <- est[(est$lhs == prods_i$y) & - (est$op == "~") & - (est$rhs %in% prods_i$prod), "est"] - pout_i$b <- est_i - names(pout_i$b) <- prods_i$prod + if (is.list(pout_i$b)) { + for (pp in pout_i$prod) { + est_i <- est[(est$lhs == prods_i$y) & + (est$op == "~") & + (est$rhs %in% pp), "est"] + pout_i$b[[pp]][] <- est_i + } + } else { + est_i <- est[(est$lhs == prods_i$y) & + (est$op == "~") & + (est$rhs %in% prods_i$prod), "est"] + pout_i$b <- est_i + names(pout_i$b) <- prods_i$prod + } return(pout_i) } } pout <- sapply(prods, tmpfct, simplify = FALSE) pout + } + +#' @noRd + +prods_group_i <- function(prods, + group_number = NULL) { + if (!is.numeric(group_number)) { + return(prods) + } + pout <- prods + for (i in seq_along(pout)) { + if (!identical(pout[[i]], NA)) { + tmp <- sapply(pout[[i]]$b, function(xx) xx[group_number]) + names(tmp) <- names(pout[[i]]$b) + pout[[i]]$b <- tmp + } + } + pout } \ No newline at end of file diff --git a/R/lav_data_used.R b/R/lav_data_used.R index d9ce127e..537e9a8d 100644 --- a/R/lav_data_used.R +++ b/R/lav_data_used.R @@ -23,7 +23,10 @@ #' lav_data_used <- function(fit, - drop_colon = TRUE) { + drop_colon = TRUE, + drop_list_single_group = TRUE) { + # Return a named list of N matrices if ngroups > 1 + # TODOs: lav_data_used_lavaan_mi() type <- NA if (inherits(fit, "lavaan")) { type <- "lavaan" @@ -36,7 +39,8 @@ lav_data_used <- function(fit, } out <- switch(type, lavaan = lav_data_used_lavaan(fit = fit, - drop_colon = drop_colon), + drop_colon = drop_colon, + drop_list_single_group = drop_list_single_group), lavaan.mi = lav_data_used_lavaan_mi(fit = fit, drop_colon = drop_colon)) out @@ -45,28 +49,40 @@ lav_data_used <- function(fit, #' @noRd lav_data_used_lavaan <- function(fit, - drop_colon = TRUE) { - dat <- lavaan::lavInspect(fit, "data") - i_excluded <- lavaan::lavInspect(fit, "empty.idx") - vnames <- colnames(dat) + drop_colon = TRUE, + drop_list_single_group = TRUE) { + # Return a named list of N matrices if ngroups > 1 + dat <- lavaan::lavInspect(fit, "data", + drop.list.single.group = FALSE) + i_excluded <- lavaan::lavInspect(fit, "empty.idx", + drop.list.single.group = FALSE) + ngp <- lavaan::lavTech(fit, "ngroups") if (drop_colon) { - vraw <- vnames[!grepl(":", colnames(dat))] - vcolon <- apply(expand.grid(vraw, vraw), 1, paste0, collapse = ":") - vkeep <- vnames[!(vnames %in% vcolon)] - dat <- dat[, vkeep] + for (k in seq_len(ngp)) { + dat_i <- dat[[k]] + vnames <- colnames(dat_i) + vraw <- vnames[!grepl(":", colnames(dat_i))] + vcolon <- apply(expand.grid(vraw, vraw), 1, paste0, collapse = ":") + vkeep <- vnames[!(vnames %in% vcolon)] + dat[[k]] <- dat_i[, vkeep] + } } - if (length(i_excluded) > 0) { - return(dat[-i_excluded, ]) - } else { - return(dat) + for (k in seq_len(ngp)) { + if (length(i_excluded[[k]]) > 0) { + dat[[k]] <- dat[[k]][-i_excluded[[k]], ] + } + } + if (drop_list_single_group && (ngp == 1)) { + dat <- dat[[1]] } + return(dat) } #' @noRd lav_data_used_lavaan_mi <- function(fit, drop_colon = TRUE) { - + # TODOs: Return a named list of N matrices if ngroups > 1 dat_list <- fit@DataList ovnames <- lavaan::lavNames(fit, "ov") fit_org <- lavaan_from_lavaam_mi(fit, data = FALSE) diff --git a/R/lavaan2lm_list.R b/R/lavaan2lm_list.R index 8bea1a89..bf918a49 100644 --- a/R/lavaan2lm_list.R +++ b/R/lavaan2lm_list.R @@ -46,9 +46,43 @@ #' #' @export - lm_from_lavaan_list <- function(fit) { + ngroups <- lavaan::lavTech(fit, "ngroups") + if (ngroups > 1) { + group_numbers <- seq_len(ngroups) + group_labels <- lavaan::lavTech(fit, "group.label") + out <- lapply(group_numbers, + function(x) { + lm_from_lavaan_list_i(fit = fit, + group_number = x) + }) + names(out) <- group_labels + class(out) <- c("lm_from_lavaan_list", class(out)) + } else { + out <- lm_from_lavaan_list_i(fit = fit) + } + out + } + +#' @noRd + +lm_from_lavaan_list_i <- function(fit, + group_number = NULL) { ptable <- lav_ptable(fit) + if ("group" %in% colnames(ptable)) { + tmp <- unique(ptable$group) + if (length(tmp) > 1) { + if (is.null(group_number)) { + stop("The model has more than one groups but group_number not set.") + } + if (!(group_number %in% tmp)) { + stop("group_number is", group_number, + "but group numbers in the model are", + paste0(tmp, collapse = ", ")) + } + ptable <- ptable[ptable$group == group_number, ] + } + } # Get all dvs (ov.nox, lv.ox) dvs <- lavaan_get_dvs(ptable) dat <- lav_data_used(fit) diff --git a/R/lavaan_helpers.R b/R/lavaan_helpers.R index e312458e..0be1b6fd 100644 --- a/R/lavaan_helpers.R +++ b/R/lavaan_helpers.R @@ -1,6 +1,7 @@ #' @noRd -lav_implied_all <- function(fit) { +lav_implied_all <- function(fit, + group_number = NULL) { type <- NA if (inherits(fit, "lavaan")) { type <- "lavaan" @@ -12,14 +13,17 @@ lav_implied_all <- function(fit) { stop("Object is not of a supported type.") } out <- switch(type, - lavaan = lav_implied_all_lavaan(fit), - lavaan.mi = lav_implied_all_lavaan_mi(fit)) + lavaan = lav_implied_all_lavaan(fit, + group_number = group_number), + lavaan.mi = lav_implied_all_lavaan_mi(fit, + group_number = group_number)) out } #' @noRd -lav_implied_all_lavaan <- function(fit) { +lav_implied_all_lavaan <- function(fit, + group_number = NULL) { ovnames <- lavaan::lavNames(fit, "ov") lvnames <- lavaan::lavNames(fit, "lv") allnames <- c(ovnames, lvnames) @@ -42,7 +46,8 @@ lav_implied_all_lavaan <- function(fit) { #' @noRd -lav_implied_all_lavaan_mi <- function(fit) { +lav_implied_all_lavaan_mi <- function(fit, + group_number = NULL) { est0 <- methods::getMethod("coef", signature = "lavaan.mi", where = asNamespace("semTools"))(fit) @@ -192,4 +197,128 @@ lav_ptable_lavaan_mi <- function(fit, ...) { out } +#' @noRd + +implied_stats_group_i <- function(object, + group_number = group_number) { + if (is.null(group_number)) { + return(object) + } + out <- lapply(object, function(x) { + if (is.list(x)) { + x[[group_number]] + } else { + x + } + }) + out + } + +#' @noRd + +group_labels_and_numbers <- function(groups = NULL, + fit) { + if (!inherits(fit, "lavaan")) { + stop("The argument 'fit' must be a lavaan object.") + } + group_labels_all <- lavaan::lavInspect(fit, "group.label") + group_numbers_all <- seq_len(lavaan::lavInspect(fit, "ngroups")) + if (is.null(groups)) { + out <- list(label = group_labels_all, + number = group_numbers_all) + return(out) + } + if (is.numeric(groups)) { + if (!all(groups %in% group_numbers_all)) { + stop("Group numbers not among the numbers in the fit object.") + } + group_numbers <- groups + group_labels <- group_labels_all[group_numbers] + } + if (is.character(groups)) { + if (!all(groups %in% group_labels_all)) { + stop("Group label(s) not among the labels in the fit object.") + } + group_labels <- groups + group_numbers <- match(groups, group_labels_all) + } + list(label = group_labels, + number = group_numbers) + } +#' @noRd + +group_labels_and_numbers_cond <- function(object, + group_label_name = "Group", + group_number_name = "Group_ID") { + if (!inherits(object, "cond_indirect_effects")) { + stop("Object must be a cond_indirect_effects-class object.") + } + group_labels <- unique(object[, group_label_name, drop = TRUE]) + group_numbers <- unique(object[, group_number_name, drop = TRUE]) + list(label = group_labels, + number = group_numbers) + } + +#' @noRd +# Check if a cond_indirect_effects-class object has wlevels. + +cond_indirect_effects_has_wlevels <- function(object) { + if (!is.null(attr(object, "wlevels"))) { + return(TRUE) + } else { + return(FALSE) + } + } + +#' @noRd +# Check if a cond_indirect_effects-class object has groups. + +cond_indirect_effects_has_groups <- function(object) { + if (isTRUE("group" %in% tolower(colnames(object)))) { + return(TRUE) + } else { + return(FALSE) + } + } + +#' @noRd +# Check if a cond_indirect_effects-class object has groups. + +indirect_list_has_groups <- function(object) { + if (!is.list(object)) { + stop("Object must be a list or an indirect_list.") + } + tmp <- sapply(object, indirect_has_groups) + if (!(!all(tmp) || all(tmp))) { + stop("Some effects are group-specific but some are not.") + } + if (all(tmp)) { + return(TRUE) + } else { + return(FALSE) + } + } + +#' @noRd +# Check if an indirect-class object has groups. + +indirect_has_groups <- function(object) { + if (isTRUE(is.numeric(object$group_number))) { + return(TRUE) + } else { + return(FALSE) + } + } + +#' @noRd + +group_labels_and_numbers_list <- function(object) { + if (!is.list(object)) { + stop("Object must be a list or an indirect_list.") + } + group_labels <- unique(sapply(object, function(xx) xx$group_label)) + group_numbers <- unique(sapply(object, function(xx) xx$group_number)) + list(label = group_labels, + number = group_numbers) + } diff --git a/R/mod_levels.R b/R/mod_levels.R index 5e717c40..a372bfa9 100644 --- a/R/mod_levels.R +++ b/R/mod_levels.R @@ -227,6 +227,11 @@ mod_levels <- function(w, reference_group_label = NULL, descending = TRUE) { fit_type <- cond_indirect_check_fit(fit) + if (fit_type == "lavaan") { + if (lavaan::lavTech(fit, "ngroups") > 1) { + stop("Multigroup models not yet supported.") + } + } w_type <- match.arg(w_type) if (w_type == "auto") { w_type <- find_w_type(w, fit) diff --git a/R/plotmod.R b/R/plotmod.R index e9b1ba1c..fa7aa1af 100644 --- a/R/plotmod.R +++ b/R/plotmod.R @@ -12,13 +12,27 @@ #' #' It plots the conditional effect from #' `x` to `y` in a model for different -#' levels of the moderators. +#' levels of the moderators. For +#' multigroup models, the group will +#' be the 'moderator' and one line is +#' drawn for each group. #' #' It does not support conditional #' indirect effects. If there is one or #' more mediators in `x`, it will raise #' an error. #' +#' ## Multigroup Models +#' +#' Since Version 0.1.14.2, support for +#' multigroup models has been added for models +#' fitted by `lavaan`. If the effect +#' for each group is drawn, the +#' `graph_type` is automatically switched +#' to `"bumble"` and the means and SDs +#' in each group will be used to determine +#' the locations of the points. +#' #' @return A [ggplot2] graph. Plotted if #' not assigned to a name. It can be #' further modified like a usual @@ -116,9 +130,11 @@ #' @param graph_type If `"default"`, the #' typical line-graph with equal #' end-points will be plotted. If -#' `"tubmle"`, then the tumble graph +#' `"tumble"`, then the tumble graph #' proposed by Bodner (2016) will be -#' plotted. Default is `"default"`. +#' plotted. Default is `"default"` +#' for single-group models, and +#' `"tumble"` for multigroup models. #' #' @param ... Additional arguments. #' Ignored. @@ -174,6 +190,26 @@ #' plot(out_2) #' plot(out_2, graph_type = "tumble") #' +#' # Multigroup models +#' +#' dat <- data_med_mg +#' mod <- +#' " +#' m ~ x + c1 + c2 +#' y ~ m + x + c1 + c2 +#' " +#' fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, +#' group = "group") +#' +#' # For a multigroup model, group will be used as +#' # a moderator +#' out <- cond_indirect_effects(x = "m", +#' y = "y", +#' fit = fit) +#' out +#' plot(out) +#' +#' #' #' @export @@ -194,11 +230,32 @@ plot.cond_indirect_effects <- function( graph_type = c("default", "tumble"), ... ) { + has_groups <- cond_indirect_effects_has_groups(x) + has_wlevels <- cond_indirect_effects_has_wlevels(x) + if (has_wlevels && has_groups) { + stop("Objects with both wlevels and groups not yet supported") + } output <- x + if (has_groups) { + tmp <- group_labels_and_numbers_cond(output) + group_labels <- tmp$label + group_numbers <- tmp$number + w_label <- ifelse(w_label == "Moderator(s)", + "Group", + w_label) + } else { + group_labels <- NULL + group_numbers <- NULL + } fit <- attr(output, "fit") fit_type <- cond_indirect_check_fit(fit) x_method <- match.arg(x_method) graph_type <- match.arg(graph_type) + if (has_groups && graph_type == "default") { + # warning("Only tumble graph is supported for multiple group models. ", + # "Changed graph_type to 'tumble'.") + graph_type <- "tumble" + } full_output <- attr(output, "full_output") full_output_1 <- full_output[[1]] x <- full_output_1$x @@ -207,12 +264,14 @@ plot.cond_indirect_effects <- function( if (!is.null(m)) { stop("The plot method does not support indirect effects.") } - wlevels <- attr(output, "wlevels") - if (fit_type == "lm") { - wlevels <- ind_to_cat(wlevels) + if (has_wlevels) { + wlevels <- attr(output, "wlevels") + if (fit_type == "lm") { + wlevels <- ind_to_cat(wlevels) + } + w_names <- colnames(wlevels) } - - w_names <- colnames(wlevels) + # mf0 is a list of datasets for multiple group models mf0 <- switch(fit_type, lavaan = lavaan::lavInspect(fit, "data"), lavaan.mi = lav_data_used(fit, drop_colon = FALSE), @@ -243,9 +302,20 @@ plot.cond_indirect_effects <- function( x_levels_list <- replicate(nrow(wlevels), x_levels, simplify = FALSE) } if (graph_type == "tumble") { - x_subsets <- lapply(split(wlevels, seq_len(nrow(wlevels))), - x_for_wlevels, - mf = mf0, x = x) + if (has_wlevels && !has_groups) { + x_subsets <- lapply(split(wlevels, seq_len(nrow(wlevels))), + x_for_wlevels, + mf = mf0, x = x) + } + if (!has_wlevels && has_groups) { + x_subsets <- lapply(mf0, function(xx) { + xx[, x, drop = TRUE] + }) + } + if (has_wlevels && has_groups) { + # TODO + # - Support objects with both wlevels and groups. + } x_levels_list <- lapply(x_subsets, gen_levels, method = x_method, @@ -258,19 +328,56 @@ plot.cond_indirect_effects <- function( x_levels_m <- do.call(rbind, x_levels_list) plot_df_xstart <- data.frame(x = x_levels_m[, 1]) plot_df_xend <- data.frame(x = x_levels_m[, 2]) - mf2 <- data.frame(lapply(as.data.frame(dat0), sum_col), - check.names = FALSE) - mf2 <- mf2[, -(which(colnames(mf2) %in% c(x, w_names)))] - plot_df_xstart <- cbind(plot_df_xstart, wlevels, mf2) - plot_df_xend <- cbind(plot_df_xend, wlevels, mf2) + if (has_groups) { + # mf2 has rows equal to ngroups + mf2 <- lapply(dat0, function(xx) { + data.frame(lapply(as.data.frame(xx), sum_col), + check.names = FALSE) + }) + mf2 <- lapply(mf2, function(xx) { + xx[, -(which(colnames(xx) %in% c(x)))] + }) + mf2 <- do.call(rbind, mf2) + } else { + # mf2 has one single row + # TO-THINK: + # - Maybe we can use different values of other variables + mf2 <- data.frame(lapply(as.data.frame(dat0), sum_col), + check.names = FALSE) + mf2 <- mf2[, -(which(colnames(mf2) %in% c(x, w_names)))] + } + if (has_groups) { + plot_df_xstart <- cbind(plot_df_xstart, mf2) + plot_df_xend <- cbind(plot_df_xend, mf2) + } else { + plot_df_xstart <- cbind(plot_df_xstart, wlevels, mf2) + plot_df_xend <- cbind(plot_df_xend, wlevels, mf2) + } colnames(plot_df_xstart)[1] <- x colnames(plot_df_xend)[1] <- x - plot_df_xstart[, y] <- stats::predict(fit_list, - x = x, y = y, m = m, - newdata = plot_df_xstart) - plot_df_xend[, y] <- stats::predict(fit_list, - x = x, y = y, m = m, - newdata = plot_df_xend) + if (has_groups) { + plot_df_xstart_i <- split(plot_df_xstart, + seq_len(nrow(plot_df_xstart)), + drop = FALSE) + plot_df_xend_i <- split(plot_df_xend, + seq_len(nrow(plot_df_xend)), + drop = FALSE) + plot_df_xstart[, y] <- mapply(stats::predict, + object = fit_list, + newdata = plot_df_xstart_i, + MoreArgs = list(x = x, y = y, m = m)) + plot_df_xend[, y] <- mapply(stats::predict, + object = fit_list, + newdata = plot_df_xend_i, + MoreArgs = list(x = x, y = y, m = m)) + } else { + plot_df_xstart[, y] <- stats::predict(fit_list, + x = x, y = y, m = m, + newdata = plot_df_xstart) + plot_df_xend[, y] <- stats::predict(fit_list, + x = x, y = y, m = m, + newdata = plot_df_xend) + } if (missing(x_label)) x_label <- x if (missing(y_label)) y_label <- y @@ -285,16 +392,44 @@ plot.cond_indirect_effects <- function( lm = lm2ptable(fit)$implied_stats) } if (x_standardized) { - x_sd <- sqrt(implied_stats$cov[x, x]) - x_mean <- implied_stats$mean[x] - if (is.null(x_mean)) x_mean <- 0 + if (has_groups) { + # x_sd and x_mean are vectors if ngroups > 1 + group_labels <- names(fit_list) + implied_stats <- implied_stats[group_labels] + x_sd <- sapply(implied_stats, function(xx) { + xx$cov[x, x] + }) + x_mean <- sapply(implied_stats, function(xx) { + out <- xx$mean[x] + out <- ifelse(is.null(out), 0, out) + out + }) + } else { + x_sd <- sqrt(implied_stats$cov[x, x]) + x_mean <- implied_stats$mean[x] + if (is.null(x_mean)) x_mean <- 0 + } plot_df_xstart[, x] <- (plot_df_xstart[, x] - x_mean) / x_sd plot_df_xend[, x] <- (plot_df_xend[, x] - x_mean) / x_sd } if (y_standardized) { - y_sd <- sqrt(implied_stats$cov[y, y]) - y_mean <- implied_stats$mean[y] - if (is.null(y_mean)) y_mean <- 0 + if (has_groups) { + # y_sd and y_mean are vectors if ngroups > 1 + group_labels <- names(fit_list) + implied_stats <- implied_stats[group_labels] + y_sd <- sapply(implied_stats, function(xx) { + xx$cov[y, y] + }) + y_mean <- sapply(implied_stats, function(xx) { + out <- xx$mean[y] + out <- ifelse(is.null(out), 0, out) + out + }) + } else { + y_sd <- sqrt(implied_stats$cov[y, y]) + y_mean <- implied_stats$mean[y] + if (is.null(y_mean)) y_mean <- 0 + } plot_df_xstart[, y] <- (plot_df_xstart[, y] - y_mean) / y_sd plot_df_xend[, y] <- (plot_df_xend[, y] - y_mean) / y_sd } @@ -326,15 +461,18 @@ plot.cond_indirect_effects <- function( title <- "Conditional Effects" } } - - plot_df_xstart$wlevels <- rownames(wlevels) - plot_df_xend$wlevels <- rownames(wlevels) + if (has_groups) { + plot_df_xstart$wlevels <- group_labels + plot_df_xend$wlevels <- group_labels + } else { + plot_df_xstart$wlevels <- rownames(wlevels) + plot_df_xend$wlevels <- rownames(wlevels) + } plot_df <- rbind(plot_df_xstart, plot_df_xend) - p <- ggplot2::ggplot() + - ggplot2::geom_point(ggplot2::aes_string(x = x, - y = y, - colour = "wlevels"), + ggplot2::geom_point(ggplot2::aes(x = .data[[x]], + y = .data[[y]], + colour = .data[["wlevels"]]), data = plot_df, size = point_size) + ggplot2::geom_segment(ggplot2::aes( @@ -343,7 +481,7 @@ plot.cond_indirect_effects <- function( y = plot_df_xstart[, y], yend = plot_df_xend[, y], colour = plot_df_xstart$wlevels - ), size = line_width) + ), linewidth = line_width) if (note_standardized & !is.null(cap_std)) { if (!is.null(cap_txt)) { @@ -352,6 +490,11 @@ plot.cond_indirect_effects <- function( cap_txt <- cap_std } } + if (has_groups && graph_type == "default") { + cap_txt <- paste0(cap_txt, + "\n", + "Graph type is set to tumble for multiple group models.") + } out <- p + ggplot2::labs(title = title, caption = cap_txt) + @@ -369,3 +512,6 @@ plot.cond_indirect_effects <- function( } out } + +utils::globalVariables(".data") + diff --git a/R/print_cond_indirect_effect.R b/R/print_cond_indirect_effect.R index e32f34d6..e822fb97 100644 --- a/R/print_cond_indirect_effect.R +++ b/R/print_cond_indirect_effect.R @@ -107,6 +107,8 @@ print.cond_indirect_effects <- function(x, digits = 3, pvalue_digits = 3, se = FALSE, ...) { + # TODO: + # - Support cases with both moderators and groups. full_output <- attr(x, "full_output") x_i <- full_output[[1]] my_call <- attr(x, "call") @@ -114,6 +116,15 @@ print.cond_indirect_effects <- function(x, digits = 3, boot_ci <- !is.null(x_i$boot_ci) has_ci <- FALSE ci_type <- NULL + has_groups <- ("group" %in% tolower(colnames(x))) + if (has_groups) { + group_labels <- unique(x$Group) + group_numbers <- unique(x$Group_ID) + } else { + group_labels <- NULL + group_numbers <- NULL + } + has_wlevels <- !is.null(attr(x, "wlevels")) if (!is.null(x_i$boot_ci)) { has_ci <- TRUE ci_type <- "boot" @@ -178,9 +189,15 @@ print.cond_indirect_effects <- function(x, digits = 3, } } out1 <- data.frame(out, check.names = FALSE) - wlevels <- attr(x, "wlevels") - w0 <- colnames(attr(wlevels, "wlevels")) - w1 <- colnames(wlevels) + if (has_wlevels) { + wlevels <- attr(x, "wlevels") + w0 <- colnames(attr(wlevels, "wlevels")) + w1 <- colnames(wlevels) + } else { + wlevels <- NULL + w0 <- NULL + w1 <- NULL + } mcond <- names(x_i$components) cond_str <- "" cond_str2 <- "" @@ -194,8 +211,15 @@ print.cond_indirect_effects <- function(x, digits = 3, cat("\n== Conditional effects ==\n") } cat("\n Path:", path) - cat("\n Conditional on moderator(s):", paste0(w0, collapse = ", ")) - cat("\n Moderator(s) represented by:", paste0(w1, collapse = ", ")) + if (has_wlevels) { + cat("\n Conditional on moderator(s):", paste0(w0, collapse = ", ")) + cat("\n Moderator(s) represented by:", paste0(w1, collapse = ", ")) + } + if (has_groups) { + tmp <- paste0(group_labels, "[", group_numbers, "]", + collapse = ", ") + cat("\n Conditional on group(s):", tmp) + } xold <- x x <- out1 cat("\n\n") @@ -243,9 +267,15 @@ print.cond_indirect_effects <- function(x, digits = 3, cat(" - The 'ind' column shows the", cond_str, "effects.", sep = " ") } cat("\n ") + tmp <- ifelse(has_wlevels && has_groups, + "the moderators and/or groups.", + ifelse(has_wlevels, + "the moderator(s).", + "the group(s).")) cat(strwrap(paste("\n -", paste(sQuote(mcond), collapse = ","), "is/are the path coefficient(s) along the path", - "conditional on the moderators."), exdent = 3), sep = "\n") + "conditional on", + tmp), exdent = 3), sep = "\n") cat("\n") } invisible(x) diff --git a/R/print_indirect.R b/R/print_indirect.R index c4483dc0..51a86d76 100644 --- a/R/print_indirect.R +++ b/R/print_indirect.R @@ -134,6 +134,11 @@ print.indirect <- function(x, standardized <- (standardized_x && standardized_y) has_ci <- FALSE ci_type <- NULL + if (is.numeric(x$group_number)) { + has_group <- TRUE + } else { + has_group <- FALSE + } if (isTRUE(!is.null(x$boot_ci))) { has_ci <- TRUE ci_type <- "boot" @@ -194,6 +199,11 @@ print.indirect <- function(x, } else { path <- paste(x0, "->", y0) } + if (has_group) { + path <- paste0(x$group_label, "[", + x$group_number, "]: ", + path) + } std_str <- "" if (standardized) { std_str <- paste0("(Both ", sQuote(x0), @@ -216,7 +226,8 @@ print.indirect <- function(x, cat("\n== Conditional", cond_str2, "Effect", std_str, " ==") } else { - cat("\n==", cond_str2, "Effect ==") + cat("\n==", cond_str2, "Effect", + std_str, "==") } cat("\n") ptable <- data.frame(Factor = "Path:", Value = path) @@ -283,7 +294,7 @@ print.indirect <- function(x, } else { if (is.null(x$op)) { ptable <- rbind(ptable, - c(ifelse(has_m, "Indirect Effect", "Effect"), + c(ifelse(has_m, "Indirect Effect:", "Effect:"), formatC(x$indirect, digits = digits, format = "f"))) } else { ptable <- rbind(ptable, @@ -292,6 +303,12 @@ print.indirect <- function(x, } if (has_ci) {ptable <- rbind(ptable, b_row, b_row2, b_row3)} } + if (has_group) { + # ptable <- rbind(ptable, + # c("Group Label:", x$group_label)) + # ptable <- rbind(ptable, + # c("Group Number:", x$group_number)) + } ptable <- data.frame(lapply(ptable, format)) colnames(ptable) <- c("", "") print(ptable, row.names = FALSE) @@ -342,6 +359,13 @@ print.indirect <- function(x, cat(strwrap(tmp1), sep = "\n") } } + print_note <- FALSE + if (standardized_x || + standardized_y || + has_group) { + print_note <- TRUE + } + note_str <- character(0) if (has_m & !is.list(mpathnames)) { if (has_w) { out <- data.frame(mpathnames, m0c, m0) @@ -356,8 +380,27 @@ print.indirect <- function(x, cat("\n") print(out, digits = digits, row.names = FALSE) if (standardized_x || standardized_y) { - cat("\nNOTE: The effects of the component paths are from the model, not standardized.") - } + note_str <- c(note_str, + strwrap("- The effects of the component paths are from the model, not standardized.", + exdent = 2)) + if (has_group) { + note_str <- c(note_str, + strwrap("- SD(s) in the selected group is/are used in standardiziation.", + exdent = 2)) + } + } + } + if (has_group) { + note_str <- c(note_str, + strwrap("- The group label is printed before each path.", + exdent = 2)) + note_str <- c(note_str, + strwrap("- The group number in square brackets is the number used internally in lavaan.", + exdent = 2)) + } + if (length(note_str) > 0) { + cat("\nNOTE:\n") + cat(note_str, sep = "\n") } cat("\n") invisible(x) diff --git a/R/total_indirect_effect_list.R b/R/total_indirect_effect_list.R index 75ebc247..22457563 100644 --- a/R/total_indirect_effect_list.R +++ b/R/total_indirect_effect_list.R @@ -70,13 +70,55 @@ total_indirect_effect <- function(object, if (!is.list(object)) { stop("object is not a list") } + has_groups <- indirect_list_has_groups(object) + if (has_groups) { + tmp <- group_labels_and_numbers_list(object) + group_labels <- tmp$label + group_numbers <- tmp$number + i1 <- sapply(object, function(xx) xx$group_label) + i2 <- split(seq_along(object), i1) + out0 <- lapply(i2, function(xx) { + object_i <- object[xx] + total_indirect_effect_i(object_i, + x = x, + y = y) + }) + out0 <- out0[group_labels] + tmp <- sapply(out0, function(xx) !identical(xx, NA)) + if (!any(tmp)) { + stop("In all groups, no paths from ", x, " to ", y, ".") + } + out0 <- out0[tmp] + } else { + # ix <- sapply(object, function(z) z$x) + # iy <- sapply(object, function(z) z$y) + # i0 <- (ix %in% x) & (iy %in% y) + # if (isFALSE(any(i0))) { + # return(NA) + # } + # p1 <- object[i0] + out0 <- total_indirect_effect_i(object = object, + x = x, + y = y) + if (identical(out0, NA)) { + stop("No paths from ", x, " to ", y, ".") + } + } + out0 + } + +#' @noRd + +total_indirect_effect_i <- function(object, + x, + y) { ix <- sapply(object, function(z) z$x) iy <- sapply(object, function(z) z$y) i0 <- (ix %in% x) & (iy %in% y) if (isFALSE(any(i0))) { - stop("No valid path was found.") + return(NA) } p1 <- object[i0] out <- Reduce(`+`, p1) out - } + } \ No newline at end of file diff --git a/README.md b/README.md index 16be47c7..e14d4d61 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ [![R-CMD-check](https://github.com/sfcheung/manymome/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/sfcheung/manymome/actions/workflows/R-CMD-check.yaml) -(Version 0.1.14.1, updated on 2024-03-24, [release history](https://sfcheung.github.io/manymome/news/index.html)) +(Version 0.1.14.5, updated on 2024-03-31, [release history](https://sfcheung.github.io/manymome/news/index.html)) # manymome @@ -35,6 +35,9 @@ by multiple regression. The package was introduced in: while Monte Carlo is supported for models fitted by `lavaan::sem()`. +- Multigroup models fitted by `lavaan::sem()` + are also supported in 0.1.14.2 and later versions. + # Advantages - **A Simpler Workflow** @@ -47,7 +50,11 @@ by multiple regression. The package was introduced in: nearly any variable, to nearly any other variables, conditional on nearly any moderators, and at any levels of the moderators. - (See `vignette("manymome")` for details.) + (See `vignette("manymome")` for details.) This is particularly + convenient for multigroup models fitted by `lavaan::sem()`, + which are supported in 0.1.14.2 and later versions + (see [this guide](https://sfcheung.github.io/manymome/articles/med_mg.html), + for an illustration). - **Supports Both SEM-Based and Regression-Based Analysis** @@ -61,6 +68,10 @@ by multiple regression. The package was introduced in: No limit on the number of predictors, mediators, and outcome variables, other than those by `lavaan::sem()` and `lm()`. + For multigroup models fitted by `lavaan::sem()`, + there is no inherent limit on the number of groups, + other than the limit due to `lavaan::sem(), if any + (supported in 0.1.14.2 and later versions). - **Supports Standardized Effects** @@ -78,7 +89,7 @@ by multiple regression. The package was introduced in: Supports datasets with missing data through `lavaan::sem()` with full information maximum likelihood (`fiml`). - Since version 0.1.9.8, it also supports missing data handled + In version 0.1.9.8 or later, it also supports missing data handled by multiple imputation if the models are fitted by `semTools::sem.mi()` or `semTools::runMI()` (see `vignette("do_mc_lavaan_mi")`). @@ -103,6 +114,14 @@ by multiple regression. The package was introduced in: latent variables for models fitted by `lavaan::sem()` (see `vignette("med_lav")`). +- **Support Treating Group As a Moderator** + + For multigroup models fitted by `lavaan::sem()`, it supports + comparing the direct or indirect effects along any path + between any two groups. That is, it uses the grouping variable + as a moderator (illustrated [here](https://sfcheung.github.io/manymome/articles/med_mg.html); + supported in 0.1.14.2 and later versions). + # Limitations Despite the aforementioned advantages, the current version of @@ -110,8 +129,6 @@ Despite the aforementioned advantages, the current version of - Does not (officially) support categorical predictors. -- Does not support multisample models (although `lavaan` does). - - Does not support multilevel models (although `lavaan` does). - For bootstrapping, only supports nonparametric bootstrapping and percentile diff --git a/_pkgdown.yml b/_pkgdown.yml index ae1b9082..d5ef4a88 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -39,6 +39,10 @@ navbar: - text: Moderated Mediation href: articles/mome_lm.html - text: ------- + - text: "" + - text: Mediation in Multigroup Models + href: articles/med_mg.html + - text: ------- - text: "" - text: Monte Carlo Confidence Intervals with Multiple Imputation href: articles/do_mc_lavaan_mi.html diff --git a/data-raw/test_data_2_med_mg.R b/data-raw/test_data_2_med_mg.R new file mode 100644 index 00000000..cc3a7b96 --- /dev/null +++ b/data-raw/test_data_2_med_mg.R @@ -0,0 +1,49 @@ +# Generate data +library(lavaan) +set.seed(3143214) +n <- 100 +ivs <- MASS::mvrnorm(n, c(10, 2, 5), diag(3)) +x <- ivs[, 1] +c1 <- ivs[, 2] +c2 <- ivs[, 3] +m <- 10 + .9 * x + .2 * c1 - .2 * c2 + rnorm(n, 0, 1) +y <- 5 + .6 * m + .2 * x + .1 * c1 - .1 * c2 + rnorm(n, 0, 2) +dat <- data.frame(x, m, y, c1, c2) +head(dat) +colMeans(dat) +summary(lm_m <- lm(m ~ x + c1 + c2, dat)) +summary(lm_y <- lm(y ~ m + x + c1 + c2, dat)) +dat1 <- dat +n <- 150 +ivs <- MASS::mvrnorm(n, c(10, 2, 5), diag(3)) +x <- ivs[, 1] +c1 <- ivs[, 2] +c2 <- ivs[, 3] +m <- 10 + .4 * x + .2 * c1 - .2 * c2 + rnorm(n, 0, 1) +y <- 5 + .9 * m + .2 * x + .1 * c1 - .1 * c2 + rnorm(n, 0, 2) +dat <- data.frame(x, m, y, c1, c2) +head(dat) +colMeans(dat) +summary(lm_m <- lm(m ~ x + c1 + c2, dat)) +summary(lm_y <- lm(y ~ m + x + c1 + c2, dat)) +dat2 <- dat +dat1$group <- "Group A" +dat2$group <- "Group B" +dat <- rbind(dat1, dat2) +mod <- +" +m ~ c(a1, a2) * x + c1 + c2 +y ~ c(b1, b2) * m + x + c1 + c2 +a1b1 := a1 * b1 +a2b2 := a2 * b2 +abdiff := a2b2 - a1b1 +adiff := a2 - a1 +bdiff := b2 - b1 +" +fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "group") +parameterEstimates(fit)[41:45, ] +head(dat) +data_med_mg <- dat +usethis::use_data(data_med_mg, overwrite = TRUE) + diff --git a/data-raw/test_data_med_mg.R b/data-raw/test_data_med_mg.R new file mode 100644 index 00000000..11bed522 --- /dev/null +++ b/data-raw/test_data_med_mg.R @@ -0,0 +1,41 @@ +# Generate data +library(lavaan) +n1 <- 100 +n2 <- 100 +mod1 <- +" +m11 ~ .4*x1 + .0*x2 + .3*c1 + .1*c2 +m12 ~ .5*m11 + .0*x1 + .0*x2 + .0*c1 + .0*c2 +m2 ~ .0*x1 + .4*x2 + .1*c1 + .1*c2 +y1 ~ .1*m11 + .4*m12 + .0*m2 + .1*x1 + .0*x2 + .0*c1 + .0*c2 +y2 ~ .0*m11 + .0*m12 + .4*m2 + .0*x1 + .2*x2 + .0*c1 + .0*c2 +" +mod2 <- +" +m11 ~ .2*x1 + .0*x2 + .3*c1 + .1*c2 +m12 ~ .5*m11 + .0*x1 + .0*x2 + .0*c1 + .0*c2 +m2 ~ .0*x1 + .4*x2 + .1*c1 + .1*c2 +y1 ~ .1*m11 + .2*m12 + .0*m2 + .1*x1 + .0*x2 + .0*c1 + .0*c2 +y2 ~ .0*m11 + .0*m12 + .4*m2 + .0*x1 + .2*x2 + .0*c1 + .0*c2 +" +dat1 <- simulateData(model = mod1, sample.nobs = n1, seed = 1234) +dat2 <- simulateData(model = mod2, sample.nobs = n2, seed = 5678) +head(dat1) +dat1 <- as.data.frame(scale(dat, + center = c(-10, -5, -5, -5, -8, -5, -10, -5, -10), scale = FALSE)) +dat1 <- as.data.frame(scale(dat, + center = c(-10, -5, -5, -5, -8, -7, -12, -5, -10), scale = FALSE)) +mod <- +" +m11 ~ x1 + x2 + c1 + c2 +m12 ~ m11 + c1 + c2 +m2 ~ x1 + x2 + c1 + c2 +y1 ~ m11 + m12 + x1 + x2 + c1 + c2 +y2 ~ m2 + x1 + x2 + c1 + c2 +" +dat <- rbind(dat1, dat2) +dat$group <- rep(c("Group A", "Group B"), times = c(n1, n2)) +fit <- sem(mod, dat, group = "group") +summary(fit, rsquare = TRUE) +data_med_complicated_mg <- dat +usethis::use_data(data_med_complicated_mg, overwrite = TRUE) diff --git a/data/data_med_complicated_mg.rda b/data/data_med_complicated_mg.rda new file mode 100644 index 00000000..ce652a64 Binary files /dev/null and b/data/data_med_complicated_mg.rda differ diff --git a/data/data_med_mg.rda b/data/data_med_mg.rda new file mode 100644 index 00000000..930e56f5 Binary files /dev/null and b/data/data_med_mg.rda differ diff --git a/man/all_indirect_paths.Rd b/man/all_indirect_paths.Rd index ad3c30cf..136ceb1b 100644 --- a/man/all_indirect_paths.Rd +++ b/man/all_indirect_paths.Rd @@ -5,7 +5,13 @@ \alias{all_paths_to_df} \title{Enumerate All Indirect Effects in a Model} \usage{ -all_indirect_paths(fit = NULL, exclude = NULL, x = NULL, y = NULL) +all_indirect_paths( + fit = NULL, + exclude = NULL, + x = NULL, + y = NULL, + group = NULL +) all_paths_to_df(all_paths) } @@ -36,6 +42,19 @@ the outcome variables in at least one regression equation will be included in the search.} +\item{group}{Either the group number +as appeared in the \code{\link[=summary]{summary()}} +or \code{\link[lavaan:parameterEstimates]{lavaan::parameterEstimates()}} +output of a \link[lavaan:lavaan-class]{lavaan::lavaan} object, +or the group label as used in +the \link[lavaan:lavaan-class]{lavaan::lavaan} object. +Used only when the number of +groups is greater than one. Default +is \code{NULL}. If not specified by the model +has more than one group, than paths +that appears in at least one group +will be included in the output.} + \item{all_paths}{An \code{all_paths}-class object. For example, the output of \code{\link[=all_indirect_paths]{all_indirect_paths()}}.} } @@ -60,6 +79,17 @@ and \code{m}, to be used by \code{indirect_effect()}. \details{ It makes use of \code{\link[igraph:all_simple_paths]{igraph::all_simple_paths()}} to identify paths in a model. +\subsection{Multigroup Models}{ + +Since Version 0.1.14.2, support for +multigroup models has been added for models +fitted by \code{lavaan}. If a model has more +than one group and \code{group} is not +specified, than paths in all groups +will be returned. If \code{group} is +specified, than only paths in the +selected group will be returned. +} } \section{Functions}{ \itemize{ @@ -100,6 +130,32 @@ out3 <- all_indirect_paths(fit, exclude = c("c1", "c2"), out3 names(out3) +# Multigroup models + +data(data_med_complicated_mg) +mod <- +" +m11 ~ x1 + x2 + c1 + c2 +m12 ~ m11 + c1 + c2 +m2 ~ x1 + x2 + c1 + c2 +y1 ~ m11 + m12 + x1 + x2 + c1 + c2 +y2 ~ m2 + x1 + x2 + c1 + c2 +" +fit <- sem(mod, data_med_complicated_mg, group = "group") +summary(fit) + +all_indirect_paths(fit, + x = "x1", + y = "y1") +all_indirect_paths(fit, + x = "x1", + y = "y1", + group = 1) +all_indirect_paths(fit, + x = "x1", + y = "y1", + group = "Group B") + } \seealso{ \code{\link[=indirect_effect]{indirect_effect()}}, \code{\link[=lm2list]{lm2list()}}. diff --git a/man/cond_indirect.Rd b/man/cond_indirect.Rd index 778b337c..26989095 100644 --- a/man/cond_indirect.Rd +++ b/man/cond_indirect.Rd @@ -38,7 +38,8 @@ cond_indirect( ci_out = NULL, save_ci_full = FALSE, save_ci_out = TRUE, - ci_type = NULL + ci_type = NULL, + group = NULL ) cond_indirect_effects( @@ -67,6 +68,7 @@ cond_indirect_effects( mc_out = NULL, ci_out = NULL, ci_type = NULL, + groups = NULL, ... ) @@ -89,6 +91,7 @@ indirect_effect( make_cluster_args = list(), progress = TRUE, save_boot_full = FALSE, + save_boot_out = TRUE, mc_ci = FALSE, mc_out = NULL, save_mc_full = FALSE, @@ -96,7 +99,8 @@ indirect_effect( ci_out = NULL, save_ci_full = FALSE, save_ci_out = TRUE, - ci_type = NULL + ci_type = NULL, + group = NULL ) many_indirect_effects(paths, ...) @@ -186,7 +190,7 @@ them from \code{fit}.} \code{TRUE}, \code{boot_out} is \code{NULL}, and bootstrap standard errors not requested if \code{fit} is a -\linkS4class{lavaan} object, this function +\link[lavaan:lavaan-class]{lavaan::lavaan} object, this function will do bootstrapping on \code{fit}. \code{R} is the number of bootstrap samples. Default is 100. For Monte Carlo @@ -304,6 +308,16 @@ other arguments supplied, will override \code{boot_ci} and \code{mc_ci}.} +\item{group}{Either the group number +as appeared in the \code{\link[=summary]{summary()}} +or \code{\link[lavaan:parameterEstimates]{lavaan::parameterEstimates()}} +output of a \link[lavaan:lavaan-class]{lavaan::lavaan} object, +or the group label as used in +the \link[lavaan:lavaan-class]{lavaan::lavaan} object. +Used only when the number of +groups is greater than one. Default +is \code{NULL}.} + \item{wlevels}{The output of \code{\link[=merge_mod_levels]{merge_mod_levels()}}, or the moderator(s) to be passed to @@ -383,6 +397,17 @@ output is a list of the outputs from for creating the levels of moderators. Default is \code{list()}.} +\item{groups}{Either a vector of +group numbers +as appeared in the \code{\link[=summary]{summary()}} +or \code{\link[lavaan:parameterEstimates]{lavaan::parameterEstimates()}} +output of a \link[lavaan:lavaan-class]{lavaan::lavaan} object, +or a vector of group labels as used in +the \link[lavaan:lavaan-class]{lavaan::lavaan} object. +Used only when the number of +groups is greater than one. Default +is \code{NULL}.} + \item{...}{For \code{\link[=many_indirect_effects]{many_indirect_effects()}}, these are arguments to be passed to \code{\link[=indirect_effect]{indirect_effect()}}.} @@ -399,8 +424,7 @@ these are arguments to be passed to These two classes of objects have their own print methods for printing -the results (see \code{\link[=print.indirect]{print.indirect()}} -and \code{\link[=print.cond_indirect_effects]{print.cond_indirect_effects()}}). +the results (see \code{\link[=print.indirect]{print.indirect()}} and \code{\link[=print.cond_indirect_effects]{print.cond_indirect_effects()}}). They also have a \code{coef} method for extracting the estimates (\code{\link[=coef.indirect]{coef.indirect()}} and @@ -515,6 +539,36 @@ If \code{boot_out} or \code{mc_out} is set, arguments such as \code{R}, \code{seed}, and \code{parallel} will be ignored. +\subsection{Multigroup Models}{ + +Since Version 0.1.14.2, support for +multigroup models has been added for models +fitted by \code{lavaan}. Both bootstrapping +and Monte Carlo confidence intervals +are supported. When used on +a multigroup model: +\itemize{ +\item For \code{\link[=cond_indirect]{cond_indirect()}} and +\code{\link[=indirect_effect]{indirect_effect()}}, users need to +specify the \code{group} argument +(by number or label). When using +\code{\link[=cond_indirect_effects]{cond_indirect_effects()}}, if +\code{group} is not set, all groups wil +be used and the indirect effect +in each group will be computed, +kind of treating group as a moderator. +\item For \code{\link[=many_indirect_effects]{many_indirect_effects()}}, +the paths can be generated from a +multigroup models. +\item Currently, \code{\link[=cond_indirect_effects]{cond_indirect_effects()}} +does not support a multigroup model +with moderators on the path selected. +The function \code{\link[=cond_indirect]{cond_indirect()}} does +not have this limitation but users +need to manually specify the desired +value of the moderator(s). +} +} } \section{Functions}{ \itemize{ @@ -583,6 +637,49 @@ cond_indirect_effects(x = "x", y = "m1", cond_indirect_effects(x = "x", y = "y", m = "m1", wlevels = w1levels, fit = fit) +# Multigroup models for cond_indirect_effects() + +dat <- data_med_mg +mod <- +" +m ~ x + c1 + c2 +y ~ m + x + c1 + c2 +" +fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, + group = "group") + +# If a model has more than one group, +# it will be used as a 'moderator'. +cond_indirect_effects(x = "x", y = "y", m = "m", + fit = fit) + + +# Multigroup model for indirect_effect() + +dat <- data_med_mg +mod <- +" +m ~ x + c1 + c2 +y ~ m + x + c1 + c2 +" +fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, + group = "group") + +# If a model has more than one group, +# the argument 'group' must be set. +ind1 <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group A") +ind1 +ind2 <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = 2) +ind2 + # Examples for many_indirect_effects(): @@ -597,18 +694,38 @@ y ~ m12 + m2 + m11 + x + c1 + c2 " fit <- sem(mod, data_serial_parallel, fixed.x = FALSE) - # All indirect paths from x to y paths <- all_indirect_paths(fit, x = "x", y = "y") paths - # Indirect effect estimates out <- many_indirect_effects(paths, fit = fit) out +# Multigroup models for many_indirect_effects() + +data(data_med_complicated_mg) +mod <- +" +m11 ~ x1 + x2 + c1 + c2 +m12 ~ m11 + c1 + c2 +m2 ~ x1 + x2 + c1 + c2 +y1 ~ m11 + m12 + x1 + x2 + c1 + c2 +y2 ~ m2 + x1 + x2 + c1 + c2 +" +fit <- sem(mod, data_med_complicated_mg, group = "group") +summary(fit) + +paths <- all_indirect_paths(fit, + x = "x1", + y = "y1") +paths +# Indirect effect estimates for all paths in all groups +out <- many_indirect_effects(paths, + fit = fit) +out } \seealso{ diff --git a/man/data_med_complicated_mg.Rd b/man/data_med_complicated_mg.Rd new file mode 100644 index 00000000..c08ca082 --- /dev/null +++ b/man/data_med_complicated_mg.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dat_5_mg.R +\docType{data} +\name{data_med_complicated_mg} +\alias{data_med_complicated_mg} +\title{Sample Dataset: A Complicated +Mediation Model With Two Groups} +\format{ +A data frame with 300 rows +and 5 variables: +\describe{ +\item{x1}{Predictor 1. Numeric.} +\item{x2}{Predictor 2. Numeric.} +\item{m11}{Mediator 1 in Path 1. Numeric.} +\item{m12}{Mediator 2 in Path 1. Numeric.} +\item{m2}{Mediator in Path 2. Numeric.} +\item{y1}{Outcome variable 1. Numeric.} +\item{y2}{Outcome variable 2. Numeric.} +\item{c1}{Control variable. Numeric.} +\item{c2}{Control variable. Numeric.} +\item{group}{Group variable. Character. 'Group A' or 'Group B'} +} +} +\usage{ +data_med_complicated_mg +} +\description{ +A mediation model with +two predictors, two pathways, and +two groups. +} +\examples{ +library(lavaan) +data(data_med_complicated_mg) +dat <- data_med_complicated_mg +mod <- +" +m11 ~ x1 + x2 + c1 + c2 +m12 ~ m11 + c1 + c2 +m2 ~ x1 + x2 + c1 + c2 +y1 ~ m11 + m12 + x1 + x2 + c1 + c2 +y2 ~ m2 + x1 + x2 + c1 + c2 +" +fit <- sem(mod, dat, group = "group") +summary(fit) +} +\keyword{datasets} diff --git a/man/data_med_mg.Rd b/man/data_med_mg.Rd new file mode 100644 index 00000000..2510f247 --- /dev/null +++ b/man/data_med_mg.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dat_2_med_mg.R +\docType{data} +\name{data_med_mg} +\alias{data_med_mg} +\title{Sample Dataset: Simple +Mediation With Two Groups} +\format{ +A data frame with 100 rows +and 5 variables: +\describe{ +\item{x}{Predictor. Numeric.} +\item{m}{Mediator. Numeric.} +\item{y}{Outcome variable. Numeric.} +\item{c1}{Control variable. Numeric.} +\item{c2}{Control variable. Numeric.} +\item{group}{Group variable. Character. "Group A" or "Group B"} +} +} +\usage{ +data_med_mg +} +\description{ +A simple mediation +model with two groups. +} +\examples{ +library(lavaan) +data(data_med_mg) +mod <- +" +m ~ c(a1, a2) * x + c1 + c2 +y ~ c(b1, b2) * m + x + c1 + c2 +a1b1 := a1 * b1 +a2b2 := a2 * b2 +abdiff := a2b2 - a1b1 +" +fit <- sem(mod, data_med_mg, fixed.x = FALSE, + group = "group") +parameterEstimates(fit) +} +\keyword{datasets} diff --git a/man/do_boot.Rd b/man/do_boot.Rd index d7a53b7c..e16294f6 100644 --- a/man/do_boot.Rd +++ b/man/do_boot.Rd @@ -105,6 +105,15 @@ It determines the type of the fit object automatically and then calls \code{\link[=lm2boot_out]{lm2boot_out()}}, \code{\link[=fit2boot_out]{fit2boot_out()}}, or \code{\link[=fit2boot_out_do_boot]{fit2boot_out_do_boot()}}. +\subsection{Multigroup Models}{ + +Since Version 0.1.14.2, support for +multigroup models has been added for models +fitted by \code{lavaan}. The implementation +of bootstrapping is identical to +that used by \code{lavaan}, with resampling +done within each group. +} } \examples{ data(data_med_mod_ab1) diff --git a/man/do_mc.Rd b/man/do_mc.Rd index b6cf824b..e1fc4320 100644 --- a/man/do_mc.Rd +++ b/man/do_mc.Rd @@ -103,6 +103,12 @@ each call to that the same set of Monte Carlo estimates is used in all subsequent analysis. +\subsection{Multigroup Models}{ + +Since Version 0.1.14.2, support for +multigroup models has been added for models +fitted by \code{lavaan}. +} } \section{Functions}{ \itemize{ diff --git a/man/indirect_i.Rd b/man/indirect_i.Rd index 9ba7fdb3..9a799360 100644 --- a/man/indirect_i.Rd +++ b/man/indirect_i.Rd @@ -21,7 +21,8 @@ indirect_i( data = NULL, expand = TRUE, warn = TRUE, - allow_mixing_lav_and_obs = TRUE + allow_mixing_lav_and_obs = TRUE, + group = NULL ) } \arguments{ @@ -122,6 +123,16 @@ omitted intentionally.} \code{TRUE}, it accepts a path with both latent variables and observed variables. Default is \code{TRUE}.} + +\item{group}{Either the group number +as appeared in the \code{\link[=summary]{summary()}} +or \code{\link[lavaan:parameterEstimates]{lavaan::parameterEstimates()}} +output of an \code{lavaan}-class object, +or the group label as used in +the \code{lavaan}-class object. +Used only when the number of +groups is greater than one. Default +is NULL.} } \value{ It returns an diff --git a/man/math_indirect.Rd b/man/math_indirect.Rd index 6f27b1b6..18a57091 100644 --- a/man/math_indirect.Rd +++ b/man/math_indirect.Rd @@ -33,8 +33,7 @@ For now, only \code{+} operator and \code{-} operator are supported. These operators can be used to estimate and test a function of effects between -the same pair of variables but along -different paths. +the same pair of variables. For example, they can be used to compute and test the total effects @@ -50,8 +49,7 @@ not valid if \item the two paths do not start from the same variable, \item the two paths do not end at the -same variable, (c) a path appears in -both objects, +same variable, \item moderators are involved but they are not set to the same values in both objects, and @@ -61,6 +59,28 @@ both objects, and estimates stored in \code{mc_out}, if any, are not identical. } +\subsection{Multigroup Models}{ + +Since Version 0.1.14.2, support for +multigroup models has been added for models +fitted by \code{lavaan}. Both bootstrapping +and Monte Carlo confidence intervals +are supported. These operators can +be used to compute and test the +difference of an indirect effect +between two groups. This can also +be used to compute and test the +difference between a function of +effects between groups, for example, +the total indirect effects between +two groups. + +The operators are flexible and allow +users to do many possible computations. +Therefore, users need to make sure +that the function of effects is +meaningful. +} } \examples{ library(lavaan) @@ -93,6 +113,35 @@ out123 <- out1 + out2 + out3 out123 coef(out1) + coef(out2) + coef(out3) +# Multigroup model with indirect effects + +dat <- data_med_mg +mod <- +" +m ~ x + c1 + c2 +y ~ m + x + c1 + c2 +" +fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, + group = "group") + +# If a model has more than one group, +# the argument 'group' must be set. +ind1 <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group A") +ind1 +ind2 <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = 2) +ind2 + +# Compute the difference in indirect effects between groups +ind2 - ind1 + } \seealso{ \code{\link[=indirect_effect]{indirect_effect()}} and diff --git a/man/plot.cond_indirect_effects.Rd b/man/plot.cond_indirect_effects.Rd index a4d0955f..84df5324 100644 --- a/man/plot.cond_indirect_effects.Rd +++ b/man/plot.cond_indirect_effects.Rd @@ -115,9 +115,11 @@ points as used in \item{graph_type}{If \code{"default"}, the typical line-graph with equal end-points will be plotted. If -\code{"tubmle"}, then the tumble graph +\code{"tumble"}, then the tumble graph proposed by Bodner (2016) will be -plotted. Default is \code{"default"}.} +plotted. Default is \code{"default"} +for single-group models, and +\code{"tumble"} for multigroup models.} \item{...}{Additional arguments. Ignored.} @@ -142,12 +144,26 @@ output. It plots the conditional effect from \code{x} to \code{y} in a model for different -levels of the moderators. +levels of the moderators. For +multigroup models, the group will +be the 'moderator' and one line is +drawn for each group. It does not support conditional indirect effects. If there is one or more mediators in \code{x}, it will raise an error. +\subsection{Multigroup Models}{ + +Since Version 0.1.14.2, support for +multigroup models has been added for models +fitted by \code{lavaan}. If the effect +for each group is drawn, the +\code{graph_type} is automatically switched +to \code{"bumble"} and the means and SDs +in each group will be used to determine +the locations of the points. +} } \examples{ library(lavaan) @@ -190,6 +206,26 @@ out_2 <- cond_indirect_effects(wlevels = out_mm_2, x = "x", y = "m3", fit = fit2 plot(out_2) plot(out_2, graph_type = "tumble") +# Multigroup models + +dat <- data_med_mg +mod <- +" +m ~ x + c1 + c2 +y ~ m + x + c1 + c2 +" +fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, se = "none", baseline = FALSE, + group = "group") + +# For a multigroup model, group will be used as +# a moderator +out <- cond_indirect_effects(x = "m", + y = "y", + fit = fit) +out +plot(out) + + } \references{ diff --git a/rebuild_vignettes.R b/rebuild_vignettes.R index 48dafb4f..4903281a 100644 --- a/rebuild_vignettes.R +++ b/rebuild_vignettes.R @@ -12,3 +12,12 @@ knitr::knit("do_mc.Rmd.original", output = "do_mc.Rmd") knitr::knit("do_mc_lavaan_mi.Rmd.original", output = "do_mc_lavaan_mi.Rmd") setwd(base_dir) + +# For articles + +base_dir <- getwd() + +setwd("vignettes/articles") +knitr::knit("med_mg.Rmd.originaL", output = "med_mg.Rmd") + +setwd(base_dir) diff --git a/tests/testthat/test_cond_indirect_math.R b/tests/testthat/test_cond_indirect_math.R index 54e7f42a..ea4d9caf 100644 --- a/tests/testthat/test_cond_indirect_math.R +++ b/tests/testthat/test_cond_indirect_math.R @@ -137,7 +137,8 @@ test_that("math for indirect: mediation", { expect_equal(out1m1m2m3boot2$indirect_raw, outm_boot$indirect_raw + outm2_boot$indirect_raw + outm3_boot$indirect_raw) expect_equal(outm1minus2boot$indirect_raw, outm_boot$indirect_raw - outm2_boot$indirect_raw) expect_equal(outm1minus3boot$indirect_raw, outm_boot$indirect_raw - outm3_boot$indirect_raw) - expect_error(outm_boot + outm_boot) + # No longer a requirement + # expect_error(outm_boot + outm_boot) expect_equal(coef(outm1minus2), outm1minus2$indirect, ignore_attr = TRUE) expect_equal(coef(outm1plus3), outm1plus3$indirect, ignore_attr = TRUE) expect_equal(coef(outm1minus2boot), outm1minus2boot$indirect, ignore_attr = TRUE) diff --git a/tests/testthat/test_get_prod_mi.R b/tests/testthat/test_get_prod_mi.R index 40140f63..3476cc61 100644 --- a/tests/testthat/test_get_prod_mi.R +++ b/tests/testthat/test_get_prod_mi.R @@ -26,7 +26,8 @@ out_mi_6 <- get_prod(x = "w2", y = "m2", fit = fit1_mi) out_mi_7 <- get_prod(x = "m2", y = "m3", fit = fit1_mi) # No need to compare b values -out_mi_3$b[] <- out_3$b +out_mi_3$b <- NULL +out_3$b <- NULL out_mi_6$b[] <- out_6$b test_that("get_prod for lavaan.mi", { diff --git a/tests/testthat/test_mg_boot.R b/tests/testthat/test_mg_boot.R new file mode 100644 index 00000000..781ac363 --- /dev/null +++ b/tests/testthat/test_mg_boot.R @@ -0,0 +1,1054 @@ +skip_on_cran() + +library(testthat) +library(manymome) +suppressMessages(library(lavaan)) + +dat <- modmed_x1m3w4y1 +n <- nrow(dat) +set.seed(860314) +dat$gp <- sample(c("gp1", "gp2", "gp3"), n, replace = TRUE) +dat$city <- sample(c("alpha", "beta", "gamma", "sigma"), n, replace = TRUE) + +dat <- cbind(dat, factor2var(dat$gp, prefix = "gp", add_rownames = FALSE)) +dat <- cbind(dat, factor2var(dat$city, prefix = "city", add_rownames = FALSE)) + +mod <- +" +m3 ~ m1 + x +y ~ m2 + m3 + x + w4 + x:w4 +" + +dat$xw4 <- dat$x * dat$w4 +dat$m3w4 <- dat$m3 * dat$w4 +mod2 <- +" +m3 ~ m1 + x +y ~ m2 + m3 + x + w4 + xw4 + w3 + m3:w3 + m3w4 +" + +# This model is not exactly identical to the previous one +# due the labelled variances +mod2_chk <- +" +m3 ~ m1 + c(a1, a2, a3)*x +y ~ m2 + c(b1, b2, b3)*m3 + x + w4 + xw4 + w3 + c(d31, d32, d33)*m3:w3 + c(d41, d42, d43)*m3w4 +ab1 := a1*b1 +ab2 := a2*b2 +ab3 := a3*b3 +ab2_d := a2*(b2 + 1*d32 + 2*d42) +ab1_d := a1*(b1 + 3*d31 + (-2)*d41) +x ~~ c(v_x1, v_x2, v_x3) * x +ab1_stdx := a1*b1*sqrt(v_x1) +ab2_stdx := a2*b2*sqrt(v_x2) +ab3_stdx := a3*b3*sqrt(v_x3) +" + +dat$m3w3 <- dat$m3 * dat$w3 +mod3 <- +" +m3 ~ m1 + x +y ~ m2 + m3 + x + w4 + xw4 + w3 + m3w3 + m3w4 +" + +mod_med <- +" +m1 ~ x +m2 ~ m1 +y ~ m2 + x +" + + +# Check against lavaan + +fit <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) +fit_boot <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, + se = "boot", bootstrap = 5, + warn = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2"), + iseed = 2345) +do_boot_out <- fit2boot_out_do_boot(fit, R = 5, + seed = 2345, + progress = FALSE, + parallel = FALSE) +lav_boot <- lavInspect(fit_boot, "boot") + +test_that("Check against lavaan boot", { + expect_equal(do_boot_out[[3]]$est$est[1:4], + unname(lav_boot[3, 1:4])) + }) + +# get_implied_i_lavaan + +fit_tmp <- sem(mod, dat[-c(1:10), ], meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) +my_implied <- get_implied_i(coef(fit), fit_tmp) +lav_implied <- lavInspect(fit, "implied") + +test_that("Check against lavaan implied", { + expect_equal(unclass(my_implied$cov$gp3), + unclass(lav_implied$gp3$cov)) + }) + +# get_prod + +fit2 <- sem(mod2, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) +fit2_chk <- sem(mod2_chk, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +fit2_ng <- sem(mod2, dat, meanstructure = TRUE, fixed.x = FALSE) +dat_tmp <- lav_data_used(fit2) +est_tmp <- lav_est(fit2, se = FALSE, ci = FALSE) +est_tmp2 <- est_tmp +est_tmp2$est <- est_tmp2$est * .5 +est_tmp_ng <- lav_est(fit2_ng, se = FALSE, ci = FALSE) +est_tmp2_ng <- est_tmp_ng +est_tmp2_ng$est <- est_tmp2_ng$est * .5 + +test_that("get_prod and friends", { + expect_true(setequal(c("x", "w4"), + find_product(dat_tmp, "xw4"))) + expect_true(setequal(names(find_all_products(dat_tmp)), + c("m3w4", "xw4"))) + tmp <- get_b(x = "xw4", + y = "y", + est = est_tmp) + tmpchk <- est_tmp[(est_tmp$rhs == "xw4") & + (est_tmp$op == "~"), "est"] + expect_equal(unname(tmp), + unname(tmpchk)) + tmp <- get_prod(x = "x", + y = "y", + fit = fit2, + expand = TRUE) + expect_true(length(tmp$b$xw4) == 3) + }) + +tmp1 <- get_prod(x = "x", + y = "y", + fit = fit2, + expand = TRUE) +tmp1_ng <- get_prod(x = "x", + y = "y", + fit = fit2_ng, + expand = TRUE) + +# indirect_i + +suppressWarnings(tmp2 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = "gp1")) +suppressWarnings(tmp3 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 1)) +# tmp2_chk <- est_tmp[(est_tmp$lhs == "m3") & +# (est_tmp$rhs == "x") & +# (est_tmp$group == 2), "est"] * +# est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3") & +# (est_tmp$group == 2), "est"] +tmp2_chk <- coef(fit2_chk, type = "user")["ab2"] +# tmp3_chk <- est_tmp[(est_tmp$lhs == "m3") & +# (est_tmp$rhs == "x") & +# (est_tmp$group == 1), "est"] * +# est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3") & +# (est_tmp$group == 1), "est"] +tmp3_chk <- coef(fit2_chk, type = "user")["ab1"] + +test_that("indirect_effect and multigrop", { + expect_equal(unname(coef(tmp2)), + unname(tmp2_chk), + tolerance = 1e-05) + # Can't just compare them. Don't know why. + expect_equal(unname(coef(tmp3)) - unname(tmp3_chk), + 0) + }) + +# mod_levels +# mod_levels_list + +test_that("mod_levels: multigroup", { + expect_error(mod_levels(fit2, w = "w3", w_method = "percentile")) + expect_error(mod_levels_list("w3", "w4", fit = fit2)) + }) + + + +# cond_indirect + +suppressWarnings(tmp2 <- cond_indirect(x = "x", + y = "y", + m = "m3", + fit = fit2, + wvalues = c(w3 = 1, w4 = 2), + group = 2)) +suppressWarnings(tmp3 <- cond_indirect(x = "x", + y = "y", + m = "m3", + fit = fit2, + wvalues = c(w3 = 3, w4 = -2), + group = "gp3")) +# tmp2_chk <- est_tmp[(est_tmp$lhs == "m3") & +# (est_tmp$rhs == "x") & +# (est_tmp$group == 2), "est"] * +# (est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3") & +# (est_tmp$group == 2), "est"] + +# est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3:w3") & +# (est_tmp$group == 2), "est"] + +# est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3w4") & +# (est_tmp$group == 2), "est"] * 2) +tmp2_chk <- coef(fit2_chk, type = "user")["ab2_d"] +# tmp3_chk <- est_tmp[(est_tmp$lhs == "m3") & +# (est_tmp$rhs == "x") & +# (est_tmp$group == 1), "est"] * +# (est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3") & +# (est_tmp$group == 1), "est"] + +# est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3:w3") & +# (est_tmp$group == 1), "est"] * 3+ +# est_tmp[(est_tmp$lhs == "y") & +# (est_tmp$rhs == "m3w4") & +# (est_tmp$group == 1), "est"] * -2) +tmp3_chk <- coef(fit2_chk, type = "user")["ab1_d"] + +test_that("indirect_effect and multigrop", { + expect_equal(unname(coef(tmp2)), + unname(tmp2_chk), + tolerance = 1e-5) + expect_equal(unname(coef(tmp3)), + unname(tmp3_chk), + tolerance = 1e-4) + }) + +# indirect_i: stdx / stdy + +suppressWarnings(tmp2 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = "gp1", + standardized_x = TRUE)) +suppressWarnings(tmp3 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 3, + standardized_y = TRUE)) +suppressWarnings(tmp4 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 3, + standardized_y = TRUE, + standardized_x = TRUE)) +sd_x_2 <- sqrt(lavInspect(fit2, "implied")$gp1$cov["x", "x"]) +sd_y_2 <- sqrt(lavInspect(fit2, "implied")$gp1$cov["y", "y"]) +sd_x_3 <- sqrt(lavInspect(fit2, "implied")[[3]]$cov["x", "x"]) +sd_y_3 <- sqrt(lavInspect(fit2, "implied")[[3]]$cov["y", "y"]) +tmp2_chk <- est_tmp[(est_tmp$lhs == "m3") & + (est_tmp$rhs == "x") & + (est_tmp$group == 2), "est"] * + est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3") & + (est_tmp$group == 2), "est"] * sd_x_2 / 1 +tmp3_chk <- est_tmp[(est_tmp$lhs == "m3") & + (est_tmp$rhs == "x") & + (est_tmp$group == 3), "est"] * + est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3") & + (est_tmp$group == 3), "est"] * 1 / sd_y_3 + +test_that("indirect_effect and multigrop", { + expect_equal(unname(coef(tmp2)), + tmp2_chk) + expect_equal(unname(coef(tmp3)), + tmp3_chk) + expect_equal(unname(coef(tmp4)), + tmp3_chk * sd_x_3) + }) + +## Math + +suppressWarnings(tmp3b <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 2, + standardized_y = TRUE)) +suppressWarnings(tmp3c <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 1, + standardized_y = TRUE)) + +tmpmath3a <- tmp3b - tmp3 +tmpmath3b <- tmp3b + tmp3 +tmpmath3c <- tmpmath3b - tmp3c +test_that("indirect_effect and multigrop: Math", { + expect_equal(coef(tmpmath3a), + coef(tmp3b) - coef(tmp3)) + expect_equal(coef(tmpmath3b), + coef(tmp3b) + coef(tmp3)) + expect_equal(coef(tmpmath3c), + coef(tmp3b) + coef(tmp3) - coef(tmp3c)) + }) + +# cond_indirect: stdx / stdy + +suppressWarnings(tmp2 <- cond_indirect(x = "x", + y = "y", + m = "m3", + fit = fit2, + wvalues = c(w3 = 1, w4 = 2), + group = 2, + standardized_x = TRUE)) +suppressWarnings(tmp3 <- cond_indirect(x = "x", + y = "y", + m = "m3", + fit = fit2, + wvalues = c(w3 = 3, w4 = -2), + group = "gp3", + standardized_y = TRUE)) +sd_x_2 <- sqrt(lavInspect(fit2, "implied")$gp1$cov["x", "x"]) +sd_y_2 <- sqrt(lavInspect(fit2, "implied")$gp1$cov["y", "y"]) +sd_x_3 <- sqrt(lavInspect(fit2, "implied")[[1]]$cov["x", "x"]) +sd_y_3 <- sqrt(lavInspect(fit2, "implied")[[1]]$cov["y", "y"]) +tmp2_chk <- est_tmp[(est_tmp$lhs == "m3") & + (est_tmp$rhs == "x") & + (est_tmp$group == 2), "est"] * + (est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3") & + (est_tmp$group == 2), "est"] + + est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3:w3") & + (est_tmp$group == 2), "est"] + + est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3w4") & + (est_tmp$group == 2), "est"] * 2) * sd_x_2 +tmp3_chk <- est_tmp[(est_tmp$lhs == "m3") & + (est_tmp$rhs == "x") & + (est_tmp$group == 1), "est"] * + (est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3") & + (est_tmp$group == 1), "est"] + + est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3:w3") & + (est_tmp$group == 1), "est"] * 3+ + est_tmp[(est_tmp$lhs == "y") & + (est_tmp$rhs == "m3w4") & + (est_tmp$group == 1), "est"] * -2) * 1 / sd_y_3 + +test_that("indirect_effect and multigrop", { + expect_equal(unname(coef(tmp2)), + tmp2_chk) + expect_equal(unname(coef(tmp3)), + tmp3_chk) + }) + +# Check do_mc + +fit_eq <- sem(mod, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2"), + group.equal = "regressions") + +fit_mc_out <- do_mc(fit, R = 4, + seed = 2345, + progress = FALSE, + parallel = FALSE) + +fit_eq_mc_out <- do_mc(fit_eq, R = 4, + seed = 2345, + progress = FALSE, + parallel = FALSE) + +fit_mc_out <- do_mc(fit2, R = 50, + seed = 2345, + progress = FALSE, + parallel = FALSE) + +get_mc_est <- function(object, lhs, op = "~", rhs, group = NA) { + out <- sapply(object, function(x) { + esti <- x$est + out1 <- esti[(esti$lhs == lhs) & (esti$op == op) & (esti$rhs == rhs), ] + if (!is.na(group)) { + out1 <- out1[out1$group == group, ] + } + out1[, "est"] + }) + out + } + +get_mc_implied <- function(object, var, group = NA) { + out <- sapply(object, function(x) { + imp <- x$implied_stats$cov + if (!is.na(group)) { + out <- imp[[group]][var, var] + } else { + out <- imp[var, var] + } + out + }) + out + } + +# indirect_i: mc + +suppressWarnings(tmp2 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = "gp1", + mc_ci = TRUE, + mc_out = fit_mc_out)) +suppressWarnings(tmp3 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 1, + mc_ci = TRUE, + mc_out = fit_mc_out)) + +tmp2a <- get_mc_est(fit_mc_out, lhs = "m3", rhs = "x", group = 2) +tmp2b <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3", group = 2) +tmp2ab <- tmp2a * tmp2b +tmp3a <- get_mc_est(fit_mc_out, lhs = "m3", rhs = "x", group = 1) +tmp3b <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3", group = 1) +tmp3ab <- tmp3a * tmp3b + +test_that("indirect_effect and multigrop", { + expect_equal(tmp2$mc_indirect, + tmp2ab) + expect_equal(tmp3$mc_indirect, + tmp3ab) + }) + +# indirect_i: mc, stdx / stdy + +suppressWarnings(tmp2 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = "gp1", + mc_ci = TRUE, + mc_out = fit_mc_out, + standardized_y = TRUE)) +suppressWarnings(tmp3 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 1, + mc_ci = TRUE, + mc_out = fit_mc_out, + standardized_x = TRUE)) +suppressWarnings(tmp4 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 1, + mc_ci = TRUE, + mc_out = fit_mc_out, + standardized_y = TRUE, + standardized_x = TRUE)) + +tmp2a <- get_mc_est(fit_mc_out, lhs = "m3", rhs = "x", group = 2) +tmp2b <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3", group = 2) +tmp2ysd <- sqrt(get_mc_implied(fit_mc_out, var = "y", group = 2)) +tmp2ab <- tmp2a * tmp2b / tmp2ysd +tmp3a <- get_mc_est(fit_mc_out, lhs = "m3", rhs = "x", group = 1) +tmp3b <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3", group = 1) +tmp3xsd <- sqrt(get_mc_implied(fit_mc_out, var = "x", group = 1)) +tmp3ysd <- sqrt(get_mc_implied(fit_mc_out, var = "y", group = 1)) +tmp3ab <- tmp3a * tmp3b * tmp3xsd +tmp4ab <- tmp3a * tmp3b * tmp3xsd / tmp3ysd + +test_that("indirect_effect and multigrop", { + expect_equal(tmp2$mc_indirect, + tmp2ab) + expect_equal(tmp3$mc_indirect, + tmp3ab) + expect_equal(tmp4$mc_indirect, + tmp4ab) + }) + +## Math + +suppressWarnings(tmp3b <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 3, + mc_ci = TRUE, + mc_out = fit_mc_out, + standardized_x = TRUE)) + +tmpmath3a <- tmp3b - tmp3 +tmpmath3b <- tmp3b + tmp3 +tmpmath3c <- tmpmath3b + tmpmath3a +test_that("indirect_effect and multigrop, MC: Math", { + expect_equal(coef(tmpmath3a), + coef(tmp3b) - coef(tmp3)) + expect_equal(coef(tmpmath3b), + coef(tmp3b) + coef(tmp3)) + }) + + +# cond_indirect, stdx / stdy: mc + +suppressWarnings(tmp2 <- cond_indirect(x = "x", + y = "y", + m = "m3", + fit = fit2, + wvalues = c(w3 = 1, w4 = 2), + group = 2, + mc_ci = TRUE, + mc_out = fit_mc_out, + standardized_y = TRUE)) +suppressWarnings(tmp3 <- cond_indirect(x = "x", + y = "y", + m = "m3", + fit = fit2, + wvalues = c(w3 = 3, w4 = -2), + group = "gp3", + mc_ci = TRUE, + mc_out = fit_mc_out, + standardized_x = TRUE)) +suppressWarnings(tmp4 <- cond_indirect(x = "x", + y = "y", + m = "m3", + fit = fit2, + wvalues = c(w3 = 3, w4 = -2), + group = "gp3", + mc_ci = TRUE, + mc_out = fit_mc_out, + standardized_y = TRUE, + standardized_x = TRUE)) + +tmp2a <- get_mc_est(fit_mc_out, lhs = "m3", rhs = "x", group = 2) +tmp2b <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3", group = 2) +tmp2d1 <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3:w3", group = 2) +tmp2d2 <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3w4", group = 2) +tmp2ysd <- sqrt(get_mc_implied(fit_mc_out, var = "y", group = 2)) +tmp2ab <- tmp2a * (tmp2b + 1 * tmp2d1 + 2 * tmp2d2) / tmp2ysd + +tmp3a <- get_mc_est(fit_mc_out, lhs = "m3", rhs = "x", group = 1) +tmp3b <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3", group = 1) +tmp3d1 <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3:w3", group = 1) +tmp3d2 <- get_mc_est(fit_mc_out, lhs = "y", rhs = "m3w4", group = 1) +tmp3xsd <- sqrt(get_mc_implied(fit_mc_out, var = "x", group = 1)) +tmp3ysd <- sqrt(get_mc_implied(fit_mc_out, var = "y", group = 1)) +tmp3ab <- tmp3a * (tmp3b + 3 * tmp3d1 + (-2) * tmp3d2) * tmp3xsd + +tmp4ab <- tmp3a * (tmp3b + 3 * tmp3d1 + (-2) * tmp3d2) * tmp3xsd / tmp3ysd + +test_that("indirect_effect and multigrop", { + expect_equal(tmp2$mc_indirect, + tmp2ab) + expect_equal(tmp3$mc_indirect, + tmp3ab) + expect_equal(tmp4$mc_indirect, + tmp4ab) + }) + +# All direct paths + +mod_tmp <- +" +m3 ~ c(NA, 0, NA)*m1 +y ~ m3 +" + +fit_tmp <- sem(mod_tmp, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +test_that("All direct path: Multiple group", { + expect_equal(all_indirect_paths(fit), + all_indirect_paths(fit2), + ignore_attr = TRUE) + expect_equal(all_indirect_paths(fit, group = 2), + all_indirect_paths(fit2, group = "gp1"), + ignore_attr = TRUE) + expect_true(length(all_indirect_paths(fit_tmp, group = 2)) == 0) + }) + +# Many direct path + +mod_tmp <- +" +m3 ~ c(NA, 0, NA)*m1 + c(NA, NA, 0)*x +m2 ~ c(0, 0, NA)*m1 + c(NA, NA, 0)*x +w3 ~ c(NA, 0, 0)*m2 +y ~ c(NA, 0, NA)*m3 + w3 +" +mod_tmp_2 <- +" +m3 ~ c(NA, 0, NA)*m1 + c(NA, NA, 0)*x +m2 ~ c(NA, NA, NA)*m1 + c(NA, NA, NA)*x +w3 ~ c(NA, NA, NA)*m2 +y ~ c(0, NA, NA)*m3 + c(NA, NA, 0)*w3 +" +mod_tmp_ng <- +" +m3 ~ m1 + x +m2 ~ m1 + x +w3 ~ m2 +y ~ m3 + w3 +" + +fit_tmp <- sem(mod_tmp, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) +fit_tmp_2 <- sem(mod_tmp_2, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) +fit_tmp_ng <- sem(mod_tmp_ng, dat, meanstructure = TRUE, fixed.x = FALSE) + +all_tmp <- all_indirect_paths(fit_tmp) +all_paths <- all_paths_to_df(all_tmp) +all_ind <- many_indirect_effects(all_tmp, fit = fit_tmp) +ind_chk <- indirect_effect(x = "x", + y = "w3", + m = "m2", + fit = fit_tmp, + group = "gp3") +all_ind_2 <- many_indirect_effects(all_indirect_paths(fit_tmp_2), fit = fit_tmp_2) +all_ind_ng <- many_indirect_effects(all_indirect_paths(fit_tmp_ng), fit = fit_tmp_ng) + +test_that("many_indirect: multiple group", { + expect_equal(coef(all_ind[[3]]), + coef(ind_chk)) + }) + +# indirect_effects_from_list + +test_that("indirect_effects_from_list: multiple group", { + expect_equal(unname(indirect_effects_from_list(all_ind)$ind), + unname(coef(all_ind))) + }) + +# total_indirect_effect + +test_that("total_indirect_effect: multiple group", { + expect_equal(length(total_indirect_effect(all_ind_2, x = "x", y = "y")), 2) + expect_error(total_indirect_effect(all_ind_2, x = "m3", y = "w3")) + }) + +# Mediation only + +fit_med <- sem(mod_med, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +tmp1 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med) +tmp1_chk1 <- indirect_effect(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + group = 1) +tmp1_chk2 <- indirect_effect(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + group = 2) +tmp1_chk3 <- indirect_effect(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + group = 3) + +tmp2 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + groups = c(2, 1)) + +tmp3 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + groups = c("gp1", "gp3")) + +test_that("cond_indirect_effects for multiple group", { + expect_equal(unname(coef(tmp1)), + unname(c(coef(tmp1_chk1), + coef(tmp1_chk2), + coef(tmp1_chk3)))) + expect_equal(unname(coef(tmp2)), + unname(c(coef(tmp1_chk2), + coef(tmp1_chk1)))) + expect_equal(unname(coef(tmp3)), + unname(c(coef(tmp1_chk2), + coef(tmp1_chk1)))) + expect_error(tmp2 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + groups = c(10, 20))) + }) + + +# Group labels helpers + +chk1 <- lavTech(fit2, "group.label") +test_that("group labels helpers", { + expect_equal(group_labels_and_numbers(c(3, 1), fit2)$label, + chk1[c(3, 1)]) + expect_equal(group_labels_and_numbers(c("gp1", "gp3"), fit2)$number, + c(2, 1)) + expect_error(group_labels_and_numbers(c("gp5", "gp3"), fit2)) + expect_error(group_labels_and_numbers(10, fit2)) + expect_error(group_labels_and_numbers(1:2, "test")) + expect_equal(group_labels_and_numbers(fit = fit2)$label, + chk1) + expect_equal(group_labels_and_numbers(fit = fit2)$number, + seq_along(chk1)) + }) + +# print.cond_indirect_effects + +fit_med <- sem(mod_med, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +tmp1 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med) + +test_that("print.cond_indirect_effects: Multiple groups", { + expect_output(print(tmp1), + "Conditional on group(s)", + fixed = TRUE) + }) + + +# coef.cond_indirect_effects + +fit_med <- sem(mod_med, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +tmp1 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med) +tmp1_3 <- indirect_effect(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + group = 3) + +test_that("coef.cond_indirect_effects with multiple groups", { + expect_equal(unname(coef(tmp1)[3]), + unname(coef(tmp1_3))) + }) + +# [.cond_indirect_effects + +fit_med <- sem(mod_med, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +tmp1 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med) + +test_that("[.cond_indirect_effects: Multiple groups", { + expect_equal(unname(coef(tmp1[c(1, 3), ])), + as.data.frame(tmp1)[c(1, 3), "ind"]) + }) + +# cond_indirect_diff + +fit_med <- sem(mod_med, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +tmp1 <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med) + +tmp1_diff <- cond_indirect_diff(tmp1, + from = 3, + to = 1) + +tmp1_1 <- indirect_effect(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + group = 1) +tmp1_3 <- indirect_effect(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + group = 3) + +test_that("cond_indirect_diff: Multiple groups", { + expect_equal(unname(coef(tmp1_diff)), + unname(coef(tmp1_1) - coef(tmp1_3))) + }) + +# plot.cond_indirect_effects + +fit_med <- sem(mod_med, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +tmp1 <- cond_indirect_effects(x = "m2", + y = "y", + fit = fit_med) + +test_that("plot.cond_indirect_effects: multiple groups", { + expect_no_error(p <- plot(tmp1)) + expect_true(setequal(unique(p$layers[[1]]$data$wlevels), + unique(tmp1$Group))) + }) + + + + + +skip("Long tests: Test in interactive sections") + +# Indirect with bootstrap + +fit3 <- sem(mod3, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) + +fit2_boot_out <- do_boot(fit2, + R = 50, + seed = 1234, + parallel = FALSE, + progress = FALSE) + +suppressWarnings(fit2_chk_boot <- sem(mod2_chk, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2"), + se = "bootstrap", + bootstrap = 50, + iseed = 1234)) + +fit2_chk_boot_out <- do_boot(fit2_chk_boot) + +fit_med <- sem(mod_med, dat, meanstructure = TRUE, fixed.x = FALSE, + group = "gp", + group.label = c("gp3", "gp1", "gp2")) +fit_med_boot_out <- do_boot(fit_med, + R = 50, + seed = 1234, + parallel = FALSE, + progress = FALSE) + +suppressWarnings(tmp2 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = "gp1", + boot_ci = TRUE, + boot_out = fit2_boot_out)) +suppressWarnings(tmp3 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 3, + boot_ci = TRUE, + boot_out = fit2_boot_out)) + +suppressWarnings(tmp2_chk_boot <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2_chk_boot, + group = "gp1", + boot_ci = TRUE)) +suppressWarnings(tmp3_chk_boot <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2_chk_boot, + group = 3, + boot_ci = TRUE)) + +est_chk <- parameterEstimates(fit2_chk_boot) + +test_that("indirect_effect and multigrop", { + i <- match("ab2", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp2))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + i <- match("ab3", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp3))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + i <- match("ab2", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp2_chk_boot))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + i <- match("ab3", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp3_chk_boot))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + }) + +## Math + +suppressWarnings(tmp3b <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 2, + boot_ci = TRUE, + boot_out = fit2_boot_out)) + +tmpmath3a <- tmp3b - tmp3 +tmpmath3b <- tmp3b + tmp3 +test_that("indirect_effect and multigrop, boot: Math", { + expect_equal(coef(tmpmath3a), + coef(tmp3b) - coef(tmp3)) + expect_equal(coef(tmpmath3b), + coef(tmp3b) + coef(tmp3)) + expect_equal(tmpmath3a$boot_indirect, + tmp3b$boot_indirect - tmp3$boot_indirect) + expect_equal(tmpmath3b$boot_indirect, + tmp3b$boot_indirect + tmp3$boot_indirect) + }) + + +# Indirect with bootstrap: stdx / stdy + +suppressWarnings(tmp2 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = "gp1", + boot_ci = TRUE, + boot_out = fit2_boot_out, + standardized_x = TRUE)) +suppressWarnings(tmp3 <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 3, + boot_ci = TRUE, + boot_out = fit2_boot_out, + standardized_x = TRUE)) + +suppressWarnings(tmp2_chk_boot <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2_chk_boot, + group = "gp1", + boot_ci = TRUE, + standardized_x = TRUE)) +suppressWarnings(tmp3_chk_boot <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2_chk_boot, + group = 3, + boot_ci = TRUE, + standardized_x = TRUE)) + +est_chk <- parameterEstimates(fit2_chk_boot) + +test_that("indirect_effect and multigrop", { + i <- match("ab2_stdx", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp2))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + i <- match("ab3_stdx", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp3))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + i <- match("ab2_stdx", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp2_chk_boot))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + i <- match("ab3_stdx", est_chk$lhs) + expect_equal(unname(as.vector(confint(tmp3_chk_boot))), + unname(unlist(est_chk[i, c("ci.lower", "ci.upper")])), + tolerance = 1e-4) + }) + +## Math + +suppressWarnings(tmp3b <- indirect_effect(x = "x", + y = "y", + m = "m3", + fit = fit2, + group = 1, + boot_ci = TRUE, + boot_out = fit2_boot_out, + standardized_x = TRUE)) + +tmpmath3a <- tmp3b - tmp3 +tmpmath3b <- tmp3b + tmp3 +test_that("indirect_effect and multigrop, boot: Math", { + expect_equal(coef(tmpmath3a), + coef(tmp3b) - coef(tmp3)) + expect_equal(coef(tmpmath3b), + coef(tmp3b) + coef(tmp3)) + expect_equal(tmpmath3a$boot_indirect, + tmp3b$boot_indirect - tmp3$boot_indirect) + expect_equal(tmpmath3b$boot_indirect, + tmp3b$boot_indirect + tmp3$boot_indirect) + }) + +# confint.cond_indirect_effects + +tmp1_boot <- cond_indirect_effects(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + boot_ci = TRUE, + boot_out = fit_med_boot_out) +suppressWarnings(tmp1_2 <- indirect_effect(x = "x", + y = "y", + m = c("m1", "m2"), + fit = fit_med, + group = 3, + boot_ci = TRUE, + boot_out = fit_med_boot_out)) + +tmp1_boot_ci <- confint(tmp1_boot) + +test_that("confint.cond_indirect_effects with multiple groups", { + expect_equal(unname(unlist(tmp1_boot_ci[3, ])), + unname(as.vector(confint(tmp1_2)))) + }) + +# indirect_effects_from_list + +all_tmp <- all_indirect_paths(fit_med) +all_ind <- many_indirect_effects(all_tmp, + fit = fit_med, + boot_ci = TRUE, + boot_out = fit_med_boot_out) + +test_that("indirect_effects_from_list: multiple group", { + expect_equal(unname(indirect_effects_from_list(all_ind)$ind), + unname(coef(all_ind))) + expect_equal(unname(indirect_effects_from_list(all_ind)$`CI.lo`), + unname(confint(all_ind)[, 1])) + }) diff --git a/vignettes/articles/manymome_draw_med_mg-1.png b/vignettes/articles/manymome_draw_med_mg-1.png new file mode 100644 index 00000000..c89b3e1a Binary files /dev/null and b/vignettes/articles/manymome_draw_med_mg-1.png differ diff --git a/vignettes/articles/manymome_draw_mod_mg-1.png b/vignettes/articles/manymome_draw_mod_mg-1.png new file mode 100644 index 00000000..8342629f Binary files /dev/null and b/vignettes/articles/manymome_draw_mod_mg-1.png differ diff --git a/vignettes/articles/med_complicated-1.png b/vignettes/articles/med_complicated-1.png new file mode 100644 index 00000000..6676119d Binary files /dev/null and b/vignettes/articles/med_complicated-1.png differ diff --git a/vignettes/articles/med_mg.Rmd b/vignettes/articles/med_mg.Rmd new file mode 100644 index 00000000..962cdc12 --- /dev/null +++ b/vignettes/articles/med_mg.Rmd @@ -0,0 +1,933 @@ +--- +title: "Multigroup Models With Mediation Effects" +author: "Shu Fai Cheung & Sing-Hang Cheung" +date: "2024-03-31" +output: + html_document: + fig.align: "center" + toc: true + number_sections: false +bibliography: references.bib +csl: apa.csl +--- + + + +# Introduction + +This article is a brief illustration of how to use +[manymome](https://sfcheung.github.io/manymome/index.html) +([Cheung & Cheung, 2023](https://doi.org/10.3758/s13428-023-02224-z)) +to compute and test indirect effects +in a multigroup model fitted by +`lavaan`. [^mgver] + +This article only focuses on issues specific +to multigroup models. Readers are assumed to have basic +understanding on using `manymome`. +Please refer to +the [Get Started](https://sfcheung.github.io/manymome/articles/manymome.html) +guide for a full introduction, and +[this section](https://sfcheung.github.io/manymome/articles/manymome.html#mediation-only) +on an illustration on a mediation model. + +[^mgver]: Support for multigroup model was introduced +in Version 0.1.14.2. + +# Model + +This is the sample data set that comes with the package: + + +```r +library(manymome) +dat <- data_med_mg +print(head(dat), digits = 3) +#> x m y c1 c2 group +#> 1 10.11 17.0 17.4 1.9864 5.90 Group A +#> 2 9.75 16.6 17.5 0.7748 4.37 Group A +#> 3 9.81 17.9 14.9 0.0973 6.96 Group A +#> 4 10.15 19.7 18.0 2.3974 5.75 Group A +#> 5 10.30 17.7 20.7 3.2225 5.84 Group A +#> 6 10.01 18.9 20.7 2.3631 4.51 Group A +``` + +Suppose this is the model being fitted, with `c1` and +`c2` the control variables. The grouping variable is `group`, +with two possible values, `"Group A"` and `"Group B"`. + +![Simple Mediation Model](manymome_draw_med_mg-1.png) + +# Fitting the Model + +We first fit this multigroup model in +`lavaan::sem()` as usual. There is +no need to label any parameters because +`manymome` will extract the parameters +automatically. + + +```r +mod_med <- +" +m ~ x + c1 + c2 +y ~ m + x + c1 + c2 +" +fit <- sem(model = mod_med, + data = dat, + fixed.x = FALSE, + group = "group") +``` + +These are the estimates of the paths: + + +```r +summary(fit, + estimates = TRUE) +#> lavaan 0.6.17 ended normally after 1 iteration +#> +#> Estimator ML +#> Optimization method NLMINB +#> Number of model parameters 40 +#> +#> Number of observations per group: +#> Group A 100 +#> Group B 150 +#> +#> Model Test User Model: +#> +#> Test statistic 0.000 +#> Degrees of freedom 0 +#> Test statistic for each group: +#> Group A 0.000 +#> Group B 0.000 +#> +#> Parameter Estimates: +#> +#> Standard errors Standard +#> Information Expected +#> Information saturated (h1) model Structured +#> +#> +#> Group 1 [Group A]: +#> +#> Regressions: +#> Estimate Std.Err z-value P(>|z|) +#> m ~ +#> x 0.880 0.093 9.507 0.000 +#> c1 0.264 0.104 2.531 0.011 +#> c2 -0.316 0.095 -3.315 0.001 +#> y ~ +#> m 0.465 0.190 2.446 0.014 +#> x 0.321 0.243 1.324 0.186 +#> c1 0.285 0.204 1.395 0.163 +#> c2 -0.228 0.191 -1.195 0.232 +#> +#> Covariances: +#> Estimate Std.Err z-value P(>|z|) +#> x ~~ +#> c1 -0.080 0.107 -0.741 0.459 +#> c2 -0.212 0.121 -1.761 0.078 +#> c1 ~~ +#> c2 -0.071 0.104 -0.677 0.499 +#> +#> Intercepts: +#> Estimate Std.Err z-value P(>|z|) +#> .m 10.647 1.156 9.211 0.000 +#> .y 6.724 2.987 2.251 0.024 +#> x 9.985 0.111 90.313 0.000 +#> c1 2.055 0.097 21.214 0.000 +#> c2 4.883 0.107 45.454 0.000 +#> +#> Variances: +#> Estimate Std.Err z-value P(>|z|) +#> .m 1.006 0.142 7.071 0.000 +#> .y 3.633 0.514 7.071 0.000 +#> x 1.222 0.173 7.071 0.000 +#> c1 0.939 0.133 7.071 0.000 +#> c2 1.154 0.163 7.071 0.000 +#> +#> +#> Group 2 [Group B]: +#> +#> Regressions: +#> Estimate Std.Err z-value P(>|z|) +#> m ~ +#> x 0.597 0.081 7.335 0.000 +#> c1 0.226 0.087 2.610 0.009 +#> c2 -0.181 0.078 -2.335 0.020 +#> y ~ +#> m 1.110 0.171 6.492 0.000 +#> x 0.264 0.199 1.330 0.183 +#> c1 -0.016 0.186 -0.088 0.930 +#> c2 -0.072 0.165 -0.437 0.662 +#> +#> Covariances: +#> Estimate Std.Err z-value P(>|z|) +#> x ~~ +#> c1 0.102 0.079 1.299 0.194 +#> c2 -0.050 0.087 -0.574 0.566 +#> c1 ~~ +#> c2 0.109 0.083 1.313 0.189 +#> +#> Intercepts: +#> Estimate Std.Err z-value P(>|z|) +#> .m 7.862 0.924 8.511 0.000 +#> .y 1.757 2.356 0.746 0.456 +#> x 10.046 0.082 121.888 0.000 +#> c1 2.138 0.078 27.515 0.000 +#> c2 5.088 0.087 58.820 0.000 +#> +#> Variances: +#> Estimate Std.Err z-value P(>|z|) +#> .m 0.998 0.115 8.660 0.000 +#> .y 4.379 0.506 8.660 0.000 +#> x 1.019 0.118 8.660 0.000 +#> c1 0.906 0.105 8.660 0.000 +#> c2 1.122 0.130 8.660 0.000 +``` + +# Generate Bootstrap estimates + +We can use `do_boot()` to generate +the bootstrap estimates first +(see [this article](https://sfcheung.github.io/manymome/articles/do_boot.html) +for an illustration on this function). +The argument `ncores` can be omitted +if the default value is acceptable. + + +```r +fit_boot_out <- do_boot(fit = fit, + R = 5000, + seed = 53253, + ncores = 8) +#> 8 processes started to run bootstrapping. +#> The expected CPU time is about 0 second(s). +``` + +# Estimate Indirect Effects + +## Estimate Each Effect by `indirect_effect()` + +The function `indirect_effect()` can be used to as usual +to estimate an indirect effect +and form its bootstrapping or Monte Carlo +confidence interval along a path in a model +that starts with any numeric variable, ends with +any numeric variable, through any numeric variable(s). +A detailed illustration can be found in +[this section](https://sfcheung.github.io/manymome/articles/manymome.html#est_indirect). + +For a multigroup model, the only +difference is that users need to specify +the group using the argument `group`. +It can be set to the group label +as used in `lavaan` (`"Group A"` +or `"Group B"` in this example) +or the group number used in `lavaan` + + +```r +ind_gpA <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group A", + boot_ci = TRUE, + boot_out = fit_boot_out) +``` + +This is the output: + + +```r +ind_gpA +#> +#> == Indirect Effect == +#> +#> Path: Group A[1]: x -> m -> y +#> Indirect Effect: 0.409 +#> 95.0% Bootstrap CI: [0.096 to 0.753] +#> +#> Computation Formula: +#> (b.m~x)*(b.y~m) +#> Computation: +#> (0.87989)*(0.46481) +#> +#> Percentile confidence interval formed by nonparametric bootstrapping +#> with 5000 bootstrap samples. +#> +#> Coefficients of Component Paths: +#> Path Coefficient +#> m~x 0.880 +#> y~m 0.465 +#> +#> NOTE: +#> - The group label is printed before each path. +#> - The group number in square brackets is the number used internally in +#> lavaan. +``` + +The indirect effect from `x` to `y` through `m` in +`"Group A"` is 0.409, +with a 95% confidence interval of +[0.096, 0.753], +significantly different from zero (*p* < .05). + +We illustrate computing the indirect effect in +`"Group B"`, using group number: + + +```r +ind_gpB <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = 2, + boot_ci = TRUE, + boot_out = fit_boot_out) +``` + +This is the output: + + +```r +ind_gpB +#> +#> == Indirect Effect == +#> +#> Path: Group B[2]: x -> m -> y +#> Indirect Effect: 0.663 +#> 95.0% Bootstrap CI: [0.411 to 0.959] +#> +#> Computation Formula: +#> (b.m~x)*(b.y~m) +#> Computation: +#> (0.59716)*(1.11040) +#> +#> Percentile confidence interval formed by nonparametric bootstrapping +#> with 5000 bootstrap samples. +#> +#> Coefficients of Component Paths: +#> Path Coefficient +#> m~x 0.597 +#> y~m 1.110 +#> +#> NOTE: +#> - The group label is printed before each path. +#> - The group number in square brackets is the number used internally in +#> lavaan. +``` + +The indirect effect from `x` to `y` through `m` in +`"Group B"` is 0.663, +with a 95% confidence interval of +[0.096, 0.753], +also significantly different from zero (*p* < .05). + +## Treating Group as a "Moderator" + +Instead of computing the indirect effects one-by-one, +we can also treat the grouping variable as +a "moderator" and use +`cond_indirect_effects()` to compute +the indirect effects along a path for +all groups. The detailed illustration +of this function can be found [here](https://sfcheung.github.io/manymome/articles/manymome.html#conditional-indirect-effects). +When use on a multigroup model, +wwe can omit the argument `wlevels`. +The function will automatically identify +all groups in a model, and compute +the indirect effect of the requested +path in each model. + + +```r +ind <- cond_indirect_effects(x = "x", + y = "y", + m = "m", + fit = fit, + boot_ci = TRUE, + boot_out = fit_boot_out) +``` + +This is the output: + + +```r +ind +#> +#> == Conditional indirect effects == +#> +#> Path: x -> m -> y +#> Conditional on group(s): Group A[1], Group B[2] +#> +#> Group Group_ID ind CI.lo CI.hi Sig m~x y~m +#> 1 Group A 1 0.409 0.096 0.753 Sig 0.880 0.465 +#> 2 Group B 2 0.663 0.411 0.959 Sig 0.597 1.110 +#> +#> - [CI.lo to CI.hi] are 95.0% percentile confidence intervals by +#> nonparametric bootstrapping with 5000 samples. +#> - The 'ind' column shows the indirect effects. +#> - 'm~x','y~m' is/are the path coefficient(s) along the path conditional +#> on the group(s). +``` + +The results are identical to those computed +individually using `indirect_effect()`. Using +`cond_indirect_effects()` is convenient when +the number of groups is more than two. + +# Compute and Test Between-Group difference + +There are several ways to compute and test +the difference in indirect effects between +two groups. + +## Using the Math Operator `-` + +The math operator `-` (described [here](https://sfcheung.github.io/manymome/reference/math_indirect.html)) +can be used if the indirect effects +have been computed individually by +`indirect_effect()`. We have already +computed the path `x->m->y` before +for the two groups. Let us compute the +differences: + + +```r +ind_diff <- ind_gpB - ind_gpA +ind_diff +#> +#> == Indirect Effect == +#> +#> Path: Group B[2]: x -> m -> y +#> Path: Group A[1]: x -> m -> y +#> Function of Effects: 0.254 +#> 95.0% Bootstrap CI: [-0.173 to 0.685] +#> +#> Computation of the Function of Effects: +#> (Group B[2]: x->m->y) +#> -(Group A[1]: x->m->y) +#> +#> +#> Percentile confidence interval formed by nonparametric bootstrapping +#> with 5000 bootstrap samples. +#> +#> NOTE: +#> - The group label is printed before each path. +#> - The group number in square brackets is the number used internally in +#> lavaan. +``` + +The difference in indirect effects from `x` to `y` through `m` +is 0.254, +with a 95% confidence interval of +[-0.173, 0.685], +not significantly different from zero (*p* < .05). Therefore, +we conclude that the two groups are not significantly +different on the indirect effects. + +## Using `cond_indirect_diff()` + +If the indirect effects are computed using +`cond_indirect_effects()`, we can use the function +`cond_indirect_diff()` to compute the difference +(described [here](https://sfcheung.github.io/manymome/reference/cond_indirect_diff.html)) +This is more convenient than using the math +operator when the number of groups is +greater than two. + +Let us use `cond_indirect_diff()` on the +output of `cond_indirect_effects()`: + + +```r +ind_diff2 <- cond_indirect_diff(ind, + from = 1, + to = 2) +ind_diff2 +#> +#> == Conditional indirect effects == +#> +#> Path: x -> m -> y +#> Conditional on group(s): Group B[2], Group A[1] +#> +#> Group Group_ID ind CI.lo CI.hi Sig m~x y~m +#> 1 Group B 2 0.663 0.411 0.959 Sig 0.597 1.110 +#> 2 Group A 1 0.409 0.096 0.753 Sig 0.880 0.465 +#> +#> == Difference in Conditional Indirect Effect == +#> +#> Levels: +#> Group +#> To: Group B [2] +#> From: Group A [1] +#> +#> Levels compared: Group B [2] - Group A [1] +#> +#> Change in Indirect Effect: +#> +#> x y Change CI.lo CI.hi +#> Change x y 0.254 -0.173 0.685 +#> +#> - [CI.lo, CI.hi]: 95% percentile confidence interval. +``` + +The convention is `to` row minus `from` row. +Though may sound not intuitive, the printout +always states clearly which group is subtracted +from which group. The results are identical +to those using the math operator. + +# Advanced Skills + +## Standardized Indirect Effects + +Standardized indirect effects can be computed +as for single-group models (described [here](https://sfcheung.github.io/manymome/articles/manymome.html#standardized-indirect-effect)), +by setting `standardized_x` and/or `standardized_y`. +This is an example: + + +```r +std_gpA <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group A", + boot_ci = TRUE, + boot_out = fit_boot_out, + standardized_x = TRUE, + standardized_y = TRUE) +std_gpA +#> +#> == Indirect Effect (Both 'x' and 'y' Standardized) == +#> +#> Path: Group A[1]: x -> m -> y +#> Indirect Effect: 0.204 +#> 95.0% Bootstrap CI: [0.049 to 0.366] +#> +#> Computation Formula: +#> (b.m~x)*(b.y~m)*sd_x/sd_y +#> Computation: +#> (0.87989)*(0.46481)*(1.10557)/(2.21581) +#> +#> Percentile confidence interval formed by nonparametric bootstrapping +#> with 5000 bootstrap samples. +#> +#> Coefficients of Component Paths: +#> Path Coefficient +#> m~x 0.880 +#> y~m 0.465 +#> +#> NOTE: +#> - The effects of the component paths are from the model, not +#> standardized. +#> - SD(s) in the selected group is/are used in standardiziation. +#> - The group label is printed before each path. +#> - The group number in square brackets is the number used internally in +#> lavaan. +``` + + +```r +std_gpB <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group B", + boot_ci = TRUE, + boot_out = fit_boot_out, + standardized_x = TRUE, + standardized_y = TRUE) +std_gpB +#> +#> == Indirect Effect (Both 'x' and 'y' Standardized) == +#> +#> Path: Group B[2]: x -> m -> y +#> Indirect Effect: 0.259 +#> 95.0% Bootstrap CI: [0.166 to 0.360] +#> +#> Computation Formula: +#> (b.m~x)*(b.y~m)*sd_x/sd_y +#> Computation: +#> (0.59716)*(1.11040)*(1.00943)/(2.58386) +#> +#> Percentile confidence interval formed by nonparametric bootstrapping +#> with 5000 bootstrap samples. +#> +#> Coefficients of Component Paths: +#> Path Coefficient +#> m~x 0.597 +#> y~m 1.110 +#> +#> NOTE: +#> - The effects of the component paths are from the model, not +#> standardized. +#> - SD(s) in the selected group is/are used in standardiziation. +#> - The group label is printed before each path. +#> - The group number in square brackets is the number used internally in +#> lavaan. +``` + +In `"Group A"`, the (completely) standardized indirect effect +from `x` to `y` through `m` is +0.204. In +`"Group B"`, this effect is +0.259. + +Note that, unlike single-group model, in multigroup models, +the standardized indirect effect in a group uses the +the standard deviations of `x`- and `y`-variables in this group +to do the standardization. Therefore, two groups can have +different unstandardized +effects on a path but similar standardized effects on the +same path, or have similar unstandardized effects on a path +but different standardized effects on this path. This is a +known phenomenon in multigroup structural equation model. + +The difference in the two completely standardized indirect +effects can computed and tested using the math operator `-`: + + +```r +std_diff <- std_gpB - std_gpA +std_diff +#> +#> == Indirect Effect (Both 'x' and 'y' Standardized) == +#> +#> Path: Group B[2]: x -> m -> y +#> Path: Group A[1]: x -> m -> y +#> Function of Effects: 0.055 +#> 95.0% Bootstrap CI: [-0.133 to 0.245] +#> +#> Computation of the Function of Effects: +#> (Group B[2]: x->m->y) +#> -(Group A[1]: x->m->y) +#> +#> +#> Percentile confidence interval formed by nonparametric bootstrapping +#> with 5000 bootstrap samples. +#> +#> NOTE: +#> - The group label is printed before each path. +#> - The group number in square brackets is the number used internally in +#> lavaan. +``` + +The difference in completely standardized indirect effects +from `x` to `y` through `m` +is 0.055, +with a 95% confidence interval of +[-0.133, 0.245], +not significantly different from zero (*p* < .05). Therefore, +we conclude that the two groups are also not significantly +different on the completely standardized indirect effects. + +The function `cond_indirect_effects()` and +`cond_indirect_diff()` can also be used with standardization: + + +```r +std <- cond_indirect_effects(x = "x", + y = "y", + m = "m", + fit = fit, + boot_ci = TRUE, + boot_out = fit_boot_out, + standardized_x = TRUE, + standardized_y = TRUE) +std +#> +#> == Conditional indirect effects == +#> +#> Path: x -> m -> y +#> Conditional on group(s): Group A[1], Group B[2] +#> +#> Group Group_ID std CI.lo CI.hi Sig m~x y~m ind +#> 1 Group A 1 0.204 0.049 0.366 Sig 0.880 0.465 0.409 +#> 2 Group B 2 0.259 0.166 0.360 Sig 0.597 1.110 0.663 +#> +#> - [CI.lo to CI.hi] are 95.0% percentile confidence intervals by +#> nonparametric bootstrapping with 5000 samples. +#> - std: The standardized indirect effects. +#> - ind: The unstandardized indirect effects. +#> - 'm~x','y~m' is/are the path coefficient(s) along the path conditional +#> on the group(s). +``` + + +```r +std_diff2 <- cond_indirect_diff(std, + from = 1, + to = 2) +std_diff2 +#> +#> == Conditional indirect effects == +#> +#> Path: x -> m -> y +#> Conditional on group(s): Group B[2], Group A[1] +#> +#> Group Group_ID std CI.lo CI.hi Sig m~x y~m ind +#> 1 Group B 2 0.259 0.166 0.360 Sig 0.597 1.110 0.663 +#> 2 Group A 1 0.204 0.049 0.366 Sig 0.880 0.465 0.409 +#> +#> == Difference in Conditional Indirect Effect == +#> +#> Levels: +#> Group +#> To: Group B [2] +#> From: Group A [1] +#> +#> Levels compared: Group B [2] - Group A [1] +#> +#> Change in Indirect Effect: +#> +#> x y Change CI.lo CI.hi +#> Change x y 0.055 -0.133 0.245 +#> +#> - [CI.lo, CI.hi]: 95% percentile confidence interval. +#> - x standardized. +#> - y standardized. +``` + +The results, again, are identical to those using +`indirect_effect()` and the math operator `-`. + +## Finding All Indirect Paths in a Multigroup Model + +Suppose a model which has more than one, or has many, +indirect paths, is fitted to this dataset: + + +```r +dat2 <- data_med_complicated_mg +print(head(dat2), digits = 2) +#> m11 m12 m2 y1 y2 x1 x2 c1 c2 group +#> 1 1.05 1.17 0.514 0.063 1.027 1.82 -0.365 0.580 -0.3221 Group A +#> 2 -0.48 0.71 0.366 -1.278 -1.442 0.18 -0.012 0.620 -0.8751 Group A +#> 3 -1.18 -2.01 -0.044 -0.177 0.152 0.32 -0.403 0.257 -0.1078 Group A +#> 4 3.64 1.47 -0.815 1.309 0.052 0.98 0.139 0.054 1.2495 Group A +#> 5 -0.41 -0.38 -1.177 -0.151 0.255 -0.36 -1.637 0.275 0.0078 Group A +#> 6 0.18 -1.00 -0.119 -0.588 0.036 -0.53 0.349 0.618 -0.4073 Group A +``` + +![A Complicated Path Model](med_complicated-1.png) + +We first fit this model in `lavaan`: + + +```r +mod2 <- +" +m11 ~ x1 + x2 +m12 ~ m11 + x1 + x2 +m2 ~ x1 + x2 +y1 ~ m2 + m12 + m11 + x1 + x2 +y2 ~ m2 + m12 + m11 + x1 + x2 +" +fit2 <- sem(mod2, data = dat2, group = "group") +``` + +The function `all_indirect_paths()` can be used on a +multigroup model to identify indirect paths. The search +can be restricted by setting arguments such as +`x`, `y`, and `exclude` (see the [help page](file:///C:/GitHub/manymome/docs/reference/all_indirect_paths.html) +for details). + +For example, the following identify all paths from `x1` +to `y1`: + + +```r +paths_x1_y1 <- all_indirect_paths(fit = fit2, + x = "x1", + y = "y1") +``` + +If the `group` argument is not specified, it will automatically +identify all paths in all groups, as shown in the printout: + + +```r +paths_x1_y1 +#> Call: +#> all_indirect_paths(fit = fit2, x = "x1", y = "y1") +#> Path(s): +#> path +#> 1 Group A.x1 -> m11 -> m12 -> y1 +#> 2 Group A.x1 -> m11 -> y1 +#> 3 Group A.x1 -> m12 -> y1 +#> 4 Group A.x1 -> m2 -> y1 +#> 5 Group B.x1 -> m11 -> m12 -> y1 +#> 6 Group B.x1 -> m11 -> y1 +#> 7 Group B.x1 -> m12 -> y1 +#> 8 Group B.x1 -> m2 -> y1 +``` + +We can then use `many_indirect_effects()` to +compute the indirect effects for all paths identified: + + +```r +all_ind_x1_y1 <- many_indirect_effects(paths_x1_y1, + fit = fit2) +all_ind_x1_y1 +#> +#> == Indirect Effect(s) == +#> ind +#> Group A.x1 -> m11 -> m12 -> y1 0.079 +#> Group A.x1 -> m11 -> y1 0.106 +#> Group A.x1 -> m12 -> y1 -0.043 +#> Group A.x1 -> m2 -> y1 -0.000 +#> Group B.x1 -> m11 -> m12 -> y1 0.000 +#> Group B.x1 -> m11 -> y1 0.024 +#> Group B.x1 -> m12 -> y1 -0.000 +#> Group B.x1 -> m2 -> y1 0.004 +#> +#> - The 'ind' column shows the indirect effects. +#> +``` + +Bootstrapping and Monte Carlo confidence intervals can +be formed in the same way they are formed for single-group +models. + +## Computing, Testing, and Plotting Conditional Effects + +Though the focus is on indirect effect, +the main functions in `manymome` can also be used for +computing and plotting the effects along the direct path +between two variables. That is, we can focus on the +moderating effect of group on a direct path. + +For example, in the simple mediation model examined +above, suppose we are interested in the between-group +difference in the path from `m` to `y`, the "b path". +We can first +compute the conditional effect using `cond_indirect_effects()`, +without setting the mediator: + + +```r +path1 <- cond_indirect_effects(x = "m", + y = "y", + fit = fit, + boot_ci = TRUE, + boot_out = fit_boot_out) +path1 +#> +#> == Conditional effects == +#> +#> Path: m -> y +#> Conditional on group(s): Group A[1], Group B[2] +#> +#> Group Group_ID ind CI.lo CI.hi Sig y~m +#> 1 Group A 1 0.465 0.110 0.819 Sig 0.465 +#> 2 Group B 2 1.110 0.765 1.475 Sig 1.110 +#> +#> - [CI.lo to CI.hi] are 95.0% percentile confidence intervals by +#> nonparametric bootstrapping with 5000 samples. +#> - The 'ind' column shows the effects. +#> - 'y~m' is/are the path coefficient(s) along the path conditional on +#> the group(s). +``` + +The difference between the two paths can be tested +using bootstrapping confidence interval using +`cond_indirect_diff()`: + + +```r +path1_diff <- cond_indirect_diff(path1, + from = 1, + to = 2) +path1_diff +#> +#> == Conditional effects == +#> +#> Path: m -> y +#> Conditional on group(s): Group B[2], Group A[1] +#> +#> Group Group_ID ind CI.lo CI.hi Sig y~m +#> 1 Group B 2 1.110 0.765 1.475 Sig 1.110 +#> 2 Group A 1 0.465 0.110 0.819 Sig 0.465 +#> +#> == Difference in Conditional Indirect Effect == +#> +#> Levels: +#> Group +#> To: Group B [2] +#> From: Group A [1] +#> +#> Levels compared: Group B [2] - Group A [1] +#> +#> Change in Indirect Effect: +#> +#> x y Change CI.lo CI.hi +#> Change m y 0.646 0.148 1.152 +#> +#> - [CI.lo, CI.hi]: 95% percentile confidence interval. +``` + +Based on bootstrapping, the effect of `m` on `y` +in `"Group B"` is significantly greater than that +in `"Group A"` (*p* < .05). (This is compatible +with the conclusion on the indirect effects because +two groups can have no difference on `ab` even if +they differ on `a` and/or `b`.) + +The `plot` method for the output +of `cond_indirect_effects()` can also be used +for multigroup models: + +![Conditional Effects](manymome_draw_mod_mg-1.png) + +Note that, for multigroup models, the *tumble* +graph proposed by @bodner_tumble_2016 will +always be used. The position of a line for +a group is determined by the descriptive statistics +of this group (e.g, means and SDs). For example, +the line segment of `"Group A"` is far to the right +because `"Group A"` has a larger mean of `m` than +`"Group B"`: + + +```r +by(dat$m, dat$group, mean) +#> dat$group: Group A +#> [1] 18.43357 +#> ------------------------------------------------------------ +#> dat$group: Group B +#> [1] 13.42327 +``` + +It would be misleading if the two lines are plotted on the +same horizontal position, assuming incorrectly that the ranges +of `m` are similar in the two groups. + +The vertical positions of the two lines are similarly determined +by the distributions of other predictors in each +group (the control variables +and `x` in this example). + +Details of the `plot` method can be found +in the [help page](https://sfcheung.github.io/manymome/reference/plot.cond_indirect_effects.html). + +# Final Remarks + +There are some limitations on the support +for multigroup models. Currently, +multiple imputation is not supported. +Moreover, most functions do not (yet) +support multigroup models with +within-group moderators, except for +`cond_indirect()`. We would appreciate +users to report issues discovered when +using [manymome](https://sfcheung.github.io/manymome/index.html) +on multigroup models at [GitHub](https://github.com/sfcheung/manymome/issues). + +# Reference(s) diff --git a/vignettes/articles/med_mg.Rmd.original b/vignettes/articles/med_mg.Rmd.original new file mode 100644 index 00000000..3ad4e5c1 --- /dev/null +++ b/vignettes/articles/med_mg.Rmd.original @@ -0,0 +1,560 @@ +--- +title: "Multigroup Models With Mediation Effects" +author: "Shu Fai Cheung & Sing-Hang Cheung" +date: "`r Sys.Date()`" +output: + html_document: + fig.align: "center" + toc: true + number_sections: false +bibliography: references.bib +csl: apa.csl +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "" +) +``` + +# Introduction + +This article is a brief illustration of how to use +[manymome](https://sfcheung.github.io/manymome/index.html) +([Cheung & Cheung, 2023](https://doi.org/10.3758/s13428-023-02224-z)) +to compute and test indirect effects +in a multigroup model fitted by +`lavaan`. [^mgver] + +This article only focuses on issues specific +to multigroup models. Readers are assumed to have basic +understanding on using `manymome`. +Please refer to +the [Get Started](https://sfcheung.github.io/manymome/articles/manymome.html) +guide for a full introduction, and +[this section](https://sfcheung.github.io/manymome/articles/manymome.html#mediation-only) +on an illustration on a mediation model. + +[^mgver]: Support for multigroup model was introduced +in Version 0.1.14.2. + +# Model + +This is the sample data set that comes with the package: + +```{r dataset_me_mg} +library(manymome) +dat <- data_med_mg +print(head(dat), digits = 3) +``` + +Suppose this is the model being fitted, with `c1` and +`c2` the control variables. The grouping variable is `group`, +with two possible values, `"Group A"` and `"Group B"`. + +```{r manymome_draw_med_mg, echo = FALSE, fig.cap = "Simple Mediation Model"} +library(semPlot) +suppressMessages(library(lavaan)) +mod <- +" +m ~ x + c1 + c2 +y ~ m + x + c1 + c2 +" +fit0 <- sem(mod, dat, do.fit = FALSE, fixed.x = FALSE) +layout_m <- matrix(c(NA, "m", NA, + "x", NA, "y", + "c1", NA, NA, + "c2", NA, NA), 4, 3, byrow = TRUE) +p <- semPaths(fit0, + layout = layout_m, + nCharNodes = 0, + exoCov = FALSE, + residuals = FALSE, + sizeMan = 10, + asize = 5, + DoNotPlot = TRUE) +plot(p) +text(label = "(Covariances excluded for readability)", + x = .25, y = -1, + adj = c(.5, .5)) +``` + +# Fitting the Model + +We first fit this multigroup model in +`lavaan::sem()` as usual. There is +no need to label any parameters because +`manymome` will extract the parameters +automatically. + +```{r} +mod_med <- +" +m ~ x + c1 + c2 +y ~ m + x + c1 + c2 +" +fit <- sem(model = mod_med, + data = dat, + fixed.x = FALSE, + group = "group") +``` + +These are the estimates of the paths: + +```{r est_med_mg} +summary(fit, + estimates = TRUE) +``` + +# Generate Bootstrap estimates + +We can use `do_boot()` to generate +the bootstrap estimates first +(see [this article](https://sfcheung.github.io/manymome/articles/do_boot.html) +for an illustration on this function). +The argument `ncores` can be omitted +if the default value is acceptable. + +```{r do_boot_mg, results = FALSE} +fit_boot_out <- do_boot(fit = fit, + R = 5000, + seed = 53253, + ncores = 8) +``` + +# Estimate Indirect Effects + +## Estimate Each Effect by `indirect_effect()` + +The function `indirect_effect()` can be used to as usual +to estimate an indirect effect +and form its bootstrapping or Monte Carlo +confidence interval along a path in a model +that starts with any numeric variable, ends with +any numeric variable, through any numeric variable(s). +A detailed illustration can be found in +[this section](https://sfcheung.github.io/manymome/articles/manymome.html#est_indirect). + +For a multigroup model, the only +difference is that users need to specify +the group using the argument `group`. +It can be set to the group label +as used in `lavaan` (`"Group A"` +or `"Group B"` in this example) +or the group number used in `lavaan` + +```{r do_indirect_mg_A} +ind_gpA <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group A", + boot_ci = TRUE, + boot_out = fit_boot_out) +``` + +This is the output: + +```{r out_med_A} +ind_gpA +``` + +The indirect effect from `x` to `y` through `m` in +`"Group A"` is `r formatC(coef(ind_gpA), digits = 3, format = "f")`, +with a 95% confidence interval of +[`r paste0(formatC(confint(ind_gpA), digits = 3, format = "f"), collapse = ", ")`], +significantly different from zero (*p* < .05). + +We illustrate computing the indirect effect in +`"Group B"`, using group number: + +```{r do_indirect_mg_B} +ind_gpB <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = 2, + boot_ci = TRUE, + boot_out = fit_boot_out) +``` + +This is the output: + +```{r out_med_B} +ind_gpB +``` + +The indirect effect from `x` to `y` through `m` in +`"Group B"` is `r formatC(coef(ind_gpB), digits = 3, format = "f")`, +with a 95% confidence interval of +[`r paste0(formatC(confint(ind_gpA), digits = 3, format = "f"), collapse = ", ")`], +also significantly different from zero (*p* < .05). + +## Treating Group as a "Moderator" + +Instead of computing the indirect effects one-by-one, +we can also treat the grouping variable as +a "moderator" and use +`cond_indirect_effects()` to compute +the indirect effects along a path for +all groups. The detailed illustration +of this function can be found [here](https://sfcheung.github.io/manymome/articles/manymome.html#conditional-indirect-effects). +When use on a multigroup model, +wwe can omit the argument `wlevels`. +The function will automatically identify +all groups in a model, and compute +the indirect effect of the requested +path in each model. + +```{r} +ind <- cond_indirect_effects(x = "x", + y = "y", + m = "m", + fit = fit, + boot_ci = TRUE, + boot_out = fit_boot_out) +``` + +This is the output: + +```{r} +ind +``` + +The results are identical to those computed +individually using `indirect_effect()`. Using +`cond_indirect_effects()` is convenient when +the number of groups is more than two. + +# Compute and Test Between-Group difference + +There are several ways to compute and test +the difference in indirect effects between +two groups. + +## Using the Math Operator `-` + +The math operator `-` (described [here](https://sfcheung.github.io/manymome/reference/math_indirect.html)) +can be used if the indirect effects +have been computed individually by +`indirect_effect()`. We have already +computed the path `x->m->y` before +for the two groups. Let us compute the +differences: + +```{r} +ind_diff <- ind_gpB - ind_gpA +ind_diff +``` + +The difference in indirect effects from `x` to `y` through `m` +is `r formatC(coef(ind_diff), digits = 3, format = "f")`, +with a 95% confidence interval of +[`r paste0(formatC(confint(ind_diff), digits = 3, format = "f"), collapse = ", ")`], +not significantly different from zero (*p* < .05). Therefore, +we conclude that the two groups are not significantly +different on the indirect effects. + +## Using `cond_indirect_diff()` + +If the indirect effects are computed using +`cond_indirect_effects()`, we can use the function +`cond_indirect_diff()` to compute the difference +(described [here](https://sfcheung.github.io/manymome/reference/cond_indirect_diff.html)) +This is more convenient than using the math +operator when the number of groups is +greater than two. + +Let us use `cond_indirect_diff()` on the +output of `cond_indirect_effects()`: + +```{r} +ind_diff2 <- cond_indirect_diff(ind, + from = 1, + to = 2) +ind_diff2 +``` + +The convention is `to` row minus `from` row. +Though may sound not intuitive, the printout +always states clearly which group is subtracted +from which group. The results are identical +to those using the math operator. + +# Advanced Skills + +## Standardized Indirect Effects + +Standardized indirect effects can be computed +as for single-group models (described [here](https://sfcheung.github.io/manymome/articles/manymome.html#standardized-indirect-effect)), +by setting `standardized_x` and/or `standardized_y`. +This is an example: + +```{r} +std_gpA <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group A", + boot_ci = TRUE, + boot_out = fit_boot_out, + standardized_x = TRUE, + standardized_y = TRUE) +std_gpA +``` + +```{r} +std_gpB <- indirect_effect(x = "x", + y = "y", + m = "m", + fit = fit, + group = "Group B", + boot_ci = TRUE, + boot_out = fit_boot_out, + standardized_x = TRUE, + standardized_y = TRUE) +std_gpB +``` + +In `"Group A"`, the (completely) standardized indirect effect +from `x` to `y` through `m` is +`r formatC(coef(std_gpA), digits = 3, format = "f")`. In +`"Group B"`, this effect is +`r formatC(coef(std_gpB), digits = 3, format = "f")`. + +Note that, unlike single-group model, in multigroup models, +the standardized indirect effect in a group uses the +the standard deviations of `x`- and `y`-variables in this group +to do the standardization. Therefore, two groups can have +different unstandardized +effects on a path but similar standardized effects on the +same path, or have similar unstandardized effects on a path +but different standardized effects on this path. This is a +known phenomenon in multigroup structural equation model. + +The difference in the two completely standardized indirect +effects can computed and tested using the math operator `-`: + +```{r} +std_diff <- std_gpB - std_gpA +std_diff +``` + +The difference in completely standardized indirect effects +from `x` to `y` through `m` +is `r formatC(coef(std_diff), digits = 3, format = "f")`, +with a 95% confidence interval of +[`r paste0(formatC(confint(std_diff), digits = 3, format = "f"), collapse = ", ")`], +not significantly different from zero (*p* < .05). Therefore, +we conclude that the two groups are also not significantly +different on the completely standardized indirect effects. + +The function `cond_indirect_effects()` and +`cond_indirect_diff()` can also be used with standardization: + +```{r} +std <- cond_indirect_effects(x = "x", + y = "y", + m = "m", + fit = fit, + boot_ci = TRUE, + boot_out = fit_boot_out, + standardized_x = TRUE, + standardized_y = TRUE) +std +``` + +```{r} +std_diff2 <- cond_indirect_diff(std, + from = 1, + to = 2) +std_diff2 +``` + +The results, again, are identical to those using +`indirect_effect()` and the math operator `-`. + +## Finding All Indirect Paths in a Multigroup Model + +Suppose a model which has more than one, or has many, +indirect paths, is fitted to this dataset: + +```{r} +dat2 <- data_med_complicated_mg +print(head(dat2), digits = 2) +``` + +```{r med_complicated, echo = FALSE, fig.cap = "A Complicated Path Model"} +library(semPlot) +suppressMessages(library(lavaan)) +mod2 <- +" +m11 ~ x1 + x2 +m12 ~ m11 + x1 + x2 +m2 ~ x1 + x2 +y1 ~ m2 + m12 + m11 + x1 + x2 +y2 ~ m2 + m12 + m11 + x1 + x2 +" +fit0 <- sem(mod2, dat2, do.fit = FALSE, fixed.x = FALSE) +layout_m <- matrix(c( NA, "m11", NA, "m12", NA, + "x1", NA, NA, NA, "y1", + "x2", NA, NA, NA, "y2", + NA, NA, "m2", NA, NA), byrow = TRUE, 4, 5) +p <- semPaths(fit0, + residuals = FALSE, + sizeMan = 8, + exoCov = FALSE, + node.width = 1, + edge.label.cex = .50, + label.cex = .75, + style = "ram", + mar = c(5, 5, 5, 5), + layout = layout_m, + DoNotPlot = TRUE) +p$graphAttributes$Edges$color[18:19] <- "white" +plot(p) +text(-1.2, -1, paste("(Covariances and\ncontrol variables", + "omitted for readability)", sep = "\n"), + adj = c(0, .5)) +``` + +We first fit this model in `lavaan`: + +```{r} +mod2 <- +" +m11 ~ x1 + x2 +m12 ~ m11 + x1 + x2 +m2 ~ x1 + x2 +y1 ~ m2 + m12 + m11 + x1 + x2 +y2 ~ m2 + m12 + m11 + x1 + x2 +" +fit2 <- sem(mod2, data = dat2, group = "group") +``` + +The function `all_indirect_paths()` can be used on a +multigroup model to identify indirect paths. The search +can be restricted by setting arguments such as +`x`, `y`, and `exclude` (see the [help page](file:///C:/GitHub/manymome/docs/reference/all_indirect_paths.html) +for details). + +For example, the following identify all paths from `x1` +to `y1`: + +```{r} +paths_x1_y1 <- all_indirect_paths(fit = fit2, + x = "x1", + y = "y1") +``` + +If the `group` argument is not specified, it will automatically +identify all paths in all groups, as shown in the printout: + +```{r} +paths_x1_y1 +``` + +We can then use `many_indirect_effects()` to +compute the indirect effects for all paths identified: + +```{r} +all_ind_x1_y1 <- many_indirect_effects(paths_x1_y1, + fit = fit2) +all_ind_x1_y1 +``` + +Bootstrapping and Monte Carlo confidence intervals can +be formed in the same way they are formed for single-group +models. + +## Computing, Testing, and Plotting Conditional Effects + +Though the focus is on indirect effect, +the main functions in `manymome` can also be used for +computing and plotting the effects along the direct path +between two variables. That is, we can focus on the +moderating effect of group on a direct path. + +For example, in the simple mediation model examined +above, suppose we are interested in the between-group +difference in the path from `m` to `y`, the "b path". +We can first +compute the conditional effect using `cond_indirect_effects()`, +without setting the mediator: + +```{r} +path1 <- cond_indirect_effects(x = "m", + y = "y", + fit = fit, + boot_ci = TRUE, + boot_out = fit_boot_out) +path1 +``` + +The difference between the two paths can be tested +using bootstrapping confidence interval using +`cond_indirect_diff()`: + +```{r} +path1_diff <- cond_indirect_diff(path1, + from = 1, + to = 2) +path1_diff +``` + +Based on bootstrapping, the effect of `m` on `y` +in `"Group B"` is significantly greater than that +in `"Group A"` (*p* < .05). (This is compatible +with the conclusion on the indirect effects because +two groups can have no difference on `ab` even if +they differ on `a` and/or `b`.) + +The `plot` method for the output +of `cond_indirect_effects()` can also be used +for multigroup models: + +```{r manymome_draw_mod_mg, echo = FALSE, fig.cap = "Conditional Effects"} +plot(path1) +``` + +Note that, for multigroup models, the *tumble* +graph proposed by @bodner_tumble_2016 will +always be used. The position of a line for +a group is determined by the descriptive statistics +of this group (e.g, means and SDs). For example, +the line segment of `"Group A"` is far to the right +because `"Group A"` has a larger mean of `m` than +`"Group B"`: + +```{r} +by(dat$m, dat$group, mean) +``` + +It would be misleading if the two lines are plotted on the +same horizontal position, assuming incorrectly that the ranges +of `m` are similar in the two groups. + +The vertical positions of the two lines are similarly determined +by the distributions of other predictors in each +group (the control variables +and `x` in this example). + +Details of the `plot` method can be found +in the [help page](https://sfcheung.github.io/manymome/reference/plot.cond_indirect_effects.html). + +# Final Remarks + +There are some limitations on the support +for multigroup models. Currently, +multiple imputation is not supported. +Moreover, most functions do not (yet) +support multigroup models with +within-group moderators, except for +`cond_indirect()`. We would appreciate +users to report issues discovered when +using [manymome](https://sfcheung.github.io/manymome/index.html) +on multigroup models at [GitHub](https://github.com/sfcheung/manymome/issues). + +# Reference(s) diff --git a/vignettes/articles/references.bib b/vignettes/articles/references.bib index e69de29b..1146363f 100644 --- a/vignettes/articles/references.bib +++ b/vignettes/articles/references.bib @@ -0,0 +1,13 @@ + +@article{bodner_tumble_2016, + title = {Tumble graphs: Avoiding misleading end point extrapolation when graphing interactions from a moderated multiple regression analysis}, + volume = {41}, + issn = {1076-9986}, + doi = {10.3102/1076998616657080}, + pages = {593--604}, + number = {6}, + journaltitle = {Journal of Educational and Behavioral Statistics}, + author = {Bodner, Todd E.}, + date = {2016}, + keywords = {Regression, Moderation (Interaction), Graphs, Tumble Graph}, +}