# Power Analysis for glmer using simr

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

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
}

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