Power calculation for mixed models

Question

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