4. Supplementary Analyses

Load Data

Code
load("../data/wrangled_data/dt_items.RData")
dt_cfa <- dt_items %>% 
  select(
    -response_id,
    -age,
    -religion,
    -education,
    -gender,
    -starts_with("intuition"),
    -starts_with("logic"),
    -starts_with("past_activism") # not really a scale
  ) %>% 
  filter(country != "Jihadist") # fit statistics get calculated differently due to binary nature of scales
Code
impute_items = FALSE

if (impute_items) {
  res <- impute_flex(
    dt_cfa,                 # your other dataset
    scope = "all",          # <- this is the new mode
    id_vars = c("response_id"),  # protect all ID-like columns
    # never_impute_extra = c("created_at"),        # optionally protect extra columns
    m = 20,
    maxit = 10,
    seed = 42,
    verbose = TRUE,
    parallel_mode = "cluster",                # <— turn on cluster backend
    cores = parallel::detectCores() - 2,
    cl_type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK"
  )
  dt_cfa <- res$completed
}
Code
calculate_descriptives_table(dt_cfa, "cfa_analysis_dataframe")
Descriptive Statistics for cfa_analysis_dataframe
n mean sd median min max range se missing
country 10327 21.13 11.85 21 1 41 40 0.12 0
violent_intent_1 10327 2.00 1.89 1 1 8 7 0.02 0
violent_intent_2 10327 2.10 1.93 1 1 8 7 0.02 0
violent_intent_3 10327 1.99 1.88 1 1 8 7 0.02 0
violent_intent_4 10327 2.50 2.17 1 1 8 7 0.02 0
violent_intent_6 10327 2.04 1.94 1 1 8 7 0.02 0
peaceful_intent_7 10327 3.08 2.21 2 1 8 7 0.02 0
peaceful_intent_8 10327 3.07 2.22 2 1 8 7 0.02 0
peaceful_intent_9 10327 3.80 2.26 4 1 8 7 0.02 0
peaceful_intent_10 10327 3.79 2.25 4 1 8 7 0.02 0
peaceful_intent_11 10327 3.59 2.30 3 1 8 7 0.02 0
peaceful_intent_12 10327 3.26 2.31 3 1 8 7 0.02 0
harmonious_passion_1 10326 4.09 1.98 4 1 7 6 0.02 1
harmonious_passion_3 10327 3.78 2.08 4 1 7 6 0.02 0
harmonious_passion_5 10325 3.71 2.07 4 1 7 6 0.02 2
harmonious_passion_6 10326 3.70 2.09 4 1 7 6 0.02 1
harmonious_passion_8 10325 3.71 2.11 4 1 7 6 0.02 2
harmonious_passion_10 10324 3.78 2.07 4 1 7 6 0.02 3
obsessive_passion_2 10327 2.59 1.84 2 1 7 6 0.02 0
obsessive_passion_4 10326 2.82 2.01 2 1 7 6 0.02 1
obsessive_passion_7 10326 3.09 2.15 3 1 7 6 0.02 1
obsessive_passion_9 10325 3.37 2.18 3 1 7 6 0.02 2
obsessive_passion_11 10326 2.63 1.96 2 1 7 6 0.02 1
obsessive_passion_12 10325 2.73 2.05 2 1 7 6 0.02 2
commitment_passion_13 10326 2.92 2.02 2 1 7 6 0.02 1
commitment_passion_14 10325 4.69 1.94 5 1 7 6 0.02 2
commitment_passion_15 10323 4.33 2.07 4 1 7 6 0.02 4
commitment_passion_16 10321 3.57 2.21 3 1 7 6 0.02 6
identity_fusion_1 10326 4.09 2.02 4 1 7 6 0.02 1
identity_fusion_2 10326 3.72 1.98 4 1 7 6 0.02 1
identity_fusion_3 10325 3.70 2.07 4 1 7 6 0.02 2
identity_fusion_4 10325 3.50 2.08 3 1 7 6 0.02 2
identity_fusion_5 10325 3.20 2.03 3 1 7 6 0.02 2
identity_fusion_6 10324 3.61 2.18 4 1 7 6 0.02 3
identity_fusion_7 10324 3.29 2.04 3 1 7 6 0.02 3
ingroup_superiority_1 10325 3.29 2.09 3 1 7 6 0.02 2
ingroup_superiority_2 10325 3.39 2.10 3 1 7 6 0.02 2
ingroup_superiority_3 10325 3.51 2.06 3 1 7 6 0.02 2
ingroup_superiority_4 10324 4.15 2.09 4 1 7 6 0.02 3
collective_relative_deprivation_1 10327 5.08 1.83 6 1 7 6 0.02 0
collective_relative_deprivation_2 10327 2.90 1.90 3 1 7 6 0.02 0
collective_relative_deprivation_3 10327 2.94 1.95 2 1 7 6 0.02 0
collective_relative_deprivation_4 10327 2.85 1.95 2 1 7 6 0.02 0
collective_relative_deprivation_5 10327 2.75 1.92 2 1 7 6 0.02 0
collective_relative_deprivation_6 10327 3.16 2.05 3 1 7 6 0.02 0
perceived_discrimination_1 10327 2.22 1.47 2 1 6 5 0.01 0
perceived_discrimination_2 10327 2.08 1.44 1 1 6 5 0.01 0
perceived_discrimination_3 10327 1.92 1.39 1 1 6 5 0.01 0
perceived_discrimination_4 10327 1.97 1.46 1 1 6 5 0.01 0
perceived_discrimination_5 10327 2.08 1.48 1 1 6 5 0.01 0
perceived_discrimination_6 10327 2.09 1.50 1 1 6 5 0.01 0
perceived_discrimination_7 10327 1.97 1.45 1 1 6 5 0.01 0
perceived_discrimination_8 10327 2.01 1.46 1 1 6 5 0.01 0
perceived_discrimination_9 10327 2.13 1.51 1 1 6 5 0.01 0
perceived_discrimination_10 10327 1.90 1.41 1 1 6 5 0.01 0
perceived_discrimination_11 10327 1.97 1.45 1 1 6 5 0.01 0
activist_intent_1 10327 3.28 1.99 3 1 7 6 0.02 0
activist_intent_2 10327 3.02 1.98 3 1 7 6 0.02 0
activist_intent_3 10327 2.99 1.96 3 1 7 6 0.02 0
activist_intent_4 10327 2.97 2.01 2 1 7 6 0.02 0
moral_neutralization_1 10318 4.11 1.91 4 1 7 6 0.02 9
moral_neutralization_2 10310 2.24 1.67 1 1 7 6 0.02 17
moral_neutralization_3 10310 2.14 1.79 1 1 7 6 0.02 17
moral_neutralization_4 10304 1.75 1.45 1 1 7 6 0.01 23
moral_neutralization_5 10302 2.09 1.62 1 1 7 6 0.02 25
moral_neutralization_6 10299 2.15 1.65 1 1 7 6 0.02 28
moral_neutralization_7 10302 2.20 1.76 1 1 7 6 0.02 25
moral_neutralization_8 10304 2.86 1.91 2 1 7 6 0.02 23
moral_neutralization_9 10299 1.90 1.58 1 1 7 6 0.02 28
moral_neutralization_10 10303 1.89 1.57 1 1 7 6 0.02 24
moral_neutralization_11 10300 1.98 1.61 1 1 7 6 0.02 27
moral_neutralization_12 10298 2.11 1.60 1 1 7 6 0.02 29
moral_neutralization_13 10297 1.73 1.47 1 1 7 6 0.01 30
moral_neutralization_14 10297 1.84 1.53 1 1 7 6 0.02 30
moral_neutralization_15 10289 1.73 1.47 1 1 7 6 0.01 38
moral_neutralization_16 10292 1.93 1.57 1 1 7 6 0.02 35
radical_attitudes_1 10327 1.89 1.77 1 1 8 7 0.02 0
radical_attitudes_2 10327 1.91 1.76 1 1 8 7 0.02 0
radical_attitudes_3 10327 1.92 1.82 1 1 8 7 0.02 0
anger_1 10321 3.58 1.38 4 1 6 5 0.01 6
anger_2 10314 3.54 1.40 4 1 6 5 0.01 13
anger_3 10309 3.50 1.51 4 1 6 5 0.01 18
anger_4 10307 2.70 1.35 2 1 6 5 0.01 20
anger_5 10307 3.19 1.41 3 1 6 5 0.01 20
anger_6 10307 1.99 1.38 1 1 6 5 0.01 20
anger_7 10305 2.62 1.42 2 1 6 5 0.01 22
anger_8 10306 2.31 1.37 2 1 6 5 0.01 21
anger_9 10307 2.49 1.39 2 1 6 5 0.01 20
anger_10 10306 2.46 1.42 2 1 6 5 0.01 21
positive_affect_1 10327 3.55 1.01 4 1 5 4 0.01 0
positive_affect_2 10327 3.73 0.97 4 1 5 4 0.01 0
positive_affect_3 10327 3.63 1.04 4 1 5 4 0.01 0
positive_affect_4 10327 3.49 1.02 4 1 5 4 0.01 0
positive_affect_5 10327 3.67 1.03 4 1 5 4 0.01 0
negative_affect_6 10327 2.38 1.12 2 1 5 4 0.01 0
negative_affect_7 10327 2.51 1.14 2 1 5 4 0.01 0
negative_affect_8 10327 2.38 1.14 2 1 5 4 0.01 0
negative_affect_9 10327 2.16 1.17 2 1 5 4 0.01 0
negative_affect_10 10327 1.84 1.11 1 1 5 4 0.01 0

