# Power calculation for mixed models

My question has to do (obviously) with how to calculate the required sample size for longitudinal data analysis with mixed models.

There is an older post here : Sample size calculation for mixed models

However it is from 2 years ago and most importantly i do not have the reputation yet to comment...!So, i thought to ask a new question.

I believe that the best way to to do a power calculation, especially for more complex models, is through simulation. And this is what i did, and you can find the (rather simplified) code below.

Sims1 <- function(n, m, sigma, rho, beta0, beta1, beta2, beta3, tau0, tau1,
tau01){

ctrl <- lmeControl(opt='optim')

# The var-cov matrix
delta2 <- c(0,1,rep(2,m-2))
t <- cumsum(delta2)
Sigma <- sigma^2*rho^abs(outer(t,t,"-"))

MU <- rep(0, m)  # The mean vector

Resds <- rmvnorm(n, MU, Sigma) # Draw n times for a multivariate normal

# Create the treatment variable
rr <- rep(c("Reference","Treatment"), each = n/2)
treatment <- sample(rr)

t1 = round(runif(n, m, m))  # a helper vector
# Time vector
time2 <- t

dat <- data.frame(id=rep(1:n, each=m), time=time2, treatment =  rep(treatment, each = m))

### set up data frame

dat$eij <- as.vector(t(Resds)) ### simulate (correlated) random effects for intercepts and slopes mu <- c(0,0) S <- matrix(c(1, tau01, tau01, 1), nrow=2) tau <- c(tau0, tau1) S <- diag(tau) %*% S %*% diag(tau) U <- mvrnorm(n, mu=mu, Sigma=S) # And finally the responses dat$yij <- (beta0 + rep(U[,1], times = t1)) + (beta1 + rep(U[,2], times=t1)) * dat$time + beta2*(as.numeric(dat$treatment)-1) + beta3*dat$time*(as.numeric(dat$treatment)-1) + dat$eij # And build the models... # Full REML model assuming independence res1 <- lme(yij ~ time*treatment, random = ~ time | id,data=dat,control=ctrl) # AR(1) model res <- lme(yij ~ time*treatment, random = ~ time | id, correlation = corCAR1(form = ~ time | id), data=dat, control=ctrl) res_phi <- summary(res)$model$corStr # Satterhwaite correction model res22 <- lmerTest::lmer( yij ~ time*treatment + (time|id), data = dat) # Return the interesting parts return(list(MLM = summary(res1)$tTable, MLM_var = as.numeric(VarCorr(summary(res1))[,2]),
MLM_S = summary(res22)$coeff, MLM_full = summary(res)$tTable, MLM_full_var =   c(as.numeric(VarCorr(summary(res))[,2]), as.numeric(coef(res_phi,      unconstrained=FALSE)))))

}

This code seems to work as expected, since i am able to "capture" the values i specify as inputs(all of them). Now, I repeat the following scenario 2000 times:

Sims1( n = 20, m = 10,beta0 = 2.39, beta1 = 0.208, beta2 = -0.26, beta3 = -0.0795, tau0 = 0.619, tau1 = 0, tau01 = 0, sigma = 0.788, rho = 0)

And i calculate how often i reject the null hypothesis for the interaction term (beta3 coefficient = Time*Treatment). That comes out to be about 95%.

Then I use exactly the same "real" values to do a power calculation with a software that it is specifically for this purpose, called SPA-ML and it is from the book by Mijam Moerbeek : https://www.crcpress.com/Power-Analysis-of-Trials-with-Multilevel-Data/Moerbeek-Teerenstra/p/book/9781498729895

To my surprise, to achieve a power of 80% based on that program requires N=950 per group!! And this is for exactly the same "real" parameters i used to do the power calculation through simulations...

And right now i am really confused...On one hand, my simulations seems to work perfect and on the other hand i have a software and a book that says rather different thing.... Do you have any idea about that huge different ? Or maybe another way to do a power calculation that is proved to work, so i can validate my results ?

Thanks, John

• I seem to recall that to obtain inference on mixed models, one should fit a null model and test nested models with ANOVA (see ?anova.merMod and especially the refit options). I think this is because the variance of the random effects changes under different parameterizations. Also, can you summarize the inputs you used to the software/book? Are you sure you're testing growth (the time by treatment interaction) and not, say, the time-averaged treatment effect? – AdamO Apr 5 '18 at 13:38
• @AdamO . As for your first remark, i understand what you mean, but i am using the nlme package which provides the p-values. And here i just want an indication of what is going on and the difference in the 2 procedures, of course cannot be explained by that! But indeed i get your point! For the second issue, the software uses the 2-level variances( within & between), the number of measurements per subject, and the parameter estimate of interest( here the interaction) but in a standardized form. This is done by dividing the estimate by the variance of the random intercept component. – GiannisZ Apr 6 '18 at 11:00
• nobody ? hmm... – GiannisZ Apr 9 '18 at 11:51