Power analysis for linear mixed effects model (lme4)

Lme4PowerAnalysis.R

# TO DO: ADJUST TO OTHER PARAMETERS (E:G: WHAT IF THERE ARE STIMULUS EFFECTS? OR MORE FIX EFFS?)
createNewData <- function (J,K,model){
  # J is the number of subjects
  # K is the number of visits
  # HP is the hyperparameter matrix from get Hyperparam
  # TO DO: LOOP THROUGH ALL FE ROWS AND AUTOMATICALLY EXTRACT NAMES OF FIXED EFFECTS AND ESTIMATES
  fe <- fixef(model)
  Intercept <- fe[1] #intercept
  bVisit <- fe[2] #visit
  bDiagnosis <- fe[3] #diagnosis
  bVisitDiagnosis <- fe[4] #visit diagnosis interaction
  # TO DO: WHAT ABOUT STANDARD ERROR? HOW SHOULD IT BE INCLUDED?
  
  # TO DO: LOOP THROUGH ALL VC COMPONENTS AND AUTOMATICALLY EXTRACT NAMES OF EFFECTS AND ESTIMATES
  vc<-VarCorr(model) # variance component
  sigmaSubject <- as.numeric(attr(vc[[1]],"stddev")[1]) # random intercept by subject
  sigmaVisit <- as.numeric(attr(vc[[1]],"stddev")[2]) # random slope of visit over subject
  sigmaResiduals <- as.numeric(attr(vc,"sc"))
  sigmaCorrelation <- as.numeric(attr(vc[[1]],"correlation")[2])
  
  d=expand.grid(Visit=1:K,Child.ID=1:J)
  condition <- sample(rep(0:1, J/2))
  d$Diagnosis<-condition[d$Child.ID]
  d$Diagnosis[is.na(d$Diagnosis)]<-1
  
  ## Define variance covariance matrices:
  Sigma.u<-matrix(c(sigmaSubject^2,
                    sigmaCorrelation*sigmaSubject*sigmaVisit,
                    sigmaCorrelation*sigmaSubject*sigmaVisit,
                    sigmaVisit^2),nrow=2)
  ## subj ranef
  u<-mvrnorm(n=J,
             mu=c(0,0),Sigma=Sigma.u)
  #u[,1]=scale(u[,1])
  #u[,2]=scale(u[,2])
  #summary(lm(u[,1]~u[,2]))
  
  
  d$CHI_MLU <- rnorm(J*K,
                     (Intercept+u[,1]) +
                     (bVisit+u[,2])*d$Visit + 
                     bDiagnosis*d$Diagnosis + 
                     ((bVisit+u[,2])*d$Visit)*(bDiagnosis*d$Diagnosis),sigmaResiduals)  
  
  return(d)
}