Meta-Analysis

We could tihnk about testing a specific model across countries in a meta-analytic manner!

Measurement Invariance

Functions

Fit Functions

Code
source("../scripts/functions/fit-cfa-model-new.R")

Viz Functions

Code
# Loop through the cfa_individual_fit_results list
extract_fit_statistics_from_results <- function(fit_results) {
  
  # Initialize empty data frame
  result_individual_fit_df <- data.frame(
    cfa_model = character(),
    country = character(),
    CFI = numeric(),
    TLI = numeric(),
    RMSEA = numeric(),
    stringsAsFactors = FALSE
  )
  
  for (cfa_model_name in names(fit_results)) {
    cfa_model <- fit_results[[cfa_model_name]]
    
    for (country_name in names(cfa_model)) {
      
      # get fit statistic
      fit_table <- cfa_model[[country_name]]$fit_table
      
      # Extract the values from the fit_table
      CFI <- ifelse(fit_table$CFI > 1, 1, fit_table$CFI)
      TLI <- ifelse(fit_table$TLI > 1, 1, fit_table$TLI)
      RMSEA <- ifelse(fit_table$RMSEA < 0, 0, fit_table$RMSEA)
      
      # Create a new row with the extracted values
      new_row <- data.frame(
        cfa_model = cfa_model_name,
        country = country_name,
        CFI = CFI,
        TLI = TLI,
        RMSEA = RMSEA,
        stringsAsFactors = FALSE
      )
      
      # Append the new row to the result_df
      result_individual_fit_df <- rbind(result_individual_fit_df, new_row)
    }
  }
  return(result_individual_fit_df)
}

