03_08_comparing-models.Rmd

---
title: Comparing models example
author: Jeff Leek
output:
  rmarkdown::html_document:
    toc: true
  vignette: >  
    %\VignetteIndexEntry{Comparing models example}
    %\VignetteEngine{knitr::rmarkdown}
    \usepackage[utf8]{inputenc}
---

```{r front, child="./../front.Rmd", echo=FALSE}
```


## Dependencies

This document depends on the following packages:

```{r load_hidden, echo=FALSE, results="hide", warning=FALSE}
suppressPackageStartupMessages({
  library(devtools)
  library(Biobase)
  library(sva)
  library(bladderbatch)
})
```

```{r load}
  library(devtools)
  library(Biobase)
  library(sva)
  library(bladderbatch)
```


To install these packages you can use the code (or if you are compiling the document, remove the `eval=FALSE` from the chunk.)

```{r install_packages, eval=FALSE}
install.packages(c("devtools"))
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("Biobase","sva","bladderbatch"))
```

## Download the data


The analyses performed in this experiment are based on gene expression measurements from a bladder cancer study: [Gene expression in the urinary bladder: a common carcinoma in situ gene expression signature exists disregarding histopathological classification.](http://cancerres.aacrjournals.org/content/64/11/4040.full.pdf) The data can be loaded from the [bladderbatch](http://bioconductor.org/packages/release/data/experiment/html/bladderbatch.html) data package. 

```{r}
library(sva)
library(bladderbatch)
data(bladderdata)
```

## Set up the data

```{r}
pheno = pData(bladderEset)
edata = exprs(bladderEset)
```

## Plot data with two different model fits

```{r}
## Setting seed so the jitter will be the same
set.seed(123)

cancerjit = jitter(as.numeric(pheno$cancer))
lm1 = lm(edata[1,] ~ 1)
lm2 = lm(edata[1,] ~ pheno$cancer)

par(mfrow=c(1,2))

plot(edata[1,] ~ cancerjit,type="p",xaxt="n",xlab="Cancer Status",ylab="Expression",pch=19,col=as.numeric(pheno$cancer))
axis(1,at=c(1,2,3),c("Biopsy","Cancer","Normal"))
abline(lm1,col="darkgrey",lwd=5)


plot(edata[1,] ~ cancerjit,type="p",xaxt="n",xlab="Cancer Status",ylab="Expression",pch=19,col=as.numeric(pheno$cancer))
axis(1,at=c(1,2,3),c("Biopsy","Cancer","Normal"))
boxplot(lm2$fitted~pheno$cancer,add=T,border=1:3)

```


## Plot the residuals

```{r}
par(mfrow=c(1,2))
plot(lm1$residuals ~ cancerjit,type="p",xaxt="n",xlab="Cancer Status",ylab="Expression",pch=19,col=as.numeric(pheno$cancer),ylim=c(-1.1,1.1))
axis(1,at=c(1,2,3),c("Biopsy","Cancer","Normal"))

plot(lm2$residuals ~ cancerjit,type="p",xaxt="n",xlab="Cancer Status",ylab="Residuals",pch=19,col=as.numeric(pheno$cancer),ylim=c(-1.1,1.1))
axis(1,at=c(1,2,3),c("Biopsy","Cancer","Normal"))
```


## Session information

Here is the session information 

```{r session_info}
devtools::session_info()
```

It is also useful to compile the time the document was processed. This document was processed on: `r Sys.Date()`.