When analyzing choice experiments or preference assays—such as behavioral tests where subjects choose between pairs of substances (A vs. B, B vs. C, and A vs. C)—you often need to run pairwise comparisons. A common challenge in R is that standard contingency tables containing structural missing values (NA) break functions like chisq.test(), yielding the error: "all entries of 'x' must be nonnegative and finite".

Furthermore, running multiple tests requires p-value adjustment (e.g., Bonferroni or FDR) to avoid Type I error inflation, and visualizing the significance bars on your charts can become cumbersome. In this tutorial, we will walk through how to conduct pairwise goodness-of-fit or binomial tests, adjust p-values for multiple comparisons, and display the results cleanly using ggplot2.

The Core Problem: Structural Missing Data

In a pairwise design where tests only present two choices at a time, each assay does not include the third option. Trying to cram these into a single 3x3 contingency table introduces structural zeros or NAs, which standard Chi-square contingency tests cannot evaluate properly.

Because each assay is an independent experiment testing whether the choice deviates from a 1:1 (50/50) split, the correct statistical approach is either a one-sample Chi-Square goodness-of-fit test or an exact Binomial test (ideal for small sample counts) evaluated separately per assay, followed by multiple testing correction.

Step 1: Prepare the Data in Tidy Format

Rather than using a matrix filled with NAs, structure your data in a tidy, "long" format using tidyverse tools:

library(tidyverse)

# Prepare assay data in long format
longdata <- tibble(
  assay = c('AB', 'AB', 'AC', 'AC', 'BC', 'BC'),
  liquid = c('A', 'B', 'A', 'C', 'B', 'C'),
  count = c(9, 10, 8, 12, 9, 6)
)

Step 2: Run Pairwise Tests and Adjust P-Values

Using dplyr and purrr, we can group the data by assay, perform a goodness-of-fit test against an expected uniform distribution (equal preference, $p = 0.5$), and apply the Benjamini-Hochberg (FDR) or Bonferroni correction:

test_results <- longdata %>%
  group_by(assay) %>%
  summarise(
    # Use chisq.test with equal expected probabilities (1:1)
    p_value_chisq = chisq.test(count)$p.value,
    # Alternatively, exact binomial test (recommended for small samples)
    p_value_binom = binom.test(count[1], sum(count), p = 0.5)$p.value,
    total_count = sum(count),
    max_count = max(count),
    .groups = 'drop'
  ) %>%
  mutate(
    # Adjust p-values for multiple comparisons across the 3 assays
    p_adj = p.adjust(p_value_chisq, method = "bonferroni"),
    # Generate significance labels
    signif_label = case_when(
      p_adj < 0.001 ~ "***",
      p_adj < 0.01  ~ "**",
      p_adj < 0.05  ~ "*",
      TRUE          ~ "ns"
    )
  )

print(test_results)

Step 3: Plot with ggplot2 and Add Statistical Annotations

To display your findings clearly, combine the column plot with significance annotations. You can position custom text labels directly using geom_text(), or leverage the ggpubr package for pre-formatted comparison brackets.

library(ggplot2)

# Create the dodge bar chart with statistical indicators
ggplot(longdata, aes(x = assay, y = count, fill = liquid)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  # Add significance indicators above each assay pair
  geom_text(
    data = test_results,
    aes(x = assay, y = max_count + 1.5, label = signif_label),
    inherit.aes = FALSE,
    size = 5,
    fontface = "bold"
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title = "Preference Assay Comparisons",
    subtitle = "Pairwise goodness-of-fit tests with Bonferroni correction",
    x = "Assay Pair",
    y = "Choice Count",
    fill = "Substance"
  ) +
  theme_classic()

Key Takeaways

  • Avoid Structural NAs: Do not try to fit incomplete designs into a single $k \times k$ contingency matrix. Each assay is an independent two-choice test.
  • Binomial vs. Chi-Square: When sample sizes per pair are small (counts under 5–10), binom.test() is statistically more reliable than chisq.test() because it calculates exact probabilities rather than asymptotic approximations.
  • Always Correct P-Values: Running three tests simultaneously inflates family-wise error rates. Using p.adjust(..., method = "bonferroni") or "fdr" ensures your statistical conclusions remain rigorous and defensible.