```knitr::opts_chunk\$set(
collapse = TRUE,
comment = "#>"
)
```

# Multiple Comparisons with MixMod Objects

In this vignette we illustrate how to correct p-values for multiple comparisons using the multcomp and emmeans packages.

We start by simulating some data for a binary longitudinal outcome:

```set.seed(1234)
n <- 300 # number of subjects
K <- 4 # number of measurements per subject
t_max <- 15 # maximum follow-up time

# we constuct a data frame with the design:
# everyone has a baseline measurment, and then measurements at K time points
DF <- data.frame(id = rep(seq_len(n), each = K),
time = gl(K, 1, n*K, labels = paste0("Time", 1:K)),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ 1, data = DF)

betas <- c(-2.13, 1, rep(c(1.2, -1.2), K-1)) # fixed effects coefficients
D11 <- 1 # variance of random intercepts

# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF\$id, ]))
# we simulate binary longitudinal data
DF\$y <- rbinom(n * K, 1, plogis(eta_y))
```

We fit a mixed effects logistic regression for `y` assuming random intercepts for the random-effects part, and the main effects of `sex` and `time` for the fixed-effects part.

```fm <- mixed_model(fixed = y ~ sex + time, random = ~ 1 | id, data = DF,
family = binomial())
```

The uncorrected p-values for the 4 time points are give by the `summary()` method:

```summary(fm)
```

To perform the pairwise comparisons and obtain corrected p-values, we load the multcomp package and use the `glht()` function. Because no specific methods exist for `MixMod` object returned by `mixed_model()`, we need to specify the `vcov.` and `coef.` arguments of `glht()`, i.e.,

```library("multcomp")
fm_mc <- glht(fm, linfct = mcp(time = "Tukey"),
vcov. = vcov(fm, "fixed"), coef. = fixef)

summary(fm_mc)
```

## Interaction Model

We continue our illustration by including the interaction term between `sex` and `time`, and we focus on the difference between males and females for the various time points. We start by fitting the model:

```gm <- mixed_model(fixed = y ~ sex * time, random = ~ 1 | id, data = DF,
family = binomial())
```

To compute the estimated log odds for males and females at the different time points we use the `emmeans()` functions from the emmeans package:

```library("emmeans")
gm_mc <- emmeans(gm, ~ sex | time)

gm_mc
```

The corresponding pairwise comparisons are performed by the `pairs()` function:

```pairs(gm_mc)
```

For additional examples in testing interactions with the emmeans package check the vignette: `vignette("interactions", package = "emmeans")`.