# Correcting for multiple pairwise comparisons with GAM objects {mgcv} in R

My question concerns pairwise comparisons of factor levels in a gam object. I have a dataframe, df, containing reaction time data (in ms) to stimuli varying in shape:

subjectID     RT     shape
001           501    square
002           722    circle
003           302    square
...


gam (from the mgcv package) offers a nice way to model this data while handling random-effects variables such as subjectID (i.e. controlling for 'random' variability between subjects):

m1 = gam(data = df, formula = lrt ~ Rp * slfreq + s(subjectID, bs = "re")


I can use summary(m1) to view pairwise contrasts of individual conditions (and corresponding p values) within the factor shape, however this method only gives uncorrected p values.

My question is this: is there a method for applying corrections to these p values (for example, a Tukey correction) when dealing with gam objects in R?

• If all you want to control error rates on those stated hypothesis tests, then you can just assign the output of summary() to an object extract the p values from the parametric table and pass them to p.adjust and choose one of the options it provides. – Gavin Simpson Nov 9 '18 at 20:46
• However, often the hypothesis tests are not the ones of interest (the default contrasts will set up the model with the intercept as referring to the reference level of shape and the other coefficients are deviations from that reference level. If you want to do all pairwise comparison and control the family-wise error rate of those comparisons, then you;re out of luck I think. You could do this with glht() from the multcomp package but the model.matrix method for gam objects doesn't contain the information glht() needs to work (although to some extent it is in the gam object). – Gavin Simpson Nov 9 '18 at 20:47
• You could easily fit this model in lme4::lmer() and then use glht() to do the Tukey contrasts and family-wise adjustment for all pairwise comparisons. – Gavin Simpson Nov 9 '18 at 20:48
• To be more precise: glht() works without problem for gam() output but the mcp() function for constructing the contrast matrix automatically does not. It is possible, though, to construct the matrix "by hand" I've described this in some detail in my new answer. – Achim Zeileis Nov 9 '18 at 22:15

The glht() function for generalized linear hypotheses from the multcomp package can be used to carry out various kinds of contrasts using a range of different p-value adjustments. The contrasts you are looking for are also called "Tukey" contrasts for all pairwise comparisons. The p-value adjustments include single-step, Shaffer, Westfall, and all p.adjust methods, see ?summary.glht.

As @GavinSimpson pointed out in the comments: For gam() objects from mgcv this does not work out of the box but requires some manual intervention. For lmer() from lme4 everything works conveniently. I illustrate below how both packages can be used with multcomp to obtain equivalent results. For illustration I use the sleepstudy data from lme4 but collapse the numeric regressor Days to a three-level factor (merely for illustration purposes):

library("lme4")
data("sleepstudy", package = "lme4")
sleepstudy$$Days <- cut(sleepstudy$$Days, breaks = c(-Inf, 2.5, 5.5, Inf),
labels = c("low", "med", "high"))
m1 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)
summary(m1)
## ...
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  262.170      9.802  26.747
## Daysmed       31.217      6.365   4.905
## Dayshigh      67.433      5.954  11.326
## ...


Then glht() can be used to set up all pairwise (aka Tukey) contrasts for the Days factor. The summary() method then applies the p-value adjustment (single-step, by default).

library("multcomp")
g1 <- glht(m1, linfct = mcp(Days = "Tukey"))
summary(g1)
##   Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
## Fit: lmer(formula = Reaction ~ Days + (1 | Subject), data = sleepstudy)
##
## Linear Hypotheses:
##                 Estimate Std. Error z value Pr(>|z|)
## med - low == 0    31.217      6.365   4.905 2.28e-06 ***
## high - low == 0   67.433      5.954  11.326  < 1e-06 ***
## high - med == 0   36.216      5.954   6.083  < 1e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)


The same model can be fitted with gam() as described in the question.

library("mgcv")
m2 <- gam(Reaction ~ Days + s(Subject, bs = "re"), data = sleepstudy)
summary(m2)
## ...
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  262.170      9.802  26.747  < 2e-16 ***
## Daysmed       31.217      6.365   4.905 2.27e-06 ***
## Dayshigh      67.433      5.954  11.326  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ...


However, the mcp(Days = "Tukey") method for describing the Tukey contrasts does not cooperate with gam() output and hence fails:

g2 <- glht(m2, linfct = mcp(Days = "Tukey"))
## Error in linfct[[nm]] %*% C :
##   requires numeric/complex matrix/vector arguments


However, it is not difficult (albeit a bit technical and tedious) to set up the contrast matrix by hand.

contr <- matrix(0, nrow = 3, ncol = length(coef(m2)))
colnames(contr) <- names(coef(m2))
rownames(contr) <- c("med - low", "high - low", "high - med")
contr[, 2:3] <- rbind(c(1, 0), c(0, 1), c(-1, 1))


The first columns of the contrast matrix show what is needed here: As the low coefficient is constrained to zero in the model, med - low is simply med and analogously for high - low. The last row then shows the contrast for high - med:

contr[, 1:5]
##            (Intercept) Daysmed Dayshigh s(Subject).1 s(Subject).2
## med - low            0       1        0            0            0
## high - low           0       0        1            0            0
## high - med           0      -1        1            0            0


And with this contrast matrix we can conduct the pairwise comparison with glht():

g2 <- glht(m2, linfct = contr)
summary(g2)
##   Simultaneous Tests for General Linear Hypotheses
##
## Fit: gam(formula = Reaction ~ Days + s(Subject, bs = "re"), data = sleepstudy)
##
## Linear Hypotheses:
##                 Estimate Std. Error z value Pr(>|z|)
## med - low == 0    31.217      6.365   4.905 2.35e-06 ***
## high - low == 0   67.433      5.954  11.326  < 1e-06 ***
## high - med == 0   36.216      5.954   6.083  < 1e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)


Another quite convenient way to indicate the contrasts to be tested is through character strings. This can set up linear functions based on the effect names from names(coef(m2)). And for factors with fewer levels (and hence fewer Tukey contrasts) this works quite nicely - but if the comparisons become more complex it's possibly easier to constract the contrast matrix as above.

g3 <- glht(m2, linfct = c("Daysmed = 0", "Dayshigh = 0", "Dayshigh - Daysmed = 0"))
summary(g3)
##   Simultaneous Tests for General Linear Hypotheses
##
## Fit: gam(formula = Reaction ~ Days + s(Subject, bs = "re"), data = sleepstudy)
##
## Linear Hypotheses:
##                         Estimate Std. Error z value Pr(>|z|)
## Daysmed == 0              31.217      6.365   4.905 2.53e-06 ***
## Dayshigh == 0             67.433      5.954  11.326  < 1e-06 ***
## Dayshigh - Daysmed == 0   36.216      5.954   6.083  < 1e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

• Achim, thank you for such a comprehensive answer! I have now implemented this solution (it worked perfectly!), but I have a couple of follow-up questions. To examine contrasts for interactions, can I just add a new contrast to 'contr' for every level of the interacting factor? Second, I cited Tukey in my question (which I assume is the default for glht), but can any other methods be used, such as Bonferoni, Holm, FDR...? Thanks! – Lyam Nov 9 '18 at 23:38
• Re: "Tukey". Be careful not to confuse the "Tukey contrasts" with the p-value adjustment method. The default method used by summary.glht is "single-step", others like Bonferroni or Holm are also available (but typically single-step would be preferred over these). For a list of options see ?summary.glht and ?p.adjust (as mentioned at the beginning of my response). Re: More levels. Yes, then you need to build a contrast matrix with more rows and adding suitable -1 and 1 entries. – Achim Zeileis Nov 10 '18 at 0:15
• One more follow-up question: does anyone know specifically which test is being used in the "single-step" method for summary.glht? My understanding is that single-step is a family of multiple comparison tests, not a single method (as opposed to Bonferroni, Holm, etc). I would like to know for reporting purposes. – Lyam Nov 12 '18 at 1:40
• Single-step computes the p-value by using the asymptotic normal distribution (or t distribution in case of lm) for the multivariate test statistic. See the details and references sections of ?summary.glht. Specifically, Hothorn et al. (2008, Biometrical Journal) and Bretz et al. (2010, CRC Press) provide the exact details. – Achim Zeileis Nov 12 '18 at 16:30