07-fallacies-correlation-and-regression.qmd

---
aliases:
  - fallacies-correlation-and-regression.html
---
# Fallacies, Correlation and Regression


```{r}
#| include: false
source("./_common.R")
```

> ... in which we hear Stories of Warplanes,
  Correlation and Regression and explore the Datasaurus Dozen.
  

{{< youtube vX1WGlFNtBo >}}

## Setup

```{r}
library(tidyverse)
library(glue)
library(broom)
```

## Data Considerations

### 1943

It is 1943. The second World War is well underway, ravaging large parts of Europe.
Military aircraft that had first entered the stage in World War I are now
reaching their peak importance as they rain fire from the skies. But the Allied
forces are facing a problem. As warplanes get better, so do anti-aircraft
systems. In an effort to improve the survival of their fleet, the US military
starts examining the planes returning from skirmishes with the opposing forces.
They characterize the pattern of bullet holes in the metal hull, meticulously
noting down each hit that the plane sustained. The resulting picture is better
summarized in the modern, redrawn version below.

![Figure from Wikipedia @SurvivorshipBias2020.](images/survivorship-bias.png)

After taking a look at the data they gathered, the military is ready to rush
into action. To improve the endurance of their aircraft, the plan is to
reinforce the parts of the plane that were most often hit by bullets. With
stronger wings and a sturdier body of the plane, they think, surely more pilots
will come back from their missions safely. They were wrong.

But the pilots where in luck. The military also consulted with the Statistics
Research Group at Columbia University. A man named Abraham Wald worked there. In
his now unclassified report "A method of estimating plane vulnerability based on
damage of survivors", he argued against the generals' conclusion
[@waldReprintMethodEstimating1980]. Instead of the most-hit parts of the planes,
the least-hit parts are to be reinforced.

![Cover of "A method of estimating plane vulnerability based on damage of survivors" @waldReprintMethodEstimating1980](images/paste-57122EF0.png)

> Instead of the most-hit parts, the least-hit parts are to be reinforced.

The reason for this seemingly counterintuitive result is what is now known as
*survivorship bias*. The data that was collected contained only survivors, those
planes that sustained damage not severe enough to hinder them from coming back
after their mission. The aircraft that where hit in other places simply didn't
make it back. Consequently, Wald advised to reinforce the engines and the fuel
tanks.

### Thinking further

This is but one of a multitude of biases, specifically a selection bias, that
will influence the quality of the inferences you can draw from available data.
Keep in mind, data is not objective and never exists in a vacuum. There is
always context to consider. The way the data was collected is just one of them.
A lot of these ideas seem obvious in hindsight, which incidentally is another
bias that social psychologists call *hindsight bias*, but they can sometimes be
hard to spot.

A common saying is that music was better back in the days, or that all the old
music still holds up while the new stuff on the radio just sounds the same.
Well, not quite. This is also survivorship bias at work. All the bad and
forgettable songs from the past just faded into oblivion, never to be mentioned
again, while the songs people generally agreed to be good survived the ravages
of time unscathed.
A similar thing happens with success in general, not just songs.
If you ask any CEO high up the corporate ladder, a millionaire, or the author of a
book that reads "How to get rich", they are sure to have a witty anecdote about
how their persistence, or their brilliance, or charisma got them to where they are
now. What we are not seeing is all the people just as witty, just as charismatic
or even just as persistent that where simply not as lucky. Very few people will
tell you this. Because it takes a whole lot of courage to admit that ones
success is based on luck and privilege.

And to take it back to the scientific context: When you are planning an
experiment for the lab, always ask whether your data collection process can in
some way be biased towards what you are trying to show.

I leave you with this:

![](./images/survivorship-tweet.png)

And from this cautionary tale we jump straight back into RStudio.

## Sidenotes

### Glue and Inline R Code

Using `paste` to create a text in which the
values of variables are inserted can be painful.

```{r}
name <- "Jannik"
age <- 26
text <- paste(name, "is", age, "years old.")
text
```

