knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library("GLMMadaptive")
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
time for the
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
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
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
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
emmeans() functions from the emmeans package:
library("emmeans") gm_mc <- emmeans(gm, ~ sex | time) gm_mc
The corresponding pairwise comparisons are performed by the
For additional examples in testing interactions with the emmeans package check the
vignette("interactions", package = "emmeans").
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.