# Power calculation in R: Multilevel model

I need to estimate the statistical power of a multilevel test in order to estimate the sample size needed. The levels of the model include country, center and individual person. I will be measuring the difference in treatment effect between two treatments.

I found the following R code on http://www.stat.columbia.edu/~gelman/stuff_for_blog/chap20.pdf

fake<-function(J,K){
time<-rep(seq(0,1,length=K),J)
person<-rep(1:J,each=K)
treatment<-sample(rep(0:1,J/2))
treatment1<-treatment[person]

mu.a.true <- 4.8
g.0.true <- -0.5
g.1.true <- 0.5
sigma.y.true <- 0.7
sigma.a.true <- 1.3
sigma.b.true <- 0.7

a.true<-rnorm(J,mu.a.true,sigma.a.true)
b.true<-rnorm(J,g.0.true+g.1.true*treatment,sigma.b.true)

y<-rnorm(J*K,a.true[person]+b.true[person]*time,sigma.y.true)
return(data.frame(y,time,person,treatment1))
}

power<-function(J,K,n.sims=1000){
signif<-rep(NA,n.sims)
for (s in 1:n.sims){
fake<-fake(J,K)
lme.power<-lmer(y~time+time:treatment1+(1+time|person),data=fake)
theta.hat<-fixef(lme.power)["time:treatment1"]
theta.se<-se.fixef(lme.power)["time:treatment1"]
signif[s]<-(theta.hat-2*theta.se)>0
}
power<-mean(signif)
return(power)
}


My intention was to edit this to my needs. They use the two functions above to estimate, for different numbers of people, the power level. They explain that this is a multilevel estimate but as the treatment is at the group level I don't think this is strictly true. Please could someone explain the multilevel part of the code and advise me on editing it as required?

-

## migrated from stackoverflow.comApr 3 at 13:54

This question came from our site for professional and enthusiast programmers.