TeachingSimulationsBasics.Rmd

---
title: "TeachingSimulationsBasics"
author: "Ian Dworkin"
date: "10/01/2022"
output: 
  html_document: 
    toc: yes
    number_sections: yes
    keep_md: yes
editor_options: 
  chunk_output_type: console
---

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

# The basics of Simulations in R

##  What do I mean by simulate? Why Monte Carlo?

## Our first example

Let us simulate 100 observations from a normally distributed variable with a mean of 5 and a standard deviation of 2  $\sim N(\mu = 5, \sigma = 2)$


We will use the `rnorm` function which generates **r**andom values from a **norm**al distribution

```{r}
norm_simulated <- rnorm(n = 100, mean  = 5, sd = 2)

par(mfrow=c(1,1))

plot(norm_simulated, pch = 20)

hist(norm_simulated, freq = F)

curve(dnorm(x, mean = 5, sd = 2), from = 0, to = 10, 
      add = T, col = "red")
```


## look at summary statistics (mean and sd) from this.

```{r}
mean(norm_simulated)
sd(norm_simulated)
```


## repeat this a few times. what do you observe?


```{r}
norm_simulated1 <- rnorm(n = 100, mean  = 5, sd = 2)
norm_simulated2 <- rnorm(n = 100, mean  = 5, sd = 2)
norm_simulated3 <- rnorm(n = 100, mean  = 5, sd = 2)
norm_simulated4 <- rnorm(n = 100, mean  = 5, sd = 2)

par(mfrow = c(2,2))
hist(norm_simulated1)
hist(norm_simulated2)
hist(norm_simulated3)
hist(norm_simulated4)
par(mfrow = c(1,1))

mean(norm_simulated1)
mean(norm_simulated2)
mean(norm_simulated3)
mean(norm_simulated4)
```

So for each (iteration) of the simulation, we have effectively sampled 100 observations from a population with a normal distribution with  $\sim N(\mu = 5, \sigma = 2)$. If we wanted to change the number of observations per iteration we simply change n. Similarly we could change the values for mean (or make it some function like a line), the standard deviation and the type of probability distribution we are sampling from.

### Please simulate 30 observations from 
$\sim N(\mu = 10, \sigma = 3)$ 5 times.


## more efficient performing the simulation 

Of course performing this one iteration at a time is not particularly useful, so instead we find a way to repeat it efficiently.

 One way to do this in R is to use the `replicate function`, which just repeats the function call n times. Notice that there is the n=4 in the replicate which refers to how many times to repeat the call to the function (in this case the rnorm). the n=100 in rnorm is how many samples are drawn from the population (distribution).
 
 
```{r}
norm_sim_all_2 <- replicate(n = 4, 
                            rnorm(n = 100, mean = 5, sd = 2)) 

par(mfrow = c(2,2))

apply(X = norm_sim_all_2, MARGIN = 2, FUN = hist)
apply(X = norm_sim_all_2, MARGIN = 2, FUN = mean)
apply(X = norm_sim_all_2, MARGIN = 2, FUN = sd)

par(mfrow=c(1,1))

```
 

let's do this again, but repeat the sampling process (number of iterations) 2000 times. $\sim N(\mu = 7, \sigma = 1.8)$

```{r}
norm_sim_all_3 <- replicate(n=2000, 
                            rnorm(n = 100, mean = 7, sd = 1.8)) 
```
 
We can ask about the spread of the means from our simulated trials. How might you code this?

```{r}
hist(apply(X = norm_sim_all_3, MARGIN = 2, FUN = mean),
     main = "distribution of means from simulation",
     xlab = "mean values")


# standard deviation among sampled means (the standard error of the mean)
sd(apply(X = norm_sim_all_3, MARGIN = 2, FUN = mean))
```


Let's repeat this experiment with lower sample size (say 25)

```{r}
norm_sim_all_4 <- replicate(n = 2000, 
                            rnorm(n = 25, mean = 5,sd = 2))

# standard deviation among sampled means (the standard error of the mean)
sd(apply(norm_sim_all_4,2, mean))

hist(apply(norm_sim_all_4,2,mean),
     main = "distribution of means from simulation",
     xlab = "mean values")
```

