---
title: "Linear Models"
author: "Statistical Computing, 36-350"
date: "Wednesday November 3, 2016"
---

Linear regression modeling
===

The linear model is arguably the **most widely used** statistical model, has a place in nearly every application domain of statistics

Given response $Y$ and predictors $X_1,\ldots,X_p$, in a **linear regression model**, we posit:

$$
Y = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p + \epsilon, 
\quad \text{where $\epsilon \sim N(0,\sigma^2)$}
$$

Goal is to estimate parameters, also called coefficients $\beta_0,\beta_1,\ldots,\beta_p$ 

Prostate cancer data, revisited
===

Recall, data set from 97 men who have prostate cancer (from the book [The Elements of Statistical Learning](http://statweb.stanford.edu/~tibs/ElemStatLearn/)):
 
```{r}
pros.df = read.table("http://www.stat.cmu.edu/~ryantibs/statcomp-F16/data/pros.dat")
dim(pros.df)
head(pros.df)
```

Fitting a linear regression model with `lm()`
===

We can use `lm()` to fit a linear regression model. The first argument is a **formula**, of the form `Y ~ X1 + X2 + ... + Xp`, where `Y` is the response and `X1`, ..., `Xp` are the predictors. These refer to column names of variables in a data frame, that we pass in through the `data` argument

E.g., for the prostate data, to regress the response variable `lpsa` onto the predictors variables `lcavol` and `lweight`:

```{r}
pros.lm = lm(lpsa ~ lcavol + lweight, data=pros.df)
class(pros.lm) # Really, a specialized list
names(pros.lm) # Here are its components
pros.lm # It has a special way of printing
```

Utility functions
===

Linear models in R come with a bunch of **utility** functions, such as `coef()`, `fitted()`, `residuals()`, `summary()`, `plot()`, `predict()`, for retrieving coefficients, fitted values, residuals, producing summaries, producing diagnostic plots, making predictions, respectively

These tasks can also be done manually, by extracting at the components of the returned object from `lm()`, and manipulating them appropriately. But this is **discouraged**, because:

- The manual strategy is more tedious and error prone
- Once you master the utility functions, you'll be able to retrieve coefficients, fitted values, make predictions, etc., in the same way for model objects returned by `glm()`, `gam()`, and many others

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

So, what were the regression coefficients that we estimated? Use the `coef()` function, to retrieve them:

```{r}
pros.coef = coef(pros.lm) # Vector of 3 coefficients
pros.coef
```

What does a linear regression coefficient mean, i.e., how do you interpret it? Note, from our linear model:

$$
\mathbb{E}(Y|X_1,\ldots,X_p) = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p
$$

So, increasing predictor $X_j$ by one unit, while holding all other predictors fixed, increases the expected response by $\beta_j$

E.g., increasing `lcavol` (log cancer volume) by one unit (one cc), while holding `lweight` (log cancer weight) fixed, increases the expected value of `lpsa` (log PSA score) by $\approx 0.65$

Retrieving fitted values with `fitted()`
===

What does our model predict for the log PSA scores of the 97 mean in question? And how do these compare to the actual log PSA scores? Use the `fitted()` function, then plot the actual values versus the fitted ones:

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

Displaying an overview with `summary()`
==

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

```{r}
summary(pros.lm) # Special way of summarizing
```

This tells us:

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

Running diagnostics with `plot()`
===

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

```{r}
plot(pros.lm) # Special way of plotting
```

The results are pretty good:

- "Residuals vs Fitted" plot: points appear randomly scattered, no particular pattern
- "Normal Q-Q" plot: points are mostly along the 45 degree line, so residuals look close to a normal distribution
- "Scale-Location" and "Residuals vs Leverage" plots: no obvious 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 the Modern Regression 36-401 course

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 linear model estimate his log PSA score would be? Use `predict()`:

```{r}
pros.new = data.frame(lcavol=4.1, lweight=4.5) # Must set up a new data frame
pros.pred = predict(pros.lm, newdata=pros.new) # Now call predict with new df
pros.pred
```

We'll learn much more about making/evaluating predictions in our week of mini-lectures on "Statistical Prediction"

Some handy shortcuts
===

Here are some handy shortcuts, for fitting linear models with `lm()` (there are also many others):

- No intercept (no $\beta_0$ in the mathematical model): use `0 +` on the right-hand side of the formula, as in:
    ```{r}
    summary(lm(lpsa ~ 0 + lcavol + lweight, data=pros.df))
    ```
- Include all predictors: use `.` on the right-hand side of the formula, as in:
    ```{r}
    summary(lm(lpsa ~ ., data=pros.df))
    ```
- Include interaction terms: use `:` between two predictors of interest, to include the interaction between them as a predictor, as in:
    ```{r}
    summary(lm(lpsa ~ lcavol + lweight + lcavol:svi, data=pros.df))
    ```