2
$\begingroup$

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

$\endgroup$
5
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ Thanks Dmiitris, I'm looking forward to going through your course material soon once I free up my schedule! $\endgroup$ – J Doe Jan 30 at 15:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.