How about if we have much larger sample size (say 1000)?

```{r}
norm_sim_all_5 <- replicate(n = 2000,
                            rnorm(1000, 5, 2))


sd(apply(norm_sim_all_5, 2, mean))

hist(apply(norm_sim_all_5, 2, mean),
     main = "distribution of means from simulation",
     xlab = "mean values")
```


# Do you see a pattern for the precision of the mean across repeated simulated samples?

```{r}

sd(apply(norm_sim_all_4, 2, mean)) # n=25

sd(apply(norm_sim_all_2, 2, mean)) # n=100

sd(apply(norm_sim_all_5, 2, mean)) # n=1000
```


Notice that as sample size (for a given simulated sample) gets larger, the variation (sd among the means) among the estimated means gets smaller!

We have just learned about monte carlo methods. When we perform simulations from that distributions, as our n (sample size) increase we will get closer to the expected value (which in this case we provided). That is the law of large numbers.
 
**This hopefully reinforces why the sample size influences the estimated the standard error so clearly!**

## Using for loops instead  

We can also run such simulations using for loops, which can be more efficient when the number of simulations (or parameters being estimated) is large

# How many iterations of the simulation do we want to do? We will call this "N"

```{r}
N <- 2000
```

 First we initialize a variable to store the mean values we are going to estimate from each iteration of the simulation
 
```{r}
simulated_means <- rep(NA, N) # Create an empty vector of length N

head(simulated_means)  # Filled with "missing data" until we fill it
```
 

Now we use the for loop to repeat this from i=1 until i=N (which is 2000 in this case )


```{r}

for (i in 1:N) {
	sim_data <- rnorm(n = 100, mean  = 5, sd = 2)
	simulated_means[i] <- mean(sim_data)
	rm(sim_data) # While this gets over-written each iteration of the loop, 
	   #we want to remove it after the last iteration.
}
```

 the vector "simulated_means" now contains all of the means from the N iterations of the simulation.
 
```{r}
hist(simulated_means)
sd(simulated_means)
```




## Simulations for a simple linear regression

Instead of simulating just a mean and sd, let's specify something more interesting... How about a line/regression!!!!!

let's have a regression $Y \sim N(\mu = a + b*x, \sigma = sd)$

Specifically with an intercept of 5 (a = 5) and a slope of 0.7 (b = 0.7).

$$Y \sim N(\mu = 5 + 7*x, \sigma = sd)$$

```{r}
a = 5
b = 0.7
x <- seq(2, 20) # our predictor variable
y_fixed <- a + b*x  # the deterministic part of our model

par(mfrow = c(2, 1))
plot(y_fixed ~ x, pch = 20, col = "red")  # The deterministic relationship.
```

# But our response is randomly sampled (conditioned on x), so we need to be able to incorporate this random variation. So let's say sd=2.5

$$Y \sim N(\mu = 5 + 7*x, \sigma = 2.5)$$

```{r}

y.sim.1 <- rnorm(length(x), mean = y_fixed, sd = 2.5) 
```

`length` just looks at the number of elements in the vector. So we have now simulated values of our response variable `y` given a deterministic and stochastic component to the model!


```{r}
plot(y.sim.1 ~ x, pch = 20, col = "red")
abline(a = 5, b = 0.7, lwd = 2) # Expected relationship based on the "known" parameters we used.
```

This is the expected relationship based on the "known" parameters we used.

But what are the actual parameter estimates for the regression based on the simulated data (for y)? We can use the `lm` function to run a linear regression.

```{r}
y.sim.1.lm <- lm(y.sim.1 ~ x)

summary(y.sim.1.lm)  # notice parameter estimates and RSE!

confint(y.sim.1.lm) # does it include our expected (TRUE) values?

plot(y.sim.1 ~ x, pch = 20, col = "red")

abline(a = 5, b = 0.7, lwd = 2) 
abline(reg=y.sim.1.lm, lty = 2, col = "purple", lwd = 2) # estimated values based on simulated data.
```

Not identical, but similar to the "true" (which is known because we chose it) relationship. The point is, we have now done a little simulation of our regression model! We have the R tools we need to really begin learning statistics and about uncertainty!