# Not sure what test to run to compare the proportions of a count outcome across a variable with 5 levels, between two groups (using R for analyses)

For a case-control study that I've conducted I've measured how many times a count outcome (count_outcome) occurred in each of 4 areas (area), between healthy controls and patients (group).

Some example data looks as follows:

library(tidyverse)
library(magrittr)

mydata <- structure(list(pat_id = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3,
4, 4, 4, 4, 5, 5, 5, 5, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30,
42, 42, 42, 43, 43, 43, 43, 44, 44, 44, 51, 51, 51, 51, 52, 52,
52, 52, 53, 53, 53, 53, 54, 54, 54, 54, 55, 55, 55, 55, 56, 56,
56, 56, 57, 57, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60,
60, 60, 61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 64, 64,
64, 64, 65, 65, 65, 65, 66, 66, 66, 66, 67, 67, 67, 67, 68, 68,
68, 68, 69, 69, 69, 69, 70, 70, 70, 70, 71, 71, 71, 71, 72, 72,
72, 72, 73, 73, 73, 73, 74, 74, 74, 74, 75, 75, 75, 75, 76, 76,
76, 76, 77, 77, 77, 77, 78, 78, 78, 78, 79, 79, 79, 79, 80, 80,
80, 80, 81, 81, 81, 81, 82), group = c("Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Patient", "Patient", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Healthy control", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Patient", "Patient", "Patient",
"Patient", "Patient", "Patient", "Healthy control", "Healthy control",
"Healthy control", "Healthy control", "Patient", "Patient", "Patient",
"Patient", "Healthy control", "Healthy control", "Healthy control",
"Healthy control", "Patient"), area = c("Area 1", "Area 2", "Area 3",
"Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2",
"Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 2", "Area 3", "Area 4", "Area 3", "Area 4", "Area 1", "Area 2",
"Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 2", "Area 3", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4",
"Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3",
"Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2",
"Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4",
"Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3",
"Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2",
"Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4",
"Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3",
"Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2",
"Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4",
"Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3",
"Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1", "Area 2",
"Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4", "Area 1",
"Area 2", "Area 3", "Area 4", "Area 1", "Area 2", "Area 3", "Area 4",
"Area 2"), count_outcome = c(1, 0, 0, 0, 3, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2,
1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 0, 1, 0,
1, 2, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 2, 0, 0, 2, 0, 0, 0,
1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,
1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0)), row.names = c(NA, -165L), class = "data.frame")

pat_id           group   area count_outcome
1      1 Healthy control Area 1             1
2      1 Healthy control Area 2             0
3      1 Healthy control Area 3             0
4      1 Healthy control Area 4             0
5      2 Healthy control Area 1             3
6      2 Healthy control Area 2             0


The total amount of the count outcome per group is:

mydata %>%
group_by(group) %>%
summarise(sum_outcome=sum(count_outcome))

> A tibble: 2 x 2
group           sum_outcome
<chr>                 <dbl>
1 Healthy control          44
2 Patient                  37


And the total amount of the outcome per group per area is:

mydata %>%
group_by(group, area) %>%
summarise(sum_outcome=sum(count_outcome))

# A tibble: 8 x 3
# Groups:   group [2]
group           area   sum_outcome
<chr>           <chr>        <dbl>
1 Healthy control Area 1          18
2 Healthy control Area 2           4
3 Healthy control Area 3           9
4 Healthy control Area 4          13
5 Patient         Area 1          21
6 Patient         Area 2           4
7 Patient         Area 3           6
8 Patient         Area 4           6


Together these make the following contingency table:

- Patients Healthy controls
Area 1 21 (56.8%) 18 (40.9%)
Area 2 4 (10.8%) 4 (9.1%)
Area 3 6 (16.2%) 9 (20.5%)
Area 4 6 (16.2%) 13 (29.5%)
Total 37 (100%) 44 (100%)

I think it is important to note here that some participants, as they had an outcome in multiple areas, contribute multiple observations to this contingency table.

What I'm trying to do is figure out if, overall, the distribution of my count outcome across the areas is different between the groups (patients and healthy controls). In other words, I want assess if, overall, the proportion of outcomes for the areas is different between the groups. I'd like to stress that I don't want to assess the difference in the outcome between groups per area.

