Huntington's disease protector alleles

hd_protector_alleles.r

# Model of control-only Huntington's disease analysis
# If you have some CAG expansion cases who are non-symptomatic then what distinguishes them from non expansion individuals?
# This script shows the analysis is akin to an interaction analysis

library(dplyr)
n <- 1000000
cag <- rpois(n, 2)
mod <- rbinom(n, 2, 0.4) - 1

d <- rep(0, n)
d[cag > 5] <- 1
d[cag > 5] <- rbinom(sum(cag > 5), 1, plogis(mod[cag > 5]))

table(d)

summary(lm(d ~ cag * mod))

tapply(d[cag>5], mod[cag>5], mean)

dat <- tibble(d, cag, mod, rep=as.numeric(cag>5))
summary(glm(rep ~ mod, data=dat, family="binomial", subset=dat$d == 0))