Skip to content

Instantly share code, notes, and snippets.

@tomcarpenter

tomcarpenter/example.R

Last active July 6, 2020 04:11
Show Gist options
  • Save tomcarpenter/6dab87e91a083e7f4c439501965961b7 to your computer and use it in GitHub Desktop.
Save tomcarpenter/6dab87e91a083e7f4c439501965961b7 to your computer and use it in GitHub Desktop.
Download ZIP
Code for graph
Raw
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment