1
$\begingroup$

There are many questions on this wonderful site about post-hoc contrasts within chi-squares, however none of the posts I have seen have quite answered the questions I have about the procedure. I will use data from this CV post, concerning disparities in gender publication rates at five different journals, to illustrate. I will use R code (relying heavily on the tidyverse packages) but this is only because I know R. Anyone, R-users or non R-users, who would can help please weigh in.

# toy data
m <- data.frame(male = c(8278, 2183, 3844, 590, 2341),
                female = c(1659, 422, 899, 137, 662)) 
row.names(m) <- c("PRx", "AXJ", "JSPRAS", "CPX", "PRX-TO")

# get row totals
m %>%
  mutate(tot = male + female,
         percF = round(female/tot*100,2)) %>%
    add_column(journal = c("PRx", "AXJ", "JSPRAS", "CPX", "PRX-TO")) %>%
      column_to_rownames(var = "journal") -> mLarge

mLarge

#        male female  tot percF
# PRx    8278   1659 9937 16.70
# AXJ    2183    422 2605 16.20
# JSPRAS 3844    899 4743 18.95
# CPX     590    137  727 18.84
# PRX-TO 2341    662 3003 22.04

Males are clearly getting more papers published than females across all journals. The reason is unknown as we only know raw numbers, but let's leave that aside. The question I would like to be able to answer with these data is: are there significant differences between journals in the proportion of articles written by males vs females? It looks like the proportion of articles written by females is slightly higher (22%) at the PRX-TO journal than the AXJ journal (16%) so it would be nice to test whether this difference in rates of female-authored manuscripts between these two journals is significant.

In the original post I linked to, the OP did not formulate their question this way and perhaps this is why I did not get the answer I wanted from that post.

The overall chi square is highly significant

# run chi-squares 
(chisq.test(m) -> chisq_m)

# Pearson's Chi-squared test
# 
# data:  m
# X-squared = 53.78, df = 4, p-value = 5.851e-11 

But, as is often the case, this does not tell us very much. To delve a bit deeper we can use the chisq.posthoc.test() function from the chisq-posthoc.test package in R to obtain p-values for post-hoc comparisons

# chi-square post hoc test
library(chisq.posthoc.test)
chisq.posthoc.test(m)

#    Dimension     Value       male     female
# 1        PRx Residuals  4.6018779 -4.6018779
# 2        PRx  p values  0.0000420  0.0000420
# 3        AXJ Residuals  2.5314042 -2.5314042
# 4        AXJ  p values  0.1136070  0.1136070
# 5     JSPRAS Residuals -1.9805914  1.9805914
# 6     JSPRAS  p values  0.4763710  0.4763710
# 7        CPX Residuals -0.6160692  0.6160692
# 8        CPX  p values  1.0000000  1.0000000
# 9     PRX-TO Residuals -6.2610684  6.2610684
# 10    PRX-TO  p values  0.0000000  0.0000000

Now from what I have read and comprehended thus far, this output is the result of an analysis of standardised residuals from the chi-square test, specifically the standardised difference between the expected value and the observed values. The p-values in each cell reflect (i) how different the observed values in each cell are from the expected values, and (ii) how much the cell contributes to the overall chi-square statistic.

What I cannot see in the output of the post-hoc test, is how to answer my question, namely is the differences in proportion of male vs female publications in the AXJ journal significantly different from the difference in proportion in the PRX-TO journal?

So I would really appreciate knowing (i) if it is possible to test the comparison I am interested in with a chi-square post-hoc test, (ii) if it is possible, how to perform that comparison.

