source("http://www.uvm.edu/~rsingle/Rdata/scripts_stat295F14.R")
 fms <- otherdata("FMS_small.tsv", sep="\t")
 attach(fms)

 x <- table(esr1_rs1042717,pre.BMI>25)
 x
 or.AA.GG <- (x[1,1]*x[3,2])/(x[1,2]*x[3,1])
 or.GA.GG <- (x[2,1]*x[3,2])/(x[2,2]*x[3,1])
 or.AA.GG  
 or.GA.GG

#swap the order of the columns (equivalent to the reciprocal of the above ORs)
 x <- table(esr1_rs1042717, !(pre.BMI>25) )
 x
 
 
 or.AA.GG <- (x[1,1]*x[3,2])/(x[1,2]*x[3,1])
 or.GA.GG <- (x[2,1]*x[3,2])/(x[2,2]*x[3,1])
 or.AA.GG  
 or.GA.GG

 chisq.test( table(esr1_rs1042717,pre.BMI>25) ,correct=F )

#-------------------------
 model1 <- glm(pre.BMI>25 ~ esr1_rs1042717, family=binomial)
 model1
 exp(model1$coeff)
 exp(model1$coeff)^-1

 or.AA.GA <- (x[1,1]*x[2,2])/(x[1,2]*x[2,1])
 or.AA.GA

 or.AA.GG <- (x[1,1]*x[3,2])/(x[1,2]*x[3,1])
 or.AA.GG

#-------------------------
 snp <- esr1_rs1042717 #make a copy with a shorter name for convenience
 attributes(snp)
 contrasts(snp)  #default contrast coding (treatment contrasts)
 
 contrasts(snp) <- contr.treatment              #default values (1,2,3) for contrast names 
 attributes(snp)                                #Not the most helpful names for interpretation

   model2 <- glm(pre.BMI>25 ~ snp, family=binomial)
   model2

 contrasts(snp) <- contr.treatment(levels(snp)) #factor levels ("AA","GA","GG") as contrast names
 attributes(snp)

 contrasts(snp) <- contr.treatment(levels(snp), base=3) #set "GG" as baseline level
 attributes(snp)
 
 model2 <- glm(pre.BMI>25 ~ snp, family=binomial)
 model2
 summary(model2)
 exp(model2$coeff)

#95% CI for OR
 m2.sum <- summary(model2) 
 ci.l <- exp( m2.sum$coeff[,1] - 1.96*(m2.sum$coeff[,2]) )
 ci.u <- exp( m2.sum$coeff[,1] + 1.96*(m2.sum$coeff[,2]) )
 rbind(ci.l, ci.u)

#-------------------------
#Recessive model via collapsing the contingency table
 x
 xx <- matrix(NA,2,2)
 xx[1,1] <- x[1,1]
 xx[1,2] <- x[1,2]
 xx[2,1] <- x[2,1] + x[3,1]
 xx[2,2] <- x[2,2] + x[3,2]
 rownames(xx) <- c("AA","OTH")
 xx
 
 (xx[1,1]*xx[2,2])/(xx[1,2]*xx[2,1])


#Dummy coding for Additive (log-odds scale), Dominant, and Recessive models
        #AA GA GG
 Add <- c(2, 1, 0) #Number of A alleles FOR THIS SNP
 Dom <- c(1, 1, 0) #With A as the 'variant' allele
 Rec <- c(1, 0, 0) #With A as the 'variant' allele

 #Additive model on the log-odds scale
  contrasts(snp,1) <- cbind(Add)
  contrasts(snp)
  model3 <- glm(pre.BMI>25 ~ snp, family=binomial)
 #model3
  exp(model3$coeff)

 #Dominant model 
  contrasts(snp,1) <- cbind(Dom)
  contrasts(snp)
  model4 <- glm(pre.BMI>25 ~ snp, family=binomial)
 #model4
  exp(model4$coeff)

 #Recessive model 
  contrasts(snp,1) <- cbind(Rec)
  contrasts(snp)
  model5 <- glm(pre.BMI>25 ~ snp, family=binomial)
 #model5
  exp(model5$coeff)

 #Note that the recessive model agrees best with the 2-parameter model results
  exp(model2$coeff)

#-------------------------
#LRT =-2*(Null deviance - Residual deviance)
#AIC = 2k - 2*log(L)
 
 model2
 with(model2, null.deviance - deviance)
 lr.obs <-  with(model2, null.deviance - deviance)
 lr.obs
 1-pchisq(lr.obs,2)