# for testing
# df <- result_individual_fit_df
# cfa_model <- unique_models[1]
Code
# for testing
# df = result_fit_individual_df
# cfa_model = unique_models[1]
# Function to plot histograms for a given cfa_model
visualize_individual_fit_statistics <- function(df_fit_statistics, cfa_model) {
  
  # Set thresholds for acceptable and good fit
  cfi_good <- 0.95
  cfi_acceptable <- 0.90
  tli_good <- 0.95
  tli_acceptable <- 0.90
  rmsea_good <- 0.06
  rmsea_acceptable <- 0.08
  
  # select relevant data
  df_filtered <- df_fit_statistics %>% filter(cfa_model == !!cfa_model)
  
  # if there is more than one row (multiple countries) make the plots
  if (nrow(df_filtered) > 1) {
    
    # make minimum and maximum for plots
    CFI_min <- min(df_filtered$CFI, na.rm = TRUE)
    CFI_max_plot <- 1.01
    if (CFI_min > 0.5) {
      CFI_min_plot = 0.5
    } else {
      CFI_min_plot = CFI_min
    }
    
    TLI_min <- min(df_filtered$TLI, na.rm = TRUE)
    TLI_max_plot <- 1.01
    if (TLI_min > 0.5) {
      TLI_min_plot = 0.5
    } else {
      TLI_min_plot = TLI_min
    }
    
    RMSEA_max <- max(df_filtered$RMSEA, na.rm = TRUE)
    RMSEA_min_plot <- 0
    if (RMSEA_max < 0.2) {
      RMSEA_max_plot = 0.2
    } else {
      RMSEA_max_plot = RMSEA_max
    }
    
    # Create a list to store the plots
    plots <- list()
    
    # CFI Histogram
    plots[['CFI']] <- ggplot(df_filtered, aes(x = CFI)) +
      geom_histogram(binwidth = 0.02, fill = "lightblue", color = "black") +
      geom_vline(xintercept = c(cfi_acceptable, cfi_good), linetype = "dashed", color = "red") +
      ggtitle(paste("CFI")) +
      coord_cartesian(xlim = c(CFI_min_plot, CFI_max_plot)) +
      # scale_x_continuous(limits = c(CFI_min_plot, CFI_max_plot)) +
      theme_minimal()
    
    # TLI Histogram
    plots[['TLI']] <- ggplot(df_filtered, aes(x = TLI)) +
      geom_histogram(binwidth = 0.02, fill = "lightgreen", color = "black") +
      geom_vline(xintercept = c(tli_acceptable, tli_good), linetype = "dashed", color = "red") +
      ggtitle(paste("TLI")) +
      coord_cartesian(xlim = c(TLI_min_plot, TLI_max_plot)) +
      # scale_x_continuous(limits = c(TLI_min_plot, TLI_max_plot)) +
      theme_minimal()
    
    # RMSEA Histogram
    plots[['RMSEA']] <- ggplot(df_filtered, aes(x = RMSEA)) +
      geom_histogram(binwidth = 0.01, fill = "lightcoral", color = "black") +
      geom_vline(xintercept = c(rmsea_acceptable, rmsea_good), linetype = "dashed", color = "red") +
      ggtitle(paste("RMSEA")) +
      coord_cartesian(xlim = c(RMSEA_min_plot, RMSEA_max_plot)) +
      # scale_x_continuous(limits = c(RMSEA_min_plot, RMSEA_max_plot)) +
      theme_minimal()
    
    # Combine the plots into a single row
    grid.arrange(plots[['CFI']], plots[['TLI']], plots[['RMSEA']], ncol = 3)
    
    # Identify and print countries with bad fit
    below_acceptable <- df_filtered %>%
      mutate(
        below_cfi = ifelse(is.na(CFI), paste(country, "(CFI not converged)"), 
                           ifelse(CFI < cfi_acceptable, country, NA)),
        below_tli = ifelse(is.na(TLI), paste(country, "(TLI not converged)"), 
                           ifelse(TLI < tli_acceptable, country, NA)),
        above_rmsea = ifelse(is.na(RMSEA), paste(country, "(RMSEA not converged)"), 
                             ifelse(RMSEA > rmsea_acceptable, country, NA))
      ) %>%
      select(country, below_cfi, below_tli, above_rmsea) %>%
      filter(!is.na(below_cfi) | !is.na(below_tli) | !is.na(above_rmsea))
    
    # Convert the below_acceptable table to a grob for printing below the plots
    if (nrow(below_acceptable) > 0) {
      below_acceptable %>% 
        kable(., caption = paste("Fit Statistics below acceptable for", cfa_model)) %>%
        kable_styling(full_width = F, latex_options = c("hold_position", "scale-down")) %>% 
        print()
    } else {
      cat("All countries are within acceptable fit for", cfa_model, "<br>")
    }
  } else {
    df_filtered %>% 
      kable(., caption = paste("Fit statistics for", cfa_model)) %>%
      kable_styling(full_width = F, latex_options = c("hold_position", "scale-down")) %>% 
      print()
    
  }
}

Compare Functions

Code
# Define the function
run_anova_with_indices <- function(unique_models, fit_configural, fit_metric, fit_scalar) {
  results_list <- list()

  pb <- progress::progress_bar$new(
    format = "  Running ANOVA [:bar] :percent in :elapsed",
    total  = length(unique_models), clear = FALSE, width = 60
  )

  get_fit <- function(x, m) x[[m]]$Grouped$fit_object
  get_tbl <- function(x, m) x[[m]]$Grouped$fit_table

  for (model in unique_models) {
    cf_fit <- get_fit(fit_configural, model)
    mt_fit <- get_fit(fit_metric,    model)
    sc_fit <- get_fit(fit_scalar,    model)

    # pull indices from the precomputed tables (keeps behavior identical)
    cf_tbl <- get_tbl(fit_configural, model)
    mt_tbl <- get_tbl(fit_metric,     model)
    sc_tbl <- get_tbl(fit_scalar,     model)

    # Build the indices columns (clamp like your original pipeline did downstream)
    TLI   <- c(cf_tbl$TLI,   mt_tbl$TLI,   sc_tbl$TLI)
    CFI   <- c(cf_tbl$CFI,   mt_tbl$CFI,   sc_tbl$CFI)
    RMSEA <- c(cf_tbl$RMSEA, mt_tbl$RMSEA, sc_tbl$RMSEA)

    # If any fit is NULL (e.g., didn’t converge), return NA LRT rows but keep indices
    if (is.null(cf_fit) || is.null(mt_fit) || is.null(sc_fit)) {
      anova_tbl <- tibble::tibble(
        name        = c("configural","metric","scalar"),
        Df          = NA_real_,
        Chisq       = NA_real_,
        `Df diff`   = NA_real_,
        `Chisq diff`= NA_real_,
        `Pr(>Chisq)`= NA_real_,
        TLI   = round(pmin(TLI, 1), 3),
        CFI   = round(pmin(CFI, 1), 3),
        rmsea = round(pmax(RMSEA, 0), 3)
      )
    } else {
      # Use lavaan's LRT explicitly (avoids the 'anova' dispatch issue)
      lrt <- as.data.frame(lavaan::lavTestLRT(cf_fit, mt_fit, sc_fit))

      # Make column names match your original cbind() output as closely as possible
      # lavTestLRT typically gives: Df, Chisq, Df.diff, Chisq.diff, Pr(>Chisq)
      # Normalize names in case of slight variations
      names(lrt) <- sub("^df\\b", "Df", names(lrt), ignore.case = TRUE)
      names(lrt) <- sub("Chisq\\.diff", "Chisq diff", names(lrt), ignore.case = TRUE)
      names(lrt) <- sub("Df\\.diff",    "Df diff",    names(lrt), ignore.case = TRUE)

      anova_tbl <-
        cbind(lrt,
              TLI   = round(pmin(TLI, 1), 3),
              CFI   = round(pmin(CFI, 1), 3),
              rmsea = round(pmax(RMSEA, 0), 3)) |>
        tibble::as_tibble() |>
        dplyr::mutate(name = c("configural","metric","scalar")) |>
        dplyr::select(name, dplyr::everything())
    }

    results_list[[model]] <- anova_tbl
    pb$tick()
  }

  results_list
}

