---
title: 'Lecture 8: Fitting Data'
author: "Statistical Computing, 36-350"
date: "Wednesday October 7, 2015"
---

In previous episodes
===

- Loading and saving data sets in R format
- Loading and saving data sets in other structured formats
- Intro to regression modeling

Today
===

- Using data frames for statistical purposes
- Basics of regression modeling

So you've got a data frame
===

Recall prostate cancer data set from last time:

```{r}
pros = read.table("http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/prostate.data")
nrow(pros)
head(pros)
```

What can we do with it?

- Plot it: examine multiple variables and distributions
- Test it: compare groups of individuals to each other
- Model it: try to predict, say, one variable from the others

---

Recall the variables that we have measured over 97 men:

- `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 

---

Last time we plotted the distribution of each varaible individually;
and computed the correlations between variables; we learned that

- `lcp`, the log amount of capsular penetration, is very skewed;
  a bunch of men with little (or none?), then a big spread
- `lcp` and `lcavol`, the log cancer volumne, are highly positive 
  correlated
- `gleason`, which records the current Gleason score, and `pgg45`, 
  which records the historical percentage of Gleason scores 4 or 5, 
  are highly positively correlated

Visualizing multiple relationships, with `pairs()`
===

Can easily look at multiple relationships at once, using 
the `pairs()` function:

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

---

Numerical correlations:

```{r}
cor(pros[,c("lcavol","lweight","lcp","lbph")])
```

Inspecting relationships over a subset of the observations
===

We see that `lcp` takes a bunch of really low values, that appear
to have little relationship with other variables.
For exploratory purposes, let's get rid of them and see what
the correlations look like

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

---

Numerical correlations:

```{r}
cor(pros.subset[,c("lcavol","lweight","lcp","lbph")])
```

(Go back to our doctor friend ... ask if this makes sense ...
he tells us that when the recorded capsular penetration value 
is that small, it basically means that it couldn't be accurately 
measured at that scale)

Testing means between two different groups
===

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

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

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

> "When the pathologist’s report following radical prostatectomy 
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 weight of the 
cancer?

Plot it first:

```{r}
pros$svi = factor(pros$svi) # Making it a factor helps with plotting
plot(pros$svi, pros$lweight, main="lweight versus svi",
     xlab="SVI? Yes or no", ylab="Log cancer weight")
```

---

Tough to tell! Let's try a simple two-sample t-test:

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

Not much significance, it appears ...

---

What about the relationship to the volume of cancer?

```{r}
plot(pros$svi, pros$lcavol, main="lcavol versus svi",
     xlab="SVI? Yes or no", ylab="Log cancer volume")
```

---

Try a two-sample t-test again:

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

Wow, looks very significant! Should go report this to 
our doctor friend ...

<small>(Warning: what's wrong with carrying out 5 more
t-tests, for the rest of the variables, and taking the
significant ones?)</small>

Basic regression modeling
===

Let's use `lm()` to fit a basic linear model. Recall that this
is a model of the form
$$
Y = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p +  \mathrm{noise}
$$
for *response* variable $Y$ and *predictor* variables $X_1, \ldots X_p$.
Goal is to estimate the coefficients $\beta_0, \beta_1, \ldots \beta_p$ 

The first argument of `lm()` is a formula, of the form
`Y ~ X1 + X2 + ... Xp`. For example:

```{r}
pros.lm.1 = lm(lpsa ~ lcavol + lweight, data=pros)
class(pros.lm.1)
names(pros.lm.1)
pros.lm.1
```

Getting an overview with `summary()`
==

The function `summary()` gives us a nice summary of the
linear model we just fit:

```{r}
summary(pros.lm.1)
```

This tells us 

- The quantiles of the residuals: ideally, this is a perfect normal dist
- The coefficient estimates
- Their standard errors
- P-values for their individual significance
- (Adjusted) R-squared value: how much variability is explained by the model?
- F-statistic for the overall model significance

Getting rid of the intercept
===

To get rid of the intercept term, we add `+0` to the right-hand
side of the formula:

```{r}
pros.lm.2 = lm(lpsa ~ lcavol + lweight+ 0, data=pros)
summary(pros.lm.2)
```

Is the intercept approriate? In general, depends on the variables that 
make up the linear model 

Getting estimated coefficients with `coef()`
===

How can we access the values of estimated regression coefficients?
Well, we could inspect the names of the fitted object:

```{r}
names(pros.lm.1)
```

Then grab the appropriate component, as in:

```{r}
pros.lm.1$coefficients
```

---

But this is not the best strategy. Instead, let's use the 
`coef()` function:

```{r}
coef(pros.lm.1)
```

Learning to use R's built-in functions like `coef()` to access
certain aspects of fitted statistical models is going to pay off
in the long run, because most object types (beyond those fit by 
`lm()`) will allow for this kind of access

Getting the fitted values with `fitted()`
===

What does our model predict for the log PSA scores of the
97 mean in question? How do these compare to the acutal log PSA
scores?  

To answer the first question, we could extract the coefficients
and then manually multiply and sum the appropriate values:

```{r}
pros.coef = coef(pros.lm.1)
pros.fits.1 = pros.coef[1] + pros.coef[2]*pros$lweight + 
  pros.coef[3]*pros$lcavol
