Correlation

]

# Week 5: .fancy[Correlation]

### <svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M243.4 2.587C251.4-.8625 260.6-.8625 268.6 2.587L492.6 98.59C506.6 104.6 514.4 119.6 511.3 134.4C508.3 149.3 495.2 159.1 479.1 160V168C479.1 181.3 469.3 192 455.1 192H55.1C42.74 192 31.1 181.3 31.1 168V160C16.81 159.1 3.708 149.3 .6528 134.4C-2.402 119.6 5.429 104.6 19.39 98.59L243.4 2.587zM256 128C273.7 128 288 113.7 288 96C288 78.33 273.7 64 256 64C238.3 64 224 78.33 224 96C224 113.7 238.3 128 256 128zM127.1 416H167.1V224H231.1V416H280V224H344V416H384V224H448V420.3C448.6 420.6 449.2 420.1 449.8 421.4L497.8 453.4C509.5 461.2 514.7 475.8 510.6 489.3C506.5 502.8 494.1 512 480 512H31.1C17.9 512 5.458 502.8 1.372 489.3C-2.715 475.8 2.515 461.2 14.25 453.4L62.25 421.4C62.82 420.1 63.41 420.6 63.1 420.3V224H127.1V416z"/></svg> EMSE 4572: Exploratory Data Analysis
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M224 256c70.7 0 128-57.31 128-128s-57.3-128-128-128C153.3 0 96 57.31 96 128S153.3 256 224 256zM274.7 304H173.3C77.61 304 0 381.6 0 477.3c0 19.14 15.52 34.67 34.66 34.67h378.7C432.5 512 448 496.5 448 477.3C448 381.6 370.4 304 274.7 304z"/></svg> John Paul Helveston
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M96 32C96 14.33 110.3 0 128 0C145.7 0 160 14.33 160 32V64H288V32C288 14.33 302.3 0 320 0C337.7 0 352 14.33 352 32V64H400C426.5 64 448 85.49 448 112V160H0V112C0 85.49 21.49 64 48 64H96V32zM448 464C448 490.5 426.5 512 400 512H48C21.49 512 0 490.5 0 464V192H448V464z"/></svg> September 28, 2022

]

---

# Quick Rmd updates

---

# .fancy[.blue[Tip of the week]]

# All data are biased

---

---

## [.center[Abraham Wald]](https://en.wikipedia.org/wiki/Abraham_Wald)

]

]

---

## Today's data

```r
msleep <- read_csv(here::here('data', 'msleep.csv'))
```

## New packages:

```r
install.packages('HistData')
install.packages('palmerpenguins')
install.packages('GGally')
```

---

# Week 5: .fancy[Correlation]

## 1. What is correlation?
## 2. Visualizing correlation

## BREAK

## 3. Linear models
## 4. Visualizing linear models

---

# Week 5: .fancy[Correlation]

## 1. .orange[What is correlation?]
## 2. Visualizing correlation

## BREAK

## 3. Linear models
## 4. Visualizing linear models

---

## Some pretty racist origins in [eugenics](https://en.wikipedia.org/wiki/Eugenics) ("well born")

### [Sir Francis Galton](https://en.wikipedia.org/wiki/Francis_Galton) (1822 - 1911)

- Charles Darwin's cousin.
- "Father" of [eugenics](https://en.wikipedia.org/wiki/Eugenics).
- Interested in heredity.

]

### [Karl Pearson](https://en.wikipedia.org/wiki/Karl_Pearson) (1857 - 1936)

- Galton's ([hero-worshiping](https://en.wikipedia.org/wiki/Apotheosis)) protégé.
- Defined correlation equation.
- "Father" of mathematical statistics.

]

???

The beautiful irony is that human genetics was also the field that conclusively demonstrated the biological falsity of race.

---

# Galton's family data