# summary(fit_configural[[model]]$Grouped_Countries$overall)

Model Structure

Code
#| label: cfa_model_lavaan
# cfa_model_intent <- ' 
#   violent_intent =~ violent_intent_1 + violent_intent_2 + violent_intent_3 + violent_intent_4 + violent_intent_6
#   peaceful_intent =~ peaceful_intent_7 + peaceful_intent_8 + peaceful_intent_9 + peaceful_intent_10 + peaceful_intent_11 + peaceful_intent_12
#   harmonious_passion =~ harmonious_passion_1 + harmonious_passion_3 + harmonious_passion_5 + harmonious_passion_6 + harmonious_passion_8 + harmonious_passion_10
#   obsessive_passion =~ obsessive_passion_2 + obsessive_passion_4 + obsessive_passion_7 + obsessive_passion_9 + obsessive_passion_11 + obsessive_passion_12
#   commitment_passion =~ commitment_passion_13 + commitment_passion_14 + commitment_passion_15 + commitment_passion_16
#   identity_fusion =~ identity_fusion_1 + identity_fusion_2 + identity_fusion_3 + identity_fusion_4 + identity_fusion_5 + identity_fusion_6 + identity_fusion_7
#   ingroup_superiority =~ ingroup_superiority_1 + ingroup_superiority_2 + ingroup_superiority_3 + ingroup_superiority_4
#   collective_relative_deprivation =~ collective_relative_deprivation_1 + collective_relative_deprivation_2 + collective_relative_deprivation_3 + collective_relative_deprivation_4 + collective_relative_deprivation_5 + collective_relative_deprivation_6
#   perceived_discrimination =~ perceived_discrimination_1 + perceived_discrimination_2 + perceived_discrimination_3 + perceived_discrimination_4 + perceived_discrimination_5 + perceived_discrimination_6 + perceived_discrimination_7 + perceived_discrimination_8 + perceived_discrimination_9 + perceived_discrimination_10 + perceived_discrimination_11
#   activist_intent =~ activist_intent_1 + activist_intent_2 + activist_intent_3 + activist_intent_4
#   moral_neutralization =~ moral_neutralization_1 + moral_neutralization_2 + moral_neutralization_3 + moral_neutralization_4 + moral_neutralization_5 + moral_neutralization_6 + moral_neutralization_7 + moral_neutralization_8 + moral_neutralization_9 + moral_neutralization_10 + moral_neutralization_11 + moral_neutralization_12 + moral_neutralization_13 + moral_neutralization_14 + moral_neutralization_15 + moral_neutralization_16
#   radical_attitudes =~ radical_attitudes_1 + radical_attitudes_2 + radical_attitudes_3
#   anger =~ anger_1 + anger_2 + anger_3 + anger_4 + anger_5 + anger_6 + anger_7 + anger_8 + anger_9 + anger_10
#   positive_affect =~ positive_affect_1 + positive_affect_2 + positive_affect_3 + positive_affect_4 + positive_affect_5
#   negative_affect =~ negative_affect_6 + negative_affect_7 + negative_affect_8 + negative_affect_9 + negative_affect_10
# '

