06_standardError.Rmd

---
title: "Intro to Inference"
output:
  bookdown::html_document2:
    includes:
      in_header: assets/06_standardError_image.html
      after_body: assets/foot.html
---

```{r setup, include = FALSE}
library(ggplot2)
library(tweenr)
library(parallel)
library(MASS)
```


:::obj
**Module learning objectives**

1. Determine how to quantify the uncertainty of an estimate
1. Describe the concept of statistical inference 
1. Interpret sampling distributions and explain how they are influenced by sample size
1. Define and calculate standard error
1. Use the standard error to construct 95% confidence intervals
:::

# How accurate is our estimate of the mean?

Let's revisit the first few days during which we collected data stored in the vector `heights_island1`. We were able to verify that the heights were normally distributed and calculated our sample mean, ${\bar{x}}$. However, we know that ${\bar{x}}$ is only an *estimate* of the true population mean, ${\mu}$, which is the true value of interest. It is unlikely that we will ever know the value of ${\mu}$, since access to all possible observations is rare. Therefore we will have to rely on ${\bar{x}}$ estimates from random samples drawn from the population as the best approximation of ${\mu}$.

Not all sample means are created equal. Some are better estimates than others. Recall the [animation](03_mean.html#mean_animation) showing the relationship between sample size and variability of the mean. As we learned from this animation, in the long-run, large samples are necessary to get an accurate estimate of ${\mu}$.


<div class= "alert alert-note">
> **A note about language:** here, words like "accuracy", "precision", and "uncertainty" are used in a rather fast and loose way. We're using the laymen's application of these terms to refer to the long-run variability of estimates produced from repeated, independent trials. There are stricter, more formal statistical uses for these words, but for right now, we're going to ignore these nuances so that we can move on with understanding these concepts in broad strokes.

</div>

One reason we care about our sample estimate's accuracy is because we want to be able to answer questions about the population by making inferences. **Statistical inference** uses math to draw conclusions about the population based on a subset of the full picture (i.e. a sample). Subsets of data are of course limited, so it's therefore important to acknowledge that the strength of the conclusions drawn about the population is dependent on the precision of the sample estimate. For example, say that we guess that the population mean value of giraffe heights on Island 1 is less than 11 cm. We can make some inferences about whether or not this is a good guess based on what we learn from our sample of giraffe heights. We'll revisit this question a few times below. 

# Creating a sampling distribution

The mean of our sample of 50 giraffes from Island 1 was:

```{r, echo=FALSE}
set.seed(12)
heights_island1 <- rnorm(50, 10, 2)
``` 

```{r}
mean(heights_island1)
``` 

How can we quantify the accuracy of this estimate, given its sample size?

In theory, one way to illustrate this is to generate data not just from a single sample but from many samples of the same size (N) drawn from the same population. 

Imagine that after you collected all 50 measurements for `heights_island1`, you wake up one morning with no memory of collecting data at all---and so you go out and collect 50 giraffe heights again and subsequently calculate the mean. Further imagine that this groundhog day (or more correctly, groundhog *week*) situation repeats itself many, many times.

When you finally return to your sanity, you find stacks of notebooks filled with mean values from each of your individual data collections. 

<center>![](images/06_standardError/Notebooks.jpg){width=600px}</center>

Instead of viewing this as a massive waste of time, you make the best out of the situation and create a histogram of all the means. In other words you create a plot showing the distribution of the sample means, also known as a **sampling distribution**. 

The animation below illustrates the process of creating the sampling distribution for 1,000 sample means.

On the left side, each histogram represents a sample (e.g. `heights_island1` would be one sample, and we're flashing through 1,000 of them in total). Correspondingly, each dot signifies an observation. After each sample histogram is completed, ${\bar{x}}$ is calculated. This ${\bar{x}}$ value is then subsequently added to the histogram of the sampling distribution on the right. As you can see below, this process is repeated, allowing the sampling distribution to build up.

<center>
```{r fig.show="animate", animation.hook = 'gifski', fig.width=7, fig.height=3, echo=FALSE, message=FALSE, warning=FALSE, results = 'hide', interval=0.08, loop=FALSE, cache=TRUE}

ppplot <- function(sub){
x <- round(rnorm(50, 9.7, 2.1))
m <- mean(x)
df <- data.frame(x = x, y = 23)
dfs <- list(df)
for(i in seq_len(nrow(df))) {
  dftemp <- tail(dfs, 1)
  dftemp[[1]]$y[i] <- sum(dftemp[[1]]$x[seq_len(i)] == dftemp[[1]]$x[i])
  dfs <- append(dfs, dftemp)
}
dfs <- append(dfs, dfs[rep(length(dfs), 3)])
dft <- tween_states(dfs, 10, 1, 'cubic-in', 50)
dft$y <- dft$y - 0.5
dft <- dft[dft$y != 23, ]

ppl <- function(frame){
  p <- ggplot(data = dft[dft$.frame==frame,], aes(x=x, y=y)) + 
    geom_point(shape=16, color="green3", size = 4) + 
    ylim(0,16) + xlim(3,17) + 
    theme_light() + 
    theme(panel.border = element_blank(), panel.grid.minor=element_blank()) + 
    labs(x="Giraffe Heights", y=NULL)
  df <- data.frame()
  p2 <- ggplot(df) + geom_point() + xlim(0, 16) + ylim(3, 17)+theme_void()
  p3 <- ggplot(df) + geom_point() + xlim(8.7, 10.7) + ylim(0, 150) + 
    theme_light() + 
    theme(panel.border=element_blank(), panel.grid.minor=element_blank()) + labs(x = "Sample means", y = NULL)
  cowplot::plot_grid(p, p2, p3, align ="h", rel_widths = c(1, 0.55, 1), ncol = 3)
  
}

ppl2 <- function(frame){
  p <- ggplot(data=dft[dft$.frame==53,], aes(x=x, y=y)) + 
    geom_point(shape = 16, color = "green3", size = 4) + 
    ylim(0,16) + xlim(3,17) + 
    theme_light() + 
    theme(panel.border = element_blank(), panel.grid.minor = element_blank()) + 
    geom_vline(xintercept = m, linetype = 2) + 
    labs(x = "Giraffe Heights", y = NULL)
  p
  df <- data.frame()
  lb1 <- paste0("bar(x)", "[", sub, "]", " == ", round(m,2))
  p2 <- ggplot(df) + 
    geom_point() + 
    xlim(0, 16) + ylim(3, 17) + 
    theme_void() + 
    annotate("text", x = 8, y=10, label = lb1, parse = TRUE, size = 7) + 
    annotate("segment", x = 1, xend = 15, y = 8, yend = 8, colour = "black", size = 1, arrow = arrow(type = "closed", length = unit(0.3,"cm")))
  p3 <- ggplot(df) + geom_point() + xlim(8.7, 10.7) + ylim(0, 150)+theme_light()+theme(panel.border=element_blank(), panel.grid.minor=element_blank())+annotate("segment", x = m, xend = m, y = 20, yend = 4, colour = "black", size=1, arrow=arrow(type = "closed", length = unit(0.3,"cm")))+labs(x="Sample means", y=NULL)
  cowplot::plot_grid(p,p2,p3, align="h", rel_widths = c(1,0.55,1), ncol = 3)
}

pf <- list(lapply(seq(1, 53, 2), function(x) ppl(x)), lapply(rep(53, 3), function(x) ppl(x)), lapply(1:40, function(x) ppl2()))
return(pf)
}
mclapply(1:3, function(x) ppplot(x), mc.cores = 8, mc.cleanup = TRUE)

circleFun <- function(center=c(0,0), diameter=1, npoints=100, start=0, end=2, filled=TRUE){
  tt <- seq(start*pi, end*pi, length.out=npoints)
  df <- data.frame(
    x = center[1] + diameter / 2 * cos(tt),
    y = center[2] + diameter / 2 * sin(tt)
  )
  if(filled==TRUE) { 
    df <- rbind(df, center)
  }
  return(df)
}
fullCircle <- circleFun(c(1, -1), 2.3, start=0, end=2, filled=FALSE)
fullCircle2 <- circleFun(c(1, -1), 2, start=0, end=2, filled=FALSE)
fullCircle3 <- circleFun(c(1, -1), 1.3, start=0, end=2, filled=FALSE)
fullCircle4 <- circleFun(c(1, -1), 0.3, start=0, end=2, filled=FALSE)
fullCircle5 <- circleFun(c(1, -1), 0.1, start=0, end=2, filled=FALSE)

tris <- circleFun(c(1, -1), 1.6, start=1.2, end=-0.2, filled=FALSE, npoints=50)
tris2 <- circleFun(c(1, -1), 0.2, start=1.4, end=0, filled=FALSE, npoints=50)
tris3 <- circleFun(c(1, -1), 0.2, start=1, end=-0.4, filled=FALSE,npoints=50)

s <- c(rep(1,10), 1:50)

trii <- lapply(s, function(x) data.frame(x=c(tris[x,1],tris2[x,1],tris3[x,1]), y=c(tris[x,2],tris2[x,2],tris3[x,2])))

quarterCircle <- circleFun(c(1,-1), diameter = 1.85, start=1, end=1.25, filled=TRUE)
quarterCircle2 <- circleFun(c(1,-1), diameter = 1.85, start=0.75, end=1, filled=TRUE)
quarterCircle3 <- circleFun(c(1,-1), diameter = 1.85, start=0.5, end=0.75, filled=TRUE)
quarterCircle4 <- circleFun(c(1,-1), diameter = 1.85, start=0.25, end=0.5, filled=TRUE)
quarterCircle5 <- circleFun(c(1,-1), diameter = 1.85, start=0.25, end=0, filled=TRUE)
quarterCircle6 <- circleFun(c(1,-1), diameter = 1.85, start=2, end=1.75, filled=TRUE)

x <- round(rnorm(50, 9.7, 2.1))
m <- mean(x)
df <- data.frame(x = x, y = 23)
dfs <- list(df)
for(i in seq_len(nrow(df))) {
  dftemp <- tail(dfs, 1)
  dftemp[[1]]$y[i] <- sum(dftemp[[1]]$x[seq_len(i)] == dftemp[[1]]$x[i])
  dfs <- append(dfs, dftemp)
}
dfs <- append(dfs, dfs[rep(length(dfs), 3)])
dft <- tween_states(dfs, 10, 1, 'cubic-in', 50)
dft$y <- dft$y - 0.5
dft <- dft[dft$y != 23, ]
dft <- dft[dft$.frame %in% c(1:26, seq(27, 53, 2)),]
dft$.frame <- rep(1:40, each=50)

plots <- function(dd){
p <- ggplot() + 
  geom_polygon(data=fullCircle, aes(x, y), color="#40596b", fill="#40596b") +
  geom_polygon(data=fullCircle2, aes(x, y), color="white", fill="white") +
  geom_polygon(data=quarterCircle, aes(x,y), color="#cdd6e0", fill="#cdd6e0") + 
  geom_polygon(data=quarterCircle2, aes(x,y), color="#acb3ba", fill="#acb3ba") + 
  geom_polygon(data=quarterCircle3, aes(x,y), color="#ffd15c", fill="#ffd15c") +
  geom_polygon(data=quarterCircle4, aes(x,y), color="#f8b64c", fill="#f8b64c") +
  geom_polygon(data=quarterCircle5, aes(x,y), color="#ff7058", fill="#ff7058") +
  geom_polygon(data=quarterCircle6, aes(x,y), color="#f1543f", fill="#f1543f") +
  geom_polygon(data=fullCircle3, aes(x,y), color="white", fill="white") +
  geom_polygon(data=fullCircle4, aes(x,y), color="#40596b", fill="#40596b") +
  geom_polygon(data=trii[[dd]], aes(x,y), color="#40596b", fill="#40596b") +
  geom_polygon(data=fullCircle5, aes(x,y), color="white", fill="white") +
  coord_equal() +
  theme_void()
  ddd <- ifelse(dd<20, 1, ifelse(dd<35, 2, ifelse(dd<45, 3, ifelse(dd<50, 4, ifelse(dd<53, 5, ifelse(dd<55, 6, ifelse(dd<57, 7, ifelse(dd<59, 8, base::sample(1:10,1)))))))))
  p2 <- ggplot(data=dft[dft$.frame==dd,],aes(x=x, y=y))+geom_point(shape=16, color="green3", size=4)+ylim(0,16)+xlim(3,17)+theme_light()+theme(panel.border=element_blank(), panel.grid.minor=element_blank())+labs(x="Giraffe Heights", y=NULL)
  df <- data.frame()
  p3 <- ggplot(df) + geom_point() + xlim(8.7, 10.7) + ylim(0, 150)+theme_light()+theme(panel.border=element_blank(), panel.grid.minor=element_blank())+labs(x="Sample means", y=NULL)
  cowplot::plot_grid(p2,p,p3, align="h",rel_widths = c(1,0.55, 1), ncol=3)
}

lapply(seq(1,40,3), function(x) plots(x))

hists <- function(x){
  x <- round(rnorm(50, 9.7, 2.1))
  m <- mean(x)
  return(m)
}  

dh <- do.call(rbind, lapply(1:1000, function(x) hists()))

 hh <- function(x){
  d <- data.frame(Height=dh[1:x])
  return(d)
}

dhh <<- lapply(1:1000, function(x) hh(x))

plots2 <- function(dd){
p <- ggplot() + 
  geom_polygon(data=fullCircle, aes(x, y), color="#40596b", fill="#40596b") +
  geom_polygon(data=fullCircle2, aes(x, y), color="white", fill="white") +
  geom_polygon(data=quarterCircle, aes(x,y), color="#cdd6e0", fill="#cdd6e0") + 
  geom_polygon(data=quarterCircle2, aes(x,y), color="#acb3ba", fill="#acb3ba") + 
  geom_polygon(data=quarterCircle3, aes(x,y), color="#ffd15c", fill="#ffd15c") +
  geom_polygon(data=quarterCircle4, aes(x,y), color="#f8b64c", fill="#f8b64c") +
  geom_polygon(data=quarterCircle5, aes(x,y), color="#ff7058", fill="#ff7058") +
  geom_polygon(data=quarterCircle6, aes(x,y), color="#f1543f", fill="#f1543f") +
  geom_polygon(data=fullCircle3, aes(x,y), color="white", fill="white") +
  geom_polygon(data=fullCircle4, aes(x,y), color="#40596b", fill="#40596b") +
  geom_polygon(data=trii[[dd]], aes(x,y), color="#40596b", fill="#40596b") +
  geom_polygon(data=fullCircle5, aes(x,y), color="white", fill="white") +
  coord_equal() +
  theme_void()
  ddd <- ifelse(dd<20, 1, ifelse(dd<35, 2, ifelse(dd<45, 3, ifelse(dd<50, 4, ifelse(dd<53, 5, ifelse(dd<55, 6, ifelse(dd<57, 7, ifelse(dd<59, 8, base::sample(1:10,1)))))))))
  set.seed(ddd)
  x <- round(rnorm(50, 9.7, 2.1))
m <- mean(x)
df <- data.frame(x = x, y = 23)
dfs <- list(df)
for(i in seq_len(nrow(df))) {
  dftemp <- tail(dfs, 1)
  dftemp[[1]]$y[i] <- sum(dftemp[[1]]$x[seq_len(i)] == dftemp[[1]]$x[i])
  dfs <- append(dfs, dftemp)
}
dfs <- append(dfs, dfs[rep(length(dfs), 3)])
dft <- tween_states(dfs, 10, 1, 'cubic-in', 50)
dft$y <- dft$y - 0.5
dft <- dft[dft$y != 23, ]
  p2 <- ggplot(data=dft[dft$.frame==53,],aes(x=x, y=y))+geom_point(shape=16, color="green3", size=4)+ylim(0,16)+xlim(3,17)+theme_light()+theme(panel.border=element_blank(), panel.grid.minor=element_blank())+geom_vline(xintercept = m, linetype=2)+labs(x="Giraffe Heights", y=NULL)
  df <- data.frame()

p3 <- ggplot(data = dhh[[dd-40]], aes(x = Height)) +
  geom_histogram(binwidth = 0.1, color = "white", fill="green3") +
  theme_light() +
  scale_y_continuous(limits = c(0,150)) +
  labs(x=NULL, y=NULL) +
  xlim(8.7, 10.7) + 
  theme(panel.border=element_blank(), panel.grid.minor=element_blank()) +
  labs(x="Sample means", y=NULL)
  cowplot::plot_grid(p2,p,p3, align="h",rel_widths = c(1,0.55, 1), ncol=3)
}
lapply(seq(41,50,2), function(x) plots2(x))
lapply(seq(51,60,1), function(x) plots2(x))


plots3 <- function(dd){
  
p <- ggplot() + 
  geom_polygon(data=fullCircle, aes(x, y), color="#40596b", fill="#40596b") +
  geom_polygon(data=fullCircle2, aes(x, y), color="white", fill="white") +
  geom_polygon(data=quarterCircle, aes(x,y), color="#cdd6e0", fill="#cdd6e0") + 
  geom_polygon(data=quarterCircle2, aes(x,y), color="#acb3ba", fill="#acb3ba") + 
  geom_polygon(data=quarterCircle3, aes(x,y), color="#ffd15c", fill="#ffd15c") +
  geom_polygon(data=quarterCircle4, aes(x,y), color="#f8b64c", fill="#f8b64c") +
  geom_polygon(data=quarterCircle5, aes(x,y), color="#ff7058", fill="#ff7058") +
  geom_polygon(data=quarterCircle6, aes(x,y), color="#f1543f", fill="#f1543f") +
  geom_polygon(data=fullCircle3, aes(x,y), color="white", fill="white") +
  geom_polygon(data=fullCircle4, aes(x,y), color="#40596b", fill="#40596b") +
  geom_polygon(data=trii[[60]], aes(x,y), color="#40596b", fill="#40596b") +
  geom_polygon(data=fullCircle5, aes(x,y), color="white", fill="white") +
  coord_equal() +
  theme_void()

x <- round(rnorm(50, 9.7, 2.1))
m <- mean(x)
df <- data.frame(x = x, y = 23)
dfs <- list(df)
for(i in seq_len(nrow(df))) {
  dftemp <- tail(dfs, 1)
  dftemp[[1]]$y[i] <- sum(dftemp[[1]]$x[seq_len(i)] == dftemp[[1]]$x[i])
  dfs <- append(dfs, dftemp)
}
dfs <- append(dfs, dfs[rep(length(dfs), 3)])
dft <- tween_states(dfs, 10, 1, 'cubic-in', 50)
dft$y <- dft$y - 0.5
dft <- dft[dft$y != 23, ]

  p2 <- ggplot(data=dft[dft$.frame==53,],aes(x=x, y=y))+geom_point(shape=16, color="green3", size=4)+ylim(0,16)+xlim(3,17)+theme_light()+theme(panel.border=element_blank(), panel.grid.minor=element_blank())+geom_vline(xintercept = m, linetype=2)+labs(x="Giraffe Heights", y=NULL)
  
hh <- function(x){
  d <- data.frame(Height=dh[1:x])
  return(d)
}

dhh <- lapply(1:1000, function(x) hh(x))
p3 <- ggplot(data = dhh[[dd+40]], aes(x = Height)) +
  geom_histogram(binwidth = 0.1, color = "white", fill="green3") +
  theme_light() +
  scale_y_continuous(limits = c(0,150)) +
  labs(x=NULL, y=NULL) +
  xlim(8.7, 10.7) + 
  theme(panel.border=element_blank(), panel.grid.minor=element_blank()) +
  labs(x="Sample means", y=NULL)

  cowplot::plot_grid(p2,p,p3, align="h",rel_widths = c(1,0.55, 1), ncol = 3)
}
mclapply(seq(1,300, 20), function(x) plots3(x), mc.cores = 8, mc.cleanup = TRUE)

plots3.2 <- function(dd){
sub <- dd+40  

x <- round(rnorm(50, 9.7, 2.1))
m <- mean(x)
df <- data.frame(x = x, y = 23)
dfs <- list(df)
for(i in seq_len(nrow(df))) {
  dftemp <- tail(dfs, 1)
  dftemp[[1]]$y[i] <- sum(dftemp[[1]]$x[seq_len(i)] == dftemp[[1]]$x[i])
  dfs <- append(dfs, dftemp)
}
dfs <- append(dfs, dfs[rep(length(dfs), 3)])
dft <- tween_states(dfs, 10, 1, 'cubic-in', 50)
dft$y <- dft$y - 0.5
dft <- dft[dft$y != 23, ]

df <- data.frame()
  lb1 <- paste0("bar(x)", "[", sub, "]", " == ", round(m,1))
  p <- ggplot(df) + geom_point() + xlim(0, 16) + ylim(3, 17)+theme_void()+annotate("text", x = 1, y=10, label=lb1, parse = TRUE, size=7,hjust = 0)+annotate("segment", x = 1, xend = 15, y = 8, yend = 8, colour = "black", size=1, arrow=arrow(type = "closed", length = unit(0.3,"cm")))

  p2 <- ggplot(data=dft[dft$.frame==53,],aes(x=x, y=y))+geom_point(shape=16, color="green3", size=4)+ylim(0,16)+xlim(3,17)+theme_light()+theme(panel.border=element_blank(), panel.grid.minor=element_blank())+geom_vline(xintercept = m, linetype=2)+labs(x="Giraffe Heights", y=NULL)
  
hh <- function(x){
  d <- data.frame(Height=dh[1:x])
  return(d)
}

dhh <- lapply(1:1000, function(x) hh(x))
p3 <- ggplot(data = dhh[[dd+40]], aes(x = Height)) +
  geom_histogram(binwidth = 0.1, color = "white", fill="green3") +
  theme_light() +
  scale_y_continuous(limits = c(0,150)) +
  labs(x=NULL, y=NULL) +
  xlim(8.7, 10.7) + 
  theme(panel.border=element_blank(), panel.grid.minor=element_blank()) +
  labs(x="Sample means", y=NULL)

  cowplot::plot_grid(p2,p,p3, align="h",rel_widths = c(1,0.55, 1), ncol = 3)
}
mclapply(seq(301,960, 20), function(x) plots3.2(x), mc.cores = 8, mc.cleanup = TRUE)

x <- round(mvrnorm(50, 9.8, 2.1^2, empirical = T))
m <- mean(x)
df <- data.frame(x = x, y = 23)
dfs <- list(df)
for(i in seq_len(nrow(df))) {
  dftemp <- tail(dfs, 1)
  dftemp[[1]]$y[i] <- sum(dftemp[[1]]$x[seq_len(i)] == dftemp[[1]]$x[i])
  dfs <- append(dfs, dftemp)
}
dfs <- append(dfs, dfs[rep(length(dfs), 3)])
dft <- tween_states(dfs, 10, 1, 'cubic-in', 50)
dft$y <- dft$y - 0.5
dft <- dft[dft$y != 23, ]

plots4 <- function(dd){
  
df <- data.frame()
  lb1 <- paste0("bar(x)", "[", 1000, "]", " == ", 9.8)
  p <- ggplot(df) + geom_point() + xlim(0, 16) + ylim(3, 17)+theme_void()+annotate("text", x = 1, y=10, label=lb1, parse = TRUE, size=7,hjust = 0)+annotate("segment", x = 1, xend = 15, y = 8, yend = 8, colour = "black", size=1, arrow=arrow(type = "closed", length = unit(0.3,"cm")))

  p2 <- ggplot(data=dft[dft$.frame==53,],aes(x=x, y=y))+geom_point(shape=16, color="green3", size=4)+ylim(0,16)+xlim(3,17)+theme_light()+theme(panel.border=element_blank(), panel.grid.minor=element_blank())+geom_vline(xintercept = 9.8, linetype=2)+labs(x="Giraffe Heights", y=NULL)
  
hh <- function(x){
  d <- data.frame(Height=dh[1:x])
  return(d)
}

dhh <<- lapply(1:1000, function(x) hh(x))
p3 <- ggplot(data = dhh[[dd+40]], aes(x = Height)) +
  geom_histogram(binwidth = 0.1, color = "white", fill="green3") +
  theme_light() +
  scale_y_continuous(limits = c(0,150)) +
  labs(x=NULL, y=NULL) +
  xlim(8.7, 10.7) + 
  theme(panel.border=element_blank(), panel.grid.minor=element_blank()) +
  labs(x="Sample means", y=NULL)

  cowplot::plot_grid(p2,p,p3, align="h",rel_widths = c(1,0.55, 1), ncol = 3)
}

mclapply(rep(960, 40), function(x) plots4(x), mc.cores = 8, mc.cleanup = TRUE)
```
</center>

<br>
<br>

A histogram of the sampling distribution is shown below. It is a histogram made up of many means.

<br>

<center>
```{r, tut=FALSE, echo=FALSE, message= FALSE, warning=FALSE, fig.height=2.7, fig.width=4.5, cache= TRUE}

samp <- function(n){
  x <- rnorm(n, 9.7, 2.1)
  m <- mean(x)
  s <- sd(x)
  return(c(m,s))
}

d2 <- as.data.frame(do.call(rbind, lapply(1:1000, function(x) samp(50))))
colnames(d2) <- c("mean", "sd")

ggplot(data = dhh[[1000]], aes(x = Height)) +
  geom_histogram(binwidth = 0.1, color = "white", fill = "green3") +
  theme_light() +
  scale_y_continuous(expand = c(0,0)) +
  labs(x = "Sample means", y = NULL) +
  theme(panel.border = element_blank(), panel.grid.minor = element_blank()) 
```
</center>

Looking at the spread of ${\bar{x}}$ values that this groundhog experience generated, we can get a sense of the range of many possible estimates of ${\mu}$ that a sample of 50 giraffes can produce. 

**The sampling distribution provides us with the first hint of the precision of our original `heights_island1` estimate**, which we'll quantify in more detail later on, but for now it's enough to notice that the range of possible ${\bar{x}}$ values are between `r round(min(d2$mean),1)` and `r round(max(d2$mean), 1)`. This means that ${\bar{x}}$ values outside of this range are essentially improbable.

Let's return to our question about whether the true mean of giraffe heights on Island 1 is less than 11 cm. Our sampling distribution suggests that ${\mu}$ *is* less than 11 cm, since values greater than that are not within the range of this sampling distribution. 


# Sample size and sampling distribution

Back to the idea that larger samples are "better", we can explore what happens if we redo the groundhog scenario, this time sampling 500 individuals (instead of 50) before taking the mean each time, repeating this until thousands of ${\bar{x}}$ values have been recorded. For completeness, let's imagine the same marathon data collection using samples that are smaller---of 5 giraffes each. We compare the resulting sampling distributions from all three scenarios below. The middle sampling distribution corresponds to the sampling distribution we already generated above.

<center>
```{r, tut=FALSE, echo=FALSE, message= FALSE, warning=FALSE, fig.height=6, fig.width=6, cache = TRUE}

samp <- function(n){
  x <- rnorm(n, 9.7, 2.1)
  m <- mean(x)
  s <- sd(x)
  return(c(m,s))
}

d <- as.data.frame(do.call(rbind, lapply(1:1000, function(x) samp(5))))
#d2 <- as.data.frame(do.call(rbind, lapply(1:1000, function(x) samp(50))))
d3<- as.data.frame(do.call(rbind, lapply(1:1000, function(x) samp(500))))
colnames(d) <- colnames(d2) <- colnames(d3) <- c("mean", "sd")

p <- ggplot(data = d, aes(x = mean)) +
  geom_histogram(binwidth = 0.06, color = "white", fill="green3") +
  theme_light() +
  scale_y_continuous(expand = c(0,0)) +
  labs(x = "Sample means N=5", y = NULL) +
  theme(panel.border = element_blank(), panel.grid.minor = element_blank(), legend.position = , legend.background = ) +
  xlim(6.7,12.7) 

p2 <- ggplot(data = dhh[[1000]], aes(x = Height)) +
  geom_histogram(binwidth = 0.06, color = "white", fill="green3") +
  theme_light() +
  scale_y_continuous(expand = c(0,0)) +
  labs(x="Sample means N=50", y="Frequency") +
  theme(panel.border = element_blank(), panel.grid.minor = element_blank(), legend.position = , legend.background = ) +
  xlim(6.7,12.7)

p3 <- ggplot(data = d3, aes(x = mean)) +
  geom_histogram(binwidth = 0.06, color = "white", fill="green3") +
  theme_light() +
  scale_y_continuous(expand = c(0,0)) +
  labs(x="Sample means N=500", y = NULL) +
  theme(panel.border = element_blank(), panel.grid.minor = element_blank(), legend.position = , legend.background = ) +
  xlim(6.7, 12.7)

cowplot::plot_grid(p3,p2,p, ncol = 1, align = "hv") 
```
</center>

What do we notice?

1) All histograms look normal. 
2) All distributions have approximately the same mean.
3) Distributions generated from larger samples are less dispersed.

