Power Analysis for glmer using simr

Question

I'm trying to run a power analysis based on pilot study. I'm using a glmer mixed effects model with a binomial logit link. The model's equation is:

success ~ group + (1|user_id)

I have about 2000 observations per each group (variable number of observations per user). I'd like to estimate the required sample size per group required to have 80% Power to detect an effect size of 0.1. I'd like to extend the power analysis to atleast 20000 observations per group to see how n observations impacts power.

model = glmer("success ~ group + (1|user_id)", df, family = binomial(link='logit'))

set_fixef = function(o, s, v){
fixef(o)[s] <- v
return(o)
}

power_model = set_fixef(model, 'groupexperiment', 0.1)
power1 = powerCurve(extend(power_model, along = "user_id", n = 20), test=fixed("groupexperiment", "z"), along = "user_id")

This code isn't really doing much for me, I'm having trouble figuring out how to really estimate Power for this type of model.

The idea is that I'm mimicking a two proportion Z test, but with a mixed effects model to account for multiple observations per user

Answer 1

Indeed, the best way to estimate the power in mixed models is using simulation. The following generic code shows how this can be done in R using the GLMMadaptive package. You can suitably adapt it to fit your needs:

simulate_binary <- function (n) {
    K <- 8 # 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 random follow-up times
    DF <- data.frame(id = rep(seq_len(n), each = K),
                     time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
                     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(~ time, data = DF)

    betas <- c(-2.13, -0.25, 0.24, -0.05) # fixed effects coefficients
    D11 <- 0.48 # variance of random intercepts
    D22 <- 0.1 # variance of random slopes

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

###################################################################

library("GLMMadaptive")
M <- 1000 # number of simulations to estimate power
p_values <- numeric(M)
for (m in seq_len(M)) {
    DF_m <- simulate_binary(n = 100)
    fm_m <- mixed_model(y ~ sex * time, random = ~ time | id, 
                        data = DF_m, family = binomial())
    p_values[m] <- coef(summary(fm_m))["sexfemale:time", "p-value"]
}
# assuming a significance level of 5%, the power will be
mean(p_values < 0.05)