The `glue` package makes it a breeze.
Everything inside of curly braces in the text inside of the `glue`
function will be evaluated as regular R code,
enabling us to write text quite naturally: 

```{r}
text <- glue("{name} is {age} years old.")
text
```

I hope you are not too confused by the package and it's main function having the
same name.

```{r}
glue("{name} is {age + 10} years old.")
```

### Inline R code

Using the a backtick followed by the letter `r` we
can add the results of code right into the text sections
of Rmarkdown reports:

1 + 1 = `r 1 + 1`.

`r text`

### Best Practices

Speaking of being careful.
There is one rule I can give you to make your data analysis more secure:

> **Your raw data is sacred!**
  Do not ever modify it or save over it.

This is even more important when,
for example, using excel to preview a csv file.
Under no circumstances should you hit the save button
in excel when you are looking at the raw data.
With approximately one-fifth of genomic research papers containing
errors in the gene lists, because excel converted genes
such as _SEPT2_ (Septin 2) into dates, you can see why [@ziemannGeneNameErrors2016].
Biologists have since given up and renamed the genes that
where commonly converted into dates... but the point still stands.
This caution is of course also necessary when analyzing data
with R, not just excel. When we read in the raw data and
save a processed version, we create a new file, or even
better, a new folder for it. A good convention for example
would be do divide your data into a `raw` and `derived` folder.

## Covariance, Correlation and Regression

