Statistical comparison of counts/proportions between multiple groups in R?

Question

I am struggling with choosing the appropriate statistical analysis for a dataset in which there are multiple groups (groupA-groupE), each having a certain number of counts in two categories (healthy or sick) that can also be presented as a proportion (sick/[sick+healthy]):

df <- data.frame(group = c("groupA","groupB","groupC","groupD","groupE"),
n_sick = c(12, 32, 99, 37, 48),
n_healthy = c(36, 250, 120, 68, 93))
df %<>% mutate(tot = n_sick + n_healthy, prop = n_sick / tot)

I would like to find out which groups have statistically different proportions from each other, akin to what one would do in e.g. ANOVA/Kruskal-Wallis+post-hoc tests when there are replicates. However, the only available data in this case is the number of sick and healthy individuals per group. I understand the I can use the Chi-square test to uncover if there is a relationship between proportion and group, but to my knowledge there are no standardized post-hoc tests for testing intergroup differences (correct me if I'm wrong!).

I found that the stats package contains a function called pairwise_prop_test that looks as though it does what I want it to do:

t <- as.table(rbind(
  c(12, 32, 99, 37, 48),
  c(36, 250, 120, 68, 93)))
dimnames(t) <- list(
  condition = c("n_sick", "n_healthy"),
  group = c("groupA","groupB","groupC","groupD","groupE"))

pairwise_prop_test(t, p.adjust.method = "bonferroni")

... which spits out a p.adj value for each comparison between groups (e.g. groupB differs from groupC). Is this the correct approach to take, or am I overlooking something? Thank you for taking the time to help!

Answer 1

It would seem like an obvious idea to use logistic regression, which you can specify in many ways:

df <- data.frame(group = c("groupA","groupB","groupC","groupD","groupE"),
                 n_sick = c(12, 32, 99, 37, 48),
                 n_healthy = c(36, 250, 120, 68, 93)) %>%
  mutate(tot = n_sick + n_healthy, 
         prop = n_sick / tot)

# Option 1 of specifying a logistic regression
glm(cbind(n_sick, n_healthy) ~ 0 + group,
    data=df,
    family=binomial(link="logit"))

# Option 2 of specifying a logistic regression
glm(prop ~ 0 + group, 
    weights=tot,
    data=df,
    family=binomial(link="logit"))

# Option 3 of specifying a logistic regression (creating individual records)
glm(y ~ 0 + group, 
    data=df %>% 
      mutate(dat=map2(n_sick, n_healthy, \(x, y) tibble(y=c(rep(1, x), rep(0, y))))) %>%
      unnest(dat),
    family=binomial(link="logit"))

Above, I showed three ways of specifying the response with different input data formats. You can also parameterize the model differently, here I directly estimate log-odds for each group, but if you write ~ 1 + group instead, one group would end up as the reference group and the coefficients for the other groups would be relative to the reference group.

Thereafter, you can form comparisons. E.g. the below gives you all pairwise comparisons without any multiplicity adjustment (which might well be needed here):

fit1 <- glm(cbind(n_sick, n_healthy) ~ 0 + group,
            data=df,
            family=binomial(link="logit"))

emmeans::emmeans(fit1, ~group) %>%
  emmeans::contrast(method="pairwise", adjust="none")

This of course all assumes that comparing the groups like this is appropriate. That would e.g. make sense, if this is a randomized controlled trial with 5 treatment groups, where each treatment might or might not prevent you from getting sick. It's less clear that this is appropriate, if this is an observational comparison, in which case this kind of simple comparison is of course difficult for saying anything other than the groups are somehow different. E.g. if you compare people that happen to be prescribed drugs, they might already differ in terms of risk factors for getting sick like age, medical background etc. before they get prescribed the drugs in ways that influence the outcome - thus interpreting these comparisons as estimating causal effects of these drugs (or whatever else defines these groups) would be wrong. In such a scenario, more sophisticated methods for observational studies (e.g. propensity score or some of the many other options) would be needed.