Galton, F. (1886). ["Regression towards mediocrity in hereditary stature"](http://www.stat.ucla.edu/~nchristo/statistics100C/history_regression.pdf). _The Journal of the Anthropological Institute of Great Britain and Ireland_ 15: 246-263.

**Galton's question**: Does marriage selection indicate a relationship between the heights of husbands and wives?<br>(He called this "assortative mating")

"midparent height" is just a scaled mean:

```r
midparentHeight =  (father + 1.08*mother)/2
```

]

```r
library(HistData)

galtonScatterplot <- ggplot(GaltonFamilies) +
    geom_point(aes(x = midparentHeight,
                   y = childHeight),
               size = 0.5, alpha = 0.7) +
    theme_classic() +
    labs(x = 'Midparent height (inches)',
         y = 'Child height (inches)')
```

]]

---
class: center, middle

# How do you measure correlation?

<br>

# Pearson came up with this:

# `$r = \frac{\text{Cov}(x, y)}{\text{sd}(x) * \text{sd}(y)}$`

---

# How do you measure correlation?

## `$r = \frac{\text{Cov}(x, y)}{\text{sd}(x) * \text{sd}(y)}$`

Assumptions:
1. Variables must be interval or ratio
2. Linear relationship

]]

]

---

# How do you _interpret_ `$r$`?

## `$r = \frac{\text{Cov}(x, y)}{\text{sd}(x) * \text{sd}(y)}$`

Interpretation:
- `$-1 \le r \le 1$`
- Closer to 1 is stronger correlation
- Closer to 0 is weaker correlation

]

```r
cor(x = GaltonFamilies$midparentHeight,
    y = GaltonFamilies$childHeight,
    method = 'pearson')
```

```
#> [1] 0.3209499
```

]

]

---

## What does `$r$` mean?

- `$\pm 0.1 - 0.3$`: Weak
- `$\pm 0.3 - 0.5$`: Moderate
- `$\pm 0.5 - 0.8$`: Strong
- `$\pm 0.8 - 1.0$`: Very strong
]]

]

---

# Visualizing correlation is...um...easy, right?

<br>

# [guessthecorrelation.com](http://guessthecorrelation.com/)

# Click [here](https://docs.google.com/presentation/d/1uN5FfQ4QiiaJ1SI1A5SRji58vWzubk2qwGi9Lykz5z0/edit?usp=sharing) to vote!

---

## The datasaurus

### (More [here](https://www.autodeskresearch.com/publications/samestats))

]

]

---

# Coefficient of determination: `$r^2$`

Percent of variance in one variable that is explained by the other variable

]]

`$r$` | `$r^2$`
----|------
0.1 | 0.01
0.2 | 0.04
0.3 | 0.09
0.4 | 0.16
0.5 | 0.25
0.6 | 0.36
0.7 | 0.49
0.8 | 0.64
0.9 | 0.81
1.0 | 1.00

]

---

## You should report both `$r$` and `$r^2$`

<br>

### Correlation between parent and child height is 0.32, therefore 10% of the variance in the child height is explained by the parent height.

---

# Correlation != Causation

### X causes Y

- Training causes improved performance

### Y causes X

- (Good / bad) performance causes people to train harder.

### Z causes both X & Y

- Commitment and motivation cause increased training and better performance.

---

## Be weary of dual axes!

## ([They can cause spurious correlations](https://www.tylervigen.com/spurious-correlations))

]

]

---

# Outliers

---

---

---

---

## Pearson correlation is highly sensitive to outliers

---

# Spearman's rank-order correlation

# `$r = \frac{\text{Cov}(x, y)}{\text{sd}(x) * \text{sd}(y)}$`

- Separately rank the values of X & Y.
- Use Pearson's correlation on the _ranks_ instead of the `$x$` & `$y$` values.

]

Assumptions:

- Variables can be ordinal, interval or ratio
- Relationship must be monotonic (i.e. does not require linearity)

]

---

## Spearman correlation more robust to outliers

---

