---
title: "Statistical Thinking: Week 2 Lab"
date: 
output:
  pdf_document: default
  html_document:
    df_print: paged
  word_document: default
---

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

# Introduction

In this week's lab, we are going to look at how to build synthetic experiments.
The main reson for doing this is it let's us see how the sample size of 
and experiment can effect the variability of the objects of interest.

We are going to learn the programming pattern for doing a simualtion study and
apply it in a few situations.

### Learning Goals for Week 1

- Learn how to simulate random numbers in R
- Learn how to simulate an experiment in tidy form
- Identify the random variation in an estimator
- Explore the randomness graphically.

# Random numbers on a computer

Computers, by and large, cannot produce a truly random sequence of numbers^[The wikipedia page on this goes into some detail https://en.wikipedia.org/wiki/Random_number_generation]. Instead, they produce things called _pseudo-random number sequences_, which are sequences of numbers that, for all intents and purposes, behave like random numbers.

While a description of how these work is well beyond the scope of this course, we do need to know that the numbers are not truly random. And this is not a bad thing. Because the numbers aren't truly random, we can make our code _reproduce_ the same results at a later time. (We would not be able to do this if the results were truly random!)

The key concept here is that of a _seed_, which is an integer that essentially starts the sequence of pseudo-random numbers. In R, we can manipulate the seed with the `set.seed()` function, which accepts an integer argument. You can use any integer you want. (In other applications, there can be security implications to using a short or common key, but that's not important here!)

So let's demonstrate this. I'm going to set my seed, and then draw some random numbers (in this case 5 poissons and one normal).

```{r}
set.seed(30127)
rpois(n = 5, lambda = 1)
rnorm(1, mean = 0, sd = 2)
```

Now if I re-set my seed to the same number and re-run those commands _in the same order_ I should get the same sequence of random numbers.

```{r}
set.seed(30127)
rpois(n = 5, lambda = 1)
rnorm(1, mean = 0, sd = 2)
```

_An important caveat:_ You need to run all of the commands in exactly the same order to get the same sequence of numbers. This is because each time a random number is needed we are taken one more step along the pseudorandom sequence, so if we miss one of the functions we may miss one of the steps!


# A simple experiment

Imagine I have 100 random numbers from a normal distribuiton with mean $3$ and standard deviation $1$. I can use R to get a realisation of this experiment.
```{r}
set.seed(30127)

draws <- rnorm(n = 100, mean = 3, sd = 1)
```

### Questions:

1. What is the mean of these draws? Is it exactly 3?

```{r, echo = FALSE}
mean(draws)
```

2. Repeat the "experiment" (in this case, the experiment is drawing 100 normal variables with mean 3 and standard devation 1) four more times and look at how much the sample mean changes.

# Repeating the experiment

If we want to explore how variable the sample mean of 100 normally distributed samples is, we are going to need to do the following two things a lot of times:

1. Simulate 100 $N(3,1)$ variables
2. Compute the sample mean

If we were only going to do this a handful of times (like above), we could do it manually. But If I want to repeat the experiment a thousand times, I'm going to need to tell R to do it.

We are going to have to do things like this a lot during the course, so I want tyou to pay close attention to the _programming pattern_ we use to  perform this task.

## Method: Enumerate the experiments

One way we could simulate 3 experiments manually would be to write code like 
```{r}
set.seed(30127)

draws1 <- rnorm(100, 3, 1)
draws2 <- rnorm(100, 3, 1)
draws3 <- rnorm(100, 3, 1)
```

If you run this code, you will see that all of the numbers are different and they will have different means.

But this is very labour intensive! The problem is the choice to use the name of the variable to encode which replicate of the experiment we are doing.

Instead of doing this, we can instead expand our data so that instead of just recording the random number we now record the pair `(experiment, random_number)`.  

```{r, message =FALSE}
library(tidyverse)
set.seed(30127)

draws_1 <- tibble(experiment = rep(1, times = 100), 
                 draw = rnorm(100, 3, 1))
draws_2 <- tibble(experiment = rep(2, times = 100), 
                 draw = rnorm(100, 3, 1))
draws_3 <- tibble(experiment = rep(3, times = 100), 
                 draw = rnorm(100, 3, 1))

```

(The function `rep(1, times = 100)` repeats the number 1 one hundred times.)

This is quite silly at the moment  because now I have the experiment number in two different places: in the actual data (good!) and in the 
variable name (bad!). 

So let's put this together in one variable. (To do this, I'm going to have to do something fancy with `rep()`)

```{r}
set.seed(30127)

draws_combined <- tibble(
  experiment = rep(1:3, each = 100),
  draw = rnorm(3 * 100, 3, 1)
)
```
The `rep()` function is pretty powerful. When you give it a sequence of numbers (`1:3` expands to `1, 2, 3`) you have two options:

- `times = n`: Repeats the whole sequence `n` times. `r rep(1:3,times =3)`
- `each = n`: Repeats each elemnt in the sequence `n` times  `r rep(1:3,each =3)`

The nice thing about this pattern is that I can use it to generate _any number_ of experiments!

I then need to compute the mean of the draws in each experiment. This is easier to do with the manual version as I could just write 
```{r}
mean(draws1)
```

So what is the equivalent here?

We need to use the `group_by`/`summarise` pattern from the `dplyr` package (part of tidyverse!). This does kinda what you'd expect: it first groups the data by a variable (in this case we would choose experiment so that each row from the same experiment is in the same group) and then we can use `summarise` to compute summaries of the groups.

```{r}
sampling_distribution <- draws_combined %>% 
         group_by(experiment) %>% 
         summarise(sample_mean = mean(draw))
```

We read this as:

- Take the tibble `draws_combined` and 
- Group it by experiment
- Make a new variable called `sample_mean` that is the sample mean of each experiment.



### Questions

3. Run the experiment 100,000 times and make a histogram of the sampling distribution^[The technical term of the distribution of a statistic (like the mean) under repeated sampling is the _sampling distribution_.] of the sample mean.

```{r, echo = FALSE, fig.show='hide'}
set.seed(30127)

N <- 100000
draws_combined <- tibble(
  experiment = rep(1:N, each = 100),
  draw = rnorm(N * 100, 3, 1)
)

sampling_distribution <- draws_combined %>% 
         group_by(experiment) %>% 
         summarise(sample_mean = mean(draw))

sampling_distribution %>% ggplot(aes(sample_mean)) +
  geom_histogram(fill = "purple", bins = 30) + theme_bw() +
  xlab("Sample mean")

```

4. It turns out that you can compute the sampling distribution of the sample mean exactly: it is normally distributed with mean `3` and variance `1/100`. Overlay the empirical histogram with the theoretical density. (See last week's bonus video for how to do this!)

```{r, echo = FALSE, fig.show='hide'}
sampling_distribution %>% ggplot(aes(sample_mean)) +
  geom_histogram(aes(y = after_stat(density)), 
                 fill = "purple",
                 bins = 50) +
  stat_function(fun = function(x) dnorm(x, 3,sqrt(1/100)), 
                colour = "red", linetype = "dashed", size = 1.2) +
  theme_bw() +
  xlab("Sample mean")

```


5. For the same experiment, plot the sampling distribution of the standard deviation.

```{r, echo = FALSE, fig.show='hide'}
sampling_distribution <- draws_combined %>% 
         group_by(experiment) %>% 
         summarise(sample_mean = sd(draw))

sampling_distribution %>% ggplot(aes(sample_mean)) +
  geom_histogram(fill = "purple", bins = 30) + theme_bw() +
  xlab("Sample standard deviation")

```

6. For the same experiment, plot the sampling distribution of the largest observation in the sample.

```{r, echo = FALSE, fig.show='hide'}
sampling_distribution <- draws_combined %>% 
         group_by(experiment) %>% 
         summarise(sample_mean = max(draw))

sampling_distribution %>% ggplot(aes(sample_mean)) +
  geom_histogram(fill = "purple", bins = 30) + theme_bw() +
  xlab("Sample standard deviation")

```


# Using the pattern to visualise more complex uncertainty

Sometimes, we want to know if a sample has a particular distribution. An example of this would be when we do linear regression and we want to check if the residuals are normally distributed. 

While there are some formal statistical tests for doing this, the most useful techniques are graphical. So we are going to look at one of those today and show how we can extend the idea of the sampling distribution that we developed in the previous example to graphs.

## Scenario: We have a sample that is supposed to be normal.

```{r}
data <- rt(20, df = 10)
```

We have some data that we have observed. (In this case, it's from a T distribution with 10 degrees of freedom, but pretend we don't know that.) We were told that this data should be $N(0,1)$. How can we check?

One choice would be to plot a histogram and compare it against the theoretical density, but it turns out to be more useful to compare the empircal cumulative density function (eCDF) of the sample to the theoretical CDF of the standard normal distribution^[This comparision is conceptually similar to the qq-plot that you may have seen in a different class.]

```{r}
tibble(sample = data) %>% ggplot(aes(x = sample)) +
  stat_ecdf() + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed")
```

This is all well and good, but is the empirical CDF close enough to the theoretical one to say that it's from the same distribution? To answer that, we need to assess how variable and emiprical CDF can be when it's computed based on 20 samples.

## Repeating the process

So. We need to assess the variability in an empirical CDF based on 20 samples. How do we do this? We follow the pattern above and make _multiple_ data sets with 20 points that come from a $N(0,1)$ distribution! And then we plot their empirical CDFs!



```{r}
set.seed(30127)
N <- 100
experiments <- tibble(experiment = rep(1:N, each = 20),
                      draw = rnorm(N * 20, 0, 1))

ggplot(experiments, aes(x = draw)) + stat_ecdf() + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed")

```

Oh no! That does not look right! Or, to be more precise, that _looks too good_! Nothing on that plot looks like the empircal CDF from 20 samples. So we've done something wrong.

What we accidentally did is plot the empircal CDF of 2000 (20 x 100) samples! This should be an identical graph:

```{r}
set.seed(30127)
tibble(draw = rnorm(2000)) %>%
  ggplot(aes(x = draw)) + stat_ecdf() + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed")

```

(Those two graphs really are the same because we used the same seed and generated the same number of normals!)

So what did we forget to do? We forgot to _group_ by experiment. Previously, we did this with the `group_by`/`summarise` pattern.  For plots, we instead add a `group` aesthetic.


```{r}
N <- 100
experiments <- tibble(experiment = rep(1:N, each = 20),
                      draw = rnorm(N * 20, 0, 1))

ggplot(experiments, aes(x = draw)) + 
  stat_ecdf(aes(group = experiment)) + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed")

```

Here we only needed `group = experiment` for the `stat_ecdf` layer, so we put it in that part of the ggplot call. It would still work in the global aesthetic but it would have thrown a warning because the `group` aesthetic doesn't make sense in `stat_function()`.

So this is ... ok. But it's very ugly. The problem is that we have plotted 101 lines on top of each other on this plot and it is visually very messy. We can fix this with two tricks: colour and transparancy.

```{r}
ggplot(experiments, aes(x = draw)) + 
  stat_ecdf(aes(group = experiment),
            colour = "lightgrey",
            alpha = 10/100) + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed") + 
  theme_classic() + xlab("x")
```

This looks a lot better:

- Removing the grey background and the grid with `theme_classic()` reduces the visual clutter.
- Making the empirical CDFs a lighter colour is useful to make them less prominent.
- The `alpha` parameter makes the lines semi-transparent. This helps reduce their intensity. This parameter is difficult to set and you should try a bunch of different values. (`alpha = 1` is the default value, and `alpha = 0` is fully transparent.)

Why try to de-emphasise the empircal CDFs? Because we don't care aboout the individual lines. We only care about their variability! Because this gives us a visual _window_ where the deviation from the theoretical CDF is purely due to the small sample size. If the empirical CDF of the sample goes outside this window, we are seeing something qualitatively different from our 100 draws.

Finally, let's add in the sample.

To do this, we need to add a new layer that is an empirical CDF of the actual data we care about. This layer needs to have _different data_ to the rest of the plot, and we can do that by adding a `data = ...` argument
to the `stat_ecdf()` layer. We also need to define an aesthetic in this layer so that ggplot knows what to put on the x-axis.

As always, ggplot build the plot in the order it's specified. So it builds the basic axes, then it adds the semi-transparent theoretical empirical CDFs, then it adds the data empirical CDF, and finally it adds the theoretical CDF.

```{r}
ggplot(experiments, aes(x = draw)) + 
  stat_ecdf(aes(group = experiment),
            colour = "lightgrey",
            alpha = 10/100) + 
  stat_ecdf(data = tibble(real_data = data), 
            aes(x = real_data)) + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed") + 
  theme_classic() + 
  xlab("x")
```

In this case, the data empirical CDF does not seem any stranger than the empirical CDFs build from 20 samples from the true distribution. So we cannot say from this plot that the data does not come from a N(0,1) distribution. 

But there is an _important thing to notice_: We also haven't proved that the does come from a N(0,1) distribution! In fact, we know that it does not! What we have shown is that _we cannot distinguish between the sample and 20 draws from a standard normal distribution_.

### Questions

7. Can you tell the difference when the experiment has 2000 observations instead of 20. (This will take a little longer to run!)

```{r, eval = FALSE, echo = FALSE}
set.seed(30127)
N <- 100
data <- rt(2000, df = 10)
experiments <- tibble(experiment = rep(1:N, each = 2000),
                      draw = rnorm(N * 2000, 0, 1))
ggplot(experiments, aes(x = draw)) + 
  stat_ecdf(aes(group = experiment),
            colour = "lightgrey",
            alpha = 10/100) + 
  stat_ecdf(data = tibble(real_data = data), 
            aes(x = real_data)) + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed") + 
  theme_classic() + 
  xlab("x")
```

8. Can you tell the difference between 20 samples from a T-distribution with 5 degrees of freedom and a N(0,1) distribution?

```{r, eval = FALSE, echo = FALSE}
set.seed(30127)
N <- 100
data <- rt(20, df = 5)
experiments <- tibble(experiment = rep(1:N, each = 20),
                      draw = rnorm(N * 20, 0, 1))
ggplot(experiments, aes(x = draw)) + 
  stat_ecdf(aes(group = experiment),
            colour = "lightgrey",
            alpha = 10/100) + 
  stat_ecdf(data = tibble(real_data = data), 
            aes(x = real_data)) + 
  stat_function(fun = pnorm, colour = "blue",
                linetype = "dashed") + 
  theme_classic() + 
  xlab("x")
```