library(lme4) # library needed for mixed modelling
simpvals <- numeric(100) # empty holding vector for P values from repeated simulations

## SIMULATIONS SET 1: effect size

effectsize <- matrix(0, 100, 10) # matrix for storing calculated power estimates - 100 power estimates from each of 10 effect sizes
colnames(effectsize) <- c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10") # 10 effect sizes
for (k in 1:10) { # k = effect size
	print(k)
	for (j in 1:100) {
	  for (i in 1:100) { # 100 power calculations, each based on 100 simulations
	      bCell <- 1 # effect size for differences between cell types
    	  bTreatment <- k # effect size for differences between treatment
    	  b0 <- 0 # model intercept
   		  VBlasto <- 10 # variance within blastocyst - from real data
  	      Verror <- 100 # error variance - from real data
    	  Blasto <- rep(1:8, each = 4) # 8 blastocysts with 4 pseudoreplicated RNAseq samples within each
    	  Cell <- rep(c("A","B","C","D"), 8) # 4 cell types, split equally between blastocysts
    	  Treatment <- rep(c("W","X","Y","Z"), each = 8) # 4 treatments, all samples from same blastocyst = same treatment
    	  Blasto.re <- rnorm(8, 0, sqrt(VBlasto)) # random effects within blastocyst
    	  eps <- rnorm(4*8, 0, sqrt(Verror)) # error variance
    	  # generate simulated data as sum of all effects and variance:
    	  GeneEx <- b0 + bCell*(Cell=="B") + bCell*(Cell=="C") + -bCell*(Cell=="D") +
    	  	bTreatment*(Treatment=="X") + -bTreatment*(Treatment=="Y") + bTreatment*(Treatment=="Z") +
     	    Blasto.re[Blasto] + eps
    	  # form this into a dataframe with cell, treatment and blastocyst factors:
   		  simdata <- data.frame(Blasto = paste("b", Blasto, sep=""), Cell = Cell, Treatment = Treatment,
                          GeneEx = GeneEx)
    	  # mixed linear model analysis:
   		  m1 <- lmer(GeneEx ~ Treatment + Cell + (1|Blasto), data = simdata, REML = F)
   		  m2 <- lmer(GeneEx ~ Cell + (1|Blasto), data = simdata, REML = F)
    	  simpvals[i] <- anova(m1, m2)[2,8]     } # extract p values and store temporarily in holding vector
      effectsize[j,k] <- length(simpvals[simpvals<0.05])/length(simpvals) }  } # calculate power from 100 simulations and store

## SIMULATIONS SET 2: blastocyst sample size

blastosamps <- matrix(0, 100, 5) # matrix for storing calculated power estimates - 100 power estimates from each of 5 sample sizes
colnames(blastosamps) <- c("X48","X96","X144","X192","X240") # 5 sample sizes - note that 8 blastocysts was modelled above
for (k in 1:5) { # k = factor for calculating number of blastocysts
  print(k)
  for (j in 1:100) {
    for (i in 1:100) { # 100 power calculations, each based on 100 simulations
      bCell <- 1 # effect size for differences between cell types
      bTreatment <- 2 # effect size for differences between treatment - fixed at small effect size
      b0 <- 0 # model intercept
      VBlasto <- 10 # variance within blastocyst - from real data
      Verror <- 100 # error variance - from real data
      b <- k*48 # number of blastocysts - needs to be divisible by 8 to create balanced design
      Blasto <- rep(1:b, each = 4) # blastocyst number varies, 4 pseudoreplicated samples within each
      Cell <- rep(c("A","B","C","D"), b) # 4 cell types, split equally between blastocysts
      Treatment <- rep(c("W","X","Y","Z"), each = b) # 4 treatments, all samples from same blastocyst = same treatment
      Blasto.re <- rnorm(b, 0, sqrt(VBlasto)) # random effects within blastocyst
      eps <- rnorm(4*b, 0, sqrt(Verror)) # error variance
      # generate simulated data as sum of all effects and variance:
      GeneEx <- b0 + bCell*(Cell=="B") + bCell*(Cell=="C") + -bCell*(Cell=="D") +
        bTreatment*(Treatment=="X") + -bTreatment*(Treatment=="Y") + bTreatment*(Treatment=="Z") +
        Blasto.re[Blasto] + eps
      # form this into a dataframe with cell, treatment and blastocyst factors:
      simdata <- data.frame(Blasto = paste("b", Blasto, sep=""), Cell = Cell, Treatment = Treatment,
                            GeneEx = GeneEx)
      # mixed linear model analysis
      m1 <- lmer(GeneEx ~ Treatment + Cell + (1|Blasto), data = simdata, REML = F)
      m2 <- lmer(GeneEx ~ Cell + (1|Blasto), data = simdata, REML = F)
      simpvals[i] <- anova(m1, m2)[2,8]     } # extract p values and store temporarily in holding vector
    blastosamps[j,k] <- length(simpvals[simpvals<0.05])/length(simpvals) }  } # calculate power from 100 simulations and store
    
## EXAMPLE PLOT CODE    
    
CI <- function(data){c(mean(data)-sd(data)*1.96, mean(data)+sd(data)*1.96)}
plot(apply(effectsize, 2, mean), pch=16, xlab="Effect size", ylab="Power", las=1, ylim=c(0,1))
for (i in 1:10) lines(c(i,i), c(CI(effectsize[,i])[1],CI(effectsize[,i])[2]), lwd=2)