cfa_models <- list(
  cfa_model_intent = 
    'violent_intent =~ violent_intent_1 + violent_intent_2 + violent_intent_3 + violent_intent_4 + violent_intent_6
  peaceful_intent =~ peaceful_intent_7 + peaceful_intent_8 + peaceful_intent_9 + peaceful_intent_10 + peaceful_intent_11 + peaceful_intent_12',
  cfa_model_passion = 
    'harmonious_passion =~ harmonious_passion_1 + harmonious_passion_3 + harmonious_passion_5 + harmonious_passion_6 + harmonious_passion_8 + harmonious_passion_10
  obsessive_passion =~ obsessive_passion_2 + obsessive_passion_4 + obsessive_passion_7 + obsessive_passion_9 + obsessive_passion_11 + obsessive_passion_12
  commitment_passion =~ commitment_passion_13 + commitment_passion_14 + commitment_passion_15 + commitment_passion_16
  commitment_passion_14 ~~ commitment_passion_15',
  cfa_model_passion_one_factor = 
    'passion =~ harmonious_passion_1 + harmonious_passion_3 + harmonious_passion_5 + harmonious_passion_6 + harmonious_passion_8 + harmonious_passion_10 + obsessive_passion_2 + obsessive_passion_4 + obsessive_passion_7 + obsessive_passion_9 + obsessive_passion_11 + obsessive_passion_12 + commitment_passion_13 + commitment_passion_14 + commitment_passion_15 + commitment_passion_16
  commitment_passion_14 ~~ commitment_passion_15',
  cfa_model_passion_core = 
    'harmonious_passion =~ harmonious_passion_1 + harmonious_passion_3 + harmonious_passion_5 + harmonious_passion_6 + harmonious_passion_8 + harmonious_passion_10
  obsessive_passion =~ obsessive_passion_2 + obsessive_passion_4 + obsessive_passion_7 + obsessive_passion_9 + obsessive_passion_11 + obsessive_passion_12',
  cfa_model_passion_core_one_factor = 
    'passion =~ harmonious_passion_1 + harmonious_passion_3 + harmonious_passion_5 + harmonious_passion_6 + harmonious_passion_8 + harmonious_passion_10 + obsessive_passion_2 + obsessive_passion_4 + obsessive_passion_7 + obsessive_passion_9 + obsessive_passion_11 + obsessive_passion_12',
  cfa_model_commitment_passion = 
    'commitment_passion =~ commitment_passion_13 + commitment_passion_14 + commitment_passion_15 + commitment_passion_16
  commitment_passion_14 ~~ commitment_passion_15',
  cfa_model_identity_fusion = 
    'identity_fusion =~ identity_fusion_1 + identity_fusion_2 + identity_fusion_3 + identity_fusion_4 + identity_fusion_5 + identity_fusion_6 + identity_fusion_7',
  cfa_model_ingroup_superiority = 
    'ingroup_superiority =~ ingroup_superiority_1 + ingroup_superiority_2 + ingroup_superiority_3 + ingroup_superiority_4',
  cfa_model_collective_relative_deprivation = 
    'collective_relative_deprivation =~ collective_relative_deprivation_1 + collective_relative_deprivation_2 + collective_relative_deprivation_3 + collective_relative_deprivation_4 + collective_relative_deprivation_5 + collective_relative_deprivation_6',
  cfa_model_perceived_discrimination = 
    'perceived_discrimination =~ perceived_discrimination_1 + perceived_discrimination_2 + perceived_discrimination_3 + perceived_discrimination_4 + perceived_discrimination_5 + perceived_discrimination_6 + perceived_discrimination_7 + perceived_discrimination_8 + perceived_discrimination_9 + perceived_discrimination_10 + perceived_discrimination_11',
  cfa_model_activist_intent = 
    'activist_intent =~ activist_intent_1 + activist_intent_2 + activist_intent_3 + activist_intent_4',
  cfa_model_moral_neutralization = 
    'moral_neutralization =~ moral_neutralization_1 + moral_neutralization_2 + moral_neutralization_3 + moral_neutralization_4 + moral_neutralization_5 + moral_neutralization_6 + moral_neutralization_7 + moral_neutralization_8 + moral_neutralization_9 + moral_neutralization_10 + moral_neutralization_11 + moral_neutralization_12 + moral_neutralization_13 + moral_neutralization_14 + moral_neutralization_15 + moral_neutralization_16',
  cfa_model_radical_attitudes = 
    'radical_attitudes =~ radical_attitudes_1 + radical_attitudes_2 + radical_attitudes_3',
  cfa_model_anger =
    'anger =~ anger_1 + anger_2 + anger_3 + anger_4 + anger_5 + anger_6 + anger_7 + anger_8 + anger_9 + anger_10
  anger_9   ~~  anger_10
  anger_1   ~~  anger_2',
  cfa_model_anger_one_fac =
    'anger =~ anger_1 + anger_2 + anger_3 + anger_5 + anger_6 + anger_7
  anger_1   ~~  anger_2
  anger_6   ~~  anger_7',
  cfa_model_affect =
    'positive_affect =~ positive_affect_1 + positive_affect_2 + positive_affect_3 + positive_affect_4 + positive_affect_5
  negative_affect =~ negative_affect_6 + negative_affect_7 + negative_affect_8 + negative_affect_9 + negative_affect_10'
)
Code
# Get the unique list of countries for all analyses
unique_countries <- unique(dt_cfa$country)
unique_models <- unique(names(cfa_models))

Invariance

Configurial Invariance

Configural invariance is the most basic level of invariance, where the same factor structure is assumed across groups, but the parameters (loadings, intercepts, etc.) are allowed to vary. Individual model fit can be found in the model comparison section.

Code
# model 
fit_configural <- fit_cfa_model(
  model_list = cfa_models,
  loop_list = NULL,
  dt = dt_cfa,
  print_output = FALSE,
  group_var = "country",
  model_fixed = NULL
)
Code
# build a tidy summary once
# result_fit_configural_df <- collect_fit_summary(fit_configural)
# 
# # Use the model names present in the results/summaries
# unique_models <- unique(result_fit_configural_df$model)
# 
# for (model in unique_models) {
# 
#   node <- get_result_node(fit_configural, model, tag = "Grouped")
#   if (is.null(node)) next
# 
#   mi_tbl <- get_mod_indices_tbl(node)
#   show_output <- should_show(node, mi_tbl)
# 
#   if (show_output) {
#     h4(paste("Model:", model))
# 
#     if (!isTRUE(node$converged)) {
#       cat("This model did <b>not converge</b> — review identification, starting values, or data quality.<br>")
#     } else if (nrow(mi_tbl) >= 1) {
#       cat("This model shows <b>inadequate fit</b> (low CFI/TLI or flagged MIs).<br>")
#     } else {
#       cat("Model converged without flagged modification indices.<br>")
#     }
# 
#     # Overall fit stats (from the tidy summary)
#     print_fit_table(result_fit_configural_df, model)
# 
#     # Modification indices (if any)
#     if (nrow(mi_tbl) >= 1) {
#       cat("<br>Top modification indices (first 5 shown):<br>")
#       # If you want only the same top-5 as computed, mi_tbl already reflects that.
#       kable(mi_tbl,
#             caption = sprintf("Modification indices for %s", model)) |>
#         kable_styling(full_width = FALSE, latex_options = c("hold_position", "scale_down")) |>
#         print()
#     }
#   }
# }