I've looked at a Chi square, but that does not seem appropriate since some participants contribute multiple observations. Then I looked at the Friedman test, but thats only for 2 variables? Anyway, I can't seem to make it work. Any help is appreciated!

• A major difficulty here is "[S]ome participants, as they had an outcome in multiple areas, contribute multiple observations to this contingency table ,,," A valid contingency table for chi-squared analysis must have disjoint outcomes. Consequently, the grand total of the table (sum or row sums or sum of column sums) must be equal to the number of subjects. // In case it is helpful, I will show a valid contingency table with a couple of tests in Answer format. Dec 15, 2021 at 3:26

Valid contingency table. In case you can resolve the issue of counting some individual subjects in more than one outcome category, here is a fictitious contingency table with disjoint outcomes for 74 subjects altogether (using R). We will use it for two tests.

H = c(15, 4, 9, 12)
P = c(17, 4, 6,  7)
TBL = rbind(H,P)
rowSums(TBL); colSums(TBL)
TBL; rowSums(TBL); colSums(TBL); sum(TBL)
[,1] [,2] [,3] [,4]
H   15    4    9   12
P   17    4    6    7
H  P
40 34
[1] 32  8 15 19
[1] 74


So the table, including totals is:

 Type\Outcome  1   2   3   4    Total
------------------------------------
H             15  4   9  12      40
P             17  4   6   7      34
------------------------------------
Total         32  8  15  19      74


Are the four outcomes equally likely? You might ask whether (regardless of subject type) the four outcomes might be equally likely. In R, that null hypothesis can be tested using the procedure prop.test, which is equivalent to a "one dimensional" chi-squared test on column totals.

prop.test (c(32, 8, 15, 19), rep(74, 4))

4-sample test for equality of proportions
without continuity correction

data:  c(32, 8, 15, 19) out of rep(74, 4)
X-squared = 21.982, df = 3, p-value = 6.58e-05
alternative hypothesis: two.sided
sample estimates:
prop 1    prop 2    prop 3    prop 4
0.4324324 0.1081081 0.2027027 0.2567568


There are highly significant differences among the four (approximate) proportions $$(0.43, 0.11, 0.20, 0.26)$$ at the 5% level, as indicated by the P-value near 0, far below $$0.05 = 5\%.$$

Do Healty and Patient subjects differ as to outcome distribution? Another possible test has the null hypothesis that there is no difference between how Healthy and Patient subjects fall among the four categories of outcomes.

In R, Pearson's chi-squared test of this null hypothesis compares observed counts with expected counts (computed on the basis of the null hypothesis of homogeneity).

chisq.test(TBL)

Pearson's Chi-squared test

data:  TBL
X-squared = 1.5646, df = 3, p-value = 0.6674

Warning message:
In chisq.test(TBL) :
Chi-squared approximation may be incorrect


The warning message notes that not all expected counts are 5 or more, possibly leading to an incorrect P-value. The table of expected counts is as follows:

chisq.test(TBL)\$exp
[,1]     [,2]     [,3]     [,4]
H 17.2973 4.324324 8.108108 10.27027
P 14.7027 3.675676 6.891892  8.72973


Although some statisticians would be happy enough to have all expected counts above 3, the implementation of chisq.test in R provides the possibility of simulating a more accurate P-value:

chisq.test(TBL, sim=T)

Pearson's Chi-squared test
with simulated p-value
(based on 2000 replicates)

data:  TBL
X-squared = 1.5646, df = NA, p-value = 0.7166


The P-value $$0.72$$ is far above 5%, so there is no evidence that the distributions of healthy and patient subjects into the four outcomes differ.

• Hi @BruceET, thanks for the elaborate answer. As the fact is that participants are actually likely to have an outcome in multiple areas (vs. just one in one area) if an outcome is measured, I'm afraid a Chi-Square test wouldn't be appropriate. I thought about comparing two mixed models with a LRT, with one just containing area, and the other area*group (with count_outcome as outcome and 1|pat_id as random effect), to see whether there is an interaction between the two, but that seems overly complicated. You wouldn't know of a more straightforward way to deal with this issue? Dec 15, 2021 at 6:33