---
title: "Exercise 1 SISG Association Mapping"
author: "null"
date: "2/6/2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```




```{r echo=FALSE}
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  error = FALSE, 
  collapse = TRUE,
  comment = "#",
  fig.height = 4,
  fig.width = 8,
  fig.align = "center",
  cache = FALSE
)
```

```{r,include=FALSE}
library(dplyr)
library(broom)
library(tidyr)
library(ggplot2)
```


# 1. Load the LHON.txt data file into your R session.  

Can read the data directly from the website if your computer is connected online 

```{r}
LHON=read.table("http://faculty.washington.edu/tathornt/sisg/LHON.txt",header=TRUE)
```

 If file is saved on your computer, could also read the data in from the directory that contains the file.
 
```{r}
LHON2=read.table("LHON.txt",header=TRUE)
```

View the first few lines of the LHON data
```{r}
head(LHON)
```


Get information about the types of variables in the LHON data frame

```{r}
str(LHON)
```




#2. Logistic regression 

First create a 0 and 1 phenotype variable indicating Case/Control Status to perform the logistic regression analysis

```{r}
LHON$newpheno=with(LHON,ifelse(PHENO=="CASE",1,0))
```

What would be the reference genotype for a logistic regression analysis?  Use the levels command in R.  The first factor will be the reference genotype.
```{r}
levels(LHON$GENO)
```
 
## 2a. Perform the logistic regression analysis from session 4 for this data with CC as the reference genotype.  


```{r}
logistmod1=glm(newpheno~GENO,family=binomial(link="logit"),data=LHON)
```

View the summary results of the logistic regression model, including parameter estimates and standard errors  
```{r}
summary(logistmod1)
```


##2b.	Obtain odds ratios and confidence intervals for the CT and TT genotypes.

Can obtain the odds ratio estimates by exponentiating the coefficient estimates from the logistic regression model.  What is the odds ratio for the CT genotype?
```{r}
exp(-1.5994)
```

Can obtain a confidence interval for the  odds ratio parameter for the CT genotype use the standard error of the coefficient estimates from the logistic regression model

```{r}
myse=1.96*(.6378)
CI=c(-1.5994-myse,-1.5994+myse)
exp(CI)
```

Similarly can obtain odds ratio estimates and 95% confidence intervals for genotype TT.

```{r}
exp(-.2654)
myse=1.96*(.5349)
CI=c(-.2654-myse,-.2654+myse)
exp(CI)
```




Alternatively, can obtain the odds ratio estimates and confidence intervale for all paramaters in the logistic regression model by using the *coef()* and *confint.default()* function  

```{r}
exp(coef(logistmod1))
exp(confint.default(logistmod1))
```

Way too many significant digits to report.  Use the *round()* function 
```{r}
round(exp(coef(logistmod1)),2)
round(exp(confint.default(logistmod1)),2)
```



#3. Logistic regression with TT as the reference genotype

Use the relevel function to create a new genotype vector with reference genotype TT

```{r}
LHON$NEWGENO=with(LHON,relevel(GENO, ref = "TT"))
levels(LHON$NEWGENO)
```


Perform the logistic regression analysis TT as the reference genotype.  


```{r}
logistmod2=glm(newpheno~NEWGENO,family=binomial(link="logit"),data=LHON)
```

View the summary results of the logistic regression model, including parameter estimates and standard errors  
```{r}
summary(logistmod2)
```

Obtain the odds ratio estimates and confidence intervale for all paramaters in the logistic regression model 

```{r}
exp(coef(logistmod2))
exp(confint.default(logistmod2))
```


Why are the odds ratios different for CT now?  

The reference genotype group has changed from CC to TT, and the TT genotype group is much larger, which increases the precision of the estimate of the effect.  Can see this more clearly with a stacked barplot

```{r}
plot(factor(PHENO)~factor(GENO), data=LHON)
```

Can also conduct a logistic regression  based on an additive logistic regression model.  First create a genotype variable with an additive coding based on the counts of the number of T alleles 


```{r}
LHON$genoadd <- with(LHON, 0 + 1*(GENO=="CT") + 2*(GENO=="TT"))
```

Now perform the logistic regression analysis with the additive genotype coding
```{r}
logistmod3 <- glm(newpheno~genoadd,family=binomial(link="logit"),data=LHON)
summary(logistmod3)
```
Obtain the odds ratio estimates and confidence intervale for all paramaters in the logistic regression model 

```{r}
round(exp(coef(logistmod3)),3)
round(exp(confint.default(logistmod3)),3)
```