# extract relevant stats
# result_fit_configural_df <- extract_fit_statistics_from_results(fit_configural)
# 
# # runng above function over models
# for (model in unique_models) {
#   
#   # if there are no modification indices, the model fits well
#   mod_indices <- fit_configural[[model]][["Grouped_Countries"]]$mod_indices %>% as_tibble() %>% mutate_if(is.numeric, round, 2)
#   # only show output if model does not fit (there are no rows)
#   show_output <- nrow(mod_indices) > 1
#   
#   if (show_output) {
#     cat("<br><h4>Model:", model, "</h4><br>")
#     cat("This model has inadequate fit.<br>")
#     cat("Overall fit statistics:<br>")
#     
#     visualize_individual_fit_statistics(
#       df = result_fit_configural_df, 
#       cfa_model = model
#     )
#     
#     cat("Modification indices:<br>")
#     mod_indices %>% 
#       kable(., caption = paste("Modification indices for", model)) %>%
#       kable_styling(full_width = F, latex_options = c("hold_position", "scale-down")) %>% 
#       print()
#   }
# }

# fit_configural <- sem(cfa_model_intent, data = dt_cfa, group = "country", optim.method = "nlminb")
# fit_configural <- sem(cfa_model_intent, data = dt_cfa, group = "country", optim.method = "em", verbose = T,
#            em.iter.max = 20000, em.fx.tol = 1e-08, em.dx.tol = 1e-04)

Metric Invariance

Metric invariance constrains the factor loadings to be equal across groups. Individual model fit can be found in the model comparison section.

Code
fit_metric <- fit_cfa_model(
  model_list = cfa_models,
  loop_list = NULL,
  dt = dt_cfa,
  print_output = FALSE,
  group_var = "country",
  model_fixed = "loadings"
)
Code
# extract relevant stats
# result_fit_metric_df <- extract_fit_statistics_from_results(fit_metric)
# 
# # runng above function over models
# for (model in unique_models) {
#   
#   # if there are no modification indices, the model fits well
#   mod_indices <- fit_metric[[model]][["Grouped_Countries"]]$mod_indices %>% as_tibble() %>% mutate_if(is.numeric, round, 2)
#   # only show output if model does not fit (there are no rows)
#   show_output <- nrow(mod_indices) > 1
#   
#   if (show_output) {
#     cat("<br><h4>Model:", model, "</h4><br>")
#     cat("This model has inadequate fit.<br>")
#     cat("Overall fit statistics:<br>")
#     
#     visualize_individual_fit_statistics(
#       df = result_fit_metric_df, 
#       cfa_model = model
#     )
#     
#     cat("Modification indices:<br>")
#     mod_indices %>% 
#       kable(., caption = paste("Modification indices for", model)) %>%
#       kable_styling(full_width = F, latex_options = c("hold_position", "scale-down")) %>% 
#       print()
#   }
# }

Scalar Invariance

Scalar invariance involves constraining both the factor loadings and the intercepts to be equal across groups. Individual model fit can be found in the model comparison section.

Code
fit_scalar <- fit_cfa_model(
  model_list = cfa_models,
  loop_list = NULL,
  dt = dt_cfa,
  print_output = FALSE,
  group_var = "country",
  model_fixed = c("loadings", "intercepts")
)
Code
# extract relevant stats
# result_fit_scalar_df <- extract_fit_statistics_from_results(fit_scalar)
# 
# # runng above function over models
# # runng above function over models
# for (model in unique_models) {
#   
#   # if there are no modification indices, the model fits well
#   mod_indices <- fit_scalar[[model]][["Grouped_Countries"]]$mod_indices %>% as_tibble() %>% mutate_if(is.numeric, round, 2)
#   # only show output if model does not fit (there are no rows)
#   show_output <- nrow(mod_indices) > 1
#   
#   if (show_output) {
#     cat("<br><h4>Model:", model, "</h4><br>")
#     cat("This model has inadequate fit.<br>")
#     cat("Overall fit statistics:<br>")
#     
#     visualize_individual_fit_statistics(
#       df = result_fit_scalar_df, 
#       cfa_model = model
#     )
#     
#     cat("Modification indices:<br>")
#     mod_indices %>% 
#       kable(., caption = paste("Modification indices for", model)) %>%
#       kable_styling(full_width = F, latex_options = c("hold_position", "scale-down")) %>% 
#       print()
#   }
# }

Model Comparison

We can compare models using a chi-square difference test or by comparing fit indices like CFI and RMSEA. Our main comparison of interest is configural to metric invariance. This is because we are mainly interested in whether the relationship between the items and the underlying latent variable as well as the relationships between the variables remain constant. Intercepts and mean comparisons are not the focus of this analysis. Our main statistic of interference are TLI and CFI. RMSEA and the Chi-quared will most likely be heavily influences by the large sample size.

Code
anova_results <- run_anova_with_indices(
  unique_models, 
  fit_configural, 
  fit_metric, 
  fit_scalar
)
Code
# Initialize the model_descriptions list with unique_models as names
model_descriptions <- setNames(vector("list", length(unique_models)), unique_models)

