# Simulation for power calculation for hypothetical new cd4 experiment

# 1.  Function to simulate hypothetical data

cd4.fake <- function (J, K){
  time <- rep (seq(0,1,length=K), J) # K measurements during the year
  person <- rep (1:J, each=K)        # person ID's
  treatment <- sample (rep (0:1, J/2))
  treatment1 <- treatment[person]
#                                    # hyperparameters
  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
#                                    # person-level parameters
  a.true <- rnorm (J, mu.a.true, sigma.a.true)
  b.true <- rnorm (J, g.0.true + g.1.true*treatment, sigma.b.true)
#                                    # data
  y <- rnorm (J*K, a.true[person] + b.true[person]*time, sigma.y.true)
  return (data.frame (y, time, person, treatment1))
}

# 2.  Function to fit the model and calculate the power

cd4.power <- function (J, K, n.sims=1000){
  signif <- rep (NA, n.sims)
  for (s in 1:n.sims){
    fake <- cd4.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] <- ifelse (theta.hat - 2*theta.se > 0, 1, 0)
  }
  power <- mean (signif)
  return (power)
}

cd4.power (J=150, K=7, n.sims=10)

# 3.  Loop for power calculations

J.values <- c(15, 60, 100, 150, 200, 225, 250, 300, 400)
n.sims.values <- rep(1000,9)
K.values <- c(3,5,7,10)
power.values <- array (NA, c(length(J.values),length(K.values)))
for (i1 in 1:length(J.values)){
  for (i2 in 1:length(K.values)){
    cat ("computing power calculation for J =", J.values[i1], ", K =", K.values[i2], "\n")
    power.values[i1,i2] <- cd4.power (J=J.values[i1], K=K.values[i2], n.sims=n.sims.values[i1])
    cat ("power =", power.values[i1,i2], "\n")
  }
}

# plot all the curves

plot (c(0,max(J.values)), c(0,1), xaxs="i", yaxs="i", xlab="number of children", ylab="power", type="n")
for (i2 in 1:length(K.values)){
  lines (c(0,J.values), c(.025,power.values[,i2]))
}

# just plot the curve for K=7

postscript ("c:/books/multilevel/power.ps", height=3.6, width=4)
plot (c(0,J.values), c(.025,power.values[,3]), ylim=c(0,1), yaxs="i", xaxs="i", type="l",
  xlab="number of children", ylab="power", yaxt="n", mgp=c(1.6,.5,0))
axis (2, seq(0,1,.1),rep("",11), mgp=c(1.6,.5,0), tcl=-.2)
axis (2, seq(0,1,.5),c("0","0.5","1"), mgp=c(1.6,.5,0))
dev.off ()