---
title: "Exploratory Data Analysis"
author: "Statistical Computing, 36-350"
date: "Monday November 1, 2016"
---

Why fit statistical (regression) models?
===

You have some data $X_1,\ldots,X_p,Y$: the variables $X_1,\ldots,X_p$ are called predictors, and $Y$ is called a response. You're interested in the relationship that governs them

So you posit that $Y|X_1,\ldots,X_p \sim P_\theta$, where $\theta$ represents some unknown parameters. This is called **regression model** for $Y$ given $X_1,\ldots,X_p$. Goal is to estimate parameters. Why?

- To assess model validity, predictor importance (**inference**)
- To predict future $Y$'s from future $X_1,\ldots,X_p$'s (**prediction**)

Prostate cancer data 
===

Data set from 97 men who have prostate cancer (from the book [The Elements of Statistical Learning](http://statweb.stanford.edu/~tibs/ElemStatLearn/)). The measured variables:

- `lpsa`: log PSA score
- `lcavol`: log cancer volume
- `lweight`: log prostate weight
- `age`: age of patient
- `lbph`: log of the amount of benign prostatic hyperplasia
- `svi`: seminal vesicle invasion
- `lcp`: log of capsular penetration
- `gleason`: Gleason score 
- ` pgg45`: percent of Gleason scores 4 or 5 

```{r}
pros.df = read.table("http://www.stat.cmu.edu/~ryantibs/statcomp-F16/data/pros.dat")
dim(pros.df)
head(pros.df)
```

Some example questions we might be interested in:

- What is the relationship between `lcavol` and `lweight`?
- What is the relationship between `svi` and `lcavol`, `lweight`?
- Can we predict `lpsa` from the other variables?
- Can we predict whether `lpsa` is high or low, from other variables?

Exploratory data analysis
===

Before pursuing a specific model, it's generally a good idea to look at your data. When done in a structured way, this is called **exploratory data analysis**. E.g., you might investigate:

- What are the distributions of the variables?
- Are there distinct subgroups of samples?
- Are there any noticeable outliers?
- Are there interesting relationship/trends to model?

Distributions of prostate cancer variables
===

```{r}
colnames(pros.df) # These are the variables
par(mfrow=c(3,3), mar=c(4,4,2,0.5)) # Setup grid, margins
for (j in 1:ncol(pros.df)) {
  hist(pros.df[,j], xlab=colnames(pros.df)[j],
       main=paste("Histogram of", colnames(pros.df)[j]),
       col="lightblue", breaks=20)
}
```

What did we learn? A bunch of things! E.g.,

- `svi`, the presence of seminal vesicle invasion, is binary
- `lcp`, the log amount of capsular penetration, is very skewed, a bunch of men with little (or none?), then a big spread; why is this?
- `gleason`, takes integer values of 6 and larger; how does it relate to `pgg45`, the percentage of Gleason scores 4 or 5?
- `lspa`, the log PSA score, is close-ish to normally distributed
  
After asking our doctor friends some questions, we learn:

- When the actual capsular penetration is very small, it can't be properly measured, so it just gets arbitrarily set to 0.25 (and we can check that `min(pros.df$lcp)` $\approx \log{0.25}$)
- The variable `pgg45` measures the percentage of 4 or 5 Gleason scores that were recorded over their visit history *before* their final current Gleason score, stored in `gleason`; a higher Gleason score is worse, so `pgg45` tells us something about the severity of their cancer in the past
  
Correlations between prostate cancer variables
===

```{r}
pros.cor = cor(pros.df)
round(pros.cor,3) 
```

Some strong correlations! Let's find the biggest (in absolute value):

```{r}
pros.cor[lower.tri(pros.cor,diag=TRUE)] = 0 # Why only upper tri part?
pros.cor.sorted = sort(abs(pros.cor),decreasing=T)
pros.cor.sorted[1]
vars.big.cor = arrayInd(which(abs(pros.cor)==pros.cor.sorted[1]), 
                        dim(pros.cor)) # Note: arrayInd() is useful
colnames(pros.df)[vars.big.cor] 
```

This is not surprising, given what we know about `pgg45` and `gleason`; essentially this is saying: if their Gleason score is high now, then they likely had a bad history of Gleason scores

Let's find the second biggest correlation (in absolute value):

```{r}
pros.cor.sorted[2]
vars.big.cor = arrayInd(which(abs(pros.cor)==pros.cor.sorted[2]), 
                        dim(pros.cor))
colnames(pros.df)[vars.big.cor] 
```

This is more interesting! If we wanted to predict `lpsa` from the other variables, then it seems like we should at least include `lcavol` as a predictor

Visualizing relationships among variables, with `pairs()`
===

Can easily look at multiple scatter plots at once, using the `pairs()` function. The first argument is written like a **formula**, with no response variable. We'll hold off on describing more about formulas until we learn `lm()`

```{r}
pairs(~ lpsa + lcavol + lweight + lcp, data=pros.df)
```

Inspecting relationships over a subset of the observations
===

As we've seen, the `lcp` takes a bunch of really low values, that don't appear to have strong relationships with other variables. Let's get rid of them and see what the relationships/correlations look like

```{r}
pros.df.subset = pros.df[pros.df$lcp > min(pros.df$lcp),]
nrow(pros.df.subset) # Beware, we've lost a half of our data! 
pairs(~ lpsa + lcavol + lweight + lcp, data=pros.df.subset)
cor(pros.df.subset[,c("lpsa","lcavol","lweight","lcp")])
```

Testing means between two different groups
===

Recall that `svi`, the presence of seminal vesicle invasion, is binary:

```{r}
table(pros.df$svi)
```

From http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1476128/:

> "When the pathologist’s report following radical pros.dftatectomy describes seminal vesicle invasion (SVI) ... prostate cancer in the areolar connective tissue around the seminal vesicles and outside the prostate ...generally the outlook for the patient is poor."

Does seminal vesicle invasion relate to the volume of cancer? Weight of cancer? Let's do some plotting first:

```{r}
pros.df$svi = factor(pros.df$svi) 
par(mfrow=c(1,2))
plot(pros.df$svi, pros.df$lcavol, main="lcavol versus svi",
     xlab="SVI (0=no, 1=yes)", ylab="Log cancer volume")
plot(pros.df$svi, pros.df$lweight, main="lweight versus svi",
     xlab="SVI (0=no, 1=yes)", ylab="Log cancer weight")
```

Visually, `lcavol` looks like it has a big difference, but `lweight` perhaps does not. Let's try simple two-sample t-tests:

```{r}
t.test(pros.df$lcavol[pros.df$svi==0],
       pros.df$lcavol[pros.df$svi==1])
t.test(pros.df$lweight[pros.df$svi==0],
       pros.df$lweight[pros.df$svi==1])
```

Confirms what we saw visually

Linear model preview
===

Say we want to form a simple linear regression model of the `lpsa` (log PSA score) on `lcavol` (log cancer volume). How? Use the `lm()` function, with the first argument being a regression formula:

```{r}
pros.lm = lm(lpsa ~ lcavol, data=pros.df)
coef(pros.lm) # These are the coefficients
summary(pros.lm) # This is a summary
# Here's a plot of the relationship with the regression line
plot(pros.df$lcavol, pros.df$lpsa, 
     xlab="Log cancer volume", ylab="Log PSA score")
abline(a=coef(pros.lm)[1], b=coef(pros.lm)[2], col="red")
```

To regress `lpsa` onto (say) `lcavol` and `lweight`, we simply change the formula to `lpsa ~ lcavol + lweight`. We'll see much more in the next mini-lecture