Genotype_LogisticRegression

# create the control data
IID <- paste("CTRL",1:20000,sep="")
ctrl <- as.data.frame(IID)
ctrl$GENO <- "AA"
ctrl[19209:19992,2] <- "AC"
ctrl[19993:20000,2] <- "CC"
ctrl$GENO <- as.factor(ctrl$GENO)
ctrl$PHENO <- "CONTROL"

# create the RA data
IID <- paste("RA",1:5500,sep="")
ra <- as.data.frame(IID)
ra$GENO <- "AA"
ra[5070:5491,2] <- "AC"
ra[5492:5500,2] <- "CC"
ra$GENO <- as.factor(ra$GENO)
ra$PHENO <- "CASE"

# merge
dat <- rbind(ctrl,ra)

# classify pheno
dat$newpheno=with(dat,ifelse(PHENO=="CASE",1,0))

# logistic regression (not additive) it calculates the OR separately for Het and Hom
model1=glm(newpheno~GENO,family=binomial(link="logit"),data=dat)
summary(model1)

# plot the data to get a better idea
plot(factor(PHENO)~factor(GENO), data=dat)

# now an additive model - more reasonable for complex disease genetics
dat$genoadd <- with(dat, 0 + 1*(GENO=="AC") + 2*(GENO=="CC"))
model3 <- glm(newpheno~genoadd,family=binomial(link="logit"),data=dat)
summary(model3)

# odds ratio
exp(model3$coefficients[2])

Was adapted from this example
http://faculty.washington.edu/tathornt/SISG2019/exercises/assoc_ex01.R