We can take the mean of the sampling distribution itself-- **the mean of the sampling distribution is a mean of means.** This mean can be interpreted to be the same as a mean that would have resulted from a single large sample, made up of all the individual observations from each of the samples whose ${\bar{x}}$ values are included in the sampling distribution.

Note that if we had only generated a sampling distribution made up of samples of 5 giraffes, we would not have been able to exclude 11 cm as a possible value for ${\mu}$. In fact, if we were to draw a vertical line in the middle of each of the sampling distributions (the mean), we can tell that the population mean is likely even less than 10 cm.

In the following window, you will test the relationship between sampling distribution and sample size. The function below (behind-the-scenes code not shown) will plot a sampling distribution made up of 1000 samples, with each sample containing `N` number of observations. Try setting `N` to a few different values. What does the resulting sampling distribution looks like? See if you can confirm for yourself that the above points are true.


<!---LEARNR EX 1-->

<iframe class="interactive" id="myIframe1" src="https://tinystats.shinyapps.io/06-standardError-ex1/" scrolling="no" frameborder="no"></iframe>

<!------------->


# Standard Error of the Mean
As we've done before, we want to summarize this spread of mean estimates with a single value. We've already learned how to quantify a measure of spread--the standard deviation. If we take the standard deviations of each of the three different sampling scenarios above, then we accept that *distributions based on smaller samples should have larger standard deviations*. 

