# Testing for pairwise proportion differences

I have survey data in which participants had to tell (yes / no) whether a number of reasons could account for their desire to cease smoking.

For instance, fear of...

                     Yes   Total  Proportion
1. Heart attack      114   137    0.8321168
2. Lung cancer       124   137    0.9051095
3. Vision problems   124   137    0.9051095
4. Heart problems    118   137    0.8613139
5. Hearing problems  106   137    0.7737226


The goal is to determine whether the observed differences can be reasonnably explained by chance, or not. My client asked for p-values for pairwise comparisons, and I started reading a little bit on the topic, not being used to testing multiple proportions.

The first problem that I'm trying to address is: what would be the preferred multiple-comparison correction for the p-values? This help page lists a good number of them, with varying leniency: http://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html. I can't find any solid argument for using one over another, though.

The second question I'm asking myself is: would there be a more sound approach to this than doing pairwise comparisons from the onset? I've seen people discuss, among other things, bootstrapping, using an overall test for equality of all proportions and then doing post-hoc tests, and so on.

On a sidenote, I am fully aware of the limitations of p-values. If some approach could combine confidence intervals with these, I'd be more than happy to try it out.

EDIT: After going with the mixed model approach (@rvl's solution), I posted a follow-up question that can be found here

• I use Holm, and I think R functions like pairwise.t.test default to Holm as well. Bonferroni is the classic p-value correction for multiple comparisons, but some say it's too conservative. en.wikipedia.org/wiki/Holm%E2%80%93Bonferroni_method
– Jon
Oct 21, 2016 at 3:27
• Yes, you're right, Holm is the default. And it seems also commonly accepted that Bonferonni is an overkill. Oct 21, 2016 at 4:01
• Is this the same 137 people for all options or is that just a coincidence? Oct 21, 2016 at 12:20
• @mdewey Yes it's the same 137 people. Oct 21, 2016 at 23:24

My suggestion would be to fit a generalized linear mixed model for data with indexes for both subjects and treatments, where a 0 or a 1 (corresponding to "no" and "yes") is recorded for each subject and treatment. Then you may use a package like lsmeans to do the rest. Here's an example.

First, a fake dataset

> junk
Subj treat resp
1     1     A    0
2     2     A    0
3     3     A    0
4     4     A    1
5     1     B    0
6     2     B    0
7     3     B    0
8     4     B    1
9     1     C    1
10    2     C    0
11    3     C    0
12    4     C    1
13    1     D    1
14    2     D    0
15    3     D    1
16    4     D    1
17    1     E    0
18    2     E    0
19    3     E    1
20    4     E    1


Now, fit a model with resp as a Bernoulli response, treat as a fixed factor, and Subj as a random effect:

> library("lme4")
> junk.glmer = glmer(resp ~ treat + (1|Subj), family = binomial, data = junk)


Obtain the estimated proportions from this model using lsmeans:

> library("lsmeans")
> lsm = lsmeans(junk.glmer, "treat", type = "response")
> lsm
treat       prob        SE df    asymp.LCL asymp.UCL
A     0.06360269 0.1866110 NA 0.0001461787 0.9692838
B     0.06360303 0.1866122 NA 0.0001461783 0.9692842
C     0.55119510 0.6632401 NA 0.0063729690 0.9957657
D     0.93950571 0.1677015 NA 0.0456321983 0.9998018
E     0.55119668 0.6632396 NA 0.0063730116 0.9957657

Confidence level used: 0.95
Intervals are back-transformed from the logit scale


... and the pairwise differences of these results:

> pairs(lsm)
contrast   odds.ratio          SE df z.ratio p.value
A - B     0.999994356  2.91218818 NA   0.000  1.0000
A - C     0.055305407  0.16164708 NA  -0.990  0.8598
A - D     0.004373512  0.01540839 NA  -1.542  0.5351
A - E     0.055305054  0.16164611 NA  -0.990  0.8598
B - C     0.055305719  0.16164803 NA  -0.990  0.8598
B - D     0.004373536  0.01540848 NA  -1.542  0.5351
B - E     0.055305367  0.16164703 NA  -0.990  0.8598
C - D     0.079079279  0.20332474 NA  -0.987  0.8614
C - E     0.999993630  2.30045748 NA   0.000  1.0000
D - E    12.645457149 32.51336617 NA   0.987  0.8614

P value adjustment: tukey method for comparing a family of 5 estimates
Tests are performed on the log odds ratio scale


By default, the comparisons are done via odds ratios. You can get differences of proportions by re-gridding the LS means to the response scale:

> pairs(regrid(lsm))
contrast      estimate        SE df z.ratio p.value
A - B    -3.361231e-07 0.1734437 NA   0.000  1.0000
A - C    -4.875924e-01 0.5916295 NA  -0.824  0.9233
A - D    -8.759030e-01 0.2055713 NA  -4.261  0.0002
A - E    -4.875940e-01 0.5916291 NA  -0.824  0.9233
B - C    -4.875921e-01 0.5916291 NA  -0.824  0.9233
B - D    -8.759027e-01 0.2055721 NA  -4.261  0.0002
B - E    -4.875936e-01 0.5916287 NA  -0.824  0.9233
C - D    -3.883106e-01 0.5809326 NA  -0.668  0.9632
C - E    -1.575777e-06 0.5690885 NA   0.000  1.0000
D - E     3.883090e-01 0.5809321 NA   0.668  0.9632

P value adjustment: tukey method for comparing a family of 5 estimates

• Thanks a lot, this seems to be working fine... Not being familiar with those packages, can I ask how you would add in a covariate (gender for instance) in such a way that the differences can be assessed for each of its values? I tried a few things but no luck so far... Oct 25, 2016 at 16:17
• Well, you'd include gender in your model, then do something like lsmeans(my.model, "treat", by = "gender"), or more compactly, lsmeans(my.model, ~ treat | gender). Enter vignette("using-lsmeans") for a tutorial and pretty extensive collection of examples with lsmeans. There is also a lot of help and examples available for the lme4 package. Oct 25, 2016 at 20:04
• Ok, so I think I got this figured out. I am somewhat puzzled though with some of the results returned by lsmeans. Whereas with a glm basic model, the proportions (prob) range between .73 and .87, with glmer the five probabilities are .99 with very low SE. New question altogether maybe? Oct 26, 2016 at 12:47
• Different question, but it seems like something is wrong. GLMMs can be numerically temperamental so check to make sure there are no warnings and that the estimated subject variance is reasonable. Oct 26, 2016 at 17:59
• For reference, the follow-up question can be found here: stats.stackexchange.com/questions/242559/… Oct 27, 2016 at 19:06

There is an inferential issue here because of the related proportions. If you set up the data as a $137 \times 5$ matrix of ones and zeroes then the natural way to compare columns is the McNemar test. Since this only uses the discordant pairs and ignores those which are 1 in both columns or 0 in both columns it is quite likely that the comparison of any given pair of columns will be based on a different subset of the data. Whether this matters is a matter for discussion with the client. An overall model using GEE or a mixed effects logistic regression could be followed by tests between the coefficients for pairs of columns although the subset issue is still there.

• Thanks for your answer. I tried the mcnemar.test` in R but it seems to apply only to square matrices (2 x 2 contingency tables). I'll try the regression route. Oct 25, 2016 at 15:25
• You do have a series of 2 by 2 matrices if you tabulate each column by another. Oct 25, 2016 at 15:36
• True... Then I'd need to apply "manually" any p-value correction I gather. Oct 25, 2016 at 15:43
• I expect there is a programming solution but goinf with the regression suggestion seems easier. Oct 25, 2016 at 15:46
• I agree... I'm currently experimenting with the solution proposed by @rvl, and it looks promising. Not sure how GEE would differ from this method though. Oct 25, 2016 at 16:22