Multiple Comparisons for MixMod Objects"

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

Multiple Comparisons with MixMod Objects

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

Additive Model

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



Try the GLMMadaptive package in your browser

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

GLMMadaptive documentation built on May 2, 2019, 2:51 p.m.