Analyzing nominal responses using the
multinomial-Poisson trick
Jacob O. Wobbrock
2025-02-09
Introduction
This vignette shows how to use the multpois package for
analyzing nominal response data. Nominal responses, sometimes called
multinomial responses, are unordered categories. In certain experiments
or surveys, the dependent variable can be one of N
categories.
For example, let’s say we ask people what their favorite ice cream
flavor is: vanilla, chocolate, or strawberry. This three-category
response would be a polytomous dependent variable. Also, let’s
say we wish to ask both adults and children about their favorite ice
cream flavors to see if there is a difference by age group. We would
then have a two-level between-subjects factor, Age. If we ask
each respondent only once, this data set would represent a one-way
between-subjects design. But perhaps we ask each participant once
each season—in fall, winter, spring, and summer—to see if their
responses vary by season. Now we would have a four-level within-subjects
factor, Season, resulting in repeated measures, since each
respondent would be asked four times.
The multpois package helps us analyze this type of data,
where the dependent variable is nominal. It does so by modeling nominal
responses as counts of category choices and uses mixed-effects Poisson
regression to analyze these counts (Baker 1994, Chen & Kuo 2001).
This technique is known as the multinomial-Poisson transformation
(Guimaraes 2004) or trick (Lee et al. 2017).
R already provides options for the following situations:
If the response is dichotomous, and the factors are only
between-subjects, we can build a model using glm with
family=binomial from the base stats package.
The Anova function from the car package can be
used to produce main effects and interactions. The emmeans
function from the emmeans package can be used to produce
post hoc pairwise comparisons.
If the response is polytomous, and the factors are only
between-subjects, we can build a model using multinom from
the nnet package. The Anova function from the
car package can be used to produce main effects and
interactions. However, we cannot use the emmeans function
from the emmeans package in the usual fashion. An approach
to this issue by emmeans package author Russ Lenth is
offered below.
If the response is dichotomous, and one or more factors is
within-subjects, we can build a model using glmer with
family=binomial from the lme4 package. The
Anova function from the car package can be
used to produce main effects and interactions. The emmeans
function from the emmeans package can be used to produce
post hoc pairwise comparisons.
If the response is polytomous, and one or more factors is
within-subjects, there is no easy option similar to the three above. The
multinom function in nnet cannot accept random
factors to handle repeated measures, and the glmer function
in lme4 does not offer a family=multinomial
option. This package was created to address this case in particular,
although it can address the above three cases, also.
The first four analyses below illustrate 2×2 designs having between-
and within-subjects factors and dichotomous and polytomous responses.
(The functions in multpois are not limited to 2×2 designs;
any number of between- and within-subjects factors can be used.) The
first three examples first use existing R solutions to which the results
from multpois functions can be compared.
The fifth example returns to our ice cream scenario, above, and
analyzes a mixed factorial design with one between-subjects factor
(Age) and one within-subjects factor
(Season).
References
Baker, S.G. (1994). The multinomial-Poisson transformation. The
Statistician 43 (4), pp. 495-504. https://doi.org/10.2307/2348134
Chen, Z. and Kuo, L. (2001). A note on the estimation of the
multinomial logit model with random effects. The American
Statistician 55 (2), pp. 89-95. https://www.jstor.org/stable/2685993
Guimaraes, P. (2004). Understanding the multinomial-Poisson
transformation. The Stata Journal 4 (3), pp. 265-273. https://www.stata-journal.com/article.html?article=st0069
Lee, J.Y.L., Green, P.J.,and Ryan, L.M. (2017). On the “Poisson
trick” and its extensions for fitting multinomial regression models.
arXiv preprint available at https://doi.org/10.48550/arXiv.1707.08538
Libraries
These are the libraries needed for running the code in this
vignette:
library(car)
library(nnet)
library(lme4)
library(lmerTest)
library(emmeans)
Let’s also load our library:
Between-subjects 2×2 design with dichotomous response
Let’s load and prepare our first data set, a 2×2 between-subjects
design with a dichotomous response. Factor X1 has levels
{a, b}, factor X2 has levels
{c, d}, and response Y has categories
{yes, no}.
data(bs2, package="multpois")
bs2$PId = factor(bs2$PId)
bs2$Y = factor(bs2$Y, levels=c("yes","no"))
bs2$X1 = factor(bs2$X1)
bs2$X2 = factor(bs2$X2)
contrasts(bs2$X1) <- "contr.sum"
contrasts(bs2$X2) <- "contr.sum"
Let’s visualize this data set using a mosaic plot:
xt = xtabs( ~ X1 + X2 + Y, data=bs2)
mosaicplot(xt, main="Y by X1, X2", las=1, col=c("lightgreen","pink"))
Figure 1. Proportions of
yes (green) and no (pink) responses in four
conditions: {a, c}, {a, d},
{b, c}, and {b, d}.
Given X1 and X2 are both between-subjects
factors, and Y is a dichotomous response, we can analyze
this data set using conventional logistic regression:
m1 = glm(Y ~ X1*X2, data=bs2, family=binomial)
Anova(m1, type=3)
#> Analysis of Deviance Table (Type III tests)
#>
#> Response: Y
#> LR Chisq Df Pr(>Chisq)
#> X1 4.4440 1 0.03502 *
#> X2 0.7513 1 0.38606
#> X1:X2 4.4440 1 0.03502 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(m1, pairwise ~ X1*X2, adjust="holm")$contrasts
#> contrast estimate SE df z.ratio p.value
#> a c - b c 2.398 0.870 Inf 2.755 0.0352
#> a c - a d 1.705 0.801 Inf 2.129 0.1661
#> a c - b d 1.705 0.801 Inf 2.129 0.1661
#> b c - a d -0.693 0.847 Inf -0.819 1.0000
#> b c - b d -0.693 0.847 Inf -0.819 1.0000
#> a d - b d 0.000 0.775 Inf 0.000 1.0000
#>
#> Results are given on the log odds ratio (not the response) scale.
#> P value adjustment: holm method for 6 tests
We can also analyze this data set using the multinomial-Poisson
trick, which converts nominal responses to category counts and analyzes
these counts using Poisson regression:
m2 = glm.mp(Y ~ X1*X2, data=bs2)
Anova.mp(m2, type=3)
#> Analysis of Deviance Table (Type III tests)
#>
#> Response: Y
#> via the multinomial-Poisson trick
#> Chisq Df N Pr(>Chisq)
#> X1 4.4440 1 60 0.03502 *
#> X2 0.7513 1 60 0.38606
#> X1:X2 4.4440 1 60 0.03502 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
glm.mp.con(m2, pairwise ~ X1*X2, adjust="holm")
#> $heading
#> [1] "Pairwise comparisons via the multinomial-Poisson trick"
#>
#> $contrasts
#> Contrast Chisq Df N p.value
#> 1 a.c - a.d 4.962518 1 30 0.129510
#> 2 a.c - b.c 9.045871 1 30 0.015798
#> 3 a.c - b.d 4.962518 1 30 0.129510
#> 4 a.d - b.c 0.687412 1 30 1.000000
#> 5 a.d - b.d 0.000000 1 30 1.000000
#> 6 b.c - b.d 0.687412 1 30 1.000000
#>
#> $notes
#> [1] "P value adjustment: holm method for 6 tests"
The omnibus results from logistic regression and from the
multinomial-Poisson trick match, and the results from the post
hoc pairwise comparisons are similar.
Between-subjects 2×2 design with polytomous response
Let’s load and prepare our second data set, a 2×2 between-subjects
design with a polytomous response. Factor X1 has levels
{a, b}, factor X2 has levels
{c, d}, and response Y has categories
{yes, no, maybe}.
data(bs3, package="multpois")
bs3$PId = factor(bs3$PId)
bs3$Y = factor(bs3$Y, levels=c("yes","no","maybe"))
bs3$X1 = factor(bs3$X1)
bs3$X2 = factor(bs3$X2)
contrasts(bs3$X1) <- "contr.sum"
contrasts(bs3$X2) <- "contr.sum"
Let’s again visualize the data using a mosaic plot:
xt = xtabs( ~ X1 + X2 + Y, data=bs3)
mosaicplot(xt, main="Y by X1, X2", las=1, col=c("lightgreen","pink","lightyellow"))
Figure 2. Proportions of
yes (green), no (pink), and maybe
(yellow) responses in four conditions: {a, c},
{a, d}, {b, c}, and {b, d}.
Given X1 and X2 are both between-subjects
factors, and Y is a polytomous response, we might wish that
glm had a family=multinomial option analogous
to its family=binomial option, but it does not.
Fortunately, we can analyze polytomous response data for (only)
between-subjects factors using the multinom function from
the nnet package:
m3 = multinom(Y ~ X1*X2, data=bs3, trace=FALSE)
Anova(m3, type=3)
#> Analysis of Deviance Table (Type III tests)
#>
#> Response: Y
#> LR Chisq Df Pr(>Chisq)
#> X1 3.5327 2 0.17096
#> X2 7.8081 2 0.02016 *
#> X1:X2 4.0039 2 0.13507
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Unfortunately, emmeans does not work straightforwardly
with multinom models. A solution to this issue from Russ
Lenth, lead author of emmeans, was posted
on StackExchange:
e0 = emmeans(m3, ~ X1*X2 | Y, mode="latent")
c0 = contrast(e0, method="pairwise", ref=1)
test(c0, joint=TRUE, by="contrast")
#> contrast df1 df2 F.ratio p.value note
#> a c - b c 2 8 3.017 0.1056 d
#> a c - a d 2 8 4.552 0.0479 d
#> a c - b d 2 8 4.610 0.0466 d
#> b c - a d 2 8 0.688 0.5298 d
#> b c - b d 2 8 0.611 0.5661 d
#> a d - b d 2 8 1.308 0.3224 d
#>
#> d: df1 reduced due to linear dependence
We can also analyze this data set using the multinomial-Poisson
trick:
m4 = glm.mp(Y ~ X1*X2, data=bs3)
Anova.mp(m4, type=3)
#> Analysis of Deviance Table (Type III tests)
#>
#> Response: Y
#> via the multinomial-Poisson trick
#> Chisq Df N Pr(>Chisq)
#> X1 3.5327 2 60 0.17096
#> X2 7.8081 2 60 0.02016 *
#> X1:X2 4.0039 2 60 0.13507
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
glm.mp.con(m4, pairwise ~ X1*X2, adjust="holm")
#> $heading
#> [1] "Pairwise comparisons via the multinomial-Poisson trick"
#>
#> $contrasts
#> Contrast Chisq Df N p.value
#> 1 a.c - a.d 12.172660 2 30 0.013644
#> 2 a.c - b.c 6.990329 2 30 0.121376
#> 3 a.c - b.d 11.647010 2 30 0.014785
#> 4 a.d - b.c 1.425017 2 30 0.980826
#> 5 a.d - b.d 2.804595 2 30 0.738093
#> 6 b.c - b.d 1.252756 2 30 0.980826
#>
#> $notes
#> [1] "P value adjustment: holm method for 6 tests"
Again, the results from multinomial logistic regression and from the
multinomial-Poisson trick match. The results from the post hoc
pairwise comparisons are similar.
Within-subjects 2×2 design with dichotomous response
Let’s load and prepare our third data set, a 2×2 within-subjects
design with a dichotomous response. Factor X1 has levels
{a, b}, factor X2 has levels
{c, d}, and response Y has categories
{yes, no}. Now the PId factor is repeated
across rows, indicating participants were measured repeatedly.
data(ws2, package="multpois")
ws2$PId = factor(ws2$PId)
ws2$Y = factor(ws2$Y, levels=c("yes","no"))
ws2$X1 = factor(ws2$X1)
ws2$X2 = factor(ws2$X2)
contrasts(ws2$X1) <- "contr.sum"
contrasts(ws2$X2) <- "contr.sum"
Let’s visualize this data set using a mosaic plot:
xt = xtabs( ~ X1 + X2 + Y, data=ws2)
mosaicplot(xt, main="Y by X1, X2", las=1, col=c("lightgreen","pink"))
Figure 3. Proportions of
yes (green) and no (pink) responses in four
conditions: {a, c}, {a, d},
{b, c}, and {b, d}.
Given X1 and X2 are both within-subjects
factors, and Y is a dichotomous response, we can analyze
this using mixed-effects logistic regression. The function
glmer from the lme4 package provides this to
us:
m5 = glmer(Y ~ X1*X2 + (1|PId), data=ws2, family=binomial)
Anova(m5, type=3)
#> Analysis of Deviance Table (Type III Wald chisquare tests)
#>
#> Response: Y
#> Chisq Df Pr(>Chisq)
#> (Intercept) 0.8553 1 0.355052
#> X1 0.8553 1 0.355052
#> X2 6.6368 1 0.009989 **
#> X1:X2 4.3758 1 0.036452 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(m5, pairwise ~ X1*X2, adjust="holm")$contrasts
#> contrast estimate SE df z.ratio p.value
#> a c - b c -0.693 0.847 Inf -0.819 0.8258
#> a c - a d -2.773 0.913 Inf -3.037 0.0143
#> a c - b d -0.981 0.833 Inf -1.177 0.7176
#> b c - a d -2.079 0.847 Inf -2.456 0.0702
#> b c - b d -0.288 0.760 Inf -0.378 0.8258
#> a d - b d 1.792 0.833 Inf 2.150 0.1262
#>
#> Results are given on the log odds ratio (not the response) scale.
#> P value adjustment: holm method for 6 tests
We can also analyze this data set using the multinomial-Poisson
trick, now with an underlying mixed-effects Poisson regression
model:
m6 = glmer.mp(Y ~ X1*X2 + (1|PId), data=ws2)
Anova.mp(m6, type=3)
#> Analysis of Deviance Table (Type III Wald chisquare tests)
#>
#> Response: Y
#> via the multinomial-Poisson trick
#> Chisq Df N Pr(>Chisq)
#> X1 0.8553 1 60 0.355052
#> X2 6.6368 1 60 0.009989 **
#> X1:X2 4.3758 1 60 0.036452 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
glmer.mp.con(m6, pairwise ~ X1*X2, adjust="holm")
#> $heading
#> [1] "Pairwise comparisons via the multinomial-Poisson trick"
#>
#> $contrasts
#> Contrast Chisq Df N p.value
#> 1 a.c - a.d 9.224701 1 30 0.014328
#> 2 a.c - b.c 0.670400 1 30 0.825824
#> 3 a.c - b.d 1.385304 1 30 0.717600
#> 4 a.d - b.c 6.033594 1 30 0.070180
#> 5 a.d - b.d 4.622979 1 30 0.126184
#> 6 b.c - b.d 0.143240 1 30 0.825824
#>
#> $notes
#> [1] "P value adjustment: holm method for 6 tests"
The results from mixed-effects logistic regression and results from
the multinomial-Poisson trick match, including the results from the
post hoc pairwise comparisons.
Within-subjects 2×2 design with polytomous response
This fourth example illustrates the main reason that the
multpois package was created. Unlike the three examples
above, there are no straightforward options for analyzing nominal
responses with repeated measures and obtaining ANOVA-style results. Some
functions do offer mixed-effects multinomial regression modeling, such
as mblogit in the mclogit package, but they do
not enable ANOVA-style output. Other advanced methods exist, such as
Markov Chain Monte Carlo (MCMC) methods in the MCMCglmm
package, which does have a family=multinomial option, but
these Bayesian methods are complex and deviate from the approaches
illustrated above. Fortunately, we can again use the multinomial-Poisson
trick.
Let’s load and prepare our fourth data set, a 2×2 within-subjects
design with a polytomous response. Factor X1 has levels
{a, b}, factor X2 has levels
{c, d}, and response Y has categories
{yes, no, maybe}. Again, the PId factor is
repeated across rows, indicating participants were measured
repeatedly.
data(ws3, package="multpois")
ws3$PId = factor(ws3$PId)
ws3$Y = factor(ws3$Y, levels=c("yes","no","maybe"))
ws3$X1 = factor(ws3$X1)
ws3$X2 = factor(ws3$X2)
contrasts(ws3$X1) <- "contr.sum"
contrasts(ws3$X2) <- "contr.sum"
Let’s visualize this data set using a mosaic plot:
xt = xtabs( ~ X1 + X2 + Y, data=ws3)
mosaicplot(xt, main="Y by X1, X2", las=1, col=c("lightgreen","pink","lightyellow"))
Figure 4. Proportions of
yes (green), no (pink), and maybe
(yellow) responses in four conditions: {a, c},
{a, d}, {b, c}, and {b, d}.
Because multinom from the nnet package
cannot accept random factors, it cannot model repeated measures. And
because glmer from the lme4 package has no
family=multinomial option, it cannot model polytomous
responses. Fortunately, with the multinomial-Poisson trick, we can
analyze polytomous responses from repeated measures:
m7 = glmer.mp(Y ~ X1*X2 + (1|PId), data=ws3)
Anova.mp(m7, type=3)
#> Analysis of Deviance Table (Type III Wald chisquare tests)
#>
#> Response: Y
#> via the multinomial-Poisson trick
#> Chisq Df N Pr(>Chisq)
#> X1 6.6734 2 60 0.03555 *
#> X2 6.6713 2 60 0.03559 *
#> X1:X2 0.5514 2 60 0.75905
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
glmer.mp.con(m7, pairwise ~ X1*X2, adjust="holm")
#> $heading
#> [1] "Pairwise comparisons via the multinomial-Poisson trick"
#>
#> $contrasts
#> Contrast Chisq Df N p.value
#> 1 a.c - a.d 6.033960 2 30 0.244735
#> 2 a.c - b.c 6.034054 2 30 0.244735
#> 3 a.c - b.d 10.676280 2 30 0.028830
#> 4 a.d - b.c 0.000000 2 30 1.000000
#> 5 a.d - b.d 1.590911 2 30 1.000000
#> 6 b.c - b.d 1.589577 2 30 1.000000
#>
#> $notes
#> [1] "P value adjustment: holm method for 6 tests"
Mixed factorial 2×2 design with polytomous response
This fifth and final example is also the reason that the
multpois package was created, since we have a polytomous
response, one between-subjects factor, and one within-subjects factor.
This mixed factorial design is also known as a split-plot design. (Note:
Do not confuse mixed factorial designs with mixed-effects models. The
former contain between- and within-subjects factors; the latter contain
fixed and random effects.)
This fictional data is based on the scenario at the beginning of this
vignette. Forty respondents, half adults and half children, were
surveyed for their favorite ice cream four times, once per season. Thus,
Age is a between-subjects factor with two levels
{adult, child}, and Season is a
within-subjects factor with four levels
{fall, winter, spring, summer}. The polytomous response,
Pref, has three categories:
{vanilla, chocolate, strawberry}. The PId
factor is repeated across rows, indicating respondents were queried four
times each, once per season.
Let’s load and prepare this data set:
data(icecream, package="multpois")
icecream$PId = factor(icecream$PId)
icecream$Pref = factor(icecream$Pref, levels=c("vanilla","chocolate","strawberry"))
icecream$Age = factor(icecream$Age, levels=c("adult","child"))
icecream$Season = factor(icecream$Season, levels=c("fall","winter","spring","summer"))
contrasts(icecream$Age) <- "contr.sum"
contrasts(icecream$Season) <- "contr.sum"
Let’s visualize this data set using a mosaic plot:
xt = xtabs( ~ Age + Season + Pref, data=icecream)
mosaicplot(xt, main="Pref by Age, Season", las=1, col=c("beige","tan","pink"))
Figure 5. Proportions of
vanilla (beige), chocolate (brown), and
strawberry (pink) responses for adults and children across
the four seasons.
As in the previous example, we can use the multinomial-Poisson trick
to analyze repeated measures data with polytomous responses:
m8 = glmer.mp(Pref ~ Age*Season + (1|PId), data=icecream)
Anova.mp(m8, type=3)
#> Analysis of Deviance Table (Type III Wald chisquare tests)
#>
#> Response: Pref
#> via the multinomial-Poisson trick
#> Chisq Df N Pr(>Chisq)
#> Age 8.8866 2 160 0.011757 *
#> Season 12.4627 6 160 0.052408 .
#> Age:Season 18.3213 6 160 0.005477 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We have a main effect of Age and an
Age×Season interaction but no main effect of
Season. We can explore this further by graphically
depicting response proportions in each age group:
xt = xtabs( ~ Age + Pref, data=icecream)
mosaicplot(xt, main="Pref by Age", las=1, col=c("beige","tan","pink"))
Figure 6. Proportions of
vanilla (beige), chocolate (brown), and
strawberry (pink) responses for adults and children. The
main effect of Age emerges, with children preferring
chocolate more and strawberry less than adults.
The different proportions by Age clearly emerge,
explaining the main effect. Let’s also graphically depict the
proportions by Season:
xt = xtabs( ~ Season + Pref, data=icecream)
mosaicplot(xt, main="Pref by Season", las=1, col=c("beige","tan","pink"))
Figure 7. Proportions of
vanilla (beige), chocolate (brown), and
strawberry (pink) responses by Season.
Although there are some differences in proportion, they are not quite
statistically significant (p = 0.052).
Finally, we can again conduct post hoc pairwise comparisons.
Note, however, there are many such possible comparisons, and best
practice would require us to only conduct those comparisons driven by
hypotheses or planned in advance. For example, we might wish to limit
our pairwise comparisons to adults vs. children within each
season, not across all seasons. In any case, we first conduct all
pairwise comparisons just for illustration:
glmer.mp.con(m8, pairwise ~ Age*Season, adjust="holm")
#> $heading
#> [1] "Pairwise comparisons via the multinomial-Poisson trick"
#>
#> $contrasts
#> Contrast Chisq Df N p.value
#> 1 adult.fall - adult.spring 9.050067 2 40 0.260016
#> 2 adult.fall - adult.summer 0.425876 2 40 1.000000
#> 3 adult.fall - adult.winter 4.153318 2 40 1.000000
#> 4 adult.fall - child.fall 8.128467 2 40 0.377872
#> 5 adult.fall - child.spring 1.222476 2 40 1.000000
#> 6 adult.fall - child.summer 2.136531 2 40 1.000000
#> 7 adult.fall - child.winter 3.642132 2 40 1.000000
#> 8 adult.spring - adult.summer 7.611634 2 40 0.447111
#> 9 adult.spring - adult.winter 14.903020 2 40 0.015820
#> 10 adult.spring - child.fall 14.957280 2 40 0.015820
#> 11 adult.spring - child.spring 13.775930 2 40 0.026520
#> 12 adult.spring - child.summer 13.630800 2 40 0.027425
#> 13 adult.spring - child.winter 8.577409 2 40 0.315629
#> 14 adult.summer - adult.winter 3.697318 2 40 1.000000
#> 15 adult.summer - child.fall 7.059671 2 40 0.556890
#> 16 adult.summer - child.spring 2.676090 2 40 1.000000
#> 17 adult.summer - child.summer 2.026498 2 40 1.000000
#> 18 adult.summer - child.winter 1.974121 2 40 1.000000
#> 19 adult.winter - child.fall 1.368363 2 40 1.000000
#> 20 adult.winter - child.spring 3.817685 2 40 1.000000
#> 21 adult.winter - child.summer 0.429746 2 40 1.000000
#> 22 adult.winter - child.winter 2.355143 2 40 1.000000
#> 23 child.fall - child.spring 7.698948 2 40 0.447111
#> 24 child.fall - child.summer 2.976095 2 40 1.000000
#> 25 child.fall - child.winter 4.253790 2 40 1.000000
#> 26 child.spring - child.summer 2.069340 2 40 1.000000
#> 27 child.spring - child.winter 6.222227 2 40 0.801918
#> 28 child.summer - child.winter 2.146742 2 40 1.000000
#>
#> $notes
#> [1] "P value adjustment: holm method for 28 tests"
If we wished to compare adults vs. children in each season (fall,
winter, spring, and summer), we would first conduct all pairwise
comparisons, leaving them uncorrected for multiple
comparisons…
glmer.mp.con(m8, pairwise ~ Age*Season, adjust="none")
#> $heading
#> [1] "Pairwise comparisons via the multinomial-Poisson trick"
#>
#> $contrasts
#> Contrast Chisq Df N p.value
#> 1 adult.fall - adult.spring 9.050067 2 40 0.010834
#> 2 adult.fall - adult.summer 0.425876 2 40 0.808206
#> 3 adult.fall - adult.winter 4.153318 2 40 0.125348
#> 4 adult.fall - child.fall 8.128467 2 40 0.017176
#> 5 adult.fall - child.spring 1.222476 2 40 0.542679
#> 6 adult.fall - child.summer 2.136531 2 40 0.343604
#> 7 adult.fall - child.winter 3.642132 2 40 0.161853
#> 8 adult.spring - adult.summer 7.611634 2 40 0.022241
#> 9 adult.spring - adult.winter 14.903020 2 40 0.000581
#> 10 adult.spring - child.fall 14.957280 2 40 0.000565
#> 11 adult.spring - child.spring 13.775930 2 40 0.001020
#> 12 adult.spring - child.summer 13.630800 2 40 0.001097
#> 13 adult.spring - child.winter 8.577409 2 40 0.013723
#> 14 adult.summer - adult.winter 3.697318 2 40 0.157448
#> 15 adult.summer - child.fall 7.059671 2 40 0.029310
#> 16 adult.summer - child.spring 2.676090 2 40 0.262358
#> 17 adult.summer - child.summer 2.026498 2 40 0.363038
#> 18 adult.summer - child.winter 1.974121 2 40 0.372670
#> 19 adult.winter - child.fall 1.368363 2 40 0.504503
#> 20 adult.winter - child.spring 3.817685 2 40 0.148252
#> 21 adult.winter - child.summer 0.429746 2 40 0.806644
#> 22 adult.winter - child.winter 2.355143 2 40 0.308026
#> 23 child.fall - child.spring 7.698948 2 40 0.021291
#> 24 child.fall - child.summer 2.976095 2 40 0.225813
#> 25 child.fall - child.winter 4.253790 2 40 0.119207
#> 26 child.spring - child.summer 2.069340 2 40 0.355344
#> 27 child.spring - child.winter 6.222227 2 40 0.044551
#> 28 child.summer - child.winter 2.146742 2 40 0.341854
#>
#> $notes
#> [1] "P value adjustment: none method for 28 tests"
…and then we would extract the relevant comparisons (rows 4, 22, 11,
and 17, respectively), and manually correct their p-values to
guard against Type I errors, like so:
p.adjust(c(0.017176, 0.308026, 0.001020, 0.363038), method="holm")
#> [1] 0.051528 0.616052 0.004080 0.616052
Thus, after correction using Holm’s sequential Bonferroni procedure
(Holm 1979), we see
that adults vs. children in spring are significantly
different (p < .05). Looking again at Figure 5 visually
confirms this result.
Copyright (C) 2024-2025 Jacob O. Wobbrock <wobbrock@uw.edu>