```

This is wrong, oops! We mixed up the order of variables; forgot
that `lcavol` was supposed to come before `lweight`

---

Try again:

```{r}
pros.fits.2 = pros.coef[1] + pros.coef[2]*pros$lcavol + 
  pros.coef[3]*pros$lweight
```

This strategy is laborious and error-prone; it's better to use
`fitted()`, as in:

```{r}
pros.fits.3 = fitted(pros.lm.1)
max(abs(pros.fits.3 - pros.fits.2)) # They are the same
```

For the same reasons as before, it's good practice to use
`fitted()`, a function that works with many common statistical 
objects in R

---

Let's plot the actual values versus the fitted ones

```{r}
plot(pros.fits.3, pros$lpsa, main="Actual versus fitted",
     xlab="Fitted values", ylab="Log PSA values")
abline(a=0, b=1, lty=2, col="red")
```

Appears to be a decent fit

Running diagnostics with `plot()`
===

We can use the `plot()` function to run a series of diagnostic
tests for our regression:

```{r}
plot(pros.lm.1)
```

---

The results are pretty good

- Residuals versus fitted plot: points appear randomly scattered,
  no particular pattern
- Normal Q-Q plot: points are more or less along the 45 degree line,
  so residuals look close to a normal distribution
- Scale-location and residuals versus leverage plots: points appear
  to be groups, no points are too far from the center
  
There is a science (and an art?) to interpreting these; you'll learn
a lot more in 36-401

Making predictions with `predict()`
===

Suppose we had a new observation from a man whose log cancer volume
was 4.1, and log cancer weight was 4.5. What would our model estimate
his log PSA score would be?

Two ways: manual way, and preferred way using `predict()`:

```{r}
pros.pred.1 = pros.coef[1] + pros.coef[2]*4.1 + 
  pros.coef[3]*4.5
pros.new = data.frame(lcavol=4.1, lweight=4.5)
pros.pred.2 = predict(pros.lm.1, newdata=pros.new)
abs(pros.pred.2 - pros.pred.1) # They are the same
```

Building a bigger model
===

We only used two variables in our regression; what happens if we 
use all 8?

```{r}
pros.lm.3 = lm(lpsa ~ lcavol + lweight + age + 
                 lbph + svi + lcp + gleason + pgg45, 
               data=pros)
summary(pros.lm.3)
```

It's unclear whether the addition variables added much, based 
on their p-values

Comparing models based on test data
===

A very robust a sensible way to pit models against each other
is to compare their *test errors*. That is, we hold out some
data, and predict the log PSA values of the unseen observations

That last column in `pros`? Looks like the designers of the data
set already defined *training* and *test* sets for us

```{r}
head(pros)
table(pros$train)
```

---

Let's build our two regression models, using only 67 men---our training data:

```{r}
pros.train = pros[pros$train==TRUE,]
pros.test = pros[pros$train==FALSE,]
pros.lm.train.1 = lm(lpsa ~ lcavol + lweight, data=pros.train)
pros.lm.train.2 = lm(lpsa ~ lcavol + lweight + age + 
                       lbph + svi + lcp + gleason + pgg45, 
                     data=pros.train)
pros.pred.1 = predict(pros.lm.train.1, newdata=pros.test)
pros.pred.2 = predict(pros.lm.train.2, newdata=pros.test)
test.err.1 = mean((pros.test$lpsa - pros.pred.1)^2)
test.err.2 = mean((pros.test$lpsa - pros.pred.2)^2)
test.err.1
test.err.2
```

We may as well go with the simpler model! Aside from being
simpler (which is a plus!) it predicts better, too

Summary
===

- Plotting, exploring, and performing basic statistical tests
  in R is easy to do with data frames
- R has powerful regression modeling tools; using them is easy,
  using them properly requires some training and care (you'll see
  a lot more in 36-401 and 36-402)