## Spearman correlation more robust to outliers

]

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Pearson </th>
   <th style="text-align:right;"> Spearman </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> -0.56 </td>
   <td style="text-align:right;"> 0.53 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 0.39 </td>
   <td style="text-align:right;"> 0.69 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 0.94 </td>
   <td style="text-align:right;"> 0.81 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 0.38 </td>
   <td style="text-align:right;"> 0.76 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 0.81 </td>
   <td style="text-align:right;"> 0.79 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 0.31 </td>
   <td style="text-align:right;"> 0.70 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 0.95 </td>
   <td style="text-align:right;"> 0.81 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 0.51 </td>
   <td style="text-align:right;"> 0.75 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> -0.56 </td>
   <td style="text-align:right;"> 0.53 </td>
  </tr>
</tbody>
</table>

]

]

---

## Summary of correlation

- **Pearson's correlation**: Described the strength of a **linear** relationship between two variables that are interval or ratio in nature.
- **Spearman's rank-order correlation**: Describes the strength of a **monotonic** relationship between two variables that are ordinal, interval, or ratio. **It is more robust to outliers**.
- The **coefficient of determination** ( `$r^2$` ) describes the amount of variance in one variable that is explained by the other variable.
- **Correlation != Causation**

]

R command (hint: add `use = "complete.obs"` to drop NA values)

```r
pearson  <- cor(x, y, method = "pearson", use = "complete.obs")
spearman <- cor(x, y, method = "spearman", use = "complete.obs")
```

---

# Week 5: .fancy[Correlation]

## 1. What is correlation?
## 2. .orange[Visualizing correlation]

## BREAK

## 3. Linear models
## 4. Visualizing linear models

---

## **Scatterplots**: The correlation workhorse

```r
scatterplot <- ggplot(mtcars) +
*   geom_point(aes(x = mpg, y = hp),
               size = 2, alpha = 0.7) +
    theme_classic(base_size = 20) +
    labs(x = 'Fuel economy (mpg)',
         y = 'Engine power (hp)')

scatterplot
```

]

]

---

## Adding a correlation label to a chart

```r
corr <- cor(
    mtcars$mpg, mtcars$hp,
    method = 'pearson')

*corrLabel <- paste('r = ', round(corr, 2))
```

Add label to the chart with `annotate()`

```r
scatterplot +
*   annotate(geom = 'text',
*            x = 25, y = 310,
*            label = corrLabel,
*            hjust = 0, size = 7)
```
]

]

---

---

## Visualize all the correlations: `ggcorr()`

```r
library('GGally')
```

```r
mtcars %>%
*   ggcorr()
```

]

]

---

## Visualizing correlations: `ggcorr()`

```r
library('GGally')
```

```r
mtcars %>%
*   ggcorr(label = TRUE,
*          label_size = 3,
*          label_round = 2)
```

]

]

---

## Visualizing correlations: `ggcorr()`

```r
ggcor_mtcars_final <- mtcars %>%
    ggcorr(label = TRUE,
           label_size = 3,
           label_round = 2,
*          label_color = 'white',
*          nbreaks = 5,
*          palette = "RdBu")
```

]

]

---

## .center[Pearson]

```r
mtcars %>%
    ggcorr(label = TRUE,
           label_size = 3,
           label_round = 2,
*          method = c("pairwise", "pearson"))
```

]

## .center[Spearman]

```r
mtcars %>%
    ggcorr(label = TRUE,
           label_size = 3,
           label_round = 2,
*          method = c("pairwise", "spearman"))
```

]

---

## Correlograms: `ggpairs()`

```r
library('GGally')
```

```r
mtcars %>%
    select(mpg, cyl, disp, hp, wt) %>%
*   ggpairs()
```

- Look for linear relationships 
- View distribution of each variable

]

]

---

## Correlograms: `ggpairs()`

```r
library('GGally')
```

