This script shows an example of how to perform pairwise chi-square tests for categorical variables in a dataset.

A pairwise chi-square test helps detect associations between distinct pairs of categorical variables. By running multiple chi-square tests on each pair, it allows you to pinpoint which specific pairs exhibit significant relationships.

Step 1: load libraries and sample data

In this example, we will use the sample pq_data dataset from the vivainsights package. We will also use dplyr and purrr for data manipulation and iteration respectively, and optionally you can just load the tidyverse package instead.

library(vivainsights)
library(dplyr)
library(purrr)

sample_data <- pq_data

Step 2: simulating a dataset with categorical variables

The next step is to simulate additional categorical variables for the sample data. In this example, we will create three fake categorical variables: Teams, Regions, and Functions. We will then merge these variables with the sample data.

# Set random seed for reproducibility
set.seed(123)

# Number of unique PersonId in data
n_personid <- length(unique(sample_data$PersonId))

# Create fake categorical variables for each PersonId
cat_data <- data.frame(
  PersonId = unique(sample_data$PersonId),
  Teams = sample(c('Team 1', 'Team 2', 'Team 3'), size = n_personid, replace = TRUE),
  Regions = sample(c('East', 'South', 'West', 'North'), size = n_personid, replace = TRUE),
  Functions = sample(c('HR', 'Finance', 'Operations', 'Sales'), size = n_personid, replace = TRUE)
)

# Merge the datasets
sample_data_merged <- merge(sample_data, cat_data, by = "PersonId")

# Assign categorical variables names to `cat_vars`, alongside existing variables
cat_vars <- c("Teams", "Regions", "Functions", "Organization", "LevelDesignation")

The unit of this test is a person, not a person-week. These sample attributes are constant within each person; repeated weeks must not increase the sample size. The following check rejects changing attributes instead of arbitrarily choosing the first value. For time-varying attributes, choose a documented common snapshot date or a repeated-observation design before using this example.

person_categories <- sample_data_merged |>
  distinct(PersonId, across(all_of(cat_vars)))
if (anyNA(person_categories$PersonId) ||
    anyDuplicated(person_categories$PersonId)) {
  stop("Expected one stable categorical record per person; choose a common snapshot.")
}

Step 3: Perform pairwise chi-square tests for all categorical variables

In the following, we first use combn() to generate all combinations of variable pairs. combn() comes from the utils package, and generates all combinations of a vector of elements of a given size. Here, we set m = 2 to yield pairs.

Next, we use map_dfr() from the purrr package to loop over each combination and perform a chi-square test. The operation is similar to a for loop, but the results are row-bound (similar to rbind() or bind_rows()) and returned as a data frame. In R, it is generally more efficient to use map() functions from the purrr package than to use for loops.

Towards the end of the code, we add a significance level to the results based on the Benjamini-Hochberg adjusted p-value. The significance level is denoted by asterisks, where *** indicates adjusted p < 0.001, ** indicates adjusted p < 0.01, and * indicates adjusted p < 0.05.

Pairs use complete cases, dropping unused levels. If any expected count is below five, we use a seeded Monte Carlo chi-square p-value (9,999 replicates), not the sparse-table asymptotic approximation. Degenerate pairs with fewer than two observed levels in either variable are reported as not testable.

set.seed(124)
# Generate all combinations of variable pairs
cat_var_combinations <- combn(x = cat_vars, m = 2, simplify = FALSE)

# Use `map_dfr()` to loop over each combination
results_df <-
  map_dfr(cat_var_combinations, ~{
  var1 <- .x[1]
  var2 <- .x[2]

  # Create a contingency table
  complete <- complete.cases(person_categories[c(var1, var2)])
  contingency_table <- table(
    droplevels(factor(person_categories[[var1]][complete])),
    droplevels(factor(person_categories[[var2]][complete]))
  )

  if (any(dim(contingency_table) < 2L)) {
    return(tibble(var1, var2, chi2 = NA_real_, p = NA_real_,
                  n = sum(contingency_table), method = "Not testable"))
  }
  expected <- outer(rowSums(contingency_table), colSums(contingency_table)) /
    sum(contingency_table)
  chi_test <- chisq.test(contingency_table,
                        simulate.p.value = any(expected < 5), B = 9999)

  # Return data frame with raw p-values
  tibble(
    var1 = var1,
    var2 = var2,
    chi2 = chi_test$statistic,
    p = chi_test$p.value,
    n = sum(contingency_table),
    method = chi_test$method
  )
})

results_df <- results_df %>%
  mutate(
    p_adjusted = p.adjust(p, method = "BH"),
    significance = case_when(
      is.na(p_adjusted) ~ "Not testable",
      p_adjusted < 0.001 ~ "***",
      p_adjusted < 0.01 ~ "**",
      p_adjusted < 0.05 ~ "*",
      TRUE ~ ""
    )
  )

print(results_df)
## # A tibble: 10 × 8
##    var1         var2             chi2     p     n method p_adjusted significance
##    <chr>        <chr>           <dbl> <dbl> <int> <chr>       <dbl> <chr>
##  1 Teams        Regions          8.08 0.232   300 "Pear…      0.763 ""
##  2 Teams        Functions        3.03 0.805   300 "Pear…      0.909 ""
##  3 Teams        Organization    14.1  0.302   300 "Pear…      0.763 ""
##  4 Teams        LevelDesignati…  5.50 0.481   300 "Pear…      0.909 ""
##  5 Regions      Functions       10.6  0.305   300 "Pear…      0.763 ""
##  6 Regions      Organization    13.3  0.774   300 "Pear…      0.909 ""
##  7 Regions      LevelDesignati…  6.32 0.707   300 "Pear…      0.909 ""
##  8 Functions    Organization    12.1  0.850   300 "Pear…      0.909 ""
##  9 Functions    LevelDesignati…  4.04 0.909   300 "Pear…      0.909 ""
## 10 Organization LevelDesignati… 36.3  0.006   300 "Pear…      0.06  ""

Finally, you can export the results to csv or clipboard using the following code:

# Copy to clipboard
results_df %>% export(method = "clipboard")

# Export to csv
results_df %>% export(method = "csv", path = "chi-square-results", timestamp = FALSE)