getHyperparam <- function(x){
  sum_obj <- summary(x)
  # TO DO: LOOP THROUGH ALL FE ROWS AND AUTOMATICALLY EXTRACT NAMES OF FIXED EFFECTS AND ESTIMATES
  fe <- fixef(x)
  Intercept <- fe[1] #intercept
  Visit <- fe[2] #visit
  Diagnosis <- fe[3] #diagnosis
  VisitDiagnosis <- 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(x) # 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])

  # TO DO: AUTOMATICALLY ADJUST THE LIST
  hp <- c(Visit, Diagnosis,VisitDiagnosis,Intercept,sigmaResiduals,sigmaSubject,sigmaVisit,sigmaCorrelation)
  names(hp) <- c("Visit", "Diagnosis","VisitDiagnosis","Intercept","sigmaResiduals","sigmaSubject","sigmaVisit", "sigmaCorrelation")
  return(hp)
}

create_fake <- function (J,K,HP){
  # J is the number of subjects
  # K is the number of visits
  # HP is the hyperparameter matrix from get Hyperparam

  time <- rep (seq(1,6,length=K),J) # Define how many visits
  person <- rep (1:J, each=K) # Define how many kids

  # Add diagnosis
  condition <- sample (rep(0:1, J/2))
  condition1<-condition[person]
  condition1[is.na(condition1)]<-1

  Visit <- as.numeric(HP['Visit'])
  Diagnosis <- as.numeric(HP['Diagnosis'])
  VisitDiagnosis <- as.numeric(HP['VisitDiagnosis'])
  Intercept <- as.numeric(HP['Intercept'])
  sigmaSubject <- as.numeric(HP['sigmaSubject']) # random intercept by subject
  sigmaVisit <- as.numeric(HP['sigmaVisit']) # random slope of visit over subject
  sigmaResiduals <- as.numeric(HP['sigmaResiduals'])
  sigmaCorrelation <- as.numeric(HP['sigmaCorrelation'])

  ## 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)

  meanLU <- rnorm(J*K,
                     (Intercept + u[,1]) +
                     (Visit + u[,2])*time +
                     Diagnosis*condition +
                     ((Visit + u[,2])*time)*(Diagnosis*condition),sigmaResiduals)

  return(data.frame(CHI_MLU=meanLU,Visit=time,Diagnosis=condition,Child.ID=person))
}