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

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

## Scatter plots

Often, we would like to know if two variables are related. A *scatter plot* is simply a plot of point pairs. This simple geometric tool can help us assess relationships


### Illustrations using actual data

#### Hurricane prediction

This timely scatter plot taken from our text illustrates the average error in hurricane prediction vs the year in which those predictions were made. The plot reveals the fact that predictions are generally improving.

```{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')))
)
#fit = lm(df$error ~ df$year)
#abline(fit, lwd=2)
```

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

#### Height vs weight

Taken from our CDC data.

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

#### Galileo's ramp

Taken from [Galileo's Gravity and Motion Experiments ](http://ww2.amstat.org/publications/jse/v3n1/datasets.dickey.html).

```{r, echo=F}
x =  c(573,534,495,451,395,337,253)
y = c(1000,800,600,450,300,200,100)
plot(x,y, xlab="elevation", ylab="distance", pch=21,  bg="blue")
title(paste('correlation = ', sprintf('%.5f', 
    cor(x,y)))
)
#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")
```

The correlation is actually quite misleading in this example because a quadratic relationship is more appropriate.


### Random illustrations

Here are a few scatter plots of randomly generated data to further illustrate the ideas.

#### A perfect linear relationship

```{r, echo = F}
x = rnorm(50)
y = 1+x/2
plot(x,y, asp=1, col='blue', pch=16)
title(paste('correlation = ', cor(x,y)))
#fit1 = lm(y~x, data = data.frame(x=x,y=y))
#abline(fit1, lwd=2)
```

#### A close to linear relationship

```{r, echo = F}
x = rnorm(50)
r = rnorm(50, 0, 0.2)
y = 1+x/2 + r
plot(x,y, asp=1, col='blue', pch=16)
title(paste('correlation = ', sprintf('%.5f', cor(x,y))))
#fit2 = lm(y~x, data = data.frame(x=x,y=y))
#abline(fit2, lwd=2)
```

#### A close to linear, but negative, relationship

```{r, echo = F}
x = rnorm(50)
r = rnorm(50, 0, 0.2)
y = 1-x/2 + r
plot(x,y, asp=1, col='blue', pch=16)
title(paste('correlation = ', sprintf('%.5f', cor(x,y))))
#fit3 = lm(y~x, data = data.frame(x=x,y=y))
#abline(fit3, lwd=2)
```

#### A nonlinear relationship


```{r, echo = F}
x = rnorm(50)
r = rnorm(50, 0, 0.1)
y = x+x^2 + r
plot(x,y, col='blue', pch=16)
title(paste('correlation = ', sprintf('%.5f', cor(x,y))))
#fit4 = lm(y~x, data = data.frame(x=x,y=y))
#abline(fit4, lwd=2)
```


#### No relationship


```{r, echo = F}
x = rnorm(80)
y = rnorm(80)
plot(x,y,col='blue', pch=16, xlim = c(-2,2))
title(paste('correlation = ', sprintf('%.5f', cor(x,y))))
#fit5 = lm(y~x, data = data.frame(x=x,y=y))
#abline(fit5, lwd=2)
```


## Quantifying the relationship

There are a few ways to quantify the relationships we see. Today, we'll learn about correlation. We'll talk about another technique, called regression, later.

When it comes to quantifying relationships, we'll generally assume that we're working with a list of data points that looks like 
$$(x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n).$$

### Correlation

The basic formula for the correlation of our list of data points is

$$
R = \frac{1}{n-1}\sum_{i=1}^{n} \frac{x_i-\bar{x}}{s_x}\frac{y_i-\bar{y}}{s_y} = 
\frac{1}{n-1}\sum_{i=1}^{n} z_{x_i}z_{y_i},
$$
where $\bar{x}$, $\bar{y}$, $s_x$, and $s_y$ are the sample means and standard deviations for each variable.

#### Comments

  * The correlation is always between $-1$ and $+1$ and
    * A number close to $+1$ indicates a strong, positive linear relationship,
    * A number close to $-1$ indicates a strong, negative linear relationship,
    * A number close $0$ indicates a weak linear relationship.
  * The simplified version, with the $z_{x_i}z_{y_i}$, emphasizes that we're really just multiplying the $z$-scores together and adding the results. 
    * The correlation tends to be larger when the signs of the $z$-scores agree.
    * The correlation tends to be smaller when the signs of the $z$-scores disagree.
  * The idea is illustrated in the figure below where
    * The green points contribute positively to the correlation and
    * The red points contribute negatively to the correlation

```{r, echo=F}
set.seed(1)
plot(c(-3,3),c(0,0), type='l', xlim=c(-3,3), ylim=c(-2,2), xlab="",ylab="")
lines(c(0,0),c(-3,3))
x = rnorm(50)
r = rnorm(50, 0, 0.5)
y = x + r
df = data.frame(x=x,y=y)
good = subset(df, x*y>0)
points(good$x,good$y, asp=1, col='green', pch=16)
bad = subset(df, x*y<0)
points(bad$x,bad$y, asp=1, col='red', pch=16)
#title(paste('correlation = ', sprintf('%.5f', cor(x,y))))
##fit2 = lm(y~x, data = data.frame(x=x,y=y))
```

#### A short computation

```{r, echo=F, eval=F}
# The data
set.seed(1)
x = round(10*runif(4))
y = round(x + runif(4,-3,3))
print(x)
print(y)
print(sd(x))
print(sd(y))
print( sqrt(((3-5.5)^2+(4-5.5)^2+(6-5.5)^2+(9-5.5)^2)/3))
print( sqrt(((1-6.5)^2+(6-6.5)^2+(9-6.5)^2+(10-6.5)^2)/3))
print(cor(x,y))
```

Let's suppose our list of $x$ values and our list of $y$ values are

$$
\begin{align}
 x &= 3,4,6,9 \\
 y &= 1,6,9,10
\end{align}
$$
Then, their averages are
$$
\begin{align}
  \bar{x} &= \frac{3+4+6+9}{4} = 5.5, \text{ and } \\
  \bar{y} &= \frac{1+6+9+10}{5} = 6.5
\end{align}
$$
and their standard deviations are
$$
\begin{align}
  s_x &= \sqrt{\frac{(3-5.5)^2+(4-5.5)^2+(6-5.5)^2+(9-5.5)^2}{3}} = 2.65, \text{ and } \\
  s_y &= \sqrt{\frac{(1-6.5)^2+(6-6.5)^2+(9-6.5)^2+(10-6.5)^2}{3}} = 4.04.
\end{align}
$$
Thus, the correlation is
$$
R = \frac{1}{3}\left(
  \left(\frac{3-5.5}{2.65}\right)\left(\frac{5-6.5}{4.04}\right) +
  \left(\frac{4-5.5}{2.65}\right)\left(\frac{7-6.5}{4.04}\right) +
  \left(\frac{6-5.5}{2.65}\right)\left(\frac{7-6.5}{4.04}\right) +
  \left(\frac{9-5.5}{2.65}\right)\left(\frac{10-6.5}{4.04}\right)\right) = 0.87
$$

## Correlation vs causation

In spite of it's importance, it should be understood that correlation does *not* imply causation. 


### Lurking variables

A *lurking* or *confounding* variable is an unspecified variable that might be related to the
variables under study.
In both of the following examples, scatter plots and correlation computations show a strong relationship between the variables:

  * Number of smartphones in a country and overall public health
  * Public funding of education in a school district and exam scores
  
Are there any lurking variables?

### Explanatory and response variables

If a researcher *suspects* genuine causation, they might designate one variable as *explanatory* and one as *response*. Let's try to identify the natural explanatory and response variables in the following situations:

  * Number of hours studying a week and GPA
  * GPA and Number of hours studying a week
  * Light level and reading comprehension
  * Per capita income and percentage of population with college degrees
  
------

------
