# Chi-square test alternatives when cell expecteds condition violated: Bayesian method

Say I have data collected from an independent multinomial sampling scheme that violates Cochran's cell expecteds condition. Instead of Pearson's $\chi^2$ test of independence, can I use the Bayesian method?

If so, I would like to make sure that I am making the right interpretations. The raceDolls example in the BayesFactor package was useful, but I want to extend this to an nx5 table and make sure I understand the output. Here's an example:

set.seed(21)
#Create fake data. Yes, I could just create a table with those probabilities, but this was more fun.
df <- data.frame(subjectType = sample(LETTERS[1:5], 150, replace=TRUE, prob=c(0.15, 0.1, 0.09, 0.22, 0.44)))
df$probs <- ifelse(df$subjectType == "A", "c(0.05, 0, 0, 0.09, 0.86)",
ifelse(df$subjectType == "B", "c(0.2, 0, 0, 0.14, 0.66)", ifelse(df$subjectType == "C", "c(0, 0.23, 0, 0, 0.77)",
ifelse(df$subjectType == "D", "c(0.1, 0.15, 0.03, 0.09, 0.63)", ifelse(df$subjectType == "E", "c(0.24, 0.14, 0.02, 0.05, 0.55)",NA)))))
splitDf <- split(df, df$subjectType) splitDf <- lapply(splitDf, function(x) cbind(x, cat = eval(parse(text = paste("sample(LETTERS[22:26], nrow(x), replace=T, prob=", x$probs[1], ")")))))
df <- do.call("rbind", splitDf)
tabDf <- table(df$cat, df$subjectType)

#Tests
chisq.test(tabDf)
#fisher.test(tabDf) # common alternative. Assumes fixed totals
friedman.test(tabDf) # another test that's used as an alternative

library(BayesFactor)
contingencyTableBF(tabDf, sampleType = "indepMulti", fixedMargin = "cols")


The output shows $BF_{01}^I = 0.07$, so it 14 times more likely under the null hypothesis of independence than the alternative hypothesis.

Yes, you can use Bayes factors to evaluate the evidence for / against independence in a R×C table (ie. a table with R rows and C columns) even if some cells are 0.

Let's look at the OP's example in more detail.

The simulation fixes the total sample size, N = 150, and sample N subjects independently with the following probabilities:

#>    A    B    C    D    E
#> 0.15 0.10 0.09 0.22 0.44


So the sampling of subjects is very uneven: (in expectation) almost half of the observations come from subject E.

And here is the contingency table (seed = 21):

#>      A  B  C  D  E
#>   V  1  3  0  4 14
#>   W  0  0  2  2  9
#>   X  0  0  0  3  1
#>   Y  2  3  0  0  3
#>   Z 21 10 22 25 25


Qualitatively, E has relatively more Vs and Ws than the other subjects (relatively = with higher frequency). However, the hypothesis being evaluated is not specifically about E vs the rest; it's the overall hypothesis about independence between subjects and categories.

So taken together the uneven subject sampling and the uneven subject-specific category probabilities mean that a sample of N = 150 subjects tends to contain little evidence against the independence hypothesis even thought it's the true hypothesis.

I demonstrate this with a simulation where I generate 1,000 random contingency tables with the given probabilities for increasing larger sample sizes: N = 150, 200, 250, 300. The histograms of the log Bayes factors illustrate how, as N increases, the evidence against the null & in favor of the alternative increases, slowly due to the lack of balance.

PS: For a fixed sample size N, it would help to take similar number of observations from all subjects.

R code to reproduce the simulation:

library("BayesFactor")
library("tidyverse")

p <- diag(
c(0.15, 0.10, 0.09, 0.22, 0.44)
)
dimnames(p) <- list(
c("A", "B", "C", "D", "E"),
c("A", "B", "C", "D", "E")
)
q <- matrix(
c(
0.05, 0.00, 0.00, 0.09, 0.86,
0.20, 0.00, 0.00, 0.14, 0.66,
0.00, 0.23, 0.00, 0.00, 0.77,
0.10, 0.15, 0.03, 0.09, 0.63,
0.24, 0.14, 0.02, 0.05, 0.55
), 5, 5,
byrow = TRUE
)
dimnames(q) <- list(
c("A", "B", "C", "D", "E"),
c("V", "W", "X", "Y", "Z")
)

simulate.data <- function(N, seed = NULL) {
probs.subject <- enframe(diag(p), name = "subject", value = "p")
probs.category_by_subject <- as_tibble(q, rownames = "subject")

probs <-
inner_join(
probs.subject,
probs.category_by_subject,
by = "subject"
) %>%
pivot_longer(
V:Z,
names_to = "category",
values_to = "q"
)
sample <-
probs %>%
slice_sample(
n = N,
weight_by = p * q,
replace = TRUE
)
table(sample$$category, sample$$subject)
}

sim.evidence <- function(N) {
tab <- simulate.data(N)
bf <- contingencyTableBF(tab, sampleType = "jointMulti")
bf@bayesFactor\$bf
}

set.seed(21)

evidence <- tibble(
N = 150 = replicate(1000, sim.evidence(150)),
N = 200 = replicate(1000, sim.evidence(200)),
N = 250 = replicate(1000, sim.evidence(250)),
N = 300 = replicate(1000, sim.evidence(300))
)