```r
mtcars %>%
    select(mpg, cyl, disp, hp, wt) %>%
    ggpairs() +
*   theme_classic()
```

- Look for linear relationships 
- View distribution of each variable

]

]

---

## Your turn

Using the `penguins` data frame:

1. Find the two variables with the largest correlation in absolute value (i.e. closest to -1 or 1).

2. Create a scatter plot of those two variables.

3. Add an annotation for the Pearson correlation coefficient.

]

### .center[[palmerpenguins library](https://allisonhorst.github.io/palmerpenguins/)]

]

---

## **Simpson's Paradox**: when correlation betrays you

]

]

---

# Quiz 2

### Link is in the #class channel

]

]

---

# Week 5: .fancy[Correlation]

## 1. What is correlation?
## 2. Visualizing correlation

## BREAK

## 3. .orange[Linear models]
## 4. Visualizing linear models

---

# Palmer Penguins

The correlation of 0.87 means that the body mass (g) explains about 75% of the variation in the flipper length (mm).

]

]

---

# Palmer Penguins

The correlation of 0.87 means that the body mass (g) explains about 75% of the variation in the flipper length (mm).

<br>

**Now let's fit a model to these points!**
]

]

---

## Modeling basics

Two parts to a model:

1. **Model family**: e.g., `$y = ax + b$`
<br>

2. **Fitted model**: e.g., `$y = 3x + 7$`

]

Here is some simulated data

]

---

## Modeling basics

Two parts to a model:
1. **Model family**: linear model: `$y = ax + b$`

]

There are an infinite number of possible models

]

---

## Modeling basics

Two parts to a model:

1. **Model family**: linear model: `$y = ax + b$`

<br>

2. **Fitted model**: How to choose the "best" `$a$` and `$b$`?

]

There are an infinite number of possible models

]

---

## Modeling basics

Two parts to a model:

1. **Model family**: linear model: `$y = ax + b$`
<br>

2. **Fitted model**: How to choose the "best" `$a$` and `$b$`?

]

Compute the **.blue["residuals"]**:<br>The distance between the model line and the data<br>

]

---

## Residual: `$y_i - y_i'$`

### **Residual**: The distance between the model line and the data

]

]

---

## Sum of squared residuals: `$\text{SSR} = \sum_{i = 1}^{n} (y_i - y_i')^2$`

### **Residual**: The distance between the model line and the data

]

]

---

## Search algorithm

1): Pick a model ( `$a$` and `$b$` ):<br>

`$y = ax + b$`

<br>

]

2): Compute the SSR:<br>

`$\text{SSR} = \sum_{i = 1}^{n} (y_i - y_i')^2$`

]

3): Repeat steps 1 & 2 until the smallest SSR is found
<br><br>

]

---

## Fitting a linear model in R

```r
model <- lm(formula = y ~ x, data = data)
```

Penguin data:

```r
model <- lm(
  formula = body_mass_g ~ flipper_length_mm,
  data    = penguins)
```

Get coefficients ( `$a$` & `$b$` in `$y = ax + b$` )

```r
coef(model)
```

```
#>       (Intercept) flipper_length_mm 
#>       -5780.83136          49.68557
```

---

## Fitting a linear model in R

```r
model <- lm(formula = y ~ x,
            data = data)
```

Penguin data:

```r
model <- lm(
  formula = body_mass_g ~ flipper_length_mm,
  data    = penguins)
```

Get coefficients

```r
coef(model)
```

```
#>       (Intercept) flipper_length_mm 
#>       -5780.83136          49.68557
```

]

]

---

## Interpreting results

---

## Example write up for Penguin data

The correlation between flipper length (mm) and body mass (g) is **0.87**. Therefore, **~75%** of the variance in body mass is explained by flipper length.

The slope of the best fitting regression line indicates that body mass increased by **49.7 g** as flipper length increased by one mm.

]]

]

---

## Making predictions