p.s. To answer this question I would ordinarily conduct a binomial regression with 'male vs female' as the binary count outcome and journal as the predictor. However a colleague mentioned that you can also answer such questions of such outcomes with chi-sqaure posthoc tests, so I came onto CV for answers, and found myself confused instead (which is fine, it's a process).

Even if the answer just directs me to a decent tutorial that will tell me what I need to know, that would help.

$\endgroup$

3 Answers 3

2
+100
$\begingroup$

One ad-hoc approach would be to conduct chi-square tests pairwise across the Journals. That is, compare PRx to AXJ, and PRx to JSPRAS, and so on, each, in this case, with a 2 x 2 contingency table.

What I don't like about this approach is that it ignores all the data that isn't in the pair being considered. Perhaps someone knows how to do something similar but use the expected counts from the whole omnibus table.

Since there are only two categories in Gender in this case, another ad-hoc approach might be to construct confidence intervals about the estimate for the proportion of Male for each Journal. This suffers from the same drawback that each confidence interval is based only on the observations for one Journal, and ignores the rest of the data. One could consider those proportions without overlapping confidence intervals to be statistically distinct, though typically this doesn't line up precisely with a hypothesis test.

The following is some R code for these approaches. Note that I changed the name of one Journal, and changed the way the data are arranged.

Also, I'm using the pairwiseNominalIndependence() function from the rcompanion package, to keep things simple for the pairwise chi-square test approach. Caveat that that I wrote this function.

The BinomCI() function from the DescTools package has options for confidence intervals for a binomial proportion, and makes it easy to compute these for each row of the data frame.

Addition:

Because there are only two categories for Gender, this specific situation lends itself to using binomial logistic regression. Code and results added below. The post-hoc result are similar to the other methods.

Data = data.frame(Count    = c(8278, 2183, 3844, 590, 2341, 1659,  422,  899, 137,  662),
               Gender  = c(rep("male",5), rep("female", 5)),
               Journal = c("PRx", "AXJ", "JSPRAS", "CPX", "PRX_TO"))

Table = xtabs(Count ~ Journal + Gender, data=Data)

Table

   ###         Gender
   ### Journal  female male
   ###  AXJ       422 2183
   ###  CPX       137  590
   ###  JSPRAS    899 3844
   ###  PRx      1659 8278
   ###  PRX_TO    662 2341

library(rcompanion)

PNI = pairwiseNominalIndependence(Table, fisher=FALSE, gtest=FALSE)

PNI

   ###         Comparison  p.Chisq p.adj.Chisq
   ### 1        AXJ : CPX 1.03e-01    1.47e-01
   ### 2     AXJ : JSPRAS 3.62e-03    7.24e-03
   ### 3        AXJ : PRx 5.65e-01    6.28e-01
   ### 4     AXJ : PRX_TO 3.92e-08    1.96e-07
   ### 5     CPX : JSPRAS 9.84e-01    9.84e-01
   ### 6        CPX : PRx 1.49e-01    1.86e-01
   ### 7     CPX : PRX_TO 6.63e-02    1.10e-01
   ### 8     JSPRAS : PRx 8.05e-04    2.65e-03
   ### 9  JSPRAS : PRX_TO 1.06e-03    2.65e-03
   ### 10    PRx : PRX_TO 2.58e-11    2.58e-10

cldList(p.adj.Chisq ~ Comparison, data=PNI)

   ###    Group Letter MonoLetter
   ### 1    AXJ      a        a  
   ### 2    CPX    abc        abc
   ### 3 JSPRAS      b         b 
   ### 4    PRx      a        a  
   ### 5 PRX_TO      c          c

library(DescTools)

m = data.frame(journal = c("PRx", "AXJ", "JSPRAS", "CPX", "PRX-TO"),
               male    = c(8278, 2183, 3844, 590, 2341),
               female  = c(1659, 422, 899, 137, 662)) 

m$total = m$male + m$female

m$proportionMale = BinomCI(m$male, m$total, method="clopper-pearson")[,1]

m$lower.ci = BinomCI(m$male, m$total, method="clopper-pearson")[,2]

m$upper.ci = BinomCI(m$male, m$total, method="clopper-pearson")[,3]

m

   ###   journal male female total proportionMale  lower.ci  upper.ci
   ### 1     PRx 8278   1659  9937      0.8330482 0.8255681 0.8403326
   ### 2     AXJ 2183    422  2605      0.8380038 0.8232849 0.8519587
   ### 3  JSPRAS 3844    899  4743      0.8104575 0.7990081 0.8215235
   ### 4     CPX  590    137   727      0.8115543 0.7811964 0.8393582
   ### 5  PRX-TO 2341    662  3003      0.7795538 0.7642920 0.7942690

library(ggplot2)

qplot(x    = journal,
      y    = proportionMale,
      data = m) +
  
  geom_errorbar(aes(
    ymin  = lower.ci,
    ymax  = upper.ci,
    width = 0.15))

model = glm(factor(Gender) ~ Journal, family=binomial(), data=Data, weights=Data$Count)

library(car)

Anova(model)

library(emmeans)

marginal = emmeans(model, ~ Journal, type="response")

marginal

    ### Journal  response      SE  df asymp.LCL asymp.UCL
    ###  AXJ        0.838 0.00722 Inf     0.823     0.852
    ###  CPX        0.812 0.01450 Inf     0.781     0.838
    ###  JSPRAS     0.810 0.00569 Inf     0.799     0.821
    ###  PRx        0.833 0.00374 Inf     0.826     0.840
    ###  PRX_TO     0.780 0.00756 Inf     0.764     0.794

marginal2 = emmeans(model, ~ Journal)

pairs(marginal2)

library(multcomp)

cld(marginal2, Letters=letters)

   ### Journal emmean     SE  df asymp.LCL asymp.UCL .group
   ### PRX_TO    1.26 0.0440 Inf      1.18      1.35  a    
   ### JSPRAS    1.45 0.0370 Inf      1.38      1.53   b   
   ### CPX       1.46 0.0948 Inf      1.27      1.65  abc  
   ### PRx       1.61 0.0269 Inf      1.55      1.66    c  
   ### AXJ       1.64 0.0532 Inf      1.54      1.75    c  
   ###    
   ### Results are given on the logit (not the response) scale. 
   ### Confidence level used: 0.95 
   ### Results are given on the log odds ratio (not the response) scale. 
   ### P value adjustment: tukey method for comparing a family of 5 estimates 
   ### significance level used: alpha = 0.05 
   ### NOTE: If two or more means share the same grouping letter,
   ### then we cannot show them to be different.
   ### But we also did not show them to be the same. 

enter image description here

$\endgroup$
5
  • $\begingroup$ Thank you @Sal Mangiafico, t I definitely got my points' worth with your (and the others') answer. Thank you for your extremely comprehensive answer. Thank you so much for alerting me to the existence of the pairwiseNominalIndependence() function. Would I be right in saying that of the three methods you outlined, this is the closest to what I asked for: i.e. a way of testing interaction contrasts (in this case between-journal differences in proportion of females) within a chi-squared/contingency table analysis? So this technique amounts to a series of 2x2 chisquares, is that right? $\endgroup$
    – llewmills
    Sep 2, 2022 at 2:12
  • $\begingroup$ I really don't know what "interaction contrasts" would mean in this context. I would probably call any of these approaches just "post-hoc analysis." ... Of the three approaches I mentioned, the pairwise chi-square tests is maybe the most straightforward, easy to explain, and general. And it conceptually matches the omnibus chi-square analysis. In this example, yes, each comparison is a 2 x 2 chi-square test. Then, if you had three genders in this example, you could follow the same procedure, determine differences between the journals, with each comparison being a 2 x 3 table. $\endgroup$ Sep 2, 2022 at 14:11
  • $\begingroup$ If you do have a dichotomous outcome, the binomial regression plus emmeans post-hoc would be my preference. It has the advantage of using all the data in the model and post-hoc. And emmeans gives you confidence intervals for each Journal, which will make a nice plot. ... If you had more than 2 genders, you could use a multinomial logistic regression. This makes for a little more complicated analysis, and I don't know if or how emmeans would handle multinomial models. $\endgroup$ Sep 2, 2022 at 14:16
  • $\begingroup$ One additional comment: there are also log-linear models, as in the loglm in the MASS package. The only post-hoc I know for these would be employing a pairwise approach of these models. $\endgroup$ Sep 2, 2022 at 14:18
  • $\begingroup$ emmeans handles multinomial regression very well @ Sal Mangiafico, through the multinom() command in the nnet package. And yes 'interaction contrasts' is probably an odd expression in this context. but it is one in a way: difference in male vs female publications and then difference in that difference across journals. $\endgroup$
    – llewmills
    Sep 3, 2022 at 16:54
