Analyzing nominal responses using the multinomial-Poisson trick

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, we might ask people what their favorite ice cream flavor is: vanilla, chocolate, or strawberry. This four-category response would be a polytomous dependent variable. Perhaps we wish to ask adults and children about their favorite ice cream to see if there is a difference by age group. We would then have a two-level between-subjects factor. 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 change. Now we would have a four-level within-subjects factor, i.e., repeated measures.

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) 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 take random factors to handle repeated measures, and the glmer function in lme4 does not offer a family=multinomial option. It is this case in particular that this package was created to address, although it can address the above three scenarios, 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: Relevant academic references for this vignette.
Libraries: External R libraries needed for this vignette.
Between-subjects 2×2 design with dichotomous response: Analysis of the bs2 data set.
Between-subjects 2×2 design with polytomous response: Analysis of the bs3 data set.
Within-subjects 2×2 design with dichotomous response: Analysis of the ws2 data set.
Within-subjects 2×2 design with polytomous response: Analysis of the ws3 data set.
Mixed factorial 2×2 design with polytomous response: Analysis of the icecream data set.

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:

library(multpois)

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)
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("pink","lightgreen"))

Figure 1. Proportions of no (pink) and yes (green) 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 quite 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)
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("lightyellow","pink","lightgreen"))

Figure 2. Proportions of maybe (yellow), no (pink), and yes (green) 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 fairly 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)
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("pink","lightgreen"))

Figure 3. Proportions of no (pink) and yes (green) 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 for 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.224694  1 30 0.014328
#> 2 a.c - b.c 0.670400  1 30 0.825824
#> 3 a.c - b.d 1.385318  1 30 0.717591
#> 4 a.d - b.c 6.033595  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 is the reason that the multpois package was created. Unlike the three examples above, there are not 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 library, which does have a family=multinomial option, but these 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)
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("lightyellow","pink","lightgreen"))

Figure 4. Proportions of maybe (yellow), no (pink), and yes (green) responses in four conditions: {a, c}, {a, d}, {b, c}, and {b, d}.

Because multinom from the nnet package cannot handle 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.6707  2 60     0.0356 *
#> X2    6.6707  2 60     0.0356 *
#> X1:X2 0.5508  2 60     0.7593  
#> ---
#> 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.033961  2 30 0.244745
#> 2 a.c - b.c  6.033962  2 30 0.244745
#> 3 a.c - b.d 10.679680  2 30 0.028782
#> 4 a.d - b.c  0.000000  2 30 1.000000
#> 5 a.d - b.d  1.589577  2 30 1.000000
#> 6 b.c - b.d  1.589578  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 factors. This mixed factorial design is also known as a split-plot design.

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.

Let’s load and prepare this data set:

data(icecream, package="multpois")
icecream$PId = factor(icecream$PId)
icecream$Pref = factor(icecream$Pref)
icecream$Age = factor(icecream$Age)
icecream$Season = factor(icecream$Season)
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("tan","pink","beige"))

Figure 5. Proportions of chocolate (brown), strawberry (pink), and vanilla (beige) responses in eight conditions: {adult, fall}, {adult, spring}, {adult, summer}, {adult, winter}, {child, fall}, {child, spring}, {child, summer}, and {child, winter}.

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.9838  2 160   0.011199 * 
#> Season     12.4522  6 160   0.052609 . 
#> Age:Season 18.3118  6 160   0.005498 **
#> ---
#> 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("tan","pink","beige"))

Figure 6. Proportions of chocolate (brown), strawberry (pink), and vanilla (beige) 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("tan","pink","beige"))

