Nothing
R/postestimate_test_MGD.R
In cSEM: Composite-Based Structural Equation Modeling
Defines functions
testMGD
Documented in
testMGD
#' Tests for multi-group comparisons
#'
#' \lifecycle{stable}
#'
#' This function performs various tests proposed in the context of multigroup analysis.
#'
#' The following tests are implemented:
#' \describe{
#' \item{`.approach_mgd = "Klesel"`: Approach suggested by \insertCite{Klesel2019;textual}{cSEM}}{
#' The model-implied variance-covariance matrix (either indicator
#' (`.type_vcv = "indicator"`) or construct (`.type_vcv = "construct"`))
#' is compared across groups. If the model-implied indicator or construct correlation
#' matrix based on a saturated structural model should be compared, set `.saturated = TRUE`.
#' To measure the distance between the model-implied variance-covariance matrices,
#' the geodesic distance (dG) and the squared Euclidean distance (dL) are used.
#' If more than two groups are compared, the average distance over all groups
#' is used.}
#' \item{`.approach_mgd = "Sarstedt"`: Approach suggested by \insertCite{Sarstedt2011;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' \insertCite{Sarstedt2011;textual}{cSEM} tests if parameter k is equal
#' across all groups. If several parameters are tested simultaneously
#' it is recommended to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value). By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details. Note: the
#' test has some severe shortcomings. Use with caution.}
#' \item{`.approach_mgd = "Chin"`: Approach suggested by \insertCite{Chin2010;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' \insertCite{Chin2010;textual}{cSEM} tests if parameter k is equal
#' between two groups. If more than two groups are tested for equality, parameter
#' k is compared between all pairs of groups. In this case, it is recommended
#' to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value) since this is essentially a multiple testing setup.
#' If several parameters are tested simultaneously, correction is by group
#' and number of parameters. By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details.}
#' \item{`.approach_mgd = "Keil"`: Approach suggested by \insertCite{Keil2000;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' \insertCite{Keil2000;textual}{cSEM} tests if parameter k is equal
#' between two groups. It is assumed, that the standard errors of the coefficients are
#' equal across groups. The calculation of the standard error of the parameter
#' difference is adjusted as proposed by \insertCite{Henseler2009;textual}{cSEM}.
#' If more than two groups are tested for equality, parameter k is compared
#' between all pairs of groups. In this case, it is recommended
#' to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value) since this is essentially a multiple testing setup.
#' If several parameters are tested simultaneously, correction
#' is by group and number of parameters. By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details.}
#' \item{`.approach_mgd = "Nitzl"`: Approach suggested by \insertCite{Nitzl2010;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' Similarly to \insertCite{Keil2000;textual}{cSEM}, a single parameter k is tested
#' for equality between two groups. In contrast to \insertCite{Keil2000;textual}{cSEM},
#' it is assumed, that the standard errors of the coefficients are
#' unequal across groups \insertCite{Sarstedt2011}{cSEM}.
#' If more than two groups are tested for equality, parameter k is compared
#' between all pairs of groups. In this case, it is recommended
#' to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value) since this is essentially a multiple testing setup.
#' If several parameters are tested simultaneously, correction
#' is by group and number of parameters. By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details.}
#' \item{`.approach_mgd = "Henseler"`: Approach suggested by \insertCite{Henseler2007a;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' In doing so, the bootstrap estimates of one parameter are compared across groups.
#' In the literature, this approach is also known as PLS-MGA.
#' Originally, this test was proposed as an one-sided test.
#' In this function we perform a left-sided and a right-sided test
#' to investigate whether a parameter differs across two groups. In doing so, the significance
#' level is divided by 2 and compared to p-value of the left and right-sided test.
#' Moreover, `.approach_p_adjust` is ignored and no overall decision
#' is returned.
#' For a more detailed description, see also \insertCite{Henseler2009;textual}{cSEM}.}
#' \item{`.approach_mgd = "CI_param"`: Approach mentioned in \insertCite{Sarstedt2011;textual}{cSEM}}{
#' This approach is based on the confidence intervals constructed around the
#' parameter estimates of the two groups. If the parameter of one group falls within
#' the confidence interval of the other group and/or vice versa, it can be concluded
#' that there is no group difference.
#' Since it is based on the confidence intervals `.approach_p_adjust` is ignored.}
#' \item{`.approach_mgd = "CI_overlap"`: Approach mentioned in \insertCite{Sarstedt2011;textual}{cSEM}}{
#' This approach is based on the confidence intervals (CIs) constructed around the
#' parameter estimates of the two groups. If the two CIs overlap, it can be concluded
#' that there is no group difference.
#' Since it is based on the confidence intervals `.approach_p_adjust` is ignored.}
#' }
#'
#' Use `.approach_mgd` to choose the approach. By default all approaches are computed
#' (`.approach_mgd = "all"`).
#'
#' For convenience, two types of output are available. See the "Value" section below.
#'
#' By default, approaches based on parameter differences across groups compare
#' all parameters (`.parameters_to_compare = NULL`). To compare only
#' a subset of parameters provide the parameters in [lavaan model syntax][lavaan::model.syntax] just like
#' the model to estimate. Take the simple model:
#'
#' \preformatted{
#' model_to_estimate <- "
#' Structural model
#' eta2 ~ eta1
#' eta3 ~ eta1 + eta2
#'
#' # Each concept os measured by 3 indicators, i.e., modeled as latent variable
#' eta1 =~ y11 + y12 + y13
#' eta2 =~ y21 + y22 + y23
#' eta3 =~ y31 + y32 + y33
#' "
#' }
#' If only the path from eta1 to eta3 and the loadings of eta1 are to be compared
#' across groups, write:
#' \preformatted{
#'to_compare <- "
#' Structural parameters to compare
#' eta3 ~ eta1
#'
#' # Loadings to compare
#' eta1 =~ y11 + y12 + y13
#' "
#' }
#' Note that the "model" provided to `.parameters_to_compare`
#' does not need to be an estimable model!
#'
#' Note also that compared to all other functions in \pkg{cSEM} using the argument,
#' `.handle_inadmissibles` defaults to `"replace"` to accommodate the Sarstedt et al. (2011) approach.
#'
#' Argument `.R_permuation` is ignored for the `"Nitzl"` and the `"Keil"` approach.
#' `.R_bootstrap` is ignored if `.object` already contains resamples,
#' i.e. has class `cSEMResults_resampled` and if only the `"Klesel"` or the `"Chin"`
#' approach are used.
#'
#' The argument `.saturated` is used by `"Klesel"` only. If `.saturated = TRUE`
#' the original structural model is ignored and replaced by a saturated model,
#' i.e. a model in which all constructs are allowed to correlate freely.
#' This is useful to test differences in the measurement models between groups
#' in isolation.
#'
#' @usage testMGD(
#' .object = NULL,
#' .alpha = 0.05,
#' .approach_p_adjust = "none",
#' .approach_mgd = c("all", "Klesel", "Chin", "Sarstedt",
#' "Keil", "Nitzl", "Henseler", "CI_para","CI_overlap"),
#' .output_type = c("complete", "structured"),
#' .parameters_to_compare = NULL,
#' .eval_plan = c("sequential", "multicore", "multisession"),
#' .handle_inadmissibles = c("replace", "drop", "ignore"),
#' .R_permutation = 499,
#' .R_bootstrap = 499,
#' .saturated = FALSE,
#' .seed = NULL,
#' .type_ci = "CI_percentile",
#' .type_vcv = c("indicator", "construct"),
#' .verbose = TRUE
#' )
#'
#' @inheritParams csem_arguments
#' @param .handle_inadmissibles Character string. How should inadmissible results
#' be treated? One of "*drop*", "*ignore*", or "*replace*". If "*drop*", all
#' replications/resamples yielding an inadmissible result will be dropped
#' (i.e. the number of results returned will potentially be less than `.R`).
#' For "*ignore*" all results are returned even if all or some of the replications
#' yielded inadmissible results (i.e. number of results returned is equal to `.R`).
#' For "*replace*" resampling continues until there are exactly `.R` admissible solutions.
#' Defaults to "*replace*" to accommodate all approaches.
#'
#' @return If `.output_type = "complete"` a list of class `cSEMTestMGD`. Technically, `cSEMTestMGD` is a
#' named list containing the following list elements:
#'
#' \describe{
#' \item{`$Information`}{Additional information.}
#' \item{`$Klesel`}{A list with elements, `Test_statistic`, `P_value`, and `Decision`}
#' \item{`$Chin`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Sarstedt`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Keil`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Nitzl`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Henseler`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$CI_para`}{A list with elements, `Decision`, and `Decision_overall`}
#' \item{`$CI_overlap`}{A list with elements, `Decision`, and `Decision_overall`}
#' }
#'
#' If `.output_type = "structured"` a tibble (data frame) with the following columns
#' is returned.
#'
#' \describe{
#' \item{`Test`}{The name of the test.}
#' \item{`Comparision`}{The parameter that was compared across groups. If "overall"
#' the overall fit of the model was compared.}
#' \item{`alpha%`}{The test decision for a given "alpha" level. If `TRUE` the null
#' hypotheses was rejected; if FALSE it was not rejected.}
#' \item{`p-value_correction`}{The p-value correction.}
#' \item{`CI_type`}{Only for the "CI_para" and the "CI_overlap" test. Which confidence interval was used.}
#' \item{`Distance_metric`}{Only for Test = "Klesel". Which distance metric was used.}
#' }
#'
#' @references
#' \insertAllCited{}
#'
#' @seealso [csem()], [cSEMResults], [testMICOM()], [testOMF()]
#'
#' @example inst/examples/example_testMGD.R
#'
#' @export
testMGD <- function(
.object = NULL,
.alpha = 0.05,
.approach_p_adjust = "none",
.approach_mgd = c("all", "Klesel", "Chin", "Sarstedt",
"Keil", "Nitzl","Henseler", "CI_para","CI_overlap"),
.output_type = c("complete", "structured"),
.parameters_to_compare = NULL,
.eval_plan = c("sequential", "multicore", "multisession"),
.handle_inadmissibles = c("replace", "drop", "ignore"),
.R_permutation = 499,
.R_bootstrap = 499,
.saturated = FALSE,
.seed = NULL,
.type_ci = "CI_percentile",
.type_vcv = c("indicator", "construct"),
.verbose = TRUE
){
### Checks and errors ========================================================
# Check .approach_mgd argument choices
diff <- setdiff(.approach_mgd, args_default(TRUE)$.approach_mgd)
if(length(diff) != 0) {
stop2(
"The following error occured in the testMGD() function:\n",
"Unknown approach: ", paste0(diff, collapse = ", "), ".",
" Possible choices are: ",
paste0(args_default(TRUE)$.approach_mgd, collapse = ", "))
}
## Match arguments
.approach_mgd <- match.arg(.approach_mgd, several.ok = TRUE)
.handle_inadmissibles <- match.arg(.handle_inadmissibles) # no reference to
# args_default, because args_default()$.handle_inadmissibles
# has "drop" as default, but testMGD hast "replace".
.type_vcv <- match.arg(.type_vcv)
.output_type <- match.arg(.output_type)
# if(!all(.type_ci %in% c("CI_standard_z", "CI_standard_t", "CI_percentile",
# "CI_basic", "CI_bc", "CI_bca") )){
if(!all(.type_ci %in% args_default(.choices = TRUE)$.type_ci )){
stop2("The specified confidence interval in .type.ci is not valid.\n",
"Please choose one of the following: CI_standard_z, CI_standard_t,\n",
"CI_percentile, CI_basic, CI_bc, CI_bca.")
}
# Check if at least two groups are present
if(!inherits(.object, "cSEMResults_multi")) {
stop2(
"The following error occured in the testMGD() function:\n",
"At least two groups required."
)
}
# Sarstedt et al. (2011) is not allowed to be used in combination with
# .handle_inadmissibles == "drop as permutation test statistic are be dropped
# only because of estimations based on the permutated dataset are dropped.
if(any(.approach_mgd %in% c("all", "Sarstedt")) & .handle_inadmissibles == "drop"){
stop2(
"The following error occured in the testMGD() function:\n",
"Approach `'Sarstedt'` not supported if `.handle_inadmissibles == 'drop'`")
}
# If a nonlinear model is used Klesel et al. approach cannot be used as
# it is not clear how to calculate the model-implied VCV
# Extract model type information
if(inherits(.object[[1]], "cSEMResults_2ndorder")){
model_type <- .object[[1]]$Second_stage$Information$Model$model_type
} else {
model_type <- .object[[1]]$Information$Model$model_type
}
if(any(.approach_mgd %in% c("all", "Klesel")) & model_type == "Nonlinear"){
stop2("The following error occured in the testMGD() function:\n",
"The approach suggested by Klesel et al. (2019) cannot be applied",
" to nonlinear models as cSEM currently can not calculate",
" the model-implied VCV matrix for such models.\n",
"Consider setting `.approach_mgd = c('Chin', 'Sarstedt')`")
}
# Check if any of the group estimates are inadmissible
if(sum(unlist(verify(.object))) != 0) {
warning2(
"The following warning occured in the testMGD() function:\n",
"Initial estimation results for at least one group are inadmissible.\n",
"See `verify(.object)` for details.")
}
# Check if data for different groups is identical
if(TRUE %in% lapply(utils::combn(.object, 2, simplify = FALSE),
function(x){ identical(x[[1]], x[[2]])})){
warning2(
"The following warning occured in the testMGD() function:\n",
"At least two groups are identical. Results may not be meaningful.")
}
# If Henseler, CI_para or CI_overlap are used with
# adjustment of p-value different than "none" return warning
if(any(.approach_mgd %in% c("all", "Henseler","CI_para","CI_overlap")) & !all(.approach_p_adjust %in% "none")){
warning2(
"The following warning occured in the testMGD() function:\n",
"Currently, there is no p-value adjustment possible for the approach suggested by\n",
"Henseler (2007), CI_para, and CI_overlap. Adjustment is ignored for these approaches."
)
}
if(.verbose) {
cat(rule2("Several tests for multi-group comparisons",
type = 3), "\n\n")
}
# Order significance levels
.alpha <- .alpha[order(.alpha)]
# Get the names of the parameters to be compared.
names_all_param <- getParameterNames(.object, .model = .parameters_to_compare)
# getParameterNames() returns also measurement error and indicator correlations
# Currently they cannot be handeled, therefore an error is returned.
# FOR FUTURE: MEasurement errors and indicator correlations can be added to the parameters allowed in the comparison
if(!is.null(names_all_param$names_cor_measurement_error)|!is.null(names_all_param$names_cor_indicator)){
stop2("The following error occured in the testMGD() function:\n",
"Currenlty it is not allowed to compare measurement error covariance",
" and/or indicator covariances across groups.")
}
names_param <- unlist(names_all_param)
## Calculation of the test statistics========================================
teststat <- list()
## Klesel et al. (2019) ------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Klesel"))) {
## Get the model-implied VCV
fit <- fit(.object = .object,
.saturated = .saturated,
.type_vcv = .type_vcv)
## Compute test statistic
temp <- c(
"dG" = calculateDistance(.matrices = fit, .distance = "geodesic"),
"dL" = calculateDistance(.matrices = fit, .distance = "squared_euclidean")
)
## Save test statistic
teststat[["Klesel"]] <- temp
}
## Chin & Dibbern (2010) -----------------------------------------------------
if(any(.approach_mgd %in% c("all", "Chin"))) {
## Compute and save test statistic
teststat[["Chin"]] <- calculateParameterDifference(.object = .object,
.model = .parameters_to_compare)
}
## Sarstedt et al. (2011), Keil et al. (2000), Nitzl (2010), Henseler (2007) -----
# All these approaches require bootstrap
if(any(.approach_mgd %in% c("all", "Sarstedt", "Keil", "Nitzl", "Henseler", "CI_para","CI_overlap"))) {
# Check if .object already contains resamples; if not, run bootstrap
if(!inherits(.object, "cSEMResults_resampled")) {
if(.verbose) {
cat("Bootstrap cSEMResults objects ...\n\n")
}
.object <- resamplecSEMResults(
.object = .object,
.resample_method = "bootstrap",
.handle_inadmissibles = .handle_inadmissibles,
.R = .R_bootstrap,
.seed = .seed,
.eval_plan = .eval_plan
)
}
## Combine bootstrap results in one matrix
bootstrap_results <- lapply(.object, function(y) {
if(inherits(.object, "cSEMResults_2ndorder")) {
x <- y$Second_stage$Information$Resamples$Estimates$Estimates1
nobs <- nrow(y$First_stage$Information$Data)
} else {
x <- y$Estimates$Estimates_resample$Estimates1
nobs <- nrow(y$Information$Data)
}
path_resamples <- x$Path_estimates$Resampled
loading_resamples <- x$Loading_estimates$Resampled
weight_resamples <- x$Weight_estimates$Resampled
# Rename
cor_exo_cons_resamples <- x$Exo_construct_correlation$Resampled
n <- nrow(path_resamples)
# Calculation of the bootstrap SEs
ses <- infer(.object=y,.quantity = "sd")
path_se <- ses$Path_estimates$sd
loading_se <- ses$Loading_estimates$sd
weight_se <- ses$Weight_estimates$sd
cor_exo_cons_se <- ses$Exo_construct_correlation$sd
# Calculation of the bias
bias <- infer(.object=y,.quantity = "bias")
path_bias <- bias$Path_estimates$bias
loading_bias <- bias$Loading_estimates$bias
weight_bias <- bias$Weight_estimates$bias
cor_exo_cons_bias <- bias$Exo_construct_correlation$bias
# Structure output
out<-list(
"n" = n,
"nObs" = nobs,
"para_all" = cbind(path_resamples,loading_resamples,weight_resamples, cor_exo_cons_resamples),
"ses_all" = c(path_se, loading_se, weight_se, cor_exo_cons_se),
"bias_all" = c(path_bias, loading_bias, weight_bias, cor_exo_cons_bias)
)
# Approaches based on CIs comparison ----
if(any(.approach_mgd %in% c("all", "CI_para", "CI_overlap"))){
diff <- setdiff(.type_ci, args_default(TRUE)$.type_ci)
if(length(diff) != 0) {
stop2(
"The following error occured in the testMGD() function:\n",
"Unknown approach: ", paste0(diff, collapse = ", "), ".",
" Possible choices are: ",
paste0(args_default(TRUE)$.type_ci, collapse = ", "))
}
# calculation of the CIs
cis <- infer(.object=y, .quantity=.type_ci, .alpha = .alpha)
cis_temp <- purrr::transpose(cis)
cis_ret<-lapply(cis_temp,function(x){
path_ci <- x$Path_estimates
loading_ci <- x$Loading_estimates
weight_ci <- x$Weight_estimates
cor_exo_cons_ci <- x$Exo_construct_correlation
# deliberately name them xyz_estimates to be able to select them later via names.
list(path_estimates = path_ci, loading_estimates = loading_ci,
weight_estimates = weight_ci, cor_exo_cons_estimates = cor_exo_cons_ci)
})
names(cis_ret) <- names(cis_temp)
out[["ci_all"]] <- cis_ret
}
# Return output
return(out)
})
## Keil, Nitzl, Henseler
if(any(.approach_mgd %in% c("all", "Nitzl", "Keil","Henseler"))) {
# Calculate the difference for one parameter between two groups
# Although Henseler approach does not compute a test statistic,
# The difference is calculated to create a dummy test statistic list
diff_para_Keil <- diff_para_Nitzl <- diff_para_Henseler<- calculateParameterDifference(
.object = .object,
.model = .parameters_to_compare
)
if(any(.approach_mgd %in% c("all","Henseler"))) {
# Create dummy test statistic list containing NAs
temp <- rep(NA,length(unlist(diff_para_Henseler)))
teststat[["Henseler"]] <-relist(flesh = temp,skeleton = diff_para_Henseler)
}
# Build list that contains pairs of objects.
object_permu <- utils::combn(.object, 2, simplify = FALSE)
names(object_permu) <- sapply(object_permu, function(x) paste0(names(x)[1], '_', names(x)[2]))
# Approach suggested by Keil
if(any(.approach_mgd %in% c("all", "Keil"))) {
teststat_Keil <- lapply(names(object_permu), function(x) {
diff <- diff_para_Keil[[x]]
ses1 <- bootstrap_results[[names(object_permu[[x]][1])]]$ses_all
ses2 <- bootstrap_results[[names(object_permu[[x]][2])]]$ses_all
# Get sample size per group
n1<-bootstrap_results[[names(object_permu[[x]][1])]]$nObs
n2<-bootstrap_results[[names(object_permu[[x]][2])]]$nObs
# Calculation of the SE of the parameter difference as proposed by
# Henseler (2007a), Henseler et al. (2009)
# Assumption of equal variances
ses_total <- sqrt((n1-1)^2/(n1+n2-2)*ses1^2 +
(n2-1)^2/(n1+n2-2)*ses2^2)*sqrt(1/n1+1/n2)
test_stat <- diff/ses_total[names(diff)]
list("teststat" = test_stat, "df" = n1 + n2 - 2)
})
names(teststat_Keil) <- names(object_permu)
teststat[["Keil"]] <- teststat_Keil
}
# Approach suggested by Nitzl (2010)
if(any(.approach_mgd %in% c("all", "Nitzl"))) {
teststat_Nitzl <- lapply(names(object_permu), function(x){
diff <- diff_para_Nitzl[[x]]
ses1 <- bootstrap_results[[names(object_permu[[x]][1])]]$ses_all
ses2 <- bootstrap_results[[names(object_permu[[x]][2])]]$ses_all
# Get sample size per group
n1<-bootstrap_results[[names(object_permu[[x]][1])]]$nObs
n2<-bootstrap_results[[names(object_permu[[x]][2])]]$nObs
# Calculation of the SE of the parameter difference as proposed by
# Henseler (2007a), Henseler et al. (2009)
# Assumption of unequal variances
ses_total <- sqrt((n1-1)/(n1)*ses1^2 +
(n2-1)/(n2)*ses2^2)
test_stat <- diff/ses_total[names(diff)]
# calculation of the degrees of freedom
numerator <- ((n1-1)/n1*ses1^2+(n2-1)/n2*ses2^2)^2
denominator <- (n1-1)/n1^2*ses1^4+(n2-1)/n2^2*ses2^4
df <- round(numerator/denominator-2)
df <- df[names(diff)]
list("teststat" = test_stat, "df" = df)
})
names(teststat_Nitzl) <- names(object_permu)
teststat[["Nitzl"]] <- teststat_Nitzl
}
}
## Sarstedt et al. approach
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
## Transpose
ll <- purrr::transpose(bootstrap_results)
group_id <- rep(1:length(.object), unlist(ll$n))
# all_comb contains all parameter estimate that could potentially be compared
# plus an id column indicating the group adherance of each row.
all_comb <- cbind(do.call(rbind,ll$para_all),
"group_id" = group_id)
## Select relevant columns
all_comb <- all_comb[, c(names_param, "group_id")]
## Add test statistic Sarstedt
teststat[["Sarstedt"]] <- calculateFR(.resample_sarstedt = all_comb)
} # END approach_mgd %in% c("all", "Sarstedt")
} # END .approach_mgd %in% c("all", "Sarstedt", "Keil", "Nitzl", "Henseler")
### Permutation ==============================================================
## Preparation
# Put data of each groups in a list and combine
if(inherits(.object, "cSEMResults_2ndorder")) {
# Data is saved in the first stage
X_all_list <- lapply(.object, function(x) x$First_stage$Information$Data)
# Collect initial arguments (from the first object, but could be any other)
arguments <- .object[[1]]$Second_stage$Information$Arguments_original
} else {
X_all_list <- lapply(.object, function(x) x$Information$Data)
# Collect initial arguments (from the first object, but could be any other)
arguments <- .object[[1]]$Information$Arguments
}
X_all <- do.call(rbind, X_all_list)
# Create a vector "id" to be used to randomly select groups (permutate) and
# set id as an argument in order to identify the groups.
id <- rep(1:length(X_all_list), sapply(X_all_list, nrow))
arguments[[".id"]] <- "id"
# Permutation is only performed for approaches that require it
if(any(.approach_mgd %in% c("all", "Klesel", "Chin", "Sarstedt"))) {
# Start progress bar if required
if(.verbose){
cat("Permutation ...\n\n")
# pb <- txtProgressBar(min = 0, max = .R_permutation, style = 3)
}
# Save old seed and restore on exit! This is important since users may have
# set a seed before, in which case the global seed would be
# overwritten if not explicitly restored
old_seed <- .Random.seed
on.exit({.Random.seed <<- old_seed})
## Create seed if not already set
if(is.null(.seed)) {
set.seed(seed = NULL)
# Note (08.12.2019): Its crucial to call set.seed(seed = NULL) before
# drawing a random seed out of .Random.seed. If set.seed(seed = NULL) is not
# called sample(.Random.seed, 1) would result in the same random seed as
# long as .Random.seed remains unchanged. By resetting the seed we make
# sure that sample draws a different element everytime it is called.
.seed <- sample(.Random.seed, 1)
}
## Set seed
set.seed(.seed)
## Calculate reference distribution
ref_dist <- list()
n_inadmissibles <- 0
counter <- 0
progressr::with_progress({
progress_bar_csem <- progressr::progressor(along = 1:.R_permutation)
repeat{
# Counter
counter <- counter + 1
progress_bar_csem(message = sprintf("Permutation run = %g", counter))
# Permutate data
X_temp <- cbind(X_all, id = sample(id))
# Replace the old dataset by the new permutated dataset
arguments[[".data"]] <- X_temp
# Estimate model
Est_temp <- do.call(csem, arguments)
# Check status
status_code <- sum(unlist(verify(Est_temp)))
# Distinguish depending on how inadmissibles should be handled
if(status_code == 0 | (status_code != 0 & .handle_inadmissibles == "ignore")) {
# Compute if status is ok or .handle inadmissibles = "ignore" AND the status is
# not ok
### Calculation of the test statistic for each resample ==================
teststat_permutation <- list()
## Klesel et al. (2019) --------------------------------------------------
if(any(.approach_mgd %in% c("all", "Klesel"))) {
## Get the model-implied VCV
fit_temp <- fit(Est_temp, .saturated = .saturated, .type_vcv = .type_vcv)
## Compute test statistic
temp <- c(
"dG" = calculateDistance(.matrices = fit_temp, .distance = "geodesic"),
"dL" = calculateDistance(.matrices = fit_temp, .distance = "squared_euclidean")
)
## Save test statistic
teststat_permutation[["Klesel"]] <- temp
}
## Chin & Dibbern (2010) -------------------------------------------------
if(any(.approach_mgd %in% c("all", "Chin"))) {
## Compute and save test statistic
teststat_permutation[["Chin"]] <- calculateParameterDifference(
.object = Est_temp,
.model = .parameters_to_compare)
}
## Sarstedt et al. (2011) ------------------------------------------------
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
# Permutation of the bootstrap parameter estimates
all_comb_permutation <- all_comb
all_comb_permutation[ , "group_id"] <- sample(group_id)
teststat_permutation[["Sarstedt"]] <- calculateFR(all_comb_permutation)
}
ref_dist[[counter]] <- teststat_permutation
} else if(status_code != 0 & .handle_inadmissibles == "drop") {
# Set list element to zero if status is not okay and .handle_inadmissibles == "drop"
ref_dist[[counter]] <- NA
} else {# status is not ok and .handle_inadmissibles == "replace"
# Reset counter and raise number of inadmissibles by 1
counter <- counter - 1
n_inadmissibles <- n_inadmissibles + 1
}
# Update progres bar
# if(.verbose){
# setTxtProgressBar(pb, counter)
# }
# Break repeat loop if .R results have been created.
if(length(ref_dist) == .R_permutation) {
break
} else if(counter + n_inadmissibles == 10000) {
# Stop if 10000 runs did not result in insufficient admissible results
stop("Not enough admissible result.", call. = FALSE)
}
} # END repeat
}) # END with_progress
# close progress bar
# if(.verbose){
# close(pb)
# }
### Postprocessing ===========================================================
# Delete potential NA's
ref_dist1 <- Filter(Negate(anyNA), ref_dist)
}
# # Order significance levels
# .alpha <- .alpha[order(.alpha)]
## Klesel et al. (2019) ------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Klesel"))) {
# Collect permutation results and combine
ref_dist_Klesel <- lapply(ref_dist1, function(x) x$Klesel)
ref_dist_matrix_Klesel <- do.call(cbind, ref_dist_Klesel)
# Extract test statistic
teststat_Klesel <- teststat$Klesel
# Calculation of p-values
pvalue_Klesel <- rowMeans(ref_dist_matrix_Klesel >= teststat_Klesel)
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Klesel <- lapply(.alpha, function(x) {
pvalue_Klesel > x
})
names(decision_Klesel) <- paste0(.alpha * 100, "%")
}
## Chin & Dibbern (2010) -----------------------------------------------------
if(any(.approach_mgd %in% c("all", "Chin"))) {
# Extract test statistic
teststat_Chin <- teststat$Chin
# Create list with matrices containing the reference distribution
# of the parameter differences
ref_dist_Chin <- lapply(ref_dist1, function(x) x$Chin)
# Transpose
ref_dist_Chin_temp <- purrr::transpose(ref_dist_Chin)
names(ref_dist_Chin_temp) <- names(teststat_Chin)
ref_dist_matrices_Chin <- lapply(ref_dist_Chin_temp, function(x) {
temp <- do.call(cbind, x)
temp_ind <- stats::complete.cases(temp)
temp[temp_ind, ,drop = FALSE]
})
# Calculation of the p-values
pvalue_Chin <- lapply(1:length(ref_dist_matrices_Chin), function(x) {
# Share of values above the positive test statistic
rowMeans(ref_dist_matrices_Chin[[x]] >= abs(teststat_Chin[[x]])) +
# share of values of the reference distribution below the negative test statistic
rowMeans(ref_dist_matrices_Chin[[x]] <= (-abs(teststat_Chin[[x]])))
})
names(pvalue_Chin) <- names(ref_dist_matrices_Chin)
# Adjust p-values
padjusted_Chin <- lapply(as.list(.approach_p_adjust), function(x){
# It is important to unlist the pvalues as pAdjust needs to now how many p-values
# there are to do a proper adjustment
pvector <- stats::p.adjust(unlist(pvalue_Chin),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Chin)
})
names(padjusted_Chin) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Chin <- lapply(padjusted_Chin, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
group_comp > alpha
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# Overall decision, i.e., was any of the test belonging to one significance level
# and one p value adjustment rejected
decision_overall_Chin <- lapply(decision_Chin, function(decision_Chin_list){
lapply(decision_Chin_list,function(x){
all(unlist(x))
})
})
}
## Sarstedt et al. (2011) ----------------------------------------------------
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
# Extract test statistic
teststat_Sarstedt <- teststat$Sarstedt
# Collect permuation results and combine to matrix
ref_dist_Sarstedt <- lapply(ref_dist1, function(x) x$Sarstedt)
ref_dist_matrix_Sarstedt <- do.call(cbind, ref_dist_Sarstedt)
# Calculation of the p-value
pvalue_Sarstedt <- rowMeans(ref_dist_matrix_Sarstedt >= teststat_Sarstedt)
# Adjust pvalues:
padjusted_Sarstedt<- lapply(as.list(.approach_p_adjust), function(x){
pvector <- stats::p.adjust(pvalue_Sarstedt, method = x)
})
names(padjusted_Sarstedt) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Sarstedt <- lapply(padjusted_Sarstedt,function(p_value){
temp <- lapply(.alpha, function(alpha){
p_value > alpha
})
names(temp) <- paste(.alpha*100,"%",sep= '')
temp
})
# Decision overall
decision_overall_Sarstedt <- lapply(decision_Sarstedt, function(x){#over p-value adjustments
lapply(x, function(xx){ #over different significant levels
all(xx)
})
})
}
## Keil et al. (2000) --------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Keil"))){
# Calculate p-values
pvalue_Keil <- lapply(teststat$Keil,function(x){
p_value <- 2*(1-pt(abs(x$teststat),df = x$df))
p_value
})
# Adjust p-values in case of more than one comparison
padjusted_Keil<- lapply(as.list(.approach_p_adjust), function(x){
pvector <- stats::p.adjust(unlist(pvalue_Keil),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Keil)
})
names(padjusted_Keil) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Keil <- lapply(padjusted_Keil, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
group_comp > alpha
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# One rejection leads to overall rejection
decision_overall_Keil <- lapply(decision_Keil, function(decision_Keil_list){
lapply(decision_Keil_list,function(x){
all(unlist(x))
})
})
}
## Nitzl et al. (2010) -------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Nitzl"))){
# Calculate p-values
pvalue_Nitzl <- lapply(teststat$Nitzl,function(x){
p_value <- 2*(1-pt(abs(x$teststat),df = x$df))
p_value
})
# Adjust p-values in case of more than one comparison
padjusted_Nitzl<- lapply(as.list(.approach_p_adjust), function(x){
pvector <- stats::p.adjust(unlist(pvalue_Nitzl),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Nitzl)
})
names(padjusted_Nitzl) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Nitzl <- lapply(padjusted_Nitzl, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
group_comp > alpha
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# One rejection leads to overall rejection
decision_overall_Nitzl <- lapply(decision_Nitzl, function(decision_Nitzl_list){
lapply(decision_Nitzl_list,function(x){
all(unlist(x))
})
})
}
## Henseler (2009) approach -------------------------------------------------
if(any(.approach_mgd %in% c("all", "Henseler"))) {
# Pseudo teststat for convenience
teststat_Henseler <- teststat$Henseler
# center the bootstrap sample Sarstedt et al. (2011) Eq. 4
ll_centered <- lapply(bootstrap_results, function(x){
# Substract from each row of para_all the corresponding bias
t(apply(x$para_all,1,function(row){
row-x$bias_all
}))
})
names(ll_centered) <- names(bootstrap_results)
# Create pairs which should be compared
pairs_centered <- utils::combn(ll_centered, 2, simplify = FALSE)
names(pairs_centered) <- sapply(pairs_centered, function(x) paste0(names(x)[1], '_', names(x)[2]))
# Calculation of the probability
pvalue_Henseler <- lapply(pairs_centered,function(x){
calculatePr(.resample_centered = x,
.parameters_to_compare = names_param)
})
# Adjust p-value in case of multiple comparisons
# Adjusting p-values is not straight forward.
# First it is a one-sided test, so we need to know the hypothesis
# Just flipping the p-value is not really an option as it might cause problems
# in situation where the order of the p-values is required for the correction
# Therefore only the "none"method is applied
padjusted_Henseler<- lapply(as.list("none"), function(x){
pvector <- stats::p.adjust(unlist(pvalue_Henseler),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Henseler)
})
names(padjusted_Henseler) <- "none"
# Decision is made:
# The probability is compared to alpha and 1-alpha
# In doing so, we assess both hypotheses theta^1 <= theta^2 and
# theta^1 => theta^2
# TRUE = p-value > alpha & p-value < 1-alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Henseler <- lapply(padjusted_Henseler, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
# The alpha is divided by two to mimic a two-sided test.
group_comp > alpha/2 & group_comp < 1- alpha/2
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# Overall decision; if one coefficient is significanty different across groups
# it is rejected.
decision_overall_Henseler <- lapply(decision_Henseler, function(decision_Henseler_list){
lapply(decision_Henseler_list,function(x){
all(unlist(x))
})
})
}
# Comparison via confidence intervals ------------------
if(any(.approach_mgd %in% c("all", "CI_para", "CI_overlap"))) {
# Select CIs from the bootstrap_results
cis <- lapply(bootstrap_results,function(x){
x$ci_all
})
# Make group pairs of CIs
cis_comp <- utils::combn(cis, 2, simplify = FALSE)
names(cis_comp) <- sapply(cis_comp, function(x) paste0(names(x)[1], '_', names(x)[2]))
# Select the relevant parameters per group and make group pairs
# This is required fo both CI_para and CI_overlap
# For the latter it is used to select the appropriate parameters
param_per_group <- getRelevantParameters(.object = .object,
.model = .parameters_to_compare)
param_per_group <- purrr::transpose(param_per_group)
param_comp <- utils::combn(param_per_group, 2, simplify = FALSE)
names(param_comp) <- sapply(param_comp,
function(x)
paste0(names(x)[1], '_', names(x)[2]))
# CI_para -------------------------------------------
if(any(.approach_mgd %in% c("all", "CI_para"))) {
# Investigate whether the estimate of one group is part of the CI of another group
decision_ci_para <- lapply(.alpha, function(alpha) {
# lapply over the comparisons
tttt <- lapply(names(cis_comp), function(comp) {
# lapply over the type of CI for the first group
# As group1 and group2 have the same structure the selection does not matter
tt = lapply(names(cis_comp[[comp]][[1]]), function(interval_type){
# lapply over the different parameters
ttt <- lapply(names(cis_comp[[comp]][[1]][[interval_type]]),
function(param) {
lb <- paste0(100 * (1 - alpha), "%L")
ub <- paste0(100 * (1 - alpha), "%U")
# Estimate first group
temp_para1 <- param_comp[[comp]][[1]][[param]]
# CIs first group
temp_cis1 <- cis_comp[[comp]][[1]][[interval_type]][[param]]
temp_cis_selected1 <- temp_cis1[, names(temp_para1), drop = FALSE]
lb_temp1 <-temp_cis_selected1[lb,]
names(lb_temp1) <- colnames(temp_cis_selected1[lb,,drop=FALSE])
ub_temp1 <-temp_cis_selected1[ub,]
names(ub_temp1) <- colnames(temp_cis_selected1[ub,,drop=FALSE])
# Estimate second group
temp_para2 <- param_comp[[comp]][[2]][[param]]
# CIs second group
temp_cis2 <- cis_comp[[comp]][[2]][[interval_type]][[param]]
temp_cis_selected2 <- temp_cis2[, names(temp_para2), drop = FALSE]
lb_temp2 <-temp_cis_selected2[lb,]
names(lb_temp2) <- colnames(temp_cis_selected2[lb,,drop=FALSE])
ub_temp2 <-temp_cis_selected2[ub,]
names(ub_temp2) <- colnames(temp_cis_selected2[ub,,drop=FALSE])
# browser()
# # check whether parameter estimates of one group
# falls within boundaries of the CIs of the other group
# TRUE parameter estimate of the other group falls within the CI
# => no group difference
# Otherwise, FALSE => group difference
decision<-
(lb_temp2 <= temp_para1 &
ub_temp2 >= temp_para1) |
(lb_temp1 <= temp_para2 &
ub_temp1 >= temp_para2)
# Structure output
out=data.frame("Estimate"=temp_para1,"lb"=lb_temp2,
"ub"=ub_temp2,"Estimate"=temp_para2,"lb"=lb_temp1,
"ub"=ub_temp1,"decision"=decision)
# that can be solved more elegant
# If there is no parameter to compare the data.frame has zero rows
# If the if is removed there are problems as the number of columns does not match
# the length of the name vector
if(nrow(out)!=0){
colnames(out)=c(paste0("Est_",names(cis_comp[[comp]])[1]),
paste0("lb_",names(cis_comp[[comp]])[2]),
paste0("ub_",names(cis_comp[[comp]])[2]),
paste0("Est_",names(cis_comp[[comp]])[2]),
paste0("lb_",names(cis_comp[[comp]])[1]),
paste0("ub_",names(cis_comp[[comp]])[1]),
"Decision")
}
# Make rownames explicit
out$Name <- rownames(out)
out <- out[, c(ncol(out), 1:(ncol(out) -1))]
rownames(out) <- NULL
out
})
do.call(rbind,ttt)
})
names(tt) = names(cis_comp[[comp]][[1]])
tt
})
names(tttt) <- names(cis_comp)
tttt
})
names(decision_ci_para) <- paste0((1 - .alpha) * 100, "%")
# Decision overall
decision_overall_ci_para<-lapply(names(decision_ci_para), function(alpha) {
t <- lapply(names(decision_ci_para[[alpha]]), function(comp) {
tt <- lapply(names(decision_ci_para[[alpha]][[comp]]),function(interval_type){
# Check whether there is one FALSE among the decisions per comparison
all(decision_ci_para[[alpha]][[comp]][[interval_type]][,"Decision"])
})
names(tt) <- names(decision_ci_para[[alpha]][[comp]])
tt
})
names(t) <- names(decision_ci_para[[alpha]])
# Check whether there is one FALSE among the comparions' decisions
temp <- purrr::transpose(t)
overall_dec<-lapply(names(temp),function(x){
all(unlist(temp[[x]]))
})
names(overall_dec)<-names(temp)
overall_dec
})
names(decision_overall_ci_para) <- names(decision_ci_para)
}#end if .approach_mgd == CI_para
# CI_overlap ---------------------------------------------
if(any(.approach_mgd %in% c("all", "CI_overlap"))) {
decision_ci_overlap <- lapply(.alpha, function(alpha) {
# lapply over the group pairs
ttt <- lapply(names(cis_comp), function(comp) {
# lapply over the type of CI for the first group
# As group1 and group2 have the same structure the selection does not matter
t <- lapply(names(cis_comp[[comp]][[1]]), function(interval_type) {
# lapply over the parameters that are compared
tt <- lapply(names(cis_comp[[comp]][[1]][[interval_type]]), function(param) {
lb <- paste0(100 * (1 - alpha), "%L")
ub <- paste0(100 * (1 - alpha), "%U")
# It does not matter which group is chosen to select the relevant parameters
# as in all groups the relevant parameters are the same
para_rel <- names(param_comp[[comp]][[1]][[param]])
# Select lower and upper bound of the CIs for the two groups
lb1 <- cis_comp[[comp]][[1]][[interval_type]][[param]][lb,para_rel ]
ub1 <- cis_comp[[comp]][[1]][[interval_type]][[param]][ub,para_rel ]
lb2 <- cis_comp[[comp]][[2]][[interval_type]][[param]][lb,para_rel ]
ub2 <- cis_comp[[comp]][[2]][[interval_type]][[param]][ub,para_rel ]
# Check whether the boundaries of the CI of the first group fall
# within the boundaries of the second group
decision <- (lb2 <= lb1 & lb1 <= ub2) |
(lb2 <= ub1 & ub1 <= ub2) |
(lb1<=lb2 & ub1>=ub2) |
(lb2<=lb1 & ub2>=ub1)
# Structure output
out <- data.frame(lb1,ub1,lb2,ub2,decision)
# that can be solved more elegant
# If there is no parameter to compare the data.frame has zero rows
# If the if is removed there are problems as the number of columns does not match
# the length of the name vector
if(nrow(out)!=0){
colnames(out)=c(paste0("lb_",names(cis_comp[[comp]])[1]),
paste0("ub_",names(cis_comp[[comp]])[1]),
paste0("lb_",names(cis_comp[[comp]])[2]),
paste0("ub_",names(cis_comp[[comp]])[2]),
"Decision"
)
out$Name <- para_rel
out <- out[, c(ncol(out), 1:(ncol(out) - 1))]
rownames(out) <- NULL
out
}
out
})
do.call(rbind,tt)
})
names(t) <- names(cis_comp[[comp]][[1]])
t
})
names(ttt) <- names(cis_comp)
ttt
})
names(decision_ci_overlap) <- paste0((1 - .alpha) * 100, "%")
# Decision overall
decision_overall_ci_overlap<-lapply(names(decision_ci_overlap), function(alpha) {
t <- lapply(names(decision_ci_overlap[[alpha]]), function(comp) {
tt <- lapply(names(decision_ci_overlap[[alpha]][[comp]]),function(interval_type){
# Check whether there is one FALSE among the decisions
all(decision_ci_overlap[[alpha]][[comp]][[interval_type]][,"Decision"])
})
names(tt) <- names(decision_ci_overlap[[alpha]][[comp]])
tt
})
names(t) <- names(decision_ci_overlap[[alpha]])
# Check whether there is one FALSE among the comparions' decisions
temp <- purrr::transpose(t)
overall_dec<-lapply(names(temp),function(x){
all(unlist(temp[[x]]))
})
names(overall_dec)<-names(temp)
overall_dec
})
names(decision_overall_ci_overlap) <- names(decision_ci_overlap)
}#end if: if one CI approach is requested
}
### Return output ============================================================
out <- list()
## General Information
out[["Information"]] <- list(
"Group_names" = names(.object),
"Number_of_observations"= sapply(X_all_list, nrow),
"Approach" = .approach_mgd,
"Approach_p_adjust" = .approach_p_adjust,
"Alpha" = .alpha,
"Parameters_to_compare" = .parameters_to_compare)
## Permutation-specific information
if(any(.approach_mgd %in% c("all", "Klesel", "Chin", "Sarstedt"))) {
out[["Information"]][["Information_permutation"]] <- list(
"Number_admissibles" = length(ref_dist1),
"Total_runs" = counter + n_inadmissibles,
"Permutation_seed" = .seed,
"Permutation_values" = list(),
"Handle_inadmissibles" = .handle_inadmissibles
)
} else {
out[["Information"]][["Information_permutation"]] <- list(
"Number_admissibles" = NA,
"Total_runs" = NA,
"Permutation_seed" = NA,
"Permutation_values" = NA,
"Handle_inadmissibles" = NA
)
}
## Bootstrap-specific information
if(any(.approach_mgd %in% c("all", "Sarstedt", "Keil", "Nitzl", "Henseler",
"CI_param","CI_overlap"))) {
# Collect bootstrap information
info_boot <-lapply(.object, function(x){
if(inherits(.object, "cSEMResults_2ndorder")) {
x <- x$Second_stage$Information$Resamples$Information_resample
} else {
x <- x$Information$Information_resample
}
list(
"Number_admissibles" = x$Number_of_admissibles,
"Total_runs" = x$Number_of_runs,
"Bootstrap_seed" = x$Seed,
"Handle_inadmissibles" = x$Handle_inadmissibles
)
})
names(info_boot) <- names(.object)
info_boot <- purrr::transpose(info_boot)
out[["Information"]][["Information_bootstrap"]] <- list(
"Number_admissibles" = info_boot$Number_admissibles,
"Total_runs" = info_boot$Total_runs,
"Bootstrap_seed" = info_boot$Bootstrap_seed,
"Handle_inadmissibles" = info_boot$Handle_inadmissibles
)
} else {
out[["Information"]][["Information_bootstrap"]] <- list(
"Number_admissibles" = NA,
"Total_runs" = NA,
"Bootstrap_seed" = NA,
"Handle_inadmissibles" = NA
)
}
# Information for Klesel et al. approach
if(any(.approach_mgd %in% c("all", "Klesel"))) {
out[["Klesel"]] <- list(
"Test_statistic" = teststat_Klesel,
"P_value" = pvalue_Klesel,
"Decision" = decision_Klesel,
"VCV_type" = .type_vcv
)
out[["Information"]][["Information_permutation"]][["Permutation_values"]][["Klesel"]] <- ref_dist_matrix_Klesel
}
# Information for Chin & Dibbern approach
if(any(.approach_mgd %in% c("all", "Chin"))) {
out[["Chin"]] <- list(
"Test_statistic" = teststat_Chin,
"P_value" = padjusted_Chin,
"Decision" = decision_Chin,
"Decision_overall" = decision_overall_Chin
)
out[["Information"]][["Information_permutation"]][["Permutation_values"]][["Chin"]] <- ref_dist_matrices_Chin
}
# Information for Sarstedt et al. approach
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
out[["Sarstedt"]] <- list(
"Test_statistic" = teststat_Sarstedt,
"P_value" = padjusted_Sarstedt,
"Decision" = decision_Sarstedt,
"Decision_overall" = decision_overall_Sarstedt
)
out[["Information"]][["Information_permutation"]][["Permutation_values"]][["Sarstedt"]] <- ref_dist_matrix_Sarstedt
}
# Information for Keil approach
if(any(.approach_mgd %in% c("all", "Keil"))) {
out[["Keil"]] <- list(
"Test_statistic" = purrr::transpose(teststat_Keil)$teststat,
"P_value" = padjusted_Keil,
"Decision" = decision_Keil,
"Decision_overall" = decision_overall_Keil,
"df" = purrr::transpose(teststat_Keil)$df[[1]]
)
}
# Information for Nitzl approach
if(any(.approach_mgd %in% c("all", "Nitzl"))) {
out[["Nitzl"]] <- list(
"Test_statistic" = purrr::transpose(teststat_Nitzl)$teststat,
"P_value" = padjusted_Nitzl,
"Decision" = decision_Nitzl,
"Decision_overall" = decision_overall_Nitzl,
"df" = purrr::transpose(teststat_Nitzl)$df[[1]]
)
}
# Information for Henseler approach
if(any(.approach_mgd %in% c("all", "Henseler"))) {
out[["Henseler"]] <- list(
"Test_statistic" = teststat_Henseler,
"P_value" = padjusted_Henseler,
"Decision" = decision_Henseler,
"Decision_overall" = decision_overall_Henseler
)
}
# Information CI_param
if(any(.approach_mgd %in% c("all", "CI_para"))) {
out[["CI_para"]] <- list(
"Decision" = decision_ci_para,
"Decision_overall" = decision_overall_ci_para
)
}
# Information CI_param
if(any(.approach_mgd %in% c("all", "CI_overlap"))) {
out[["CI_overlap"]] <- list(
"Decision" = decision_ci_overlap,
"Decision_overall" = decision_overall_ci_overlap
)
}
class(out) <- "cSEMTestMGD"
if(.output_type == "complete") {
return(out)
} else {
structureTestMGDDecisions(out)
}
}
Try the cSEM package in your browser
Any scripts or data that you put into this service are public.
cSEM documentation built on Nov. 25, 2022, 1:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions testMGD
Documented in testMGD
#' Tests for multi-group comparisons
#'
#' \lifecycle{stable}
#'
#' This function performs various tests proposed in the context of multigroup analysis.
#'
#' The following tests are implemented:
#' \describe{
#' \item{`.approach_mgd = "Klesel"`: Approach suggested by \insertCite{Klesel2019;textual}{cSEM}}{
#' The model-implied variance-covariance matrix (either indicator
#' (`.type_vcv = "indicator"`) or construct (`.type_vcv = "construct"`))
#' is compared across groups. If the model-implied indicator or construct correlation
#' matrix based on a saturated structural model should be compared, set `.saturated = TRUE`.
#' To measure the distance between the model-implied variance-covariance matrices,
#' the geodesic distance (dG) and the squared Euclidean distance (dL) are used.
#' If more than two groups are compared, the average distance over all groups
#' is used.}
#' \item{`.approach_mgd = "Sarstedt"`: Approach suggested by \insertCite{Sarstedt2011;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' \insertCite{Sarstedt2011;textual}{cSEM} tests if parameter k is equal
#' across all groups. If several parameters are tested simultaneously
#' it is recommended to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value). By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details. Note: the
#' test has some severe shortcomings. Use with caution.}
#' \item{`.approach_mgd = "Chin"`: Approach suggested by \insertCite{Chin2010;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' \insertCite{Chin2010;textual}{cSEM} tests if parameter k is equal
#' between two groups. If more than two groups are tested for equality, parameter
#' k is compared between all pairs of groups. In this case, it is recommended
#' to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value) since this is essentially a multiple testing setup.
#' If several parameters are tested simultaneously, correction is by group
#' and number of parameters. By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details.}
#' \item{`.approach_mgd = "Keil"`: Approach suggested by \insertCite{Keil2000;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' \insertCite{Keil2000;textual}{cSEM} tests if parameter k is equal
#' between two groups. It is assumed, that the standard errors of the coefficients are
#' equal across groups. The calculation of the standard error of the parameter
#' difference is adjusted as proposed by \insertCite{Henseler2009;textual}{cSEM}.
#' If more than two groups are tested for equality, parameter k is compared
#' between all pairs of groups. In this case, it is recommended
#' to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value) since this is essentially a multiple testing setup.
#' If several parameters are tested simultaneously, correction
#' is by group and number of parameters. By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details.}
#' \item{`.approach_mgd = "Nitzl"`: Approach suggested by \insertCite{Nitzl2010;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' Similarly to \insertCite{Keil2000;textual}{cSEM}, a single parameter k is tested
#' for equality between two groups. In contrast to \insertCite{Keil2000;textual}{cSEM},
#' it is assumed, that the standard errors of the coefficients are
#' unequal across groups \insertCite{Sarstedt2011}{cSEM}.
#' If more than two groups are tested for equality, parameter k is compared
#' between all pairs of groups. In this case, it is recommended
#' to adjust the significance level or the p-values (in \pkg{cSEM} correction is
#' done by p-value) since this is essentially a multiple testing setup.
#' If several parameters are tested simultaneously, correction
#' is by group and number of parameters. By default
#' no multiple testing correction is done, however, several common
#' adjustments are available via `.approach_p_adjust`. See
#' \code{\link[stats:p.adjust]{stats::p.adjust()}} for details.}
#' \item{`.approach_mgd = "Henseler"`: Approach suggested by \insertCite{Henseler2007a;textual}{cSEM}}{
#' Groups are compared in terms of parameter differences across groups.
#' In doing so, the bootstrap estimates of one parameter are compared across groups.
#' In the literature, this approach is also known as PLS-MGA.
#' Originally, this test was proposed as an one-sided test.
#' In this function we perform a left-sided and a right-sided test
#' to investigate whether a parameter differs across two groups. In doing so, the significance
#' level is divided by 2 and compared to p-value of the left and right-sided test.
#' Moreover, `.approach_p_adjust` is ignored and no overall decision
#' is returned.
#' For a more detailed description, see also \insertCite{Henseler2009;textual}{cSEM}.}
#' \item{`.approach_mgd = "CI_param"`: Approach mentioned in \insertCite{Sarstedt2011;textual}{cSEM}}{
#' This approach is based on the confidence intervals constructed around the
#' parameter estimates of the two groups. If the parameter of one group falls within
#' the confidence interval of the other group and/or vice versa, it can be concluded
#' that there is no group difference.
#' Since it is based on the confidence intervals `.approach_p_adjust` is ignored.}
#' \item{`.approach_mgd = "CI_overlap"`: Approach mentioned in \insertCite{Sarstedt2011;textual}{cSEM}}{
#' This approach is based on the confidence intervals (CIs) constructed around the
#' parameter estimates of the two groups. If the two CIs overlap, it can be concluded
#' that there is no group difference.
#' Since it is based on the confidence intervals `.approach_p_adjust` is ignored.}
#' }
#'
#' Use `.approach_mgd` to choose the approach. By default all approaches are computed
#' (`.approach_mgd = "all"`).
#'
#' For convenience, two types of output are available. See the "Value" section below.
#'
#' By default, approaches based on parameter differences across groups compare
#' all parameters (`.parameters_to_compare = NULL`). To compare only
#' a subset of parameters provide the parameters in [lavaan model syntax][lavaan::model.syntax] just like
#' the model to estimate. Take the simple model:
#'
#' \preformatted{
#' model_to_estimate <- "
#' Structural model
#' eta2 ~ eta1
#' eta3 ~ eta1 + eta2
#'
#' # Each concept os measured by 3 indicators, i.e., modeled as latent variable
#' eta1 =~ y11 + y12 + y13
#' eta2 =~ y21 + y22 + y23
#' eta3 =~ y31 + y32 + y33
#' "
#' }
#' If only the path from eta1 to eta3 and the loadings of eta1 are to be compared
#' across groups, write:
#' \preformatted{
#'to_compare <- "
#' Structural parameters to compare
#' eta3 ~ eta1
#'
#' # Loadings to compare
#' eta1 =~ y11 + y12 + y13
#' "
#' }
#' Note that the "model" provided to `.parameters_to_compare`
#' does not need to be an estimable model!
#'
#' Note also that compared to all other functions in \pkg{cSEM} using the argument,
#' `.handle_inadmissibles` defaults to `"replace"` to accommodate the Sarstedt et al. (2011) approach.
#'
#' Argument `.R_permuation` is ignored for the `"Nitzl"` and the `"Keil"` approach.
#' `.R_bootstrap` is ignored if `.object` already contains resamples,
#' i.e. has class `cSEMResults_resampled` and if only the `"Klesel"` or the `"Chin"`
#' approach are used.
#'
#' The argument `.saturated` is used by `"Klesel"` only. If `.saturated = TRUE`
#' the original structural model is ignored and replaced by a saturated model,
#' i.e. a model in which all constructs are allowed to correlate freely.
#' This is useful to test differences in the measurement models between groups
#' in isolation.
#'
#' @usage testMGD(
#' .object = NULL,
#' .alpha = 0.05,
#' .approach_p_adjust = "none",
#' .approach_mgd = c("all", "Klesel", "Chin", "Sarstedt",
#' "Keil", "Nitzl", "Henseler", "CI_para","CI_overlap"),
#' .output_type = c("complete", "structured"),
#' .parameters_to_compare = NULL,
#' .eval_plan = c("sequential", "multicore", "multisession"),
#' .handle_inadmissibles = c("replace", "drop", "ignore"),
#' .R_permutation = 499,
#' .R_bootstrap = 499,
#' .saturated = FALSE,
#' .seed = NULL,
#' .type_ci = "CI_percentile",
#' .type_vcv = c("indicator", "construct"),
#' .verbose = TRUE
#' )
#'
#' @inheritParams csem_arguments
#' @param .handle_inadmissibles Character string. How should inadmissible results
#' be treated? One of "*drop*", "*ignore*", or "*replace*". If "*drop*", all
#' replications/resamples yielding an inadmissible result will be dropped
#' (i.e. the number of results returned will potentially be less than `.R`).
#' For "*ignore*" all results are returned even if all or some of the replications
#' yielded inadmissible results (i.e. number of results returned is equal to `.R`).
#' For "*replace*" resampling continues until there are exactly `.R` admissible solutions.
#' Defaults to "*replace*" to accommodate all approaches.
#'
#' @return If `.output_type = "complete"` a list of class `cSEMTestMGD`. Technically, `cSEMTestMGD` is a
#' named list containing the following list elements:
#'
#' \describe{
#' \item{`$Information`}{Additional information.}
#' \item{`$Klesel`}{A list with elements, `Test_statistic`, `P_value`, and `Decision`}
#' \item{`$Chin`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Sarstedt`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Keil`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Nitzl`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$Henseler`}{A list with elements, `Test_statistic`, `P_value`, `Decision`, and `Decision_overall`}
#' \item{`$CI_para`}{A list with elements, `Decision`, and `Decision_overall`}
#' \item{`$CI_overlap`}{A list with elements, `Decision`, and `Decision_overall`}
#' }
#'
#' If `.output_type = "structured"` a tibble (data frame) with the following columns
#' is returned.
#'
#' \describe{
#' \item{`Test`}{The name of the test.}
#' \item{`Comparision`}{The parameter that was compared across groups. If "overall"
#' the overall fit of the model was compared.}
#' \item{`alpha%`}{The test decision for a given "alpha" level. If `TRUE` the null
#' hypotheses was rejected; if FALSE it was not rejected.}
#' \item{`p-value_correction`}{The p-value correction.}
#' \item{`CI_type`}{Only for the "CI_para" and the "CI_overlap" test. Which confidence interval was used.}
#' \item{`Distance_metric`}{Only for Test = "Klesel". Which distance metric was used.}
#' }
#'
#' @references
#' \insertAllCited{}
#'
#' @seealso [csem()], [cSEMResults], [testMICOM()], [testOMF()]
#'
#' @example inst/examples/example_testMGD.R
#'
#' @export
testMGD <- function(
.object = NULL,
.alpha = 0.05,
.approach_p_adjust = "none",
.approach_mgd = c("all", "Klesel", "Chin", "Sarstedt",
"Keil", "Nitzl","Henseler", "CI_para","CI_overlap"),
.output_type = c("complete", "structured"),
.parameters_to_compare = NULL,
.eval_plan = c("sequential", "multicore", "multisession"),
.handle_inadmissibles = c("replace", "drop", "ignore"),
.R_permutation = 499,
.R_bootstrap = 499,
.saturated = FALSE,
.seed = NULL,
.type_ci = "CI_percentile",
.type_vcv = c("indicator", "construct"),
.verbose = TRUE
){
### Checks and errors ========================================================
# Check .approach_mgd argument choices
diff <- setdiff(.approach_mgd, args_default(TRUE)$.approach_mgd)
if(length(diff) != 0) {
stop2(
"The following error occured in the testMGD() function:\n",
"Unknown approach: ", paste0(diff, collapse = ", "), ".",
" Possible choices are: ",
paste0(args_default(TRUE)$.approach_mgd, collapse = ", "))
}
## Match arguments
.approach_mgd <- match.arg(.approach_mgd, several.ok = TRUE)
.handle_inadmissibles <- match.arg(.handle_inadmissibles) # no reference to
# args_default, because args_default()$.handle_inadmissibles
# has "drop" as default, but testMGD hast "replace".
.type_vcv <- match.arg(.type_vcv)
.output_type <- match.arg(.output_type)
# if(!all(.type_ci %in% c("CI_standard_z", "CI_standard_t", "CI_percentile",
# "CI_basic", "CI_bc", "CI_bca") )){
if(!all(.type_ci %in% args_default(.choices = TRUE)$.type_ci )){
stop2("The specified confidence interval in .type.ci is not valid.\n",
"Please choose one of the following: CI_standard_z, CI_standard_t,\n",
"CI_percentile, CI_basic, CI_bc, CI_bca.")
}
# Check if at least two groups are present
if(!inherits(.object, "cSEMResults_multi")) {
stop2(
"The following error occured in the testMGD() function:\n",
"At least two groups required."
)
}
# Sarstedt et al. (2011) is not allowed to be used in combination with
# .handle_inadmissibles == "drop as permutation test statistic are be dropped
# only because of estimations based on the permutated dataset are dropped.
if(any(.approach_mgd %in% c("all", "Sarstedt")) & .handle_inadmissibles == "drop"){
stop2(
"The following error occured in the testMGD() function:\n",
"Approach `'Sarstedt'` not supported if `.handle_inadmissibles == 'drop'`")
}
# If a nonlinear model is used Klesel et al. approach cannot be used as
# it is not clear how to calculate the model-implied VCV
# Extract model type information
if(inherits(.object[[1]], "cSEMResults_2ndorder")){
model_type <- .object[[1]]$Second_stage$Information$Model$model_type
} else {
model_type <- .object[[1]]$Information$Model$model_type
}
if(any(.approach_mgd %in% c("all", "Klesel")) & model_type == "Nonlinear"){
stop2("The following error occured in the testMGD() function:\n",
"The approach suggested by Klesel et al. (2019) cannot be applied",
" to nonlinear models as cSEM currently can not calculate",
" the model-implied VCV matrix for such models.\n",
"Consider setting `.approach_mgd = c('Chin', 'Sarstedt')`")
}
# Check if any of the group estimates are inadmissible
if(sum(unlist(verify(.object))) != 0) {
warning2(
"The following warning occured in the testMGD() function:\n",
"Initial estimation results for at least one group are inadmissible.\n",
"See `verify(.object)` for details.")
}
# Check if data for different groups is identical
if(TRUE %in% lapply(utils::combn(.object, 2, simplify = FALSE),
function(x){ identical(x[[1]], x[[2]])})){
warning2(
"The following warning occured in the testMGD() function:\n",
"At least two groups are identical. Results may not be meaningful.")
}
# If Henseler, CI_para or CI_overlap are used with
# adjustment of p-value different than "none" return warning
if(any(.approach_mgd %in% c("all", "Henseler","CI_para","CI_overlap")) & !all(.approach_p_adjust %in% "none")){
warning2(
"The following warning occured in the testMGD() function:\n",
"Currently, there is no p-value adjustment possible for the approach suggested by\n",
"Henseler (2007), CI_para, and CI_overlap. Adjustment is ignored for these approaches."
)
}
if(.verbose) {
cat(rule2("Several tests for multi-group comparisons",
type = 3), "\n\n")
}
# Order significance levels
.alpha <- .alpha[order(.alpha)]
# Get the names of the parameters to be compared.
names_all_param <- getParameterNames(.object, .model = .parameters_to_compare)
# getParameterNames() returns also measurement error and indicator correlations
# Currently they cannot be handeled, therefore an error is returned.
# FOR FUTURE: MEasurement errors and indicator correlations can be added to the parameters allowed in the comparison
if(!is.null(names_all_param$names_cor_measurement_error)|!is.null(names_all_param$names_cor_indicator)){
stop2("The following error occured in the testMGD() function:\n",
"Currenlty it is not allowed to compare measurement error covariance",
" and/or indicator covariances across groups.")
}
names_param <- unlist(names_all_param)
## Calculation of the test statistics========================================
teststat <- list()
## Klesel et al. (2019) ------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Klesel"))) {
## Get the model-implied VCV
fit <- fit(.object = .object,
.saturated = .saturated,
.type_vcv = .type_vcv)
## Compute test statistic
temp <- c(
"dG" = calculateDistance(.matrices = fit, .distance = "geodesic"),
"dL" = calculateDistance(.matrices = fit, .distance = "squared_euclidean")
)
## Save test statistic
teststat[["Klesel"]] <- temp
}
## Chin & Dibbern (2010) -----------------------------------------------------
if(any(.approach_mgd %in% c("all", "Chin"))) {
## Compute and save test statistic
teststat[["Chin"]] <- calculateParameterDifference(.object = .object,
.model = .parameters_to_compare)
}
## Sarstedt et al. (2011), Keil et al. (2000), Nitzl (2010), Henseler (2007) -----
# All these approaches require bootstrap
if(any(.approach_mgd %in% c("all", "Sarstedt", "Keil", "Nitzl", "Henseler", "CI_para","CI_overlap"))) {
# Check if .object already contains resamples; if not, run bootstrap
if(!inherits(.object, "cSEMResults_resampled")) {
if(.verbose) {
cat("Bootstrap cSEMResults objects ...\n\n")
}
.object <- resamplecSEMResults(
.object = .object,
.resample_method = "bootstrap",
.handle_inadmissibles = .handle_inadmissibles,
.R = .R_bootstrap,
.seed = .seed,
.eval_plan = .eval_plan
)
}
## Combine bootstrap results in one matrix
bootstrap_results <- lapply(.object, function(y) {
if(inherits(.object, "cSEMResults_2ndorder")) {
x <- y$Second_stage$Information$Resamples$Estimates$Estimates1
nobs <- nrow(y$First_stage$Information$Data)
} else {
x <- y$Estimates$Estimates_resample$Estimates1
nobs <- nrow(y$Information$Data)
}
path_resamples <- x$Path_estimates$Resampled
loading_resamples <- x$Loading_estimates$Resampled
weight_resamples <- x$Weight_estimates$Resampled
# Rename
cor_exo_cons_resamples <- x$Exo_construct_correlation$Resampled
n <- nrow(path_resamples)
# Calculation of the bootstrap SEs
ses <- infer(.object=y,.quantity = "sd")
path_se <- ses$Path_estimates$sd
loading_se <- ses$Loading_estimates$sd
weight_se <- ses$Weight_estimates$sd
cor_exo_cons_se <- ses$Exo_construct_correlation$sd
# Calculation of the bias
bias <- infer(.object=y,.quantity = "bias")
path_bias <- bias$Path_estimates$bias
loading_bias <- bias$Loading_estimates$bias
weight_bias <- bias$Weight_estimates$bias
cor_exo_cons_bias <- bias$Exo_construct_correlation$bias
# Structure output
out<-list(
"n" = n,
"nObs" = nobs,
"para_all" = cbind(path_resamples,loading_resamples,weight_resamples, cor_exo_cons_resamples),
"ses_all" = c(path_se, loading_se, weight_se, cor_exo_cons_se),
"bias_all" = c(path_bias, loading_bias, weight_bias, cor_exo_cons_bias)
)
# Approaches based on CIs comparison ----
if(any(.approach_mgd %in% c("all", "CI_para", "CI_overlap"))){
diff <- setdiff(.type_ci, args_default(TRUE)$.type_ci)
if(length(diff) != 0) {
stop2(
"The following error occured in the testMGD() function:\n",
"Unknown approach: ", paste0(diff, collapse = ", "), ".",
" Possible choices are: ",
paste0(args_default(TRUE)$.type_ci, collapse = ", "))
}
# calculation of the CIs
cis <- infer(.object=y, .quantity=.type_ci, .alpha = .alpha)
cis_temp <- purrr::transpose(cis)
cis_ret<-lapply(cis_temp,function(x){
path_ci <- x$Path_estimates
loading_ci <- x$Loading_estimates
weight_ci <- x$Weight_estimates
cor_exo_cons_ci <- x$Exo_construct_correlation
# deliberately name them xyz_estimates to be able to select them later via names.
list(path_estimates = path_ci, loading_estimates = loading_ci,
weight_estimates = weight_ci, cor_exo_cons_estimates = cor_exo_cons_ci)
})
names(cis_ret) <- names(cis_temp)
out[["ci_all"]] <- cis_ret
}
# Return output
return(out)
})
## Keil, Nitzl, Henseler
if(any(.approach_mgd %in% c("all", "Nitzl", "Keil","Henseler"))) {
# Calculate the difference for one parameter between two groups
# Although Henseler approach does not compute a test statistic,
# The difference is calculated to create a dummy test statistic list
diff_para_Keil <- diff_para_Nitzl <- diff_para_Henseler<- calculateParameterDifference(
.object = .object,
.model = .parameters_to_compare
)
if(any(.approach_mgd %in% c("all","Henseler"))) {
# Create dummy test statistic list containing NAs
temp <- rep(NA,length(unlist(diff_para_Henseler)))
teststat[["Henseler"]] <-relist(flesh = temp,skeleton = diff_para_Henseler)
}
# Build list that contains pairs of objects.
object_permu <- utils::combn(.object, 2, simplify = FALSE)
names(object_permu) <- sapply(object_permu, function(x) paste0(names(x)[1], '_', names(x)[2]))
# Approach suggested by Keil
if(any(.approach_mgd %in% c("all", "Keil"))) {
teststat_Keil <- lapply(names(object_permu), function(x) {
diff <- diff_para_Keil[[x]]
ses1 <- bootstrap_results[[names(object_permu[[x]][1])]]$ses_all
ses2 <- bootstrap_results[[names(object_permu[[x]][2])]]$ses_all
# Get sample size per group
n1<-bootstrap_results[[names(object_permu[[x]][1])]]$nObs
n2<-bootstrap_results[[names(object_permu[[x]][2])]]$nObs
# Calculation of the SE of the parameter difference as proposed by
# Henseler (2007a), Henseler et al. (2009)
# Assumption of equal variances
ses_total <- sqrt((n1-1)^2/(n1+n2-2)*ses1^2 +
(n2-1)^2/(n1+n2-2)*ses2^2)*sqrt(1/n1+1/n2)
test_stat <- diff/ses_total[names(diff)]
list("teststat" = test_stat, "df" = n1 + n2 - 2)
})
names(teststat_Keil) <- names(object_permu)
teststat[["Keil"]] <- teststat_Keil
}
# Approach suggested by Nitzl (2010)
if(any(.approach_mgd %in% c("all", "Nitzl"))) {
teststat_Nitzl <- lapply(names(object_permu), function(x){
diff <- diff_para_Nitzl[[x]]
ses1 <- bootstrap_results[[names(object_permu[[x]][1])]]$ses_all
ses2 <- bootstrap_results[[names(object_permu[[x]][2])]]$ses_all
# Get sample size per group
n1<-bootstrap_results[[names(object_permu[[x]][1])]]$nObs
n2<-bootstrap_results[[names(object_permu[[x]][2])]]$nObs
# Calculation of the SE of the parameter difference as proposed by
# Henseler (2007a), Henseler et al. (2009)
# Assumption of unequal variances
ses_total <- sqrt((n1-1)/(n1)*ses1^2 +
(n2-1)/(n2)*ses2^2)
test_stat <- diff/ses_total[names(diff)]
# calculation of the degrees of freedom
numerator <- ((n1-1)/n1*ses1^2+(n2-1)/n2*ses2^2)^2
denominator <- (n1-1)/n1^2*ses1^4+(n2-1)/n2^2*ses2^4
df <- round(numerator/denominator-2)
df <- df[names(diff)]
list("teststat" = test_stat, "df" = df)
})
names(teststat_Nitzl) <- names(object_permu)
teststat[["Nitzl"]] <- teststat_Nitzl
}
}
## Sarstedt et al. approach
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
## Transpose
ll <- purrr::transpose(bootstrap_results)
group_id <- rep(1:length(.object), unlist(ll$n))
# all_comb contains all parameter estimate that could potentially be compared
# plus an id column indicating the group adherance of each row.
all_comb <- cbind(do.call(rbind,ll$para_all),
"group_id" = group_id)
## Select relevant columns
all_comb <- all_comb[, c(names_param, "group_id")]
## Add test statistic Sarstedt
teststat[["Sarstedt"]] <- calculateFR(.resample_sarstedt = all_comb)
} # END approach_mgd %in% c("all", "Sarstedt")
} # END .approach_mgd %in% c("all", "Sarstedt", "Keil", "Nitzl", "Henseler")
### Permutation ==============================================================
## Preparation
# Put data of each groups in a list and combine
if(inherits(.object, "cSEMResults_2ndorder")) {
# Data is saved in the first stage
X_all_list <- lapply(.object, function(x) x$First_stage$Information$Data)
# Collect initial arguments (from the first object, but could be any other)
arguments <- .object[[1]]$Second_stage$Information$Arguments_original
} else {
X_all_list <- lapply(.object, function(x) x$Information$Data)
# Collect initial arguments (from the first object, but could be any other)
arguments <- .object[[1]]$Information$Arguments
}
X_all <- do.call(rbind, X_all_list)
# Create a vector "id" to be used to randomly select groups (permutate) and
# set id as an argument in order to identify the groups.
id <- rep(1:length(X_all_list), sapply(X_all_list, nrow))
arguments[[".id"]] <- "id"
# Permutation is only performed for approaches that require it
if(any(.approach_mgd %in% c("all", "Klesel", "Chin", "Sarstedt"))) {
# Start progress bar if required
if(.verbose){
cat("Permutation ...\n\n")
# pb <- txtProgressBar(min = 0, max = .R_permutation, style = 3)
}
# Save old seed and restore on exit! This is important since users may have
# set a seed before, in which case the global seed would be
# overwritten if not explicitly restored
old_seed <- .Random.seed
on.exit({.Random.seed <<- old_seed})
## Create seed if not already set
if(is.null(.seed)) {
set.seed(seed = NULL)
# Note (08.12.2019): Its crucial to call set.seed(seed = NULL) before
# drawing a random seed out of .Random.seed. If set.seed(seed = NULL) is not
# called sample(.Random.seed, 1) would result in the same random seed as
# long as .Random.seed remains unchanged. By resetting the seed we make
# sure that sample draws a different element everytime it is called.
.seed <- sample(.Random.seed, 1)
}
## Set seed
set.seed(.seed)
## Calculate reference distribution
ref_dist <- list()
n_inadmissibles <- 0
counter <- 0
progressr::with_progress({
progress_bar_csem <- progressr::progressor(along = 1:.R_permutation)
repeat{
# Counter
counter <- counter + 1
progress_bar_csem(message = sprintf("Permutation run = %g", counter))
# Permutate data
X_temp <- cbind(X_all, id = sample(id))
# Replace the old dataset by the new permutated dataset
arguments[[".data"]] <- X_temp
# Estimate model
Est_temp <- do.call(csem, arguments)
# Check status
status_code <- sum(unlist(verify(Est_temp)))
# Distinguish depending on how inadmissibles should be handled
if(status_code == 0 | (status_code != 0 & .handle_inadmissibles == "ignore")) {
# Compute if status is ok or .handle inadmissibles = "ignore" AND the status is
# not ok
### Calculation of the test statistic for each resample ==================
teststat_permutation <- list()
## Klesel et al. (2019) --------------------------------------------------
if(any(.approach_mgd %in% c("all", "Klesel"))) {
## Get the model-implied VCV
fit_temp <- fit(Est_temp, .saturated = .saturated, .type_vcv = .type_vcv)
## Compute test statistic
temp <- c(
"dG" = calculateDistance(.matrices = fit_temp, .distance = "geodesic"),
"dL" = calculateDistance(.matrices = fit_temp, .distance = "squared_euclidean")
)
## Save test statistic
teststat_permutation[["Klesel"]] <- temp
}
## Chin & Dibbern (2010) -------------------------------------------------
if(any(.approach_mgd %in% c("all", "Chin"))) {
## Compute and save test statistic
teststat_permutation[["Chin"]] <- calculateParameterDifference(
.object = Est_temp,
.model = .parameters_to_compare)
}
## Sarstedt et al. (2011) ------------------------------------------------
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
# Permutation of the bootstrap parameter estimates
all_comb_permutation <- all_comb
all_comb_permutation[ , "group_id"] <- sample(group_id)
teststat_permutation[["Sarstedt"]] <- calculateFR(all_comb_permutation)
}
ref_dist[[counter]] <- teststat_permutation
} else if(status_code != 0 & .handle_inadmissibles == "drop") {
# Set list element to zero if status is not okay and .handle_inadmissibles == "drop"
ref_dist[[counter]] <- NA
} else {# status is not ok and .handle_inadmissibles == "replace"
# Reset counter and raise number of inadmissibles by 1
counter <- counter - 1
n_inadmissibles <- n_inadmissibles + 1
}
# Update progres bar
# if(.verbose){
# setTxtProgressBar(pb, counter)
# }
# Break repeat loop if .R results have been created.
if(length(ref_dist) == .R_permutation) {
break
} else if(counter + n_inadmissibles == 10000) {
# Stop if 10000 runs did not result in insufficient admissible results
stop("Not enough admissible result.", call. = FALSE)
}
} # END repeat
}) # END with_progress
# close progress bar
# if(.verbose){
# close(pb)
# }
### Postprocessing ===========================================================
# Delete potential NA's
ref_dist1 <- Filter(Negate(anyNA), ref_dist)
}
# # Order significance levels
# .alpha <- .alpha[order(.alpha)]
## Klesel et al. (2019) ------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Klesel"))) {
# Collect permutation results and combine
ref_dist_Klesel <- lapply(ref_dist1, function(x) x$Klesel)
ref_dist_matrix_Klesel <- do.call(cbind, ref_dist_Klesel)
# Extract test statistic
teststat_Klesel <- teststat$Klesel
# Calculation of p-values
pvalue_Klesel <- rowMeans(ref_dist_matrix_Klesel >= teststat_Klesel)
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Klesel <- lapply(.alpha, function(x) {
pvalue_Klesel > x
})
names(decision_Klesel) <- paste0(.alpha * 100, "%")
}
## Chin & Dibbern (2010) -----------------------------------------------------
if(any(.approach_mgd %in% c("all", "Chin"))) {
# Extract test statistic
teststat_Chin <- teststat$Chin
# Create list with matrices containing the reference distribution
# of the parameter differences
ref_dist_Chin <- lapply(ref_dist1, function(x) x$Chin)
# Transpose
ref_dist_Chin_temp <- purrr::transpose(ref_dist_Chin)
names(ref_dist_Chin_temp) <- names(teststat_Chin)
ref_dist_matrices_Chin <- lapply(ref_dist_Chin_temp, function(x) {
temp <- do.call(cbind, x)
temp_ind <- stats::complete.cases(temp)
temp[temp_ind, ,drop = FALSE]
})
# Calculation of the p-values
pvalue_Chin <- lapply(1:length(ref_dist_matrices_Chin), function(x) {
# Share of values above the positive test statistic
rowMeans(ref_dist_matrices_Chin[[x]] >= abs(teststat_Chin[[x]])) +
# share of values of the reference distribution below the negative test statistic
rowMeans(ref_dist_matrices_Chin[[x]] <= (-abs(teststat_Chin[[x]])))
})
names(pvalue_Chin) <- names(ref_dist_matrices_Chin)
# Adjust p-values
padjusted_Chin <- lapply(as.list(.approach_p_adjust), function(x){
# It is important to unlist the pvalues as pAdjust needs to now how many p-values
# there are to do a proper adjustment
pvector <- stats::p.adjust(unlist(pvalue_Chin),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Chin)
})
names(padjusted_Chin) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Chin <- lapply(padjusted_Chin, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
group_comp > alpha
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# Overall decision, i.e., was any of the test belonging to one significance level
# and one p value adjustment rejected
decision_overall_Chin <- lapply(decision_Chin, function(decision_Chin_list){
lapply(decision_Chin_list,function(x){
all(unlist(x))
})
})
}
## Sarstedt et al. (2011) ----------------------------------------------------
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
# Extract test statistic
teststat_Sarstedt <- teststat$Sarstedt
# Collect permuation results and combine to matrix
ref_dist_Sarstedt <- lapply(ref_dist1, function(x) x$Sarstedt)
ref_dist_matrix_Sarstedt <- do.call(cbind, ref_dist_Sarstedt)
# Calculation of the p-value
pvalue_Sarstedt <- rowMeans(ref_dist_matrix_Sarstedt >= teststat_Sarstedt)
# Adjust pvalues:
padjusted_Sarstedt<- lapply(as.list(.approach_p_adjust), function(x){
pvector <- stats::p.adjust(pvalue_Sarstedt, method = x)
})
names(padjusted_Sarstedt) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Sarstedt <- lapply(padjusted_Sarstedt,function(p_value){
temp <- lapply(.alpha, function(alpha){
p_value > alpha
})
names(temp) <- paste(.alpha*100,"%",sep= '')
temp
})
# Decision overall
decision_overall_Sarstedt <- lapply(decision_Sarstedt, function(x){#over p-value adjustments
lapply(x, function(xx){ #over different significant levels
all(xx)
})
})
}
## Keil et al. (2000) --------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Keil"))){
# Calculate p-values
pvalue_Keil <- lapply(teststat$Keil,function(x){
p_value <- 2*(1-pt(abs(x$teststat),df = x$df))
p_value
})
# Adjust p-values in case of more than one comparison
padjusted_Keil<- lapply(as.list(.approach_p_adjust), function(x){
pvector <- stats::p.adjust(unlist(pvalue_Keil),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Keil)
})
names(padjusted_Keil) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Keil <- lapply(padjusted_Keil, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
group_comp > alpha
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# One rejection leads to overall rejection
decision_overall_Keil <- lapply(decision_Keil, function(decision_Keil_list){
lapply(decision_Keil_list,function(x){
all(unlist(x))
})
})
}
## Nitzl et al. (2010) -------------------------------------------------------
if(any(.approach_mgd %in% c("all", "Nitzl"))){
# Calculate p-values
pvalue_Nitzl <- lapply(teststat$Nitzl,function(x){
p_value <- 2*(1-pt(abs(x$teststat),df = x$df))
p_value
})
# Adjust p-values in case of more than one comparison
padjusted_Nitzl<- lapply(as.list(.approach_p_adjust), function(x){
pvector <- stats::p.adjust(unlist(pvalue_Nitzl),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Nitzl)
})
names(padjusted_Nitzl) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Nitzl <- lapply(padjusted_Nitzl, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
group_comp > alpha
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# One rejection leads to overall rejection
decision_overall_Nitzl <- lapply(decision_Nitzl, function(decision_Nitzl_list){
lapply(decision_Nitzl_list,function(x){
all(unlist(x))
})
})
}
## Henseler (2009) approach -------------------------------------------------
if(any(.approach_mgd %in% c("all", "Henseler"))) {
# Pseudo teststat for convenience
teststat_Henseler <- teststat$Henseler
# center the bootstrap sample Sarstedt et al. (2011) Eq. 4
ll_centered <- lapply(bootstrap_results, function(x){
# Substract from each row of para_all the corresponding bias
t(apply(x$para_all,1,function(row){
row-x$bias_all
}))
})
names(ll_centered) <- names(bootstrap_results)
# Create pairs which should be compared
pairs_centered <- utils::combn(ll_centered, 2, simplify = FALSE)
names(pairs_centered) <- sapply(pairs_centered, function(x) paste0(names(x)[1], '_', names(x)[2]))
# Calculation of the probability
pvalue_Henseler <- lapply(pairs_centered,function(x){
calculatePr(.resample_centered = x,
.parameters_to_compare = names_param)
})
# Adjust p-value in case of multiple comparisons
# Adjusting p-values is not straight forward.
# First it is a one-sided test, so we need to know the hypothesis
# Just flipping the p-value is not really an option as it might cause problems
# in situation where the order of the p-values is required for the correction
# Therefore only the "none"method is applied
padjusted_Henseler<- lapply(as.list("none"), function(x){
pvector <- stats::p.adjust(unlist(pvalue_Henseler),method = x)
# Sort them back into list
relist(flesh = pvector,skeleton = pvalue_Henseler)
})
names(padjusted_Henseler) <- "none"
# Decision is made:
# The probability is compared to alpha and 1-alpha
# In doing so, we assess both hypotheses theta^1 <= theta^2 and
# theta^1 => theta^2
# TRUE = p-value > alpha & p-value < 1-alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
decision_Henseler <- lapply(padjusted_Henseler, function(adjust_approach){ # over the different p adjustments
temp <- lapply(.alpha, function(alpha){# over the different significance levels
lapply(adjust_approach,function(group_comp){# over the different group comparisons
# check whether the p values are larger than a certain alpha
# The alpha is divided by two to mimic a two-sided test.
group_comp > alpha/2 & group_comp < 1- alpha/2
})
})
names(temp) <- paste0(.alpha*100, "%")
temp
})
# Overall decision; if one coefficient is significanty different across groups
# it is rejected.
decision_overall_Henseler <- lapply(decision_Henseler, function(decision_Henseler_list){
lapply(decision_Henseler_list,function(x){
all(unlist(x))
})
})
}
# Comparison via confidence intervals ------------------
if(any(.approach_mgd %in% c("all", "CI_para", "CI_overlap"))) {
# Select CIs from the bootstrap_results
cis <- lapply(bootstrap_results,function(x){
x$ci_all
})
# Make group pairs of CIs
cis_comp <- utils::combn(cis, 2, simplify = FALSE)
names(cis_comp) <- sapply(cis_comp, function(x) paste0(names(x)[1], '_', names(x)[2]))
# Select the relevant parameters per group and make group pairs
# This is required fo both CI_para and CI_overlap
# For the latter it is used to select the appropriate parameters
param_per_group <- getRelevantParameters(.object = .object,
.model = .parameters_to_compare)
param_per_group <- purrr::transpose(param_per_group)
param_comp <- utils::combn(param_per_group, 2, simplify = FALSE)
names(param_comp) <- sapply(param_comp,
function(x)
paste0(names(x)[1], '_', names(x)[2]))
# CI_para -------------------------------------------
if(any(.approach_mgd %in% c("all", "CI_para"))) {
# Investigate whether the estimate of one group is part of the CI of another group
decision_ci_para <- lapply(.alpha, function(alpha) {
# lapply over the comparisons
tttt <- lapply(names(cis_comp), function(comp) {
# lapply over the type of CI for the first group
# As group1 and group2 have the same structure the selection does not matter
tt = lapply(names(cis_comp[[comp]][[1]]), function(interval_type){
# lapply over the different parameters
ttt <- lapply(names(cis_comp[[comp]][[1]][[interval_type]]),
function(param) {
lb <- paste0(100 * (1 - alpha), "%L")
ub <- paste0(100 * (1 - alpha), "%U")
# Estimate first group
temp_para1 <- param_comp[[comp]][[1]][[param]]
# CIs first group
temp_cis1 <- cis_comp[[comp]][[1]][[interval_type]][[param]]
temp_cis_selected1 <- temp_cis1[, names(temp_para1), drop = FALSE]
lb_temp1 <-temp_cis_selected1[lb,]
names(lb_temp1) <- colnames(temp_cis_selected1[lb,,drop=FALSE])
ub_temp1 <-temp_cis_selected1[ub,]
names(ub_temp1) <- colnames(temp_cis_selected1[ub,,drop=FALSE])
# Estimate second group
temp_para2 <- param_comp[[comp]][[2]][[param]]
# CIs second group
temp_cis2 <- cis_comp[[comp]][[2]][[interval_type]][[param]]
temp_cis_selected2 <- temp_cis2[, names(temp_para2), drop = FALSE]
lb_temp2 <-temp_cis_selected2[lb,]
names(lb_temp2) <- colnames(temp_cis_selected2[lb,,drop=FALSE])
ub_temp2 <-temp_cis_selected2[ub,]
names(ub_temp2) <- colnames(temp_cis_selected2[ub,,drop=FALSE])
# browser()
# # check whether parameter estimates of one group
# falls within boundaries of the CIs of the other group
# TRUE parameter estimate of the other group falls within the CI
# => no group difference
# Otherwise, FALSE => group difference
decision<-
(lb_temp2 <= temp_para1 &
ub_temp2 >= temp_para1) |
(lb_temp1 <= temp_para2 &
ub_temp1 >= temp_para2)
# Structure output
out=data.frame("Estimate"=temp_para1,"lb"=lb_temp2,
"ub"=ub_temp2,"Estimate"=temp_para2,"lb"=lb_temp1,
"ub"=ub_temp1,"decision"=decision)
# that can be solved more elegant
# If there is no parameter to compare the data.frame has zero rows
# If the if is removed there are problems as the number of columns does not match
# the length of the name vector
if(nrow(out)!=0){
colnames(out)=c(paste0("Est_",names(cis_comp[[comp]])[1]),
paste0("lb_",names(cis_comp[[comp]])[2]),
paste0("ub_",names(cis_comp[[comp]])[2]),
paste0("Est_",names(cis_comp[[comp]])[2]),
paste0("lb_",names(cis_comp[[comp]])[1]),
paste0("ub_",names(cis_comp[[comp]])[1]),
"Decision")
}
# Make rownames explicit
out$Name <- rownames(out)
out <- out[, c(ncol(out), 1:(ncol(out) -1))]
rownames(out) <- NULL
out
})
do.call(rbind,ttt)
})
names(tt) = names(cis_comp[[comp]][[1]])
tt
})
names(tttt) <- names(cis_comp)
tttt
})
names(decision_ci_para) <- paste0((1 - .alpha) * 100, "%")
# Decision overall
decision_overall_ci_para<-lapply(names(decision_ci_para), function(alpha) {
t <- lapply(names(decision_ci_para[[alpha]]), function(comp) {
tt <- lapply(names(decision_ci_para[[alpha]][[comp]]),function(interval_type){
# Check whether there is one FALSE among the decisions per comparison
all(decision_ci_para[[alpha]][[comp]][[interval_type]][,"Decision"])
})
names(tt) <- names(decision_ci_para[[alpha]][[comp]])
tt
})
names(t) <- names(decision_ci_para[[alpha]])
# Check whether there is one FALSE among the comparions' decisions
temp <- purrr::transpose(t)
overall_dec<-lapply(names(temp),function(x){
all(unlist(temp[[x]]))
})
names(overall_dec)<-names(temp)
overall_dec
})
names(decision_overall_ci_para) <- names(decision_ci_para)
}#end if .approach_mgd == CI_para
# CI_overlap ---------------------------------------------
if(any(.approach_mgd %in% c("all", "CI_overlap"))) {
decision_ci_overlap <- lapply(.alpha, function(alpha) {
# lapply over the group pairs
ttt <- lapply(names(cis_comp), function(comp) {
# lapply over the type of CI for the first group
# As group1 and group2 have the same structure the selection does not matter
t <- lapply(names(cis_comp[[comp]][[1]]), function(interval_type) {
# lapply over the parameters that are compared
tt <- lapply(names(cis_comp[[comp]][[1]][[interval_type]]), function(param) {
lb <- paste0(100 * (1 - alpha), "%L")
ub <- paste0(100 * (1 - alpha), "%U")
# It does not matter which group is chosen to select the relevant parameters
# as in all groups the relevant parameters are the same
para_rel <- names(param_comp[[comp]][[1]][[param]])
# Select lower and upper bound of the CIs for the two groups
lb1 <- cis_comp[[comp]][[1]][[interval_type]][[param]][lb,para_rel ]
ub1 <- cis_comp[[comp]][[1]][[interval_type]][[param]][ub,para_rel ]
lb2 <- cis_comp[[comp]][[2]][[interval_type]][[param]][lb,para_rel ]
ub2 <- cis_comp[[comp]][[2]][[interval_type]][[param]][ub,para_rel ]
# Check whether the boundaries of the CI of the first group fall
# within the boundaries of the second group
decision <- (lb2 <= lb1 & lb1 <= ub2) |
(lb2 <= ub1 & ub1 <= ub2) |
(lb1<=lb2 & ub1>=ub2) |
(lb2<=lb1 & ub2>=ub1)
# Structure output
out <- data.frame(lb1,ub1,lb2,ub2,decision)
# that can be solved more elegant
# If there is no parameter to compare the data.frame has zero rows
# If the if is removed there are problems as the number of columns does not match
# the length of the name vector
if(nrow(out)!=0){
colnames(out)=c(paste0("lb_",names(cis_comp[[comp]])[1]),
paste0("ub_",names(cis_comp[[comp]])[1]),
paste0("lb_",names(cis_comp[[comp]])[2]),
paste0("ub_",names(cis_comp[[comp]])[2]),
"Decision"
)
out$Name <- para_rel
out <- out[, c(ncol(out), 1:(ncol(out) - 1))]
rownames(out) <- NULL
out
}
out
})
do.call(rbind,tt)
})
names(t) <- names(cis_comp[[comp]][[1]])
t
})
names(ttt) <- names(cis_comp)
ttt
})
names(decision_ci_overlap) <- paste0((1 - .alpha) * 100, "%")
# Decision overall
decision_overall_ci_overlap<-lapply(names(decision_ci_overlap), function(alpha) {
t <- lapply(names(decision_ci_overlap[[alpha]]), function(comp) {
tt <- lapply(names(decision_ci_overlap[[alpha]][[comp]]),function(interval_type){
# Check whether there is one FALSE among the decisions
all(decision_ci_overlap[[alpha]][[comp]][[interval_type]][,"Decision"])
})
names(tt) <- names(decision_ci_overlap[[alpha]][[comp]])
tt
})
names(t) <- names(decision_ci_overlap[[alpha]])
# Check whether there is one FALSE among the comparions' decisions
temp <- purrr::transpose(t)
overall_dec<-lapply(names(temp),function(x){
all(unlist(temp[[x]]))
})
names(overall_dec)<-names(temp)
overall_dec
})
names(decision_overall_ci_overlap) <- names(decision_ci_overlap)
}#end if: if one CI approach is requested
}
### Return output ============================================================
out <- list()
## General Information
out[["Information"]] <- list(
"Group_names" = names(.object),
"Number_of_observations"= sapply(X_all_list, nrow),
"Approach" = .approach_mgd,
"Approach_p_adjust" = .approach_p_adjust,
"Alpha" = .alpha,
"Parameters_to_compare" = .parameters_to_compare)
## Permutation-specific information
if(any(.approach_mgd %in% c("all", "Klesel", "Chin", "Sarstedt"))) {
out[["Information"]][["Information_permutation"]] <- list(
"Number_admissibles" = length(ref_dist1),
"Total_runs" = counter + n_inadmissibles,
"Permutation_seed" = .seed,
"Permutation_values" = list(),
"Handle_inadmissibles" = .handle_inadmissibles
)
} else {
out[["Information"]][["Information_permutation"]] <- list(
"Number_admissibles" = NA,
"Total_runs" = NA,
"Permutation_seed" = NA,
"Permutation_values" = NA,
"Handle_inadmissibles" = NA
)
}
## Bootstrap-specific information
if(any(.approach_mgd %in% c("all", "Sarstedt", "Keil", "Nitzl", "Henseler",
"CI_param","CI_overlap"))) {
# Collect bootstrap information
info_boot <-lapply(.object, function(x){
if(inherits(.object, "cSEMResults_2ndorder")) {
x <- x$Second_stage$Information$Resamples$Information_resample
} else {
x <- x$Information$Information_resample
}
list(
"Number_admissibles" = x$Number_of_admissibles,
"Total_runs" = x$Number_of_runs,
"Bootstrap_seed" = x$Seed,
"Handle_inadmissibles" = x$Handle_inadmissibles
)
})
names(info_boot) <- names(.object)
info_boot <- purrr::transpose(info_boot)
out[["Information"]][["Information_bootstrap"]] <- list(
"Number_admissibles" = info_boot$Number_admissibles,
"Total_runs" = info_boot$Total_runs,
"Bootstrap_seed" = info_boot$Bootstrap_seed,
"Handle_inadmissibles" = info_boot$Handle_inadmissibles
)
} else {
out[["Information"]][["Information_bootstrap"]] <- list(
"Number_admissibles" = NA,
"Total_runs" = NA,
"Bootstrap_seed" = NA,
"Handle_inadmissibles" = NA
)
}
# Information for Klesel et al. approach
if(any(.approach_mgd %in% c("all", "Klesel"))) {
out[["Klesel"]] <- list(
"Test_statistic" = teststat_Klesel,
"P_value" = pvalue_Klesel,
"Decision" = decision_Klesel,
"VCV_type" = .type_vcv
)
out[["Information"]][["Information_permutation"]][["Permutation_values"]][["Klesel"]] <- ref_dist_matrix_Klesel
}
# Information for Chin & Dibbern approach
if(any(.approach_mgd %in% c("all", "Chin"))) {
out[["Chin"]] <- list(
"Test_statistic" = teststat_Chin,
"P_value" = padjusted_Chin,
"Decision" = decision_Chin,
"Decision_overall" = decision_overall_Chin
)
out[["Information"]][["Information_permutation"]][["Permutation_values"]][["Chin"]] <- ref_dist_matrices_Chin
}
# Information for Sarstedt et al. approach
if(any(.approach_mgd %in% c("all", "Sarstedt"))) {
out[["Sarstedt"]] <- list(
"Test_statistic" = teststat_Sarstedt,
"P_value" = padjusted_Sarstedt,
"Decision" = decision_Sarstedt,
"Decision_overall" = decision_overall_Sarstedt
)
out[["Information"]][["Information_permutation"]][["Permutation_values"]][["Sarstedt"]] <- ref_dist_matrix_Sarstedt
}
# Information for Keil approach
if(any(.approach_mgd %in% c("all", "Keil"))) {
out[["Keil"]] <- list(
"Test_statistic" = purrr::transpose(teststat_Keil)$teststat,
"P_value" = padjusted_Keil,
"Decision" = decision_Keil,
"Decision_overall" = decision_overall_Keil,
"df" = purrr::transpose(teststat_Keil)$df[[1]]
)
}
# Information for Nitzl approach
if(any(.approach_mgd %in% c("all", "Nitzl"))) {
out[["Nitzl"]] <- list(
"Test_statistic" = purrr::transpose(teststat_Nitzl)$teststat,
"P_value" = padjusted_Nitzl,
"Decision" = decision_Nitzl,
"Decision_overall" = decision_overall_Nitzl,
"df" = purrr::transpose(teststat_Nitzl)$df[[1]]
)
}
# Information for Henseler approach
if(any(.approach_mgd %in% c("all", "Henseler"))) {
out[["Henseler"]] <- list(
"Test_statistic" = teststat_Henseler,
"P_value" = padjusted_Henseler,
"Decision" = decision_Henseler,
"Decision_overall" = decision_overall_Henseler
)
}
# Information CI_param
if(any(.approach_mgd %in% c("all", "CI_para"))) {
out[["CI_para"]] <- list(
"Decision" = decision_ci_para,
"Decision_overall" = decision_overall_ci_para
)
}
# Information CI_param
if(any(.approach_mgd %in% c("all", "CI_overlap"))) {
out[["CI_overlap"]] <- list(
"Decision" = decision_ci_overlap,
"Decision_overall" = decision_overall_ci_overlap
)
}
class(out) <- "cSEMTestMGD"
if(.output_type == "complete") {
return(out)
} else {
structureTestMGDDecisions(out)
}
}
Try the cSEM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.