---
title: "Scatter plots and regression"
#author: "Mark"
date: "9/13/2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

[Last time](09.11.17.ScatterCorrelation.html), we explored scatter plots and correlation. Today, we'll use a numerical tool called linear regression to quantify the patterns we see in a scatter plot. While correlation can help us see *if* there's a linear relationship, regression can help us postulate a formula for that relationship

## Our real world examples revisited
### Hurricane prediction

This scatter plot relating error in hurricane prediction to the year the prediction was made
uses the same data we used
before. In addition to the scatter plot and the correlation coefficient $R$,  we see a formula for a line and a graph of that line superimposed on the data.

```{r, echo=F}
df = data.frame(year=1970:2010, error=c(c(270,380,695,370,360,410,420,490,415,265,400,410,280,400,390,370,387,373,265,290,305,303,300,292,360,210,180,200,198,195,200,180,170,150,153,155,155,153,150,145,150)))
plot(df$year, df$error, ylim=c(100,700), xlab = 'Year', ylab = 'Prediction error', 
     pch=21,  bg="blue")
#title(paste('correlation = ', sprintf('%.5f', 
#    cor(df$year, df$error, use='na.or.complete')))
#)
hurricane_fit = lm(df$error ~ df$year)
abline(hurricane_fit, lwd=2)
title(
  sprintf('R = %.4f\n E = %.2f Y + %.2f', 
    cor(df$year, df$error, use='na.or.complete'),
    hurricane_fit$coefficients['df$year'], hurricane_fit$coefficients['(Intercept)'])
)
```

The correlation of about $-0.83$ is a quantitative assessment of the relationship between the variables. 

The formula $E=-8.33Y+16870.37$ yields an estimate of the error $E$ in terms of the year $Y$.

#### Height vs weight

Here's another look at our CDC data relating height and weight.

```{r, echo=F}
set.seed(6)
cdc = read.csv("../../../data/cdc.csv")
men = subset(cdc, gender=='m')
subset = men[sample(1:20000,2000),]
cdc_fit = lm(subset$weight~subset$height)
plot(subset$height, subset$weight, pch=21,  bg="blue", xlab='height', ylab='weight')
abline(cdc_fit, lwd=2)
#title(paste('correlation = ', sprintf('%.5f', 
#    cor(subset$weight,subset$height, use='na.or.complete')))
#)
title(
  sprintf('R = %.4f\n W = %.2f H%.2f', 
    cor(subset$height, subset$weight, use='na.or.complete'),
    cdc_fit$coefficients['subset$height'], cdc_fit$coefficients['(Intercept)'])
)
```

Again, the correlation of about $0.47$ is a quantitative assessment of the relationship between the variables. 

The formula $W=5.38H-190$ yields an estimate of the weight $W$ in terms of the height $H$.

#### Galileo's ramp

Even though the correlation is high for Galileo's ramp, a quadratic fit is much better 
than a linear fit.

```{r, echo=F}
x =  c(573,534,495,451,395,337,253)
y = c(1000,800,600,450,300,200,100)
f = function(x) {return(600.054 - 4.00484*x + 0.00818074*x^2)}
xs = 230:600
plot(xs,f(xs), type='l', lwd=2, xlab="elevation", ylab="distance")
points(x,y, pch=21,  bg="blue")
fit = lm(y~x)
abline(fit)
title(paste('R = ', sprintf('%.5f', 
    cor(x,y)))
)
```

## A little mathematics

There are a few mathematical details you should know to understand linear regression.

### Equations of lines

The graph of an equation of the form $y=ax+b$ is a line with slope $a$ and $y$-intercept $b$. Here's a simple example.

```{r, echo=F}
plot(function(x) x*2+1, xlim=c(-2,2), ylab='',xlab='', asp=1, main='y=2x+1')
```

We typically say that $x$ is the independent variable and that $y$ is the dependent variable. 
In stat speak, this is ... ?

Given an equation like this, it's easy to compute the dependent variable in terms of the independent. If $y=2x+1$ and $x=3$, then $y = 2\times3+1 = 7$.

### Prediction

In the context of working with data, the regression line can be used to predict values for the response variable given a value of the explanatory variable that is missing or outside the range of the given data.

#### Example 

The heights in our CDC data are all rounded to inches. What does the model predict should be the height of a man who is 5'10 1/2'' tall?

**Solution**: The formula for the regression line is $w = 5.38h - 190.21$. When the height $h$ is $70.5$, we get a weight of $w=189$ pounds.

Note that the formula is reasonably good over a *range* of values but not good for all values.

#### Question: 

What does our hurricane forecast model predict for our error in 2020?

### Least squares

What makes the regression line better than some other line? It minimizes the total squared error.

More specifically, if we have observations 
$$(x_1, y_1),(x_2, y_2),\ldots,(x_n, y_n)$$
and we assume that our model is $Y=aX+b$, then the regression line chooses $a$ and $b$ so that
$$(y_1-(ax_1+b))^2 + (y_2-(ax_2+b))^2 + \cdots + (y_n-(ax_n+b))^2$$
is as small as possible.