2
$\begingroup$

For your specific hypothesis, all you need to do is enter the subset of the data you care about into the chi-squared test.

subset_m = m[c('AXJ', 'PRX-TO'),]
subset_m
#        male female
# AXJ    2183    422
# PRX-TO 2341    662

chisq.test(subset_m)

# Pearson's Chi-squared test with Yates' continuity correction
# 
# data:  subset_m
# X-squared = 30.189, df = 1, p-value = 3.92e-08

If you hadn't planned on running this specific comparison beforehand (it is really a post-hoc test) you should also adjust your p-value for multiple comparisons. The simplest way of doing this is using Bonferroni correction. Rather than considering a result as significant if $p < .05$, only do so if $p < \frac{.05}{m}$, where $m$ is the number of post-hoc tests you could have run, which for 5 journals is 5 * (5 -1) / 2 = 10 pairs, so a criterion of $p < .005$ here.

Finally, the chisq.posthoc.test() function seems to test, individually for each cell in the data, whether the value is significantly different from the expected value. This isn't really relevant to your hypothesis.


Bonus

If you wanted to do this for every pair of journals...

pairs_of_journals = combn(rownames(m), 2, simplify = F)
map_df(pairs_of_journals, function(pair){
  pair_data = m[pair,]
  pair_label = paste(pair, collapse = ' vs. ')
  broom::tidy(chisq.test(pair_data)) %>%
    mutate(contrast = pair_label)
}) %>%
  select(contrast, everything())