# describe the individual models
model_descriptions[["cfa_model_intent"]] <- "Cross-cultural performance of the violent and peaceful intent scales is satisfactory. TLI and CFI remain stable despite increasing model restrictions."
model_descriptions[["cfa_model_radical_attitudes"]] <- "A just-identified three-item, one-factor CFA shows strong item loadings across countries. Due to the just-identified nature, TLI and CFI at 1. We are using the default lavaan <i>marker method</i> which fixes the first loading to 1."
model_descriptions[["cfa_model_passion"]] <- "Initial model fit is suboptimal; cross-country data is not ideal.<br><b>After checking:</b><br>We adjusted one residual covariance (commitment_passion_14 ~~ commitment_passion_15) which improved model fit to an acceptable level."
model_descriptions[["cfa_model_passion_core"]] <- "Model fit is acceptable. The commitment scale is often seen as a seperate entity (measuring the presence of passion) which is why we also tested it seperately."
model_descriptions[["cfa_model_commitment_passion"]] <- "Seperating this from the original three factor solution improved model fit substantially. Adjusting one parameter increased model fit dramatically (commitment_passion_14 ~~ commitment_passion_15). Either the original scale or the seperation suggest an acceptable model solution."
model_descriptions[["cfa_model_identity_fusion"]] <- "Configural and metric invariance are supported, but scalar invariance is weak. Avoid comparing intercepts or means across countries."
model_descriptions[["cfa_model_ingroup_superiority"]] <- "The ingroup superiority scale demonstrates good model fit even under scalar invariance."
model_descriptions[["cfa_model_collective_relative_deprivation"]] <- "Configural and metric invariance are defensible; scalar invariance is uncertain. Avoid comparing intercepts or means across countries."
model_descriptions[["cfa_model_perceived_discrimination"]] <- "Good model fit. TLI and CFI fit indices remain above thresholds across restrictive models."
model_descriptions[["cfa_model_activist_intent"]] <- "Good model fit across countries, including under more restrictive conditions."
model_descriptions[["cfa_model_moral_neutralization"]] <- "Model fit is poor; data from Argentina, Denmark, and Colombia did not converge.<br><b>After checking:</b> We cannot assume scalar invariance. However, both configural and metric invariance are in line with [previous implementations of invariance](https://www.researchgate.net/publication/279472182_Are_Moral_Disengagement_Neutralization_Techniques_and_Self-Serving_Cognitive_Distortions_the_Same_Developing_a_Unified_Scale_of_Moral_Neutralization_of_Aggression). As not our core variable, we use it as a variable in the random forest but cross-cultural conclusions regarding the variable should be drawn cautiously."
model_descriptions[["cfa_model_anger"]] <- "Model fit is bad. While all countries converge, the fit within each country is not satisfactory.<br><b>After checking:</b><br>Adjusting the following residual covariances improved model fit (anger_9  ~~  anger_10; anger_1   ~~  anger_2)."
model_descriptions[["cfa_model_anger_one_fac"]] <- "This is the original scale without the adjusted items. Adjusting two residual covariances resulted in satisfactory model fit (anger_1   ~~  anger_2; anger_6    ~~  anger_7)."
model_descriptions[["cfa_model_affect"]] <- "Configural and metric invariance are marginal; scalar invariance is not supported. Fit within individual countries is generally satisfactory, but conclusions should be cautious."


# Print results for each model
for (model in unique_models) {
  cat("<br><h4>Model:", model, "</h4><br>")
  
  # Print the description for the current model
  description <- model_descriptions[[model]]
  if (!is.null(description) && description != "") {
    cat("Description:<br>", description, "<br>")
  } else {
    cat("<p>No description available for this model.</p><br>")
  }
  
  anova_results[[model]] %>% 
    kable(., caption = paste("ANOVA Model comparison for", model)) %>%
    kable_styling(full_width = F, latex_options = c("hold_position", "scale-down")) %>% 
    print()
}

Model: cfa_model_intent


Description:
Cross-cultural performance of the violent and peaceful intent scales is satisfactory. TLI and CFI remain stable despite increasing model restrictions.
ANOVA Model comparison for cfa_model_intent
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 1763 369641.4 379737.4 10275.95 NA NA NA NA 0.904 0.925 0.138
metric 2123 369861.6 377350.4 11216.19 940.2403 0.0799940 360 0 0.915 0.920 0.130
scalar 2483 372778.0 377659.4 14852.55 3636.3641 0.1900858 360 0 0.901 0.891 0.141

Model: cfa_model_passion


Description:
Initial model fit is suboptimal; cross-country data is not ideal.
After checking:
We adjusted one residual covariance (commitment_passion_14 ~~ commitment_passion_15) which improved model fit to an acceptable level.
ANOVA Model comparison for cfa_model_passion
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 4100 521264.4 536705.4 19911.66 NA NA NA NA 0.878 0.898 0.124
metric 4620 522073.0 533748.0 21760.31 1848.654 0.1007185 520 0 0.882 0.890 0.121
scalar 5140 528092.8 536001.6 28820.04 7059.734 0.2234515 520 0 0.854 0.848 0.135

Model: cfa_model_passion_one_factor


No description available for this model.


ANOVA Model comparison for cfa_model_passion_one_factor
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 4223 530564.0 545114.2 29457.29 NA NA NA NA 0.811 0.837 0.154
metric 4823 532143.3 542348.0 32236.57 2779.279 0.1200841 600 0 0.820 0.823 0.150
scalar 5423 541164.5 547023.7 42457.76 10221.188 0.2523153 600 0 0.784 0.761 0.165

Model: cfa_model_passion_core


Description:
Model fit is acceptable. The commitment scale is often seen as a seperate entity (measuring the presence of passion) which is why we also tested it seperately.
ANOVA Model comparison for cfa_model_passion_core
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 2173 397469.3 408456.2 10887.57 NA NA NA NA 0.900 0.920 0.126
metric 2573 397758.6 405848.5 11976.86 1089.290 0.0827134 400 0 0.909 0.913 0.120
scalar 2973 401469.3 406662.2 16487.57 4510.703 0.2019914 400 0 0.887 0.875 0.134

Model: cfa_model_passion_core_one_factor


No description available for this model.


ANOVA Model comparison for cfa_model_passion_core_one_factor
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 2214 406104.2 416794.2 19604.46 NA NA NA NA 0.804 0.840 0.177
metric 2654 407268.3 414771.6 21648.56 2044.095 0.1203077 440 0 0.821 0.825 0.169
scalar 3094 413581.0 417897.5 28841.22 7192.660 0.2468402 440 0 0.792 0.763 0.182

Model: cfa_model_commitment_passion


Description:
Seperating this from the original three factor solution improved model fit substantially. Adjusting one parameter increased model fit dramatically (commitment_passion_14 ~~ commitment_passion_15). Either the original scale or the seperation suggest an acceptable model solution.
ANOVA Model comparison for cfa_model_commitment_passion
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 41 138249.3 142109.5 73.28548 NA NA NA NA 0.992 0.999 0.056
metric 161 138819.9 141811.0 883.91621 810.6307 0.1511674 120 0 0.956 0.971 0.134
scalar 281 140828.9 142951.0 3132.96336 2249.0471 0.2654165 120 0 0.901 0.887 0.201

Model: cfa_model_identity_fusion