### Residuals

The terms $y_i-(ax_i+b)$ are the called the *residuals*; they measure how far off our regression model is at the data points.

Often, examining the residuals of a scatter plot can be illuminating. Here are several scatter plots together with scatter plots of their residuals. The scatter plot on the left below looks slightly non-linear; the scatter plot of the residuals make this much more clear.

```{r, message=F, echo=F}
library(openintro)
GenerateLmPlot <- function(x, y, xlim, ylim1, ylim2.mult) {
  plot(x, y,
       axes = FALSE,
       pch = 20,
       col = COL[1, 2],
       cex = 1.202,
       xlim = xlim,
       ylim = ylim1)
  box()
  g <- lm(y ~ x)
  abline(g, col = COL[5], lwd=2)
  plot(x, g$residuals,
       pch = 20,
       col = COL[1, 2],
       cex = 1.202,
       xlim = xlim,
       axes = FALSE,
       ylim = ylim2.mult * c(-1, 1) * max(abs(g$residuals)))
  box()
  abline(h = 0, col = COL[5], lty = 2, lwd=2)
}

par(mfrow = c(1,2), mar = rep(0.5, 4))

MyLayOut <- matrix(1:2, 1)
layout(mat = MyLayOut,
       widths = rep(2, 2),
       heights = c(1.9),
       respect = TRUE)

# n <- 25
# x <- runif(n[1])
# y <- -8 * x + rnorm(n[1])
# GenerateLmPlot(x, y,
#                xlim = c(-0.03, 1.03),
#                ylim1 = c(-10, 1),
#                ylim2.mult = 2.5)
# 
n <- 30
x <- c(runif(n[1] - 2, 0, 4), 2, 2.1)
y <- -2 * x^2 + rnorm(n[1])
GenerateLmPlot(x, y,
               xlim = c(-0.2, 4.2),
               ylim1 = c(-35, 2),
               ylim2.mult = 1.8)

# n <- 40
# x <- runif(n[1])
# y <- 0.2 * x + rnorm(n[1])
# y[y < -2] <- -1.5
# GenerateLmPlot(x, y,
#                xlim = c(-0.03, 1.03),
#                ylim1 = c(-2, 2),
# ylim2.mult = 1.2)
title("A scatter plot and its residual plot", outer = T, line=-2)
```

It can also be illuminating to look at a histogram of the residuals. The histogram below shows the distribution of the residuals for our weight estimate example. Since it appears to be normal, we can apply our knowledge of the normal distribution to obtain confidencde in our estimates.

```{r, echo=F}
hist(cdc_fit$residuals, col='gray', xlab='', ylab='', 
  main = 'Histogram of our model estimate residuals'
)
```

## When should we use linear regression?

We will learn quite a few techniques to analyze data this semester. It's always good idea to think carefully about which technique you want to use and whether it's genuinely applicable.


### Conditions for linear regression:

First and foremost, linear regression might be useful if you want to understand the relationship between two numerical variables. Furthemore, the data should satisfy a few more conditions:

 * **Linearity**:  
 The data should *look* approximately linear. If there is a nonlinear trend, a more advanced regression method should be applied.
 * **Nearly normal residuals**:  
 Generally the residuals must be nearly normal. 
 * **Constant variability**:  
 The variability of points around the least squares line remains roughly constant.
 * **Independent observations**:  
 Be particularly cautious with time series data, which are sequential observations. Such data are rarely independent.
 
## Connection with correlation
 
It seems that correlation should be connected with slope - and it is, via the formula
$$m = r\frac{s_y}{s_x}.$$
Furthermore, the regression line should go through the point $(\bar{x},\bar{y})$. These formulae make it relatively easy to compute a regression line for given data.

### Example

Suppose that $X$ is a data set with mean $90$ and standard deviation $5$; $Y$ is a data set with mean $74$ and standard deviation $4$. Furthermore, $X$ and $Y$ have a tight correlation of $r=0.85$.

 * What is the regression line connecting $Y$ to $X$?
 * What value of $Y$ does the regression line predict if $X=85$?
 
**Solution**

We know the regression line has the form $y=mx+b$. The slope is
$$m = r\frac{s_y}{s_x} = 0.85 \times \frac{4}{5} = 0.68.$$
We now know the regression line has the more specific form $y=0.68x+b$. Once we know this, we can plug the means for $X$ and $Y$ to get $b$:
$$74 = 0.68 \times 90 + b, $$
so
$$b = 74 - 0.68 \times 90 = 12.8.$$
Thus, the final equation of the line is 
$$y = 0.68x +12.8.$$

To solve the second part, we simply plug $x=85$ into our formula for the line to find the corresponding value of $y$.

 ------
 
 ------
 
 