**Interpolation is OK**: You may predict values of `$y$` for values of `$x$` that were not observed but are within the range of the observed values of `$x$`.

]

**Extrapolation is DANGEROUS**: You generally should NOT predict values of `$y$` using values of `$x$` that are outside the observed range of `$x$`.

]

---

## Repeat: Extrapolation is **DANGEROUS**

"Extrapolation of these trends to the 2008 Olympiad indicates that the women’s 100-metre race could be won in a time of 10.57±0.232 seconds and the men’s event in 9.73±0.144 seconds. **Should these trends continue, the projections will intersect at the 2156 Olympics, when — for the first time ever — the winning women's 100-metre sprint time of 8.079 seconds will be lower than that of the men's winning time of 8.098 seconds (Fig. 1).**"

]

.font80[Tatem, A. J., Guerra, C. A., Atkinson, P. M., & Hay, S. I. (2004). Momentous sprint at the 2156 Olympics? _Nature_, 431(7008), 525-525. [View online](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3173856/)]

]

---

# Symantics

These all mean the same thing:

- "Use X to predict Y"
- "Regress Y _on_ X"
- "Regression of Y _on_ X"

```r
model <- lm(formula = y ~ x,
            data = data)
```

]]

]

---

# Symantics

```r
model <- lm(formula = y ~ x,
            data = data)
```

**Y: Dependent variable**
- Outcome variable
- Response variable
- Regressand
- Left-hand variable

]]

**X: Independent variable**
- Predictor variable
- Explanatory variable
- Regressor
- Right-hand variable

]]

---

# Week 5: .fancy[Correlation]

## 1. What is correlation?
## 2. Visualizing correlation

## BREAK

## 3. Linear models
## 4. .orange[Visualizing linear models]

---

## Adding the correlation annotation

```r
# Make the correlation label
corr <- cor(
    penguins$body_mass_g, 
    penguins$flipper_length_mm,
    method = 'pearson', 
    use = "complete.obs"
)

*corrLabel <- paste("r = ", round(corr, 2))

# Make the chart!
penguins %>% 
    ggplot(aes(x = flipper_length_mm, y = body_mass_g)) +
    geom_point(size = 1, alpha = 0.7) +
    theme_classic(base_size = 20) +
    labs(x = "Flipper length (mm)",
         y = "Body mass (g)") + 
*   annotate(
*       geom = 'text', x = 175, y = 6000,
*       label = corrLabel,
*       hjust = 0, size = 5)
```

]

]

---

```r
# Make model label
model <- lm(
  formula = body_mass_g ~ flipper_length_mm,
  data    = penguins)
coefs <- round(coef(model), 2)
*modelLabel <- paste('y = ', coefs[1], ' + ', coefs[2], 'x')

# Make the chart!
penguins %>% 
    ggplot(aes(x = flipper_length_mm, y = body_mass_g)) +
    geom_point(size = 1, alpha = 0.7) +
*   geom_smooth(method = 'lm', se = FALSE) +
    theme_classic(base_size = 20) +
    labs(x = "Flipper length (mm)",
         y = "Body mass (g)") + 
    annotate(geom = 'text', x = 175, y = 6000,
             label = corrLabel,
             hjust = 0, size = 5) +
*   annotate(
*       geom = 'text', x = 175, y = 5700,
*       label = modelLabel, color = "blue",
*       hjust = 0, size = 5)
```

]

## .center[Add model]

]

---

## Your turn

Using the `msleep` data frame:

1. Create a scatter plot of `brainwt` versus `bodywt`.
2. Include an annotation for the Pearson correlation coefficient.
3. Include an annotation for the best fit line.

Bonus: Compare your results to a log-linear relationship by converting the x and y variables to the log of x and y, like this:

```r
model <- lm(log(x) ~ log(y), data = msleep)
```

You can also convert your plot to log axes by adding these layers:

```r
plot +    
    scale_x_log10() +
    scale_y_log10()  
```