CWY2016_power_analysis.R

if(!require(simr)){install.packages("simr"); library(simr)}else{library(simr)}
#library(lme4)

## Means and Variance of Behavior data
CWL2016RT <- rbind(M=c(.793,.815,.893,.866),VAR=c(.011^2,.010^2,.016^2,.015^2))

## Got the parameters of RT distributions
## Referring to https://stats.stackexchange.com/questions/12232/calculating-the-parameters-of-a-beta-distribution-using-the-mean-and-variance

estBetaParams <-function(mu, var) {
        alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
        beta <- alpha * (1 / mu - 1)
        return(params = list(alpha = alpha, beta = beta))
}

## Setup parameters for RT, consistency, neighborhood

sim_RT_params <- estBetaParams(CWL2016RT["M",],CWL2016RT["VAR",])

sim_RT <- NULL
for(i in 1:2){
        for(j in 1:4){
                sim_RT <- c(sim_RT, 1000*rbeta(1, sim_RT_params$alpha[j], sim_RT_params$beta[j]))
        }
}

sim_CONSISTENCY <- rep(rep(c(1,2),each=2),2)
sim_NEIGHBORHOOD <- rep(rep(c(1,2),2),2)

sim_SUBJ <- rep(c("s1","s2"), each = 4)
sim_ITEM <- rep(paste0("i",1:4),2)

## Simulated average RT of two particiapnts
sim_df <- data.frame(sim_RT, sim_CONSISTENCY, sim_NEIGHBORHOOD, sim_SUBJ, sim_ITEM)


## Model of all main effects and interaction
sim_alleffects.lmer = lmer( sim_RT ~ sim_CONSISTENCY*sim_NEIGHBORHOOD + (1|sim_SUBJ) + (1|sim_ITEM), data=sim_df )

## Model without interaction
sim_nointeraction.lmer = lmer( sim_RT ~ sim_CONSISTENCY + sim_NEIGHBORHOOD + (1|sim_SUBJ) + (1|sim_ITEM), data=sim_df )

## Effect size of simulated data
fixef(sim_alleffects.lmer)
fixef(sim_nointeraction.lmer)

## Change the effect size of factor and interaction close to the original study
fixef(sim_alleffects.lmer)["sim_CONSISTENCY"] <- 75
fixef(sim_alleffects.lmer)["sim_NEIGHBORHOOD"] <- 27
fixef(sim_alleffects.lmer)["sim_CONSISTENCY:sim_NEIGHBORHOOD"] <- -48

fixef(sim_nointeraction.lmer)["sim_CONSISTENCY"] <- 75
fixef(sim_nointeraction.lmer)["sim_NEIGHBORHOOD"] <- 27

set.seed(523)
## Check the initial powers
powerSim(sim_alleffects.lmer, test=fixed("sim_CONSISTENCY"), nsim=10)
powerSim(sim_alleffects.lmer, test=fixed("sim_NEIGHBORHOOD"), nsim=10)
powerSim(sim_alleffects.lmer, test=fixed("sim_CONSISTENCY:sim_NEIGHBORHOOD"), nsim=10)

## Extend the model to 120 items
sim_alleffects.item = extend(sim_alleffects.lmer, along = "sim_ITEM", n = 120)

sim_nointeraction.item = extend(sim_nointeraction.lmer, along = "sim_ITEM", n = 120)

## Check the updated powers
powerSim(sim_alleffects.item, test=fixed("sim_CONSISTENCY"), nsim=10)
powerSim(sim_alleffects.item, test=fixed("sim_NEIGHBORHOOD"), nsim=10)
powerSim(sim_alleffects.item, test=fixed("sim_CONSISTENCY:sim_NEIGHBORHOOD"), nsim=10)

## Extend the model to 27 subjects. As equal as the original study
sim_alleffects.subj = extend(sim_alleffects.item, along = "sim_SUBJ", n = 27)

sim_nointeraction.subj = extend(sim_nointeraction.item, along = "sim_SUBJ", n = 27)

## Check the significance of the interaction
anova(sim_alleffects.subj, sim_nointeraction.subj)

## Check the updated powers
sim_alleffects.CONSISTENCY.pwr <- powerSim(sim_alleffects.subj, test=fixed("sim_CONSISTENCY", method = "t"), nsim=1000)
sim_alleffects.NEIGHBORHOOD.pwr <- powerSim(sim_alleffects.subj, test=fixed("sim_NEIGHBORHOOD", method = "t"), nsim=1000)
sim_alleffects.INTERACTION.pwr <- powerSim(sim_alleffects.subj, test=fixed("sim_CONSISTENCY:sim_NEIGHBORHOOD", method = "t"), nsim=1000)

## Power Curve of Interaction
sim_alleffects.samples = extend(sim_alleffects.item, along = "sim_SUBJ", n = 100)
pc.interaction.samples <- powerCurve(sim_alleffects.samples, test=fixed("sim_CONSISTENCY:sim_NEIGHBORHOOD", method = "t"), along = "sim_SUBJ", nsim=100, breaks = c(20,40,60,80,100))
pc.consistency.samples <- powerCurve(sim_alleffects.samples, test=fixed("sim_CONSISTENCY", method = "t"), along = "sim_SUBJ", nsim=100, breaks = c(20,40,60,80,100))