knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library("GLMMadaptive")

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)

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")`

.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.