3
$\begingroup$

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.

$\endgroup$
0

1 Answer 1

3
$\begingroup$

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.

enter image description here

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))
)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.