## Building (Vlad) Valid models

![vlad valid](http://fc06.deviantart.net/fs70/i/2013/337/0/8/vlad_valid_model_sheet_by_gabbyvee-d6wnzyi.jpg)


## Building valid models {.build}

A **valid model** is one that is a good fit to the data; i.e. where the 
assumptions made by the model are reasonable for the data set.

This not the same thing as:

- A model with high explanatory power (e.g. high $R^2$)
- A statistically significant model.

Both of these only matter if we're dealing with a valid model.


## Building valid models
Today we'll talk about:

1. Leverage
2. Standardized residuals
3. Influence


# Leverage and Influence


## Recall... the outlier

```{r, echo=FALSE, warning=FALSE, message=FALSE}
# this chunk sets the chunk options for the whole document
require(knitr)
opts_chunk$set(warning=FALSE, message=FALSE)
```

```{r, echo=FALSE}
library(openintro)
COL <- c('#55000088','#225588')
set.seed(238)
n <- c(50, 25, 78, 55, 70, 150)
m <- c(12, -4, 7, -19, 0, 40)
xr <- list(0.3, c(2), 1.42, runif(4,1.45,1.55), 5.78, -0.6)
yr <- list(-4, c(-8), 19, c(-17,-20,-21,-19), 12, -23.2)
i <- 1
x <- runif(n[i])
y <- m[i]*x + rnorm(n[i])
x <- c(x,xr[[i]])
y <- c(y,yr[[i]])
par(mar=c(4,4,1,1), las=1, mgp=c(2.5,0.5,0), cex.lab = 1.25, cex.axis = 1.25, mfrow = c(2,1))
lmPlot(x, y, col = COL[2], lCol = COL[1], lwd = 3)
```


## ...the leverage

```{r, echo=FALSE}
i <- 2
x <- runif(n[i])
y <- m[i]*x + rnorm(n[i])
x <- c(x,xr[[i]])
y <- c(y,yr[[i]])
par(mar=c(4,4,1,1), las=1, mgp=c(2.5,0.5,0), cex.lab = 1.25, cex.axis = 1.25, mfrow = c(2,1))
lmPlot(x, y, col = COL[2], lCol = COL[1], lwd = 3)
```


## ...the influence

```{r, echo=FALSE}
i <- 3
x <- runif(n[i])
y <- m[i]*x + rnorm(n[i])
x <- c(x,xr[[i]])
y <- c(y,yr[[i]])
par(mar=c(4,4,1,1), las=1, mgp=c(2.5,0.5,0), cex.lab = 1.25, cex.axis = 1.25, mfrow = c(2,1))
lmPlot(x, y, col = COL[2], lCol = COL[1], lwd = 3)
```


## ...the influence

```{r, echo=FALSE}
i <- 5
x <- runif(n[i])
y <- m[i]*x + rnorm(n[i])
x <- c(x,xr[[i]])
y <- c(y,yr[[i]])
par(mar=c(4,4,1,1), las=1, mgp=c(2.5,0.5,0), cex.lab = 1.25, cex.axis = 1.25, mfrow = c(2,1))
lmPlot(x, y, col = COL[2], lCol = COL[1], lwd = 3)
```


## ...the influence

```{r, echo=FALSE}
par(mar=c(4,4,1,1), las=1, mgp=c(2.5,0.5,0), cex.lab = 1.25, cex.axis = 1.25, mfrow = c(2,1))
lmPlot(x[1:70], y[1:70], col = COL[2], lCol = COL[1], lwd = 3, xlim = range(x), ylim = range(y))
```


## Outliers, leverage, influence {.build}

**Outliers** are points that don't fit the trend in the rest of the data.

**High leverage points** have the potential to have an unusually large influence 
on the fitted model.

**Influential points** are high leverage points that cause a very different
line to be fit than would be with that point removed.


## Quantifying leverage: $h_{ii}$ {.build}

We need a metric for the leverage of $x_i$ that incorporates

1. The distance $x_i$ is away from the bulk of the $x$'s.
2. The extent to which the fitted regression line is attracted by the given point.

For simple regression:

\[ h_{ii} = \frac{1}{n} + \frac{(x_i - \bar{x})^2}{\sum_{j = 1}^n(x_j - \bar{x})^2} \]

A general solution for multiple regression: calculate the hat matrix, pull out the
diagonal elements.


# Standardized residuals


## Why standardize?

We will need to assess whether the assumption of errors with constant variance 
is reasonable.

What if the non-constant variance that we observe isn't a property of the *errors*
but the non-constant variance that goes into making the prediction from which the
*residuals* are calculated?


## Properties of the residuals

We can write the residuals as

\[ \hat{e} = (I - H) Y \]

and calcuate the expected value and variance as

\[ E(\hat{e} | X) = 0 \]

\[ Var(\hat{e} | X) = \sigma^2 (I - H) \]

and for a particular $\hat{e}_i$:

\[ Var(\hat{e}_i) = \sigma^2 (I - h_{ii}) \]


## Standardized residuals {.build}

To be sure all of our residuals are assessed on equal footing, we divide each one
by our estimate of it's standard deviation.

\[ r_i = \frac{\hat{e}}{s \sqrt{1 - h_{ii}}} \]

Where $s$ is our usual estimate of $\sigma$.

Observations with high standardized residuals can be considered *outliers*. Rule
of thumb: $|r_i| > 2$ for small data, $|r_i| > 4$ for large.


## We've already been using these

```{r, echo=FALSE,fig.align='center', fig.height=4.5, fig.width=4.5}
LA <- read.csv("http://andrewpbray.github.io/data/LA.csv")
m1 <- lm(log(price) ~ log(sqft) + bed + city, data = LA)
plot(m1, 2)
```


# Influence


## Cook's Distance {.build}

An alternate form:

\[ D_i = \frac{r_i^2}{p + 1} \frac{h_{ii}}{1 - h_{ii}} \]

To be influential a point must:

1. Have high leverage $h_{ii}$ *and*
2. Have a high standardized residual $r_i$


# Activity 

## Activity 9 (Part I)

```{r echo = FALSE, fig.align='center'}
LA <- read.csv("http://andrewpbray.github.io/data/LA.csv")
plot(price ~ sqft, data = LA, col = "steelblue")
m1 <- lm(price ~ sqft + bed + city, data = LA)
```

In the data set LA, this scatterplot suggests two influential points but are
they influential in a MLR model?


## Activity 9 (Part I)

In the data set LA, this scatterplot suggests two influential points but are
they influential in a MLR model?

1. Fit the model $\hat{price} \sim sqrt + bed + city$.
2. By the rules of thumb, are those two points high leverage?  Outliers?
(you can extract the hat values using `influence(m1)$hat`.)
3. Calculate the Cook's distance of those two observations using the equation
provided on the last slide.
4. Generate the Cook's distance plot to double check that the values that you
calculated in 3 seem correct.








