Code for graph

example.R

library(mvtnorm)
library(TruncatedNormal)
library(rockchalk)
library(tidyverse)


# Basic information about IATs
# IAT ~ mean = 0, sd = .45, trunc @ +/-2
# preference ~ mean @ 4, SD =?, trun @ 1-7

#for hedonic, cor should be around .50, 
#for util, let's go with upper estimate of .45

# first make data generating function; r1 is the larger one
n=1000; r1=.5; r2=.45
gendat <- function(n, r1, r2){
  require(TruncatedNormal); require(mvtnorm); require(rockchalk); require(tidyverse)
  s1 <- rockchalk::lazyCor(r1, 2)  
  s2 <- rockchalk::lazyCor(r2, 2)  
  dat1 <- TruncatedNormal::rtmvnorm(n, mu=c(0, 4), sigma=s1, lb=c(-2, 1), ub=c(2, 7)) %>% as.data.frame()
  dat2 <- TruncatedNormal::rtmvnorm(n, mu=c(0, 4), sigma=s2, lb=c(-2, 1), ub=c(2, 7)) %>% as.data.frame()
  dat1$cond <- "hedon"
  dat2$cond <- "util"
  dat <- bind_rows(dat1, dat2)
  names(dat) <- c("D", "xp", "cond")
  dat$xp <- round(dat$xp)
  dat$cond <- factor(dat$cond)
  #scamble 
  dat <- dat[sample(1:nrow(dat), replace=FALSE),]
  return(dat)
}

#try it out once
gendat(1000, r1=.5, r2=.45)

#function runs it many times with given values, then returns p-value for interaction
simpower <- function(n, r1, r2){
  #generate the data
  simdata <- replicate(1000, expr=gendat(n, r1=r1, r2=r2), simplify=FALSE)
  
  #run the models
  simmodels <- sapply(simdata, function(x){
    mod <- lm(D ~ xp*cond, data=x)
    p <- summary(mod)$coef[4,4]
  })
  print(paste0("Completed for n = ", n))
  return(data.frame(p=simmodels, n=n))
}


#run function over many values for n; collapse list into data frame
out <- lapply(seq(1000, 10000, by=1000), function(x){simpower(x, r1=.5, r2=.40)}) %>% do.call(rbind, .) %>% data.frame(.)

# get power by sample size, graph
out %>% group_by(n) %>% summarise(power = sum(p < .05)/n()) %>% 
  ggplot(data=., aes(x=n, y=power))+
  geom_point()+
  geom_smooth(se=FALSE)+
  geom_line()+
  geom_hline(yintercept=.80, color="red")+
  scale_x_continuous(breaks=seq(1000, 10000, 250))+
  scale_y_continuous(breaks=seq(0,1,.05))+
  theme_light()+
  theme(axis.text.x = element_text(angle = 90))+
  ggtitle("Power for Interaction", subtitle="n = n per group (two groups total)" )
  
#run that analysis over possible lower bound correlations and try again
rbind(
  lapply(seq(1000, 10000, by=1000), function(x){simpower(x, r1=.5, r2=.40)}) %>% do.call(rbind, .) %>% data.frame(., r2=.40),
  lapply(seq(1000, 10000, by=1000), function(x){simpower(x, r1=.5, r2=.41)}) %>% do.call(rbind, .) %>% data.frame(., r2=.41),
  lapply(seq(1000, 10000, by=1000), function(x){simpower(x, r1=.5, r2=.42)}) %>% do.call(rbind, .) %>% data.frame(., r2=.42),
  lapply(seq(1000, 10000, by=1000), function(x){simpower(x, r1=.5, r2=.43)}) %>% do.call(rbind, .) %>% data.frame(., r2=.43),
  lapply(seq(1000, 10000, by=1000), function(x){simpower(x, r1=.5, r2=.44)}) %>% do.call(rbind, .) %>% data.frame(., r2=.44),
  lapply(seq(1000, 10000, by=1000), function(x){simpower(x, r1=.5, r2=.45)}) %>% do.call(rbind, .) %>% data.frame(., r2=.45)
) -> out
out$r2 <- factor(out$r2)

# get power by sample size and r2, graph
out %>% group_by(n, r2) %>% summarise(power = sum(p < .05)/n()) %>% 
  ggplot(data=., aes(x=n, y=power, color=r2))+
  geom_point()+
  geom_line()+
  geom_hline(yintercept=.80, color="red")+
  scale_x_continuous(breaks=seq(1000, 10000, 250))+
  scale_y_continuous(breaks=seq(0,1,.05))+
  theme_light()+
  theme(axis.text.x = element_text(angle = 90))+
  ggtitle("Power for Interaction: r1 = .50, r2 = ....", subtitle="n = n per group (two groups total)")
# saveRDS(out, "Power.rds")