I have asked people which food they prefer:

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:

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:

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 
  • 2
    $\begingroup$ (Following your Meta post.) I think both methods that you suggest are fine and should be basically equivalent. You do get a small difference in p-values (that I cannot explain) but it looks negligible. $\endgroup$
    – amoeba
    Aug 16, 2017 at 19:45
  • $\begingroup$ Unlike @amoeba, I can't follow the logic of your second procedure at all. It's certainly not equivalent to the first. $\endgroup$ Aug 16, 2017 at 20:45
  • 1
    $\begingroup$ You may want to look at residuals in each cell and their significance (for data meeting Chi-sq test assumptions adj. st. residual >2 in magnitude means p<.05. stats.stackexchange.com/q/74708/3277 $\endgroup$
    – ttnphns
    Aug 22, 2017 at 16:34
  • 1
    $\begingroup$ It is said in my answer there: because adj. residuals are normally distributed as standard z-scores, so 1.96 (approx. 2) is what cuts off 2.5% tail of the curve, and 2*2.5=5% (i.e. alpha .05) with 2-sided null hypothesis. $\endgroup$
    – ttnphns
    Aug 22, 2017 at 19:57
  • 1
    $\begingroup$ Yes, considering Bonferroni correction is a way to go. If you didn't plan a specific cell or two beforehand, i.e. if you are doing post hoc. $\endgroup$
    – ttnphns
    Aug 22, 2017 at 20:46

1 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.


#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) +
    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).

  • 1
    $\begingroup$ The code for the small-sample test you describe raises some questions: (1) E1 - typo for E? (2) Shouldn't probabilities (under the null) of tables where the test statistic is greater than or equal to that observed be summed to give the p-value? (3) The sampling scheme you've used has fixed row margins (a product-binomial rather than a multinomial distribution). Why not fix the column margins instead, or all margins, or just the overall total? (4) Why do "possible" values for i & j range from 60 to 120 rather than from 0 to 754 and 0 to 776? ... $\endgroup$ Aug 29, 2017 at 9:47
  • 1
    $\begingroup$ ... (5) The null hypothesis under which probabilities of each table is calculated is a simple one derived from the observed table rather than the composite one of homogeneity specified up front for Pearson's chi-square test. What's the justification for that? (Exact tests for 2-by-2 tables typically maximize the p-value over the nuisance parameter or remove the dependence of the test statistic on it by conditioning on all margins.) $\endgroup$ Aug 29, 2017 at 9:49
  • 1
    $\begingroup$ I've corrected (1), (2) and (4). The reduction of margins (4) came from a test-code which calculated the binomial for all values (taking lot's of time). (3) and (5) are though to discuss (and I do not get 5 entirely). What I did is trying to mimic the $\chi^2$-test but using the binomial distribution in place of the, $N(\mu=E,\sigma=\sqrt{E})$ distribution. I fixed the row totals assuming this was the experiment design (one time 754 people from group A are asked, one time 776 people from group B are asked). Indeed, this could also have been a test in which 884 apple eaters, etc, were gathered. $\endgroup$ Aug 30, 2017 at 12:00
  • $\begingroup$ .... I am not sure if it matters how you fix the columns totals or row totals. At least in the representation of the normal distribution it all boils down to four parameters/dimensions that have three linear restrictions (being restricted to a total, or being restricted due to a modeled expected value) and it becomes the same. But maybe if you use the binomial distribution it is different. (edit: tested it you get a different result for transpose of A) $\endgroup$ Aug 30, 2017 at 12:04
  • 1
    $\begingroup$ Wouldn't call it self-study! Though both approaches result in valid tests, which is better, or what "better" ought to mean, is a rather tricky question that's been argued about for decades. $\endgroup$ Aug 31, 2017 at 8:00

Your Answer

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