![Source: https://xkcd.com/552/](images/correlation.png)

Last week, we talked about a measure of the spread of a
random variable called the **variance**.

$$var(X) = \frac{\sum_{i=0}^{n}{(x_i-\bar x)^2}}{(n-1)}$$

Today, we are extending this idea to 2 random variables.
Because the normal distribution is so common, we are using
two normally distributed variables.


```{r}
N <- 50
df <- tibble(
  x = rnorm(N),
  y = rnorm(N)
)

m_x <- mean(df$x)
m_y <- mean(df$y)

ggplot(df, aes(x, y)) +
  geom_vline(xintercept = m_x, alpha = 0.8, color = "midnightblue") +
  geom_hline(yintercept = m_y, alpha = 0.8, color = "midnightblue") +
  geom_point(fill = "white", color = "black")
```

We also added lines for the
means of the two random
variables. Maybe I should have mentioned
this more clearly earlier on,
but the general convention in statistics is that random variables
are uppercase and concrete values from the distribution have the
same letter but lowercase.

We now get the **covariance** of X and Y as:

$$cov(X,Y)=\text{E}\left[(X-\text{E}\left[X\right])(Y-\text{E}\left[Y\right])\right]$$

The expected value $E[X]$ is just a fancy way of saying
the mean of X.
If we asses the contribution of individual points towards the
covariance, we can understand it quite intuitively.
A point that has a higher x than the mean of X and a higher
y than the mean of Y (top right quadrant) will push the covariance towards
positive values. Likewise, a point in the bottom left quadrant
will have negative differences with the X and Y mean, which cancel
each other out to result in a positive covariance.
The bottom right and top left quadrants push towards a negative
covariance. A mix of positive and negative contributions will
result in a covariance with a small absolute value.

The covariance has one problem: It will have weird units
(X times Y) and the scale is different depending on the random
variables.
So what we do is standardize it by dividing by both standard
deviations and get the **correlation coefficient**:

$$cor(X,Y)=\frac{cov(X,Y)}{\sigma_{X}\sigma_{Y}}$$

It can assume values between -1 and 1. It's full name is
_Pearson product-moment correlation coefficient_, or
_pearsons R_. We can square it to get $R^2$ (obviously),
which indicates the strength of the correlation with
values between 0 and 1 independent of the direction.
We will meet it again later.

Let us apply our knowledge to a new dataset.

### Introducing the Dataset

The `dplyr` package includes and example dataset of Star Wars
characters. Unfortunately, it was created a while ago,
so the is no baby yoda, but 87 other characters are present.

![I guess it is the baby yoda show now.](images/baby_yoda.gif)

```{r}
starwars
```

Let's look at some correlations:

### Pearson vs. Spearman (not a Boxing Match)

To compute pearsons correlation, we use the `cor` function in R.
Instead of filtering out `NA`, we can use
`use = "complete.obs"` to ignore rows with missing values in the computation.

```{r}
pearson <- cor(starwars$height, starwars$mass, use = "complete.obs")
pearson
```

When I first did this I was surprised that the correlation was so
low. We are after all talking about height and mass, which
I assumed to be highly correlated.
Let us look at the data to see what is going on.

```{r}
label_text <- glue("Pearson correlation: {round(pearson, 2)}")

jabba <- filter(starwars, str_detect(name, "Jabba"))
jabba_text <- list(x = 1100, y = 120)

starwars %>% 
  ggplot(aes(mass, height)) +
  geom_point() +
  annotate(geom = "text", x = 500, y = 75, label = label_text,
           hjust = 0) +
  annotate(geom = "curve",
           x = jabba_text$x, y = jabba_text$y,
           xend = jabba$mass, yend = jabba$height,
           curvature = .3,
           arrow = arrow(length = unit(2, "mm"))) +
  annotate(geom = "text",
           x = jabba_text$x,
           y = jabba_text$y, label = "Jabba the Hutt",
           hjust = 1.1) +
  xlim(0, 1500) +
  labs(x = "mass [kg]",
       y = "height [cm]")
```

This is the culprit! We have a massive outlier,
in all senses of the word "massive".
Luckily, there is another method to asses correlation.
Spearman's method is more resistant to outliers,
because the data is transformed into ranks first,
which negates the massive effect of outliers.


```{r}
spearman <- cor(starwars$height, starwars$mass,
                use = "complete.obs", method = "spearman")
spearman
```

Visually, this is what the points look like
after rank transformation:

```{r}
label_text <- glue("Spearman rank correlation: {round(spearman, 2)}")

starwars %>% 
  mutate(mass = rank(mass),
         height = rank(height)) %>% 
  ggplot(aes(mass, height)) +
  geom_point() +
  annotate(geom = "text", x = 0, y = 75, label = label_text,
           hjust = 0) +
  labs(x = "rank(mass)",
       y = "rank(height)")
```

Apart from `cor`, there is also `cor.test`, which gives more information. 
If we so fancy, we can use `broom` to turn the test output into
a tidy format as well.

```{r}
cortest <- cor.test(starwars$mass, starwars$height,
                # method = "spearman",
                use = "complete.obs")

cortest
```

```{r}
tidy(cortest)
```

There is another way we can specify which features to correlate.
`corr` also takes a matrix or data frame as it's x argument instead
of x and y.
We then end up with the pairwise correlation coefficients
for all columns of the dataframe.

This is known as a correlation matrix, and we can create it for
more than two features, as long as all features are numeric
(after all, what is the correlation between 1,4 and "cat" "dog"?).
Unfortunately there are only three numeric columns
in the `starwars` dataset, which makes for a pretty boring
correlation matrix.

```{r}
starwars %>%
  select(where(is.numeric)) %>% 
  head()
```

So let's look at another built-in dataset instead.
`mtcars` has some data about cars, like their
engine displacement or miles per gallon.

```{r}
mtcars %>% 
  ggplot(aes(disp, mpg)) +
  geom_point()
```

This makes for a much more interesting correlation matrix:

```{r}
cor(mtcars) %>% 
  as_tibble(rownames = "feature") %>% 
  pivot_longer(-feature) %>% 
  ggplot(aes(feature, name, fill = value)) +
  geom_raster() +
  geom_text(aes(label = round(value, 2))) +
  scale_fill_gradient2(low = "blue", high = "red",
                       mid = "white", midpoint = 0)
```

If you are working a lot with correlations, it is certainly
worth checking out the `corrr` package from the tidymodels framework:

<aside>
<a href="https://corrr.tidymodels.org/">
![](images/corrr.png){width=200}
</a>
</aside>

Its functions make these steps easier.

```{r}
corrr::correlate(mtcars) %>% 
  corrr::stretch()
```

And give use access to two different types of plots out of the box.

```{r}
corrr::correlate(mtcars) %>% 
  corrr::rplot()
```


```{r}
corrr::correlate(mtcars) %>% 
  corrr::network_plot()
```

### Difference to Linear Regression

Finally, linear regression is
a related concept, because both correlation and
linear regression quantify the strength of a linear
relationship. However, there are key differences.
When we fit a linear model like:

$$y \sim a + x * b$$

there is no error in x. We assume x is something that
is fixed, like the temperature we set for an experiment
or the dosage we used. Y on the other hand is a random
variable. In `cov(X,Y)` and `cor(X,Y)`, X and Y are both random variables,
usually things we observed, not set ourselves.

While the correlation coefficient is symmetrical and translation-scale-invariant:

$$cor(X,Y)=cor(Y,X)$$



$$cor(X,Y)=cor(X * a +b,Y * c + d)$$

The same is **not** true for linear models!

Let us look at an example where linear regression
is more appropriate than correlation.
In the `data` folder we find the IMDB ratings for 10
Star Wars movies (plus more features).


```{r}
ratings <- read_rds("data/07/starwars_movies.rds")
ratings
```

We can fit a linear model to see if the production year
has an effect on the rating.

```{r}
model <- lm(imdbRating ~ year, data = ratings)

augment(model) %>% 
  ggplot(aes(year, imdbRating)) +
  geom_smooth(method = "lm", alpha = 0.3, color = "midnightblue") +
  geom_segment(aes(x = year, y = .fitted,
                   xend = year, yend = imdbRating),
               alpha = 0.4) +
  geom_point()
```

What I added here as gray segments are the so called **residuals**.
They are what makes linear regression work.
It's full name is **Ordinary Least Squares** and the squares in
question are the squares of these residuals, the word _least_
indicates that these squares are minimized in order to find the
best fit line.

```{r}
broom::tidy(model)
```

Looks like every year decreases the estimated rating by 0.03.

One thing however is the same between correlation and
linear regression, and that is the $R^2$ value we get
from both calculations:

```{r}
summary(model)
```

We can interpret $R^2$ as the fraction of the variance of
the response variable y that can be explained by the
predictor x.

## Non-linear Least Squares

So far, we only properly dealt with linear relationships
and now it is time to get non-linear.
We will be creating a mechanistically driven
predictive model, so we have a formula of which
we want to adjust the parameters so that it fits our data.

Let's take classical Michaelis-Menten-Kinetics
There is a dataset for enzyme reaction rates included in R.
But we convert it from a 
dataframe to a tibble so that it prints
nicer:

```{r}
puromycin <- as_tibble(Puromycin)
puromycin
```

The initial rate $v_0$ of the an
enzymatic reaction was measured
for a control and a sample treated
with puromycin at different substrate
concentrations.
For every concentration we
have two replicates except for
one missing replicate.

```{r}
puromycin %>% 
  ggplot(aes(conc, rate, color = state)) +
  geom_point()
```


From our Biochemistry studies, we know
that we can express the rate depending
on the concentration with the following
formula:

$$rate=\frac{(Vm * conc)}{(K + conc)}$$

To make it easier to work with, let's
turn it into a function.

```{r}
rate <- function(conc, Vm, K) {
  Vm * conc / (K + conc)
}
```

Let's pick some arbitrary starting values.
For example, we see that the maximal velocity
could be around 200.
We also know that K is the concentration at which the half-maximal
velocity is reached.

```{r}
puromycin %>% 
  ggplot(aes(conc, rate, color = state)) +
  geom_point() +
  geom_function(fun = ~ rate(conc = .x, Vm = 200, K = 0.2),
                color = "black")
```

`geom_function` expects a function of x or an anonymous function
where the first argument is the values on the x-axis,
so this is what we did.
Well, I bet we can do better than guessing the function!
What R can do for us is the same it did for linear least squares
and that is minimizing the distance of our curve to the
datapoints.
This is the job of the `nls` function, which stands for
**Nonlinear Least Squares**.

### One model

Let's look at just the "treated" data first.

```{r}
treated <- filter(puromycin, state == "treated")
model <- nls(rate ~ rate(conc, Vm, K),
             data = treated,
             start = list(Vm = 200, K = 0.3)
             )

model
```

NlS needs starting values, so we use any guess that isn't too far off.
If it is completely wrong, the model doesn't know in which direction it
should move the parameters to improve the fit and we get an error like this:
`Error in nls(rate ~ rate(conc, Vm, K), data = puro, subset = state ==  : singular gradient`

> For this special case, R also has a self-starting model. I won't go
  into it because it is not as useful as the general concept of fitting
  arbitry functions, but you can check out `SSmicmen` for a model that
  estimes the starting values automatically.

Additionally, `nls` takes an argument `subset`, which works
like the `dplyr` verb `filter` so that we can fit
the model on a subset of the data without having to create it beforehand.

We use the `broom` package to display our model parameters in a tidy tibble.

```{r}
tidy(model)
```

With the base-R function `predict` we can make new predictions
based on a model and new data:

```{r}
head(predict(model, newdata = list(conc = seq(0, 1, 0.01))))
```

We can use the same function inside of `geom_function`:

```{r}
treated %>% 
  ggplot(aes(conc, rate, color = state)) +
  geom_point() +
  geom_function(fun = ~ predict(model, newdata = list(conc = .x)),
                color = "black")
```

Or alternatively create a new dataset of predictions
beforehand and use that with `geom_line`:

```{r}
predictions <- tibble(
  conc = seq(0, 1, 0.01),
  rate = predict(model, newdata = list(conc = conc))
)

treated %>% 
  ggplot(aes(conc, rate, color = state)) +
  geom_point() +
  geom_line(data = predictions, color = "black")
```

```{r}
augment(model) %>% 
  ggplot(aes(conc, .resid)) +
  geom_point()
```

### Multiple models

Now, what if we want to fit the model for both states?
We can resort back to our trusty `purrr` package like
we did in an earlier lecture.

We start out by creating a function that takes
a dataframe and fits our model:

```{r}
fit_micmen <- function(data) {
  nls(rate ~ rate(conc, Vm, K),
      data = data,
      start = list(Vm = 200, K = 0.3)
  ) 
}
```

And by nesting the data (grouped by state)
into a list column we can map
this function over each dataset.
And to get the fitted parameters
we map the `tidy` function from `broom`
over the fitted models.

```{r}
models <- puromycin %>% 
  group_by(state) %>% 
  nest() %>% 
  mutate(
    model = map(data, fit_micmen),
    params = map(model, tidy)
  )
```

Let's inspect the fitted parameters.

```{r}
models %>% 
  select(state, params) %>% 
  unnest(params) %>% 
  select(state, term, estimate) %>% 
  pivot_wider(names_from = term, values_from = estimate)
```

To plot our fitted models we have two options.
Firstly, we could generate the predicted values
for a number of concentrations beforehand
and then plot these:

```{r}
make_predicions <- function(model) {
  tibble(
    conc = seq(0, 1.2, 0.01),
    rate = predict(model, newdata = list(conc = conc))
  )
}

predictions <- models %>% 
  mutate(
    preds = map(model, make_predicions)
  ) %>% 
  select(state, preds) %>% 
  unnest(preds)

puromycin %>% 
  ggplot(aes(conc, rate, color = state)) +
  geom_point() +
  geom_line(data = predictions)
```

Or we use `geom_smooth`, which can take "nls" as a method as well.
We just need to make sure to pass the correct arguments.
And it can be confusing, because when we are specifying
the formula in `geom_smooth`, it always needs to be 
a formula of `y ~ x`, whereas in the normal `nls` we did
earlier, we specified the variables in terms of their
actual names (`rate` and `conc`).

```{r}
puromycin %>% 
  ggplot(aes(conc, rate, color = state)) +
  geom_point() +
  geom_smooth(method = "nls",
              formula = y ~ rate(conc = x, Vm, K),
              method.args = list(start = list(Vm = 200, K = 0.3)),
              se = FALSE
              )
```

We also need `se = FALSE`, because by default R would
try to plot a confidence interval around the fit-line
like it did for the linear model, but `nls` doesn't return one,
so we would get an error.

The unfortunate thing about this method is that we end up
fitting the model twice, once to get the estimated parameters
and the likes for ourselves and a second time in ggplot
to display the fitted lines. But in most cases this is not
a problem, because the model is not very computationally expensive.

### Excursion: A weird error message

Finally, I want to take minute to mention another approach
which we took earlier in the series when we where fitting
many linear models and show you, why it unfortunately
does not work here.

```{r}
newmodels <- puromycin %>% 
  group_by(state) %>% 
  summarise(
    model = list(nls(rate ~ rate(conc, Vm, K),
                start = list(Vm = 200, K = 0.3))
    )
  )
```

At first it looks like everything is fine.
Because we are inside of a dplyr verb `nls` know
where to look for the columns `rate` and `conc` that it should
fit, so we are not specifying its `data` argument.
However, this fails in an unexpected way when we later try
to make predictions with one of the models:

```{r, error=TRUE}
make_predicions(newmodels$model[[1]])
```

The reason for this as follows:
When `nls` fit the model it didn't remember the actual values
of `rate` and `conc`, it just made a note that these are columns
available in the data.
And because the data was not passed to it explicitly it
just wrote down that the columns are available in the
environment in which it was called, which at that
time was inside of `summarise`.
Check out the `data` argument here:

```{r}
newmodels$model[[1]]
```

It just says `parent.frame`, meaning "the environment around me".
But once it has left the context of `summarise`,
this is no longer available, so it can't find the `rate` column.
This is why is is always safer to pass the data
explicitly like we did in the approach that worked.

## Exercises

### The Datasaurus Dozen

The Datasaurus Dozen [@matejka2017] is a dataset
crafted to illustrate certain concepts.
It can be accessed from R via the `datasauRus` package.

```{r}
#| eval: false
datasauRus::datasaurus_dozen
```

- Explore the dataset before looking at the publication
  above (it contains spoilers...):
  - It actually contains 13 different datasets,
    denoted by the column `dataset`, in one tibble.
    What are the means for x and y for the different datasets?
    What are the standard deviations for x and y for the different datasets?
    What are the correlations coefficients for the different datasets?
    I bet you notice a pattern by now.
  - Now create one (or multiple) scatterplots of the data.
    What do you notice? what conclusions do you draw from this observation?
    
There is another dataset in the package to illustrate a different
point:

```{r}
#| eval: false
datasauRus::box_plots
```

- First, turn it into a tidy format, much like the `datasaurus_dozen`
  tibble.
- Now, visualize the distributions of the values for the
  5 different groups. Try out different versions of your plot
  until you are satisfied, but be sure to also include a boxplot
  and compare it to your approaches.
  What do you find?

### Fit a non-linear model

I found this gloriously 2000s website for "Statistical Reference Datasets":
<https://www.itl.nist.gov/div898/strd/index.html>
by the Information Technology Laboratory.
Not only has this official website of the United Stats Government
amazing unapologetic Word-Art,
it also features some handy datasets to
practice fitting non-linear models (<https://itl.nist.gov/div898/strd/nls/nls_main.shtml>)!

Of these I chose one for you to explore:
<https://itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Chwirut2.dat>

Because you might come across some challenges, I am leaving 
some tips below, hidden behind `details` panels,
so you can choose if and when you need them:

<details>
<summary>
Tip 1
</summary>
You can read in data that is separated by whitespace with `readr`s function
`read_table`.
</details>

<details>
<summary>
Tip 2
</summary>
You have to skip the first 60 lines and set column names manually.
</details>


<details>
<summary>
Tip 3
</summary>
The description in the dataset header also contains
the function to fit and potential starting values to try out.
Note, the `e` in the function refers to the remaining
error of the fit, so you don't need it in your function.
</details>