Description:
Configural and metric invariance are supported, but scalar invariance is weak. Avoid comparing intercepts or means across countries.
ANOVA Model comparison for cfa_model_identity_fusion
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 574 228001.9 234237.6 5158.894 NA NA NA NA 0.896 0.931 0.178
metric 814 228286.0 232783.6 5923.020 764.1262 0.0931190 240 0 0.919 0.923 0.158
scalar 1054 231891.8 234651.1 10008.741 4085.7204 0.2522372 240 0 0.890 0.865 0.184

Model: cfa_model_ingroup_superiority


Description:
The ingroup superiority scale demonstrates good model fit even under scalar invariance.
ANOVA Model comparison for cfa_model_ingroup_superiority
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 82 144400.8 147964.1 445.8304 NA NA NA NA 0.957 0.986 0.133
metric 202 144379.4 147073.5 664.3759 218.5456 0.0571023 120 1e-07 0.978 0.982 0.095
scalar 322 145517.8 147342.9 2042.8221 1378.4462 0.2040575 120 0e+00 0.948 0.932 0.146

Model: cfa_model_collective_relative_deprivation


Description:
Configural and metric invariance are defensible; scalar invariance is uncertain. Avoid comparing intercepts or means across countries.
ANOVA Model comparison for cfa_model_collective_relative_deprivation
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 369 216139.8 221484.7 1861.817 NA NA NA NA 0.934 0.961 0.127
metric 569 216291.6 220188.0 2413.636 551.8185 0.0835698 200 0 0.947 0.951 0.113
scalar 769 218198.5 220646.5 4720.559 2306.9235 0.2045098 200 0 0.917 0.896 0.143

Model: cfa_model_perceived_discrimination


Description:
Good model fit. TLI and CFI fit indices remain above thresholds across restrictive models.
ANOVA Model comparison for cfa_model_perceived_discrimination
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 1804 271619.0 281418.1 10591.81 NA NA NA NA 0.917 0.933 0.139
metric 2204 272045.5 278947.6 11818.31 1226.501 0.0905725 400 0 0.926 0.927 0.132
scalar 2604 273087.7 277092.8 13660.51 1842.194 0.1196429 400 0 0.927 0.916 0.130

Model: cfa_model_activist_intent


Description:
Good model fit across countries, including under more restrictive conditions.
ANOVA Model comparison for cfa_model_activist_intent
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 82 137553.4 141116.7 491.9483 NA NA NA NA 0.961 0.987 0.141
metric 202 137585.1 140279.3 763.6385 271.6902 0.0708424 120 0 0.979 0.982 0.105
scalar 322 138430.2 140255.4 1848.7896 1085.1511 0.1786948 120 0 0.963 0.952 0.137

Model: cfa_model_moral_neutralization


Description:
Model fit is poor; data from Argentina, Denmark, and Colombia did not converge.
After checking: We cannot assume scalar invariance. However, both configural and metric invariance are in line with previous implementations of invariance. As not our core variable, we use it as a variable in the random forest but cross-cultural conclusions regarding the variable should be drawn cautiously.
ANOVA Model comparison for cfa_model_moral_neutralization
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 4264 490548.2 504800.7 17793.51 NA NA NA NA 0.877 0.894 0.112
metric 4864 491448.5 501355.7 19893.79 2100.279 0.0996551 600 0 0.881 0.882 0.111
scalar 5464 497480.8 503042.8 27126.10 7232.310 0.2095298 600 0 0.847 0.830 0.125

Model: cfa_model_radical_attitudes


Description:
A just-identified three-item, one-factor CFA shows strong item loadings across countries. Due to the just-identified nature, TLI and CFI at 1. We are using the default lavaan marker method which fixes the first loading to 1.
ANOVA Model comparison for cfa_model_radical_attitudes
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 0 90231.46 92903.94 0.0000 NA NA NA NA 1.000 1.000 0.000
metric 80 90380.36 92473.45 308.9032 308.9032 0.1065825 80 0.0e+00 0.988 0.992 0.107
scalar 160 90369.62 91883.30 458.1618 149.2587 0.0586269 80 4.3e-06 0.992 0.990 0.086

Model: cfa_model_anger


Description:
Model fit is bad. While all countries converge, the fit within each country is not satisfactory.
After checking:
Adjusting the following residual covariances improved model fit (anger_9 ~~ anger_10; anger_1 ~~ anger_2).
ANOVA Model comparison for cfa_model_anger
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 1353 304745.0 314246.4 6692.273 NA NA NA NA 0.878 0.910 0.125
metric 1713 304819.8 311714.1 7487.110 794.8372 0.0692697 360 0 0.895 0.903 0.116
scalar 2073 307090.0 311377.3 10477.335 2990.2245 0.1703634 360 0 0.874 0.859 0.127

Model: cfa_model_anger_one_fac


Description:
This is the original scale without the adjusted items. Adjusting two residual covariances resulted in satisfactory model fit (anger_1 ~~ anger_2; anger_6 ~~ anger_7).
ANOVA Model comparison for cfa_model_anger_one_fac
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 287 195772.9 201711.3 1005.769 NA NA NA NA 0.929 0.967 0.100
metric 487 195778.1 200268.1 1410.923 405.1532 0.0638345 200 0 0.946 0.957 0.087
scalar 687 197240.5 200282.2 3273.396 1862.4739 0.1817161 200 0 0.893 0.881 0.122

Model: cfa_model_affect


Description:
Configural and metric invariance are marginal; scalar invariance is not supported. Fit within individual countries is generally satisfactory, but conclusions should be cautious.
ANOVA Model comparison for cfa_model_affect
name Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) TLI CFI rmsea
configural 1394 267266.0 276471.2 4192.295 NA NA NA NA 0.900 0.925 0.089
metric 1714 267729.2 274616.8 5295.511 1103.216 0.0985759 320 0 0.896 0.904 0.091
scalar 2034 270242.5 274812.6 8448.848 3153.337 0.1874904 320 0 0.844 0.828 0.112

Fit Seperately

Run Analyses

Code
fit_individual_results <- fit_cfa_model(
  model_list = cfa_models,
  loop_list = unique_countries,
  dt = dt_cfa,
  print_output = FALSE,
  group_var = NULL,
  model_fixed = NULL
)

References