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)