PowerAnalysis_final.R

library(lme4) library(dplyr) #This script has been adapted from https://osf.io/5476a/ and https://github.com/RLadiesNijmegen/Power/blob/master/RLadies_SimplePower.R

powerFuncLMB <-function(s, i){ si=s*i c=4 ic=i/c cr=s/c g <- s*i/4

ss <- rep(rep(1:s,each=ic),c)

#responses by condition group1 <- rnorm(g,m=679,sd=183) group2 <- rnorm(g,625,183) group3 <- rnorm(g,625,183) group4 <- rnorm(g,597,169)

sran <- cbind(1:s, rnorm(s, 0, 95.17)) colnames(sran) <- c("ss", "sran") iran <- cbind(1:i, rnorm(i, 0, 32.48)) colnames(iran) <- c("ii", "iran")

# for latin square within items design: repeat subject numbers cr times per condition. rotate which people get which items ss <- rep(rep(1:s,each=ic),c) ii1 <- rep(1:i, cr) ii2 <- rep(c((ic+1):i, 1:ic),cr) ii3 <- rep(c((2*ic+1):i, 1:(2*ic)),cr) ii4 <- rep(c((3*ic+1):i, 1:(3*ic)),cr) ii <- c(ii1, ii2, ii3, ii4)

#build dataframe Resp <- c(group1,group2,group3,group4) class(Resp) Condition <- c(rep("Unrelated", si/c), rep("LongDelay", si/c), rep("ShortDelay", si/c), rep("NoDelay", si/c)) ds <- cbind(Resp, Condition, ss, ii) ds <- as.data.frame(ds) ds <- ds[order(ss),]

#Helmert contrast: unrelated vs everything else ds$uvall <- -1/4 ds[ds$Condition == 'Unrelated',]$uvall <- 3/4

#Helmert contrast: LongDelay vs ShortDelay and NoDelay ds$lvsn <- -1/3 ds[ds$Condition == 'Unrelated',]$lvsn <- 0 ds[ds$Condition == 'LongDelay',]$lvsn <- 2/3

#Helmert contrast: ShortDelay vs NoDelay ds$svn <- -1/2 ds[ds$Condition == 'Unrelated',]$svn <- 0 ds[ds$Condition == 'LongDelay',]$svn <- 0 ds[ds$Condition == 'ShortDelay',]$svn <- 1/2

summary(ds) ds <- merge(ds, sran) ds <- merge(ds, iran) ds$Resp <- as.numeric(as.character(ds$Resp)) + ds$sran + ds$iran

#run models fit1 <- lmer(Resp ~ 1+ (uvall+lvsn+svn) + (1|ss) + (1|ii),data=ds, REML=F) summary(fit1) fit2 <- lmer(Resp ~ 1+ (uvall+lvsn+svn) - uvall + (1|ss) + (1|ii),data=ds, REML=F) fit3 <- lmer(Resp ~ 1+ (uvall+lvsn+svn) - lvsn + (1|ss) + (1|ii),data=ds, REML=F) fit4 <- lmer(Resp ~ 1+ (uvall+lvsn+svn) - svn + (1|ss) + (1|ii),data=ds, REML=F)

c(round(anova(fit2,fit1)[2,8],3),round(anova(fit3,fit1)[2,8],3), round(anova(fit4,fit1)[2,8],3)) }

out <- replicate(1000,powerFuncLMB(s=40,i=160))

power<- mean(out<0.05)