# A tibble: 10 × 5
#   contrast          statistic  p.value parameter method    
#   <chr>                 <dbl>    <dbl>     <int> <chr>     
# 1 PRx vs. AXJ        0.331    5.65e- 1         1 Pearson's…
# 2 PRx vs. JSPRAS    11.2      8.05e- 4         1 Pearson's…
# 3 PRx vs. CPX        2.08     1.49e- 1         1 Pearson's…
# 4 PRx vs. PRX-TO    44.5      2.58e-11         1 Pearson's…
# 5 AXJ vs. JSPRAS     8.47     3.62e- 3         1 Pearson's…
# 6 AXJ vs. CPX        2.66     1.03e- 1         1 Pearson's…
# 7 AXJ vs. PRX-TO    30.2      3.92e- 8         1 Pearson's…
# 8 JSPRAS vs. CPX     0.000379 9.84e- 1         1 Pearson's…
# 9 JSPRAS vs. PRX-TO 10.7      1.06e- 3         1 Pearson's…
#10 CPX vs. PRX-TO     3.37     6.63e- 2         1 Pearson's…
$\endgroup$
1
  • $\begingroup$ Thank you @Eoin. I appreciate your thoughtful response. $\endgroup$
    – llewmills
    Sep 2, 2022 at 2:13
2
$\begingroup$

One approach is to recognize the close relationship between this type of contingency table and binomial regression. I added a journal column to the m data frame

m
#   male female journal
# 1 8278   1659     PRx
# 2 2183    422     AXJ
# 3 3844    899  JSPRAS
# 4  590    137     CPX
# 5 2341    662  PRX-TO

and performed logistic binomial regression with journal as predictor.

mod1 <- glm(cbind(female,male)~journal,data=m,family=binomial)

That works in the log-odds (logit) scale.

There are many tools for doing comparisons after building a model. I used the emmeans package to evaluate all pairwise differences among the journals, with the default Tukey correction for multiple comparisons, and return values to the probability scale.

library(emmeans)
emmeans(mod1,pairwise~journal,type="response")
# $emmeans
#  journal  prob      SE  df asymp.LCL asymp.UCL
#  AXJ     0.162 0.00722 Inf     0.148     0.177
#  CPX     0.188 0.01450 Inf     0.162     0.219
#  JSPRAS  0.190 0.00569 Inf     0.179     0.201
#  PRx     0.167 0.00374 Inf     0.160     0.174
#  PRX-TO  0.220 0.00756 Inf     0.206     0.236
# 
# Confidence level used: 0.95 
# Intervals are back-transformed from the logit scale 
# 
# $contrasts
#  contrast          odds.ratio     SE  df null z.ratio p.value
#  AXJ / CPX              0.833 0.0905 Inf    1  -1.686  0.4427
#  AXJ / JSPRAS           0.827 0.0536 Inf    1  -2.939  0.0273
#  AXJ / PRx              0.965 0.0575 Inf    1  -0.605  0.9743
#  AXJ / (PRX-TO)         0.684 0.0472 Inf    1  -5.510  <.0001
#  CPX / JSPRAS           0.993 0.1011 Inf    1  -0.070  1.0000
#  CPX / PRx              1.159 0.1142 Inf    1   1.494  0.5666
#  CPX / (PRX-TO)         0.821 0.0859 Inf    1  -1.885  0.3254
#  JSPRAS / PRx           1.167 0.0534 Inf    1   3.372  0.0067
#  JSPRAS / (PRX-TO)      0.827 0.0476 Inf    1  -3.301  0.0086
#  PRx / (PRX-TO)         0.709 0.0366 Inf    1  -6.674  <.0001
# 
# P value adjustment: tukey method for comparing a family of 5 estimates 
# Tests are performed on the log odds ratio scale 

That displays odds ratios and p-values for all pairwise comparisons. For more complicated models you might be better off to pre-define the reference grid of predictor values for your post-modeling analysis via emmeans, but this model can be handled simply as shown.

$\endgroup$
1
  • $\begingroup$ Thanks @EdM. A binomial model in glm() and then emmeans() is how I would usually answer my question. I was just curious whether and equivalent analysis could be done inside a chi-square post-hoc test. It looks to me like the answer is no, at least not in the same way you would do it with a binomial regression. $\endgroup$
    – llewmills
    Sep 2, 2022 at 1:13

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.