power_interactions

power1.pdf

power2.pdf

power_interactions.r

library(dplyr)
library(parallel)
library(ggplot2)

simulate_data <- function(n, prop, target_var, b_int)
{
	x1 <- rnorm(n)
	x2 <- rbinom(n, 1, prop)

	score <- x1 + x2 + x1 * x2 * b_int
	y <- score + rnorm(n, 0, sqrt(target_var-var(score)))
	coefficients(summary(lm(y ~ x1 + x2 + x1 * x2)))[4,4]
}

conv_b_to_rsq <- function(prop, vary, b_int, varx)
{
	b_int^2 * prop * (1-prop) * varx / vary
}

interpolate_power <- function(x, y, target)
{
	approx(x, y, target)$y
}

# Set parameters to calculate
params <- expand.grid(
	n = 80000,
	prop = seq(0.1, 0.5,  by=0.05),
	target_var = 10,
	b_int = seq(0.01, 0.14, by=0.01),
	reps = 1:1000,
	p = NA
) %>% as_tibble()
params$count <- 1:nrow(params)

# Run simulations
# takes about 5 mintues with 16 cores
o <- mclapply(split(params, seq(nrow(params))), function(x) {
	message(x$count)
	x$p <- simulate_data(x$n, x$prop, x$target_var, x$b_int)
	return(x)
}, mc.cores=16) %>% bind_rows

# Convert interaction term to difference in slopes of sd / sd
o$b_int <- o$b_int / sqrt(target_var)

# Summarise
res <- o %>% 
	group_by(n, prop, b_int) %>%
	summarise(power = sum(p < 0.01) / n(), nsim = n())

# plot power curves
p1 <- ggplot(res, aes(y=power, x=prop, group=as.factor(b_int))) +
geom_point(aes(colour=as.factor(b_int))) +
geom_line(aes(colour=as.factor(b_int))) +
labs(x="proportion low mat edu", colour="interaction\neffect size")
ggsave(p1, file="power1.pdf") 


# interpolate to identify the b_int value for each prop value that achieves 80% power
const_power <- res %>% group_by(prop) %>%
	summarise(b_int = interpolate_power(power, b_int, 0.8))

# plot interpolation
p2 <- ggplot(const_power, aes(x=prop, y=b_int)) +
geom_point() +
geom_line() +
labs(x="Proportion of sample with low maternal education", y="Difference in slopes detectable (SD change in mental health per SD change in air pollution)\nbased on maternal education given 80% power and 0.01 p-value")
ggsave(p2, file="power2.pdf")

# save simulation results
save(o, res, const_power, file="powersims.rdata")

# calculate for prop=0.25
subset(const_power, prop==0.25)

powersims.rdata