In the window below, calculate the standard deviation of each of the three sampling distributions (i.e. for N = 500, N = 50, and N = 5), and confirm that the italicized point above is true. (If you're working in R locally, use your "homemade" standard deviation function from the [Variance](04_variance.html) module.)

To complete this exercise, you will need to use the objects `sampling_distribution_N500`, `sampling_distribution_N50`, `sampling_distribution_N5`, which are vectors storing the thousands of ${\bar{x}}$ values from the corresponding groundhog sampling distributions.

<!---LEARNR EX 2-->

<iframe class="interactive" id="myIframe2" src="https://tinystats.shinyapps.io/06-standardError-ex2/" scrolling="no" frameborder="no"></iframe>

<!------------->


When you calculate the standard deviation of a sampling distribution of ${\bar{x}}$ values, you are calculating the **standard error of the mean (SEM)**, or just "standard error". The SEM is the value that we use to capture the level of precision of our sample estimate. But, we need a better and more efficient way to arrive at this value without relying on a groundhog day situation. Keep reading to learn more.

<div class= "alert alert-note">
> **A note about SEM:** Here "standard error" will imply standard error of the *mean*. But we can technically calculate the standard error of any sample statistic, not just the mean. We'll talk about that more in future modules.

</div>

# Time for a tea break!
<center>![](images/06_standardError/Slingshot.jpg){width=800px}</center>

# Standard error in practice
Deriving the equation used for calculating the standard error of the mean using theory (i.e. without going out and resampling MANY times) is a bit complicated, but if you're interested, you can learn more about it [here](https://stats.stackexchange.com/questions/89154/general-method-for-deriving-the-standard-error). Instead, we can capture the relationship between **standard deviation**, **sample size**, and **standard error** with the plot below.

<center>
```{r, tut=FALSE, echo=FALSE, message= FALSE, warning=FALSE, out.height=2, cache = TRUE, fig.align="center"}

d <- data.frame(N = seq(5, 1000, 5), sd = 2.1)
d$sterr <- d$sd/sqrt(d$N)

d %>%
  plot_ly(
    width = 600, 
    height = 350,
    type = 'scatter',
    mode ='markers',
    x = ~N,
    y = ~sterr,
    marker = list(size = 10, opacity = 0.75),
    hoverinfo = "text",
    text=~paste("Standard Error:", round(sterr, 3), "\nN:", round(N, 1))
  )%>%
  config(displayModeBar = F) %>%
layout(autosize = F, width = 600, height = 350,
  xaxis = list(zeroline = F),
  yaxis = list(title = "Standard Error", zeroline = F) 
  )

```
</center>

The standard deviation in this plot is `2.1`, which represents ${\sigma}$ for giraffe heights on Island 1. This population value is technically still unknown but can be deduced in theory by repeating the groundhog day example for the standard deviation instead of for the mean. It's important to note that the plot would have the same *shape* regardless of what scenario or standard deviation we were using.

**Can you figure out what the equation is for the SEM?** Look at the plot above, hover over the points, and see if you can gather how standard error of the mean, standard deviation, and sample size are related. Here are some hints:

* SEM will be on one side of the equation, standard devation, and N will be on the other.
* The equation will involve division.
* There is one more missing piece of the puzzle: When you look at the shape of the plot above. What type of function does this remind you of? We haven't covered this explicitly, but take a look [here](https://www.mathsisfun.com/sets/functions-common.html) and see if you get any ideas.

Use the window below as a calculator to see if you can figure out the equation for the SEM.

<!---LEARNR EX 3-->

<iframe class="interactive" id="myIframe3" src="https://tinystats.shinyapps.io/06-standardError-ex3/" scrolling="no" frameborder="no"></iframe>

<!------------->

  
In case you weren't able to figure it out, remember to check the `Solutions` tab in the exercise window or take a look at this [link](https://en.wikipedia.org/wiki/Standard_error) for the equation for calculating the SEM. Recall that we're working with the sample (and not population) standard deviation ($s$), so make sure you find the correct equation.

# Confirming that the SEM equation works
Let's test out the SEM equation on our original sample of `heights_island1` and compare it to what we would have gotten by taking the standard deviation of the sampling distribution example with the N= 50 case. **Does the SEM seem like a good approximation of the standard deviation of the sampling distribution?**

Below, you will use the object `heights_island1`, which contains our single sample of N=50, and the object `sampling_distribution_N50`, which contains the data from the corresponding groundhog sampling distribution.


<!---LEARNR EX 4-->

<iframe class="interactive" id="myIframe4" src="https://tinystats.shinyapps.io/06-standardError-ex4/" scrolling="no" frameborder="no"></iframe>

<!------------->

Close enough! We wouldn't expect these to be *exactly* the same because of sampling variability.
  

# How do we apply the SEM?
Now that we have a better understanding of how to gauge the precision of our sample estimates, we can test our question about the ${\mu}$ being less than 11 cm once and for all.

To formally make inferences, we need to revisit the principles of the [empirical rule](04_variance.html#empirical) to construct confidence intervals. (Confidence intervals are just one way to make inferences-- we'll discuss other ways later.)

Remember, that the SEM is just the standard deviation of the sampling distribution, so we can apply the empirical rule. As a result, ± 2 SEM from a point estimate will capture ~95% of the sampling distribution. Actually, we were a little bit sloppy earlier when we said 2 standard deviations captures 95% of a normal distribution; this will actually give you 95.45% of the data. The true value is 1.96 standard deviations--and this is what we use to construct a 95% confidence interval (CI).

Loosely speaking, a 95% CI is the range of values that we are 95% confident contains the true mean of the population. We want to know whether our guess of 11 cm falls outside of this range of certainty. If it does -- we can be sure enough that the true ${\mu}$ of giraffe heights on Island 1 is less than 11 cm.

Use the window below to find out and make your first inference by constructing the 95% CI for the `heights_island1` mean estimate!


<!---LEARNR EX 5-->

<iframe class="interactive" id="myIframe5" src="https://tinystats.shinyapps.io/06-standardError-ex5/" scrolling="no" frameborder="no"></iframe>

<!------------->

The upper limit of our 95% CI is less than 11 cm, so the population mean of heights on island 1 is likely less than 11 cm. In the scientific community, this is a bonafide way of drawing this conclusion.

<center> ![](images/06_standardError/Babyinference.jpg){width=600px} </center>


# Things to think about

We've been a little fast and loose with our words. The formal definition of CIs is the following: 

**If we were to sample over and over again, then 95% of the time the CIs would contain the true mean.**


Importantly, some examples of what the 95% CI does NOT mean are:

* A 95% CI does **not** mean that it contains 95% of the sample data.
* A CI is not a definitive range of likely values for the sample statistic, but you can think of it as estimate of likely values for the population parameter.
* It does not mean that values outside of the 95% CI have a 5% chance of being the true mean.


The precise interpretation of CIs is quite a nuanced and rather hotly debated topic [see here](https://featuredcontent.psychonomic.org/confidence-intervals-more-like-confusion-intervals/) and becomes somewhat philosophical-- so if these definition subtleties seem confusing, don't feel bad. As mentioned in the blog post linked above, one recent paper reported that 97% of surveyed researchers endorsed at least one misconception (out of 6) about CIs. 

<script>
  iFrameResize({}, ".interactive");
</script>