Figure 7. Proportions of chocolate (brown), strawberry (pink), and vanilla (beige) responses by season. Although there are some differences in proportion, they are not quite statistically significant (p = 0.053).

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 in each season, not across all seasons. In any case, we first conduct all pairwise comparisons 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.050033  2 40 0.260040
#> 2    adult.fall - adult.summer  0.425905  2 40 1.000000
#> 3    adult.fall - adult.winter  4.152565  2 40 1.000000
#> 4      adult.fall - child.fall  8.128467  2 40 0.377872
#> 5    adult.fall - child.spring  1.222474  2 40 1.000000
#> 6    adult.fall - child.summer  2.136526  2 40 1.000000
#> 7    adult.fall - child.winter  3.642131  2 40 1.000000
#> 8  adult.spring - adult.summer  7.606836  2 40 0.447111
#> 9  adult.spring - adult.winter 14.904790  2 40 0.015820
#> 10   adult.spring - child.fall 14.957210  2 40 0.015820
#> 11 adult.spring - child.spring 13.775910  2 40 0.026520
#> 12 adult.spring - child.summer 13.630850  2 40 0.027425
#> 13 adult.spring - child.winter  8.577378  2 40 0.315629
#> 14 adult.summer - adult.winter  3.697246  2 40 1.000000
#> 15   adult.summer - child.fall  7.059707  2 40 0.556871
#> 16 adult.summer - child.spring  2.676089  2 40 1.000000
#> 17 adult.summer - child.summer  2.026498  2 40 1.000000
#> 18 adult.summer - child.winter  1.974120  2 40 1.000000
#> 19   adult.winter - child.fall  1.368365  2 40 1.000000
#> 20 adult.winter - child.spring  3.817680  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.698949  2 40 0.447111
#> 24   child.fall - child.summer  2.978168  2 40 1.000000
#> 25   child.fall - child.winter  4.252011  2 40 1.000000
#> 26 child.spring - child.summer  2.069326  2 40 1.000000
#> 27 child.spring - child.winter  6.222224  2 40 0.801918
#> 28 child.summer - child.winter  2.136527  2 40 1.000000
#> 
#> $notes
#> [1] "P value adjustment: holm method for 28 tests"

If we did wish to compare adults vs. children in each season (fall, winter, spring, and summer), we would first conduct all pairwise comparisons, leaving them uncorrected…

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.050033  2 40 0.010835
#> 2    adult.fall - adult.summer  0.425905  2 40 0.808195
#> 3    adult.fall - adult.winter  4.152565  2 40 0.125396
#> 4      adult.fall - child.fall  8.128467  2 40 0.017176
#> 5    adult.fall - child.spring  1.222474  2 40 0.542679
#> 6    adult.fall - child.summer  2.136526  2 40 0.343605
#> 7    adult.fall - child.winter  3.642131  2 40 0.161853
#> 8  adult.spring - adult.summer  7.606836  2 40 0.022294
#> 9  adult.spring - adult.winter 14.904790  2 40 0.000580
#> 10   adult.spring - child.fall 14.957210  2 40 0.000565
#> 11 adult.spring - child.spring 13.775910  2 40 0.001020
#> 12 adult.spring - child.summer 13.630850  2 40 0.001097
#> 13 adult.spring - child.winter  8.577378  2 40 0.013723
#> 14 adult.summer - adult.winter  3.697246  2 40 0.157454
#> 15   adult.summer - child.fall  7.059707  2 40 0.029309
#> 16 adult.summer - child.spring  2.676089  2 40 0.262358
#> 17 adult.summer - child.summer  2.026498  2 40 0.363038
#> 18 adult.summer - child.winter  1.974120  2 40 0.372671
#> 19   adult.winter - child.fall  1.368365  2 40 0.504502
#> 20 adult.winter - child.spring  3.817680  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.698949  2 40 0.021291
#> 24   child.fall - child.summer  2.978168  2 40 0.225579
#> 25   child.fall - child.winter  4.252011  2 40 0.119313
#> 26 child.spring - child.summer  2.069326  2 40 0.355346
#> 27 child.spring - child.winter  6.222224  2 40 0.044551
#> 28 child.summer - child.winter  2.136527  2 40 0.343605
#> 
#> $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.

Analyzing nominal responses using the multinomial-Poisson trick

Introduction

Contents

References

Libraries

Between-subjects 2×2 design with dichotomous response

Between-subjects 2×2 design with polytomous response

Within-subjects 2×2 design with dichotomous response

Within-subjects 2×2 design with polytomous response

Mixed factorial 2×2 design with polytomous response