#' Test measurement invariance of composites
#'
#' \lifecycle{stable}
#'
#' The functions performs the permutation-based test for measurement invariance
#' of composites across groups proposed by \insertCite{Henseler2016;textual}{cSEM}.
#' According to the authors assessing measurement invariance in composite
#' models can be assessed by a three-step procedure. The first two steps
#' involve an assessment of configural and compositional invariance.
#' The third steps involves mean and variance comparisons across groups.
#' Assessment of configural invariance is qualitative in nature and hence
#' not assessed by the [testMICOM()] function.
#'
#' As [testMICOM()] requires at least two groups, `.object` must be of
#' class `cSEMResults_multi`. As of version 0.2.0 of the package, [testMICOM()]
#' does not support models containing second-order constructs.
#'
#' It is possible to compare more than two groups, however, multiple-testing
#' issues arise in this case. To adjust p-values in this case several p-value
#' adjustments are available via the `approach_p_adjust` argument.
#'
#' The remaining arguments set the number of permutation runs to conduct
#' (`.R`), the random number seed (`.seed`),
#' instructions how inadmissible results are to be handled (`handle_inadmissibles`),
#' and whether the function should be verbose in a sense that progress is printed
#' to the console.
#'
#' The number of permutation runs defaults to `args_default()$.R` for
#' performance reasons. According to \insertCite{Henseler2016;textual}{cSEM}
#' the number of permutations should be at least 5000 for assessment to be
#' sufficiently reliable.
#'
#' @usage testMICOM(
#' .object = NULL,
#' .approach_p_adjust = "none",
#' .handle_inadmissibles = c("drop", "ignore", "replace"),
#' .R = 499,
#' .seed = NULL,
#' .verbose = TRUE
#' )
#'
#' @inheritParams csem_arguments
#'
#' @return
#' A named list of class `cSEMTestMICOM` containing the following list element:
#' \describe{
#' \item{`$Step2`}{A list containing the results of the test for compositional invariance (Step 2).}
#' \item{`$Step3`}{A list containing the results of the test for mean and variance equality (Step 3).}
#' \item{`$Information`}{A list of additional information on the test.}
#' }
#'
#' @references
#' \insertAllCited{}
#'
#' @example inst/examples/example_testMICOM.R
#'
#' @seealso [csem()], [cSEMResults], [testOMF()], [testMGD()]
#'
#' @export
#'
testMICOM <- function(
.object = NULL,
.approach_p_adjust = "none",
.handle_inadmissibles = c("drop", "ignore", "replace"),
.R = 499,
.seed = NULL,
.verbose = TRUE
) {
# Match arguments
.handle_inadmissibles <- match.arg(.handle_inadmissibles)
# Implementation is based on:
# Henseler et al. (2016) - Testing measurement invariance of composites using
# partial least squares
if(.verbose) {
cat(rule(center = "Test for measurement invariance based on Henseler et al (2016)",
line = "bar3"), "\n\n")
}
## If second-order
if(inherits(.object, "cSEMResults_2ndorder")) {
stop2("Currently, second-order models are not supported by `testMICOM()`.")
} else {
### Checks and errors ========================================================
## Check if at least two groups are present
if(!inherits(.object, "cSEMResults_multi")) {
stop2(
"The following error occured in the `testMICOM()` function:\n",
"At least two groups required."
)
}
if(sum(unlist(verify(.object))) != 0) {
stop2(
"The following error occured in the `testMICOM()` function:\n",
"Initial estimation results for at least one group are inadmissible.\n",
"See `verify(.object)` for details.")
}
if(.verbose & all(.object[[1]]$Information$Model$construct_type == "Common factor")) {
warning2(
"The following warning occured in the `testMICOM()` function:\n",
"All constructs are modelled as common factors.\n",
"Test results are only meaningful for composite models!")
}
# if(.verbose & length(.object) > 2) {
# warning2(
# "The following warning occured in the `testMICOM()` function:\n",
# "Comparing more than two groups inflates the familywise error rate.\n",
# "Interpret statistical significance with caution.\n",
# "Future versions of the package will likely include appropriate correction options.")
# }
# if(.approach_p_adjust != 'none'){
# stop2("P-value adjustment to control the familywise error rate not supported yet.")
# }
### Preparation ==============================================================
## Get pooled data (potentially unstandardized)
X <- .object[[1]]$Information$Data_pooled
X <- processData(X, .model = .object[[1]]$Information$Model)
## Remove id column:
# If .id has been supplied, delete column with the id name otherwise skip
# if(!is.null(.object[[1]]$Information$Arguments$.id)) {
# X <- X[, -which(colnames(X) == .object[[1]]$Information$Arguments$.id)]
# }
X <- as.matrix(X)
# Collect initial arguments (from the first object, but could be any other)
arguments <- .object[[1]]$Information$Arguments
# Create a vector "id" to be used to randomly select groups (permutate) and
# set id as an argument in order to identify the groups.
X_list <- lapply(.object, function(x) x$Information$Arguments$.data)
id <- rep(1:length(X_list), sapply(X_list, nrow))
arguments[[".id"]] <- "id"
### Step 1 - Configural invariance ===========================================
# Has to be assessed by the user prior to using the testMICOM function. See
# the original paper for details.
### Step 2 - Compositional invariance ========================================
# Procedure: see page 414 and 415 of the paper
## Compute proxies/scores using the pooled data
# (it does not matter if the data is scaled or unscaled as this does not
# affect the correlation)
H <- lapply(.object, function(x) X %*% t(x$Estimates$Weight_estimates))
## Compute the correlation of the scores for all group combinations
# Get the scores for all group combinations
H_combn <- utils::combn(H, 2, simplify = FALSE)
# Compute the correlation c for each group combination
c <- lapply(H_combn, function(x) diag(cor(x[[1]], x[[2]])))
# Set the names for each group combination
names(c) <- utils::combn(names(.object), 2, FUN = paste0,
collapse = "_", simplify = FALSE)
## Permutation ---------------------------------------------------------------
# Start progress bar
# if(.verbose){
# pb <- txtProgressBar(min = 0, max = .R, 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)
repeat{
# Counter
counter <- counter + 1
progress_bar_csem(message = sprintf("Permutation run = %g", counter))
# Permutate data
X_temp <- cbind(X, 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
## Compute weights for each group and use these to compute proxies/scores using
# the pooled data (= the original combined data). Note that these
# scores are unstandardized, however since we consider the correlation
# it does not matter whether we consider standardized or unstandardized
# proxies
H_temp <- lapply(Est_temp, function(x) X %*% t(x$Estimates$Weight_estimates))
## Compute the correlation of the scores for all group combinations
# Get the scores for all group combinations
H_combn_temp <- utils::combn(H_temp, 2, simplify = FALSE)
# Compute the correlation c for each group combination
c_temp <- lapply(H_combn_temp, function(x) diag(cor(x[[1]], x[[2]])))
# Set the names for each group combination
names(c_temp) <- utils::combn(names(H_temp), 2, FUN = paste0,
collapse = "_", simplify = FALSE)
ref_dist[[counter]] <- c_temp
} 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 progress bar
# if(.verbose){
# setTxtProgressBar(pb, counter)
# }
# Break repeat loop if .R results have been created.
if(length(ref_dist) == .R) {
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)
# }
# Delete potential NA's
ref_dist <- Filter(Negate(anyNA), ref_dist)
# Bind
temp <- do.call(rbind, lapply(ref_dist, function(x) do.call(rbind, x)))
temp <- split(as.data.frame(temp), rownames(temp))
# Order alphas (decreasing order)
# .alpha <- .alpha[order(.alpha)]
# critical_values_step2 <- lapply(lapply(temp, as.matrix), matrixStats::colQuantiles,
# probs = .alpha, drop = FALSE) # lower quantile needed, hence
# alpha and not 1 - alpha
# Calculate the p-value for the second step of MICOM
pvalue_step2<- lapply(1:length(temp), function(x) {
# Share of values above the positive test statistic
# temp contains a list of the reference distributions
rowMeans(t(temp[[x]]) <= c[[x]])
})
names(pvalue_step2) <- names(temp)
padjusted_step2 <- lapply(as.list(.approach_p_adjust), function(x){
# Select p_values per composite and only adjust those
temp <- purrr::transpose(pvalue_step2)
temp1 <- lapply(temp,function(comp){
stats::p.adjust(unlist(comp),method = x)
})
# sort them back
lapply(purrr::transpose(temp1),unlist)
})
names(padjusted_step2) <- .approach_p_adjust
# Decision
# TRUE = p-value > alpha --> not reject
# FALSE = sufficient evidence against the H0 --> reject
# decision_step2 <- lapply(padjusted_step2, 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
# as.matrix(group_comp > alpha,ncol=1)
# })
# })
# names(temp) <- paste0(.alpha*100, "%")
# temp
# })
# TRUE do not reject; FALSE: reject
### Step 3 - Equal mean values and variances==================================
# Update arguments
arguments[[".data"]] <- X
arguments[[".id"]] <- NULL
# Estimate model using pooled data set
Est <- do.call(csem, arguments)
# Procedure below:
# 1. Create list of original group id's + .R permutated id's
# 2. Extract construct scores, attach ids and convert to data.frame (for split)
# 3. Split construct scores by its id column
# --- Now the structure is a list of length 1 + .R each containing a list
# of the same length as there are numbers of groups (often just 2).
# 4. Convert group data set to a matrix and delete id column
# 5. Compute the mean and the variance for each replication (+ the original data of course)
# and data set for each group.
# 6. Transpose list
# --- Now we have a list of length 1 + .R containing two list elements
# "Mean" and "Var" which in turn contain as many list elements as there are
# groups
meanvar <- c(list(id), replicate(.R, sample(id), simplify = FALSE)) %>%
lapply(function(x) as.data.frame(cbind(Est$Estimates$Construct_scores, id = x))) %>%
lapply(function(x) split(x, f = x$id)) %>%
lapply(function(x) lapply(x, function(y) as.matrix(y[, -ncol(y), drop = FALSE]))) %>%
lapply(function(x) lapply(x, function(y) list(
"Mean" = colMeans(y),
"Var" = {tt <- matrixStats::colVars(y); names(tt) <- colnames(y); tt})
)) %>%
lapply(function(x) purrr::transpose(x))
# Compute the difference of the means for each possible group combination
# and attach a name
m <- meanvar %>%
lapply(function(x) {
tt <- utils::combn(x[["Mean"]], m = 2,
FUN = function(y) y[[1]] - y[[2]],
simplify = FALSE)
names(tt) <- utils::combn(names(.object), m = 2, FUN = paste0,
collapse = "_", simplify = FALSE)
tt
})
# Compute the log difference of the variances for each possible group combination
# and attach a name
v <- meanvar %>%
lapply(function(x) {
tt <- utils::combn(x[["Var"]], m = 2,
FUN = function(y) log(y[[1]]) - log(y[[2]]),
simplify = FALSE)
names(tt) <- utils::combn(names(.object), m = 2, FUN = paste0,
collapse = "_", simplify = FALSE)
tt
})
# Combine in a list and bind (log) differences for each replication in a list
mv <- list("Mean" = m, "Var"= v) %>%
lapply(function(x) lapply(x, function(y) do.call(rbind, y))) %>%
lapply(function(x) cbind(as.data.frame(do.call(rbind, x), row.names = FALSE),
id = rownames(do.call(rbind, x)))) %>%
lapply(function(x) split(x, f = x$id)) %>%
lapply(function(x) lapply(x, function(y) as.matrix(y[, -ncol(y), drop = FALSE])))
# Extract the original group differences. This is the first row of each group
# dataset
mv_o <- lapply(mv, function(x) lapply(x, function(y) y[1, ]))
# Test the means
teststat_mean <- mv_o$Mean
# Collect reference distribution
ref_dist_mean <- lapply(mv$Mean,function(x){x[-1,,drop=FALSE]})
# Calculate p-value
pvalue_mean <- lapply(1:length(ref_dist_mean), function(x) {
# Share of values above the positive test statistic
rowMeans(t(ref_dist_mean[[x]]) > abs(teststat_mean[[x]])) +
# share of values of the reference distribution below the negative test statistic
rowMeans(t(ref_dist_mean[[x]]) < (-abs(teststat_mean[[x]])))
})
names(pvalue_mean) <- names(ref_dist_mean)
# Adjust p-values: Correct pvalues by the number of groups
padjusted_mean <- lapply(as.list(.approach_p_adjust), function(x){
# Select p_values per composite and only adjust those
temp <- purrr::transpose(pvalue_mean)
# pAdjust needs to now how many p-values there are to do a proper adjustment
temp1 <- lapply(temp,function(comp){
stats::p.adjust(unlist(comp),method = x)
})
# sort them back
lapply(purrr::transpose(temp1),unlist)
})
names(padjusted_mean) <- .approach_p_adjust
# Make decision
# decision_mean <- lapply(padjusted_mean, 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
# as.matrix(group_comp > alpha,ncol=1)
# })
# })
# names(temp) <- paste0(.alpha*100, "%")
# temp
# })
# TRUE: do not reject; FALSE: reject
# Test the variance
teststat_var <- mv_o$Var
# Collect reference distribution
ref_dist_var <- lapply(mv$Var,function(x){x[-1,,drop=FALSE]})
# Calculate p-value
pvalue_var <- lapply(1:length(ref_dist_var), function(x) {
# Share of values above the positive test statistic
rowMeans(t(ref_dist_var[[x]]) > abs(teststat_var[[x]])) +
# share of values of the reference distribution below the negative test statistic
rowMeans(t(ref_dist_var[[x]]) < (-abs(teststat_var[[x]])))
})
names(pvalue_var) <- names(ref_dist_var)
# Adjust p-values, e.g., Bonferroni: Multiply all p-values by the number of comparisons
padjusted_var <- lapply(as.list(.approach_p_adjust), function(x){
# Select p_values per composite and only adjust those
temp <- purrr::transpose(pvalue_var)
# pAdjust needs to now how many p-values there are to do a proper adjustment
temp1 <- lapply(temp,function(comp){
stats::p.adjust(unlist(comp),method = x)
})
# sort them back
lapply(purrr::transpose(temp1),unlist)
})
names(padjusted_var) <- .approach_p_adjust
# Make decision
# decision_var <- lapply(padjusted_var, 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
# as.matrix(group_comp > alpha,ncol=1)
# })
# })
# names(temp) <- paste0(.alpha*100, "%")
# temp
# })
## Compute quantiles/critical values
# probs <- c()
# for(i in seq_along(.alpha)) {
# probs <- c(probs, .alpha[i]/2, 1 - .alpha[i]/2)
# }
#
# critical_values_step3 <- lapply(mv, function(x) lapply(x, function(y) y[-1, ])) %>%
# lapply(function(x) lapply(x, function(y) matrixStats::colQuantiles(y, probs = probs, drop = FALSE)))
#
# ## Compare critical value and test statistic
# # For Mean
# decision_m <- mapply(function(x, y) abs(x) < y[, seq(2, length(.alpha)*2, by = 2), drop = FALSE],
# x = mv_o[[1]],
# y = critical_values_step3[[1]],
# SIMPLIFY = FALSE)
# # For Var
# decision_v <- mapply(function(x, y) abs(x) < y[, seq(2, length(.alpha)*2, by = 2), drop = FALSE],
# x = mv_o[[2]],
# y = critical_values_step3[[2]],
# SIMPLIFY = FALSE)
### Return output ==========================================================
out <- list(
"Step2" = list(
"Test_statistic" = c,
"P_value" = padjusted_step2,
# "Decision" = decision_step2,
"Bootstrap_values" = ref_dist
),
"Step3" = list(
"Mean" = list(
"Test_statistic" = teststat_mean,
"P_value" = padjusted_mean
# "Decision" = decision_mean
),
"Var" = list(
"Test_statistic" = teststat_var,
"P_value" = padjusted_var
# "Decision" = decision_var
)
),
"Information" = list(
"Group_names" = names(.object),
"Number_admissibles" = length(ref_dist),
"Number_of_observations" = sapply(X_list, nrow),
"Total_runs" = counter + n_inadmissibles,
"Seed" = .seed
)
)
}
class(out) <- "cSEMTestMICOM"
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.