Test individual categories in a contingency table for significance

Question

I have asked people which food they prefer:

     choice
group apple orange pizza beer
    A   374     63   216  101
    B   510     65   125   76

Apparently group B prefers fruit and group A prefers pizza and beer, and a chi-square test shows that the overall differences between groups are significant. But how can I test for which individual choice there is a significant difference between groups?

For example, I want to know whether there is a significant difference in the preference for oranges. But I cannot, I believe, just subset the orange choices, because that way I wouldn't consider the total number of participants per group. I mean, a difference between 1 from A and 2 from B will be significant if I have only sampled three people, but not if those are three in a million.

Participants were asked to choose one from the four foods. They could not select multiple answers.

How can I test this?

---

My hunch would be to either add up the non-orange answers and test the resulting 2×2 table with a chi-square test:

     choice
group orange not orange
    A     63        691
    B     65        711

orange <- matrix(c(63, 691, 65, 711), 2, 2, TRUE,
                 list(group = c("A", "B"), choice = c("orange", "not orange"))
                )

chisq.test(orange, correct = FALSE)
# p = .9883

or to calculate the percentage of orange answers in each group, consider the two numbers as counts in a binomial distribution and test that with a binomial test:

a <- 63 / (63 + 691)
b <- 65 / (65 + 711)
all <- 63 + 691 + 65 + 711

binom.test(c(round(a * all / (a + b)), round(b * all / (a + b))))
# p = .9796

# just checkin'
all == sum(c(round(a * all / (a + b)), round(b * all / (a + b))))
[1] TRUE

Or is there a better, maybe more common way?

---

Sample data

food <- c("apple", "orange", "pizza", "beer")
dat <- data.frame(
                  group  = rep(c("A", "B"), c(754, 776)),
                  choice = c(
                             rep(food, c(374, 63, 216, 101)),
                             rep(food, c(510, 65, 125, 76))
                            )
                 )
tab <- table(dat)

---

Explanation of second procedure

We want to compare the orange answers between groups. But if we only look at the orange answers themselves, we disregard the fact that other answers could be given. So instead of comparing the absolute numbers of orange answers, what we do is weigh the absolute number of orange answers by their proportion within all the answers in each group. Or in other words, we test if there is a significant difference between the percentages of orange answers in both groups.

Given this contingency table:

     choice
group orange not orange
    A     63        691
    B     65        711

for group A, the percentage of orange answers is:

a <- 63 / (63 + 691)  # 0.08355438 * 100 = 8.36%

and for group B it is:

b <- 65 / (65 + 711)  # 0.08376289 * 100 = 8.38%

We can already tell that the difference in percentages is minimal, but this is only an example, so let's continue.

To compare the percentages, we are going to consider them as two categories (A and B) in a binomial distribution. For a binomial test, we need a vector of the same length as the overall number of answers. The overall number of answers in my study is:

all <- 63 + 691 + 65 + 711

To calculate the proportion of the binomial distribution that corresponds to the percentages of orange answers in each group, we simply "scale" (i.e. multiply by the same factor) both percentages to add up to 100% (of all observations); that is, we resolve the calculation:

a * x + b * x = all

The resolution, of course is:

x = all / (a + b)

Now we can calculate the number of observations for each category:

# for A:
a * all / (a + b)

# for B:
b * all / (a + b)

Finally we round the possibly fractional numbers to integers and perform the binomial test:

binom.test(c(round(a * all / (a + b)), round(b * all / (a + b))))

which returns:

number of successes = 764, number of trials = 1530, p-value = 0.9796
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4739876 0.5247077
sample estimates:
probability of success 
             0.4993464 

Answer 1

Your 2nd procedure with the binomial test does not make sense. You are calculating the probability to obtain 764 vs 766 heads and tails for a fair coin. Which relates to some scaling to compare groups A and B as if there was a total of 1530 oranges (instead of 128) and the probability for each group A or B is 0.5

---

You could apply the binomial/multinomial test in place of the Chi-Square test but it involves annoying integration to obtain the distribution of $\chi$, the sum of squares of the error terms.

The Chi-square test is normally resolving this annoying integration by using the Gaussian distribution as an approximation for the distribution of the measured variables. This allows the use of a transform to spherical coordinates and a single integral to obtain the result (since the multivariate Gaussian probability distribution has a elliptical/spherical shape and $\chi$ is the only variable in spherical coordinates for which the probability varies).

If you do not express $\chi$ as approximated by Gaussian distributed error terms then you would have to integrate over all the different dimensions. So say you have your table matrix(c(63, 691, 65, 711), 2, 2) then you could calculate the difference observed minus expected and use the multinomial distribution to calculate the probability for the size of this difference, by integrating over the probabilities for each possible observation with smaller $\chi$, given the expected value.

This calculation only makes sense if the numbers are small and the approximation by a Gaussian distribution is creating large discrepancies.

Example:

#matrix for the beer case 
#(since uncorrected chi^2 is zero for the orange and not interesting to model)
A = matrix(c(101, 374+63+216, 76, 510+65+125), 2, 2,TRUE) 

#A = t(A)

#modeling expected values
pt <- (A[1,1] + A[2,1]) / sum(A)
E <- c(A[1,1]+A[1,2] , A[2,1]+A[2,2]) * pt

r1 <- A[1,1]+A[1,2]
r2 <- A[2,1]+A[2,2]

#calculating chi^2
D <- E - A[,1]
chi2 <- sum(D^2/E)+sum(D^2/(c(r1,r2)-E))

#initialization of integration variable
pchi <- 0

# assuming two binomial distribution with equal p=pt based on observed data
# integrating over the possible observed values for which chi^2 is lower than observed chi^2
#calculate probability for particular set i,j and add to the integration 
for (i in 2:(r1-1)) {
  for (j in 2:(r2-1)) {
    pt2 <- (i+j)/(r1+r2)
    chi22 <- (i-pt2*r1)^2/(pt2*r1) +
      (r1-i-(1-pt2)*r1)^2/((1-pt2)*r1) +
      (j-pt2*r2)^2/(pt2*r2) +
      (r2-j-(1-pt2)*r2)^2/((1-pt2)*r2)
    if(chi22<chi2) { 
      pchi <- pchi + dbinom(i,r1,pt)*dbinom(j,r2,pt)
    }
  }
}

#p-value by multinomial
1 - pchi

#p-value by chi-square distribution
chisq.test(A, correct=0)[3]

These results are close: 0.02742869 and 0.02767172

Notes on this example (see discussion in comments):

• The nuisance parameter problem is solved by using the estimate $\hat{p}$. This is another issue that the chi-square test resolves since the approximation by the Gaussian distribution does not involve this parameter.
• The assumption was made that the row totals are fixed by the experimental design.

---

The above is about your 2nd procedure and the $\chi^2$-test in your 1st procedure.

Regarding your 2nd question (are there other procedures?). You could express estimates for the value and variance/error, or confidence intervals, of the coefficients in the multinomial distribution and use these estimates to express your thoughts about the variations between the two groups for specific individual categories. (see more here: Confidence interval and sample size multinomial probabilities)

I believe that such approach, using confidence intervals, would be better (provide a better idea of the situation) than just calculating the p-values.

(The p-value is more like a combined statistic balancing the difference between groups versus the error in the estimates, and only gives insight about the balance but not about the two separately. Imagine the following thought experiments: 1) you could obtain a high p-value and significant difference for your orange category if you just have a sufficiently large sample size, even with the small 8.36% vs 8.38% difference, and also 2) you could have a non significant difference for your pizza category, even with a factor 2 difference, if your sample size is too small).