Comparing interaction test against variance heterogeneity test in gene x smoking interaction example

interaction_vs_varhet.R

# Initial parameters
sample_size <- 120000
bmi_mean <- 25
bmi_sd <- 4
allele_freq <- 0.19
effect_size <- -0.23
proportion_smokers <- 0.5


# The transformation of BMI to z^2 allows detection of variance heterogeneity, e.g.:
bmi <- rnorm(n=sample_size, mean=bmi_mean, sd=bmi_sd)
bmi[snp==0] <- rnorm(n=sum(snp==0), mean=bmi_mean, sd=1)
bmi[snp==1] <- rnorm(n=sum(snp==1), mean=bmi_mean, sd=2)
bmi[snp==2] <- rnorm(n=sum(snp==2), mean=bmi_mean, sd=3)
bmi_z2 <- ((bmi - mean(bmi)) / sd(bmi))^2

summary(lm(bmi ~ snp)) # There is no main effect so this is not significant
summary(lm(bmi_z2 ~ snp)) # There is a big var het effect so this is highly significant


# Simulate data such that BMI is influenced by a gene x smoking interaction
# Test for interaction
# Perform variance heterogeneity test
simulateData <- function(sample_size, bmi_mean, bmi_sd, allele_freq, effect_size, proportion_smokers)
{
	smoker <- sample(c(0,1), sample_size, replace=T, prob=c(1 - proportion_smokers, proportion_smokers))
	snp <- rbinom(sample_size, 2, allele_freq)

	bmi <- rnorm(n=sample_size, mean=bmi_mean, sd=bmi_sd)
	index <- smoker == 1
	bmi[index] <- bmi[index] + snp[index] * effect_size

	p_int <- anova(lm(bmi ~ snp : smoker))$P[1]

	bmi_z2 <- ((bmi - mean(bmi)) / sd(bmi))^2
	p_vhet <- anova(lm(bmi_z2 ~ snp))$P[1]
	return(c(p_int, p_vhet))
}


# Perform 20 simulations 

nsim <- 20
results <- array(0, c(nsim, 2))
for(i in 1:nsim)
{
	cat(i, "\n")
	results[i,] <- simulateData(sample_size, bmi_mean, bmi_sd, allele_freq, effect_size, proportion_smokers)
}

results <- data.frame(Interaction = results[,1], VarHet = results[,2])
boxplot(-log10(results), ylab=expression(-log[10]*p))
abline(h=-log10(0.05))