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)
summary()to an object extract the p values from the parametric table and pass them top.adjustand choose one of the options it provides. $\endgroup$shapeand 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 withglht()from the multcomp package but themodel.matrixmethod forgamobjects doesn't contain the informationglht()needs to work (although to some extent it is in the gam object). $\endgroup$lme4::lmer()and then useglht()to do the Tukey contrasts and family-wise adjustment for all pairwise comparisons. $\endgroup$glht()works without problem forgam()output but themcp()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. $\endgroup$