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, 

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

  • $\begingroup$ 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? $\endgroup$ – AdamO Apr 5 '18 at 13:38
  • $\begingroup$ @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. $\endgroup$ – GiannisZ Apr 6 '18 at 11:00
  • $\begingroup$ nobody ? hmm... $\endgroup$ – GiannisZ Apr 9 '18 at 11:51

Your Answer

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

Browse other questions tagged or ask your own question.