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:

 Sims2( 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, max = 2000)

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