# Estimation in practice - Confidence Intervals

Now that we understand the implications of using $\bar{X}$ as an estimator, let's explore what this means in practice.

## The standard error of the mean

The standard error of the mean ($SEM$) is widely used and very important as it tells us how precisely we can estimate $\mu$ using $\bar{X}$. We want the $SEM$ to be small, since this means greater precision. 

The good new is that we've covered the $SEM$ already in passing, so let's briefly summarise what we learned so far: The sample mean $\bar{X}$ has a distribution with mean $\mu_{\bar{X}} = \mu$ and variance $\sigma^2_{\bar{X}} = \sigma^2/n$.

Recall that the standard deviation is just the square of the variance. In this case, the standard deviation of the estimator $\bar{X}$ can be written as
$$ \sigma_{\bar{X}} = \sqrt{\sigma^2_{\bar{X}}} = \sqrt{\sigma^2 /n} = \sigma / \sqrt{n} $$
which is the $SEM$! That's all there is to it! 

$$SEM = \sigma / \sqrt{n}$$

Recall that $\sigma^2$ denotes the _population_ variance, so $\sigma$ is the _population_ standard deviation. Just like $\mu$, it usually is unknown so we need to replace it with the _sample_ standard deviation, $s$, which we discussed in section \@ref(var-sq). The $SEM$ that you are dealing with in practice is

$$SEM = s/ \sqrt{n}.$$
Many scientists confuse the standard devation of the sample $s$ (which is a descriptive statistic) with the $SEM$ (which is an inferential statistic). This is an important distinction, so don't mix them up!

<dl>
  <dt>s, the standard deviation</dt>
  <dd>A measure of spread of the observed data, i.e. the sample. It's a descriptive statistic.</dd>

  <dt>SEM, the standard error of the mean</dt>
  <dd>A measure of precision in estimation. i.e. It's an inferential statistic.</dd>
</dl>

## The Central Limit Theorem {#CLT}

The term _Central Limit Theorem_ may sound scary but the underlying concept is very cool and very, very useful!^[Depending on who you read, the CLT is _the_ best thing to ever happen to Statistics. It's like what PCR is to molecular biology, or mayonnaise to french fries. Life just isn't the same after you discover these things.] It states that not only do we know the mean and variance of the estimator $\bar X$, we actually know its entire distribution!

Recall from the above videos and Figures \@ref(fig:3010PopSam) and \@ref(fig:3010PopSamLarge) that $\bar{X}$ has a pleasant-looking symmetrical distribution. As you may have guessed, the central limit theorem (CLT) links the distribution of the sample mean to the Normal distribution. This is very good news, because we spent a lot of time learning stuff about the Normal distribution, and now we can apply everything to the problem of estimating $\mu$ by $\bar X$! 

As a formula the CLT looks like this:

$$ \bar{X} \stackrel{n \rightarrow \infty}{\sim} N(\mu, \sigma^2/n).$$
When we say it in words it simply means "if the sample size $n$ approaches infinity, the sample mean $\bar{X}$ is Normally distributed with mean $\mu$ and variance $\sigma^2/n$".

Actually the only scary part that we don't understand right now is the meaning of "as $n$ approaches infinity". If we substitute this part with "if your sample size is big enough" the whole thing is not scary anymore at all. It basically says that, for large sample sizes, the sample mean is Normally distributed.

AMAZING!! The only question that remains is, "how large is large enough?" Many text books answer this question with $n=30$, but 30 is not a God-given magic number. The truth is that "large enough" depends on the distribution of the variable of interest. To explain this, let's explore the CLT in our web app.

We'll explore the Central Limit Theorem interactively in the workshop and discuss factors that influence the distribution of $\bar{X}$.


```{exercise, name="The CLT"}
According to the CLT, what are the benefits of having a large sample?
```

The app shows us that large samples not only ensure that the precision of your estimates increases (by decreasing the variance, and thus the $SEM$), but also that $\bar{X}$ follows a Normal distribution which is very easy to deal with since we know everything about it. This is the power of the CLT!

<dl>
  <dt>The Central Limit Theorem, CLT</dt>
  <dd>states that, if the sample size is large enough, $\bar{X}$ will be Normally distributed _regardless of the underlying distribution of your variable $X$!</dd>
</dl>


## Confidence intervals {#conf-int}

The CLT is so important because it provides the basis for many results in inferential statistics. Actually most of what we'll cover from this point on will have the CLT acting in the background. To start with, the CLT is behind the well-known formula for confidence intervals. To see this, let's return to our Martians.

For Martian height at site I, we saw that $\bar{x} = `r H.Means[1]`$ cm. What can we say about $\mu$? We now know that our $\bar{x}$ depends on the specific sample we've taken so it's very unlikely that the $\bar{x}$ of our sample of 10 Martians _exactly_ matches the population $\mu$. Therefore we can't just report $\bar{x}$ on its own. We should report a range around $\bar{x}$ that somehow reflects this uncertainty. This is the purpose of a confidence interval (CI).

Suppose we have a large sample, and thus the CLT applies. When we say that the "CLT applies", it means that our sample is large enough that we can assume a Normal distribution of the sample means. In this case an $n$ of 10 is pretty small, but we've seen that Normal assumption for the variable height seems justified, so the sample mean will also be Normally distributed. That means that the distribution of all possible sample means, $\bar{X}$, can be described by

$$ \bar{X} {\sim} N(\mu, \sigma^2/n).$$

```{r eval = FALSE, echo = FALSE}

# Draw simple normal curve:
muCI <- 0
sigmaCI <- 1
xRangeCI <- 4 # Set the upper limit of the x-axis

xCI <- seq(muCI-(xRangeCI*2), muCI+(xRangeCI*2), length=5000) # Define an x series based on inputed mu and xRange
hxCI <- dnorm(xCI, mean = muCI, sd = sigmaCI) # Define the density of the normal curve at each x position

dfCI <- data.frame(value = xCI, dist = hxCI)

# Base plot:
q0 <- ggplot(dfCI, aes (x= value, y = dist)) +
  geom_line() +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  coord_cartesian(xlim = c(-xRangeCI, xRangeCI), ylim = c(0, 0.5)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.margin = unit(c(.2,.8,0,0), "cm"),
        panel.spacing = unit(c(.2,0,0,0), "cm")
  )

q_NoAxis0 <- q0 +
  scale_x_continuous("$\bar{X}$", breaks = seq(-xRangeCI, xRangeCI, 1), labels = c(
    "$\\mu-4\\frac{\\sigma}{\\sqrt{n}}$", # -4
    "", # -3
    "$\\mu-2\\frac{\\sigma}{\\sqrt{n}}$", # -2
    "", # -1
    "$\\mu$", # 0
    "", # +1
    "$\\mu+2\\frac{\\sigma}{\\sqrt{n}}$", # +2
    "", # +3
    "$\\mu+4\\frac{\\sigma}{\\sqrt{n}}$" # +4
  ), expand = c(0,0)) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_line(colour = "black"),
        # axis.line = element_line(colour = "black"),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_text(colour = "black"),
        axis.title.x = element_blank(),
        plot.margin = unit(c(.2,.8,-0.3,0), "cm"),
        panel.spacing = unit(c(.2,0,-0.3,0), "cm"))

q_WithAxis <- q0 +
  scale_x_continuous("$\\bar{X}$", breaks = seq(-xRangeCI, xRangeCI, 1), labels = c(
    expression(paste(mu, "-4", SEM  )), #"$-4\\sigma$", # -4
    "", #"$-3\\sigma$", # -3
    expression(paste(mu, "-2", SEM)), #"$-2\\sigma$", # -2
    "", #"$-1\\sigma$", # -1
    expression(mu), #"$\\mu$", # 0
    "", #"$1\\sigma$", # -1
    expression(paste(mu, "+2", SEM)), #"$2\\sigma$", # -2
    "", #"$3\\sigma$", # -3
    expression(paste(mu, "+4", SEM)) #"$4\\sigma$" # -4
  # )
  # c(
  #   "$\\mu-4\\frac{\\sigma}{\\sqrt{n}}$", # -4
  #   "", # -3
  #   "$\\mu-2\\frac{\\sigma}{\\sqrt{n}}$", # -2
  #   "", # -1
  #   "$\\mu$", # 0
  #   "", # +1
  #   "$\\mu+2\\frac{\\sigma}{\\sqrt{n}}$", # +2
  #   "", # +3
  #   "$\\mu+4\\frac{\\sigma}{\\sqrt{n}}$" # +4
  ), expand = c(0,0))


# q_NoAxis + geom_vline(xintercept = muCI, colour = MidRed)

q_NoAxis <- q_NoAxis0 +
  geom_ribbon(data=subset(dfCI, value > qnorm(0.5 - 0.95/2) & value < qnorm(0.5 + 0.95/2)   ),
              aes(x=value,ymax=dist),ymin=0,fill=MidBlue) +
  geom_vline(xintercept = 0, colour = DarkRed, linetype = 1) +
  geom_vline(xintercept = muCI + qnorm(0.5 + 0.95/2), colour = MidRed, linetype = 2) +
  geom_vline(xintercept = muCI - qnorm(0.5 + 0.95/2), colour = MidRed, linetype = 2) +
  annotate("text", 
           x = 1.96+0.05, 
           y = 0.46, 
           label = expression(mu~1.96*SEM), #"$\\mu\\pm1.96\\frac{\\sigma}{\\sqrt{n}}$",
           # label = "$\\mu\\pm1.96\\frac{\\sigma}{\\sqrt{n}}$", 
           hjust = 0, 
           size = 3.5) +
  annotate("text", 
           x = 1.96+0.05, 
           y = 0.42, 
           label = expression(95*"% of"*means), #"$=95.0\\%$ of $\\bar{X}$", 
           hjust = 0, 
           size = 3.5)



q_WithAxis <- q_WithAxis +
  geom_ribbon(data=subset(dfCI, value > qnorm(0.5 - 0.95/2) & value < qnorm(0.5 + 0.95/2)   ),
              aes(x=value,ymax=dist),ymin=0,fill=MidBlue) +
  geom_vline(xintercept = 0, colour = DarkRed, linetype = 1) +
  geom_vline(xintercept = muCI + qnorm(0.5 + 0.95/2), colour = MidRed, linetype = 2) +
  geom_vline(xintercept = muCI - qnorm(0.5 + 0.95/2), colour = MidRed, linetype = 2) +
  annotate("text", 
           x = 1.96+0.05, 
           y = 0.46, 
           label = expression(mu~1.96*SEM), #"$\\mu\\pm1.96\\frac{\\sigma}{\\sqrt{n}}$",
           hjust = 0, 
           size = 3.5) +
  annotate("text", 
           x = 1.96+0.05, 
           y = 0.42, 
           label = expression(95*"% of"*means), #"$=95.0\\%$ of $\\bar{X}$", 
           hjust = 0, size = 3.5)

q_WithAxis

# 20 samples
value <- c(-0.17, 0.41, -1.07, 0.86, -1.19, -0.28, -0.66, 0.83, -0.61, -1.17, -0.23, -2.45, 1.51, 0.00, -0.23, 1.12, -0.42, -1.89, 0.10, -0.86)
dfSAM <- data.frame(position = 1:length(value), value)

# Modified:
base <- ggplot(dfSAM, aes(y = position, x = value)) +
  geom_rect(xmin = -1.96, xmax = 1.96, ymin = 0, ymax = 23, fill = MidBlue, col = NA) +
  geom_vline(xintercept = 0, colour = DarkRed, linetype = 1) +
  geom_vline(xintercept = muCI + qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  geom_vline(xintercept = muCI - qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  geom_point() +
  scale_x_continuous("$\\bar{X}$", breaks = seq(-xRangeCI, xRangeCI, 1), labels = c(
    expression(paste(mu, "-4", SEM  )), #"$-4\\sigma$", # -4
    "", #"$-3\\sigma$", # -3
    expression(paste(mu, "-2", SEM)), #"$-2\\sigma$", # -2
    "", #"$-1\\sigma$", # -1
    expression(mu), #"$\\mu$", # 0
    "", #"$1\\sigma$", # -1
    expression(paste(mu, "+2", SEM)), #"$2\\sigma$", # -2
    "", #"$3\\sigma$", # -3
    expression(paste(mu, "+4", SEM)) #"$4\\sigma$" # -4
    # "$\\mu-4\\frac{\\sigma}{\\sqrt{n}}$", # -4
    # "", # -3
    # "$\\mu-2\\frac{\\sigma}{\\sqrt{n}}$", # -2
    # "", # -1
    # "$\\mu$", # 0
    # "", # +1
    # "$\\mu+2\\frac{\\sigma}{\\sqrt{n}}$", # +2
    # "", # +3
    # "$\\mu+4\\frac{\\sigma}{\\sqrt{n}}$" # +4
  ), expand = c(0,0)) +
  coord_cartesian(xlim = c(-xRangeCI, xRangeCI), ylim = c(0, 23)) +
  scale_y_discrete("sample", labels = c("1","","","","5","","","","","10","","","","","15","","","","","20"), breaks = 1:20) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(colour = "black"),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text = element_text(colour = "black"),
        # axis.text.x = element_blank(),
        plot.margin = unit(c(-.2,.8,0,.11), "cm"), # corresponding to the top, right, bottom and left sides of the plot
        panel.spacing = unit(c(-.2,0,0,.31), "cm")
        )

```

```{r DistSEM, fig.cap='Distributions', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/CLT.png")
```
According to the CLT, the distribution of $\bar{X}$ is the Normal distribution with a mean $\mu$ and variance $\sigma^2/n$. The standard error, $SEM$, is $\sigma/\sqrt{n}$, which means that 95\% of all possible sample means will fall within the interval $\mu\pm1.96\frac{\sigma}{\sqrt{n}}$.

Now comes the cool stuff: since $\bar X$ follows a Normal distribution, we can apply the results from Table _Normal_perc_. For example we know that 68\% of all values fall within the interval $[\mu - \sigma/\sqrt{n} \; , \; \mu + \sigma/\sqrt{n}]$, i.e. 68\% of all the samples we can draw from the population will have a mean value $\bar{x}$ that falls in this interval. ^[This is exactly what we saw in Figure \@ref(fig:DifferentNormalCurves) and Table _Normalperc_ on page _DifferentNormalCurves_. The only difference is that the standard deviation that is associated with the mean is $\sigma/n$ rather then just $\sigma$.] Similarly,

<dl>
  <dt>The distribution of sample means</dt>
  <dd>if the sample size is large enough, there is a 95\% chance that the sample mean $\bar{x}$ lies within $[\mu - 1.96 \cdot \sigma/\sqrt{n} \; , \; \mu + 1.96 \cdot \sigma/\sqrt{n}] = \mu \pm 1.96 \cdot SEM$ (Figure _DistSEM_).</dd>
</dl>

Let's imagine that we took 20 independent samples from a population, as shown in Figure \@ref(fig:DistwithSamples). There will probably be only one of them (5\%) that results in a sample mean that falls outside the interval $\mu \pm 1.96 \cdot SEM$. 95\% of all potential sample means will fall in that interval. All of this follows from the CLT. 

```{r eval = FALSE}
grid.arrange(q_NoAxis, base, nrow=2)
```

```{r DistwithSamples, fig.cap='Sampling', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/SamplingBlue.png")
```
_Top_: The distribution of all sample means, $\bar{X}$, as shown in Figure \@ref(fig:DistSEM). _Bottom_: 95\% (here 19/20) Of all sample means $\bar{x}$ are expected to fall within, and 5\% (1/20) are expected to _not_ fall within, the interval $\mu\pm1.96\frac{\sigma}{\sqrt{n}}$.

As you might have realised the interval $\mu \pm 1.96 \cdot SEM$ looks a bit like the well--known formula for CI 95\% -- except it's centered around $\mu$. Obviuosly this is rather unpractical as $\mu$ is the unknown quantity that we are trying to estimate. What we have in front of us from the sample is $\bar x$!

The trick here is to realise that whenever $\bar x$ is close enough to $\mu$ to fall into $\mu \pm 1.96 \cdot SEM$, the interval $\bar x \pm 1.96 \cdot SEM$ that is centered around $\bar x$ will cover $\mu$. We know that there is a 95\% chance for the former to happen so this 95\% chance also holds for the latter. 95\% of all intervals $\bar x \pm 1.96 \cdot SEM$ will cover $\mu$! This is what we call the 95\% confidence interval:

\begin{eqnarray}\label{eq:95CI}
95\% CI = \bar{x} \pm 1.96 \cdot SEM.
\end{eqnarray}

<dl>
  <dt>95% Confidence Intervals</dt>
  <dd>If the sample size is large enough, there is a 95\% chance that the interval $\bar{x} \pm 1.96 \cdot SEM$ covers the true population mean $\mu$. The width of the interval decreases with larger sample sizes.</dd>
</dl>

An analogy I like to use is like digging for buried treasure in the desert. As an archaeologist, we want to dig as small a hole as possible and be as close as possible to the treasure. Lucky for us it's not moving! What is moving is the location and size of the hole we dig. The 95\% CI tells us that we'll find the treasure 95% of the time, but that means that there's still a 5% chance we're going to miss the true $\mu$!

While this all sounds great, it is important to remember that

> We never know which $\bar{x}$ we have, cf. Figure _BlueRedFig_. We may have the bad luck to get one of the 5\% of $\bar{x}$ values that fall outside the region $\mu \pm 1.96 \cdot SEM$ in which case our 95\% CI, $\bar{x} \pm 1.96 \cdot SEM$ will miss $\mu$.

```{r eval = FALSE, echo = FALSE}
BlueValue <- value[16]

q_NoAxisBlue <- q_NoAxis0 +
  geom_ribbon(data=subset(dfCI, value > BlueValue-1.96 & value < BlueValue+1.96),
              aes(x=value,ymax=dist),ymin=0,fill=MidBlue) +
  geom_vline(xintercept = 0, colour = DarkRed, linetype = 1) +
  geom_vline(xintercept = BlueValue + qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  geom_vline(xintercept = BlueValue - qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  annotate("text", x = BlueValue - 1.96-0.05, y = 0.46, label = "$95.0\\%$ CI $=\\bar{x}\\pm1.96\\frac{\\sigma}{\\sqrt{n}}$", hjust = 1, size = 3.5) +
  annotate("text", x = BlueValue - 1.96-0.05, y = 0.42, label = "$\\mu$ is inside the $95.0\\%$ CI", hjust = 1, size = 3.5)



baseBlue <- ggplot(dfSAM, aes(y = position, x = value)) +
  geom_rect(xmin = value[16]-1.96, xmax = value[16]+1.96, ymin = 0, ymax = 23, fill = MidBlue, col = NA) +
  geom_errorbarh(aes(xmax = BlueValue + 1.96, xmin = BlueValue - 1.96, y = 16)) +
  geom_vline(xintercept = 0, colour = DarkRed, linetype = 1) +
  geom_vline(xintercept = value[16] + qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  geom_vline(xintercept = value[16] - qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  geom_point() +
  scale_x_continuous("$\\bar{X}$", breaks = seq(-xRangeCI, xRangeCI, 1), labels = c(
    "$\\mu-4\\frac{\\sigma}{\\sqrt{n}}$", # -4
    "", # -3
    "$\\mu-2\\frac{\\sigma}{\\sqrt{n}}$", # -2
    "", # -1
    "$\\mu$", # 0
    "", # +1
    "$\\mu+2\\frac{\\sigma}{\\sqrt{n}}$", # +2
    "", # +3
    "$\\mu+4\\frac{\\sigma}{\\sqrt{n}}$" # +4
  ), expand = c(0,0)) +
  coord_cartesian(xlim = c(-xRangeCI, xRangeCI), ylim = c(0, 23)) +
  scale_y_discrete("sample", labels = c("1","","","","5","","","","","10","","","","","15","","","","","20"), breaks = 1:20) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(colour = "black"),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text = element_text(colour = "black"),
        # axis.text.x = element_blank(),
        plot.margin = unit(c(-.17,.8,0,.11), "cm"), # corresponding to the top, right, bottom and left sides of the plot
        panel.spacing = unit(c(-.17,0,0,.11), "cm")
        )

RedValue <- min(value)

q_NoAxisRed <- q_NoAxis0 +
  geom_ribbon(data=subset(dfCI, value > RedValue-1.96 & value < RedValue+1.96),
              aes(x=value,ymax=dist),ymin=0,fill=MidRed) +
  geom_vline(xintercept = 0, colour = DarkRed, linetype = 1) +
  geom_vline(xintercept = RedValue + 1.96, colour = DarkRed, linetype = 2) +
  geom_vline(xintercept = RedValue - 1.96, colour = DarkRed, linetype = 2) +
  annotate("text", x = RedValue + 1.96+0.05, y = 0.46, label = "$95.0\\%$ CI $=\\bar{x}\\pm1.96\\frac{\\sigma}{\\sqrt{n}}$", hjust = 0, size = 3.5) +
  annotate("text", x = RedValue + 1.96+0.05, y = 0.42, label = "$\\mu$ is outside the $95.0\\%$ CI", hjust = 0, size = 3.5)
#   annotate("text", x = RedValue+0.05, y = 0.45, label = "$\\mu\\pm1.96\\sigma$", hjust = 0, size = 3.5) +
#   annotate("text", x = RedValue+0.05, y = 0.42, label = "$=95.0\\%$", hjust = 0, size = 3.5)


baseRed <- ggplot(dfSAM, aes(y = position, x = value)) +
  geom_rect(xmin = min(value)-1.96, xmax = min(value)+1.96, ymin = 0, ymax = 23, fill = MidRed, col = NA) +
  geom_errorbarh(aes(xmax = RedValue + 1.96, xmin = RedValue - 1.96, y = 12)) +
  geom_vline(xintercept = 0, colour = DarkRed, linetype = 1) +
  geom_vline(xintercept = min(value) + qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  geom_vline(xintercept = min(value) - qnorm(0.5 + 0.95/2), colour = DarkRed, linetype = 2) +
  geom_point() +
  scale_x_continuous("$\\bar{X}$", breaks = seq(-xRangeCI, xRangeCI, 1), labels = c(
    "$\\mu-4\\frac{\\sigma}{\\sqrt{n}}$", # -4
    "", # -3
    "$\\mu-2\\frac{\\sigma}{\\sqrt{n}}$", # -2
    "", # -1
    "$\\mu$", # 0
    "", # +1
    "$\\mu+2\\frac{\\sigma}{\\sqrt{n}}$", # +2
    "", # +3
    "$\\mu+4\\frac{\\sigma}{\\sqrt{n}}$" # +4
  ), expand = c(0,0)) +
  coord_cartesian(xlim = c(-xRangeCI, xRangeCI), ylim = c(0, 23)) +
  scale_y_discrete("sample", labels = c("1","","","","5","","","","","10","","","","","15","","","","","20"), breaks = 1:20) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(colour = "black"),
        axis.line.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text = element_text(colour = "black"),
        # axis.text.x = element_blank(),
        plot.margin = unit(c(-.17,.8,0,.11), "cm"), # corresponding to the top, right, bottom and left sides of the plot
        panel.spacing = unit(c(-.17,0,0,.11), "cm")
        )
grid.arrange(q_NoAxisBlue, baseBlue, nrow=2)
grid.arrange(q_NoAxisRed, baseRed, nrow=2)

```

```{r BlueRedFig, fig.cap='Successful and Unsuccessully covering the population mean.', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/SamplingBlueRed.png")
```

_Left_: When centering the 95\% CI over an $\bar{x}$, we will cover $\mu$ 95\% of the time. _Right_: In 5\% of the cases, we will fail to cover $\mu$.

From equation (_eq:95CI_) we see that the smaller the $SEM$ is, the narrower is the CI. Conceptually, this should make sense, since you've seen that the $SEM$ reflects the uncertainty that is associated with using $\bar{X}$ as an estimator for the unknown $\mu$.

```{exercise, name="The size of confidence intervals"}
What can you do to decrease the width of your confidence interval? Hint: think about the $SEM$!
```

Aside from the two parameters that influence the size of $SEM$ (and thus the CI), there is another factor to consider - the desired level of confidence. So far we have only talked about the 95\% CI, but you may also want to be 99\% certain that the interval you state covers the true population mean.

```{exercise, name="The relationship between confidence and precision"}
If we increased the CI from 95\% to 99\%, would the interval get smaller or larger?
```

```{exercise, name="Reporting confidence intervals"}
Look at the properties of the normal curve in Table \ref{Normal_perc} on page _DifferentNormalCurves_. What CI do scientists refer to when they report $\bar{x} \pm SEM$ in publications? Hint: Give the percent.
```

This is pretty much all you need to know about confidence intervals, but there are two points to keep in mind.

First, we ignored that fact that in any practical situation we don't know the population variance, $\sigma^2$, or standard deviation $\sigma$. That means we also don't _really_ know the $SEM$! Remeber, we're using equation XX, with $s$, not equation XX, with $\sigma$. Nonetheless, since we're talking about situations where the sample size is large, otherwise the CLT would not apply, we can quite safely replace $\sigma$ with the sample estimate $s$. We'll discuss the case when we need to account for the uncertaintly in using $s$ in section _sec:the_t_distributiuon_.

Second, we deliberately repeated "if your sample size is large enough" many times because all of this holds true only in that case. If your sample size is small, it's still possible to remedy the situation and come up with similar statements _if_ your variable of interest (or a transformed version) can safely be assumed to be Normally distributed in the population. If you have a small sample size and are not sure if the variable that you are looking at follows a Normal distribution in the population, you know so little about your system that you should be content with sample description or use non--parametric methods some of which we'll cover later on.

## Sample size calculation for beginners {#width}

The previous section tells us that we can be 95\% certain that the true mean height of all Martians is covered by the interval $[`r Height.95CI`]$. Are you happy with this result? In some cases, the calculated CIs are really unsatisfactorily large, so what can we do to add precision and obtain small confidence intervals? The bad news is that once the study is done, all we can do is accept what we have and aim for more precision next time we do a study. There is, however, a way to influence the width of the confidence interval when you start thinking about it BEFORE you even start conducting the study.

As we have seen, the width of the CI depends on the $SEM$ which in turn depends on the sample size. Large samples give small $SEM$s which give small CIs and therefore high precision. If you know that in your particular case, it's absolutely important to measure some quantity with a precision that results in a 95\% CI of a certain width $w$, you can just look at formula (_eq:95CI_) again to see that

\begin{eqnarray*}
w & = & 2 \times 1.96 \times SEM \\
 & = & 3.92 \times SEM \\
 & = & 3.92 \times \frac{\sigma}{\sqrt{n}}
\end{eqnarray*}
which means that achieving a width of $w$ would need a sample size of
\begin{eqnarray*}
\sqrt{n} & = &  \frac{3.92 \times \sigma}{w} \\
n & = &  \frac{15.37 \times \sigma^2}{w^2}
\end{eqnarray*}

would do the trick! As expected, this formula tells us that we need large samples if the desired width $w$ is small (something small in the denominator makes the fraction large) and/or if the variance is large.^[Caution: in order to use the formula, you need to have an idea of the magnitude of the variance BEFORE you do the study! You could use values from previous similar studies or make some other informed guess.]

```{exercise, name="Basic sample size calculations"}
How can you calculate the sample size required for a 99\% CI to have a certain width $w$?
```

## Using error bars in figures

To reiterate some of the things we've learned so far about uncertainty and precision, let us consider error bars which put things into a practical context. We are now equipped to critically examine the use of error bars, but before we do so, let's remind ourselves how they are defined:

Error bars are a graphical representation of the variability of data and are used on graphs to indicate the error, or uncertainty in a reported measurement. They give a general idea of how accurate a measurement is, or conversely, how far from the reported value the true (error free) value might be. Error bars often represent one standard deviation of uncertainty, one standard error, or a certain confidence interval (e.g., a 95\% interval).^[Wikipedia, The Free Encyclopedia, s.v. "Error bar," (accessed Januar 04, 2018), `https://en.wikipedia.org/wiki/Error\_bar`.]

```{exercise, name="Types of error bars"}
Often you see plots with error bars around the sample mean of the data. What is the true (error free) value in this case?
```

<!-- %the population mean $\mu$ which is unknown -->

```{r eval = FALSE, echo = FALSE}

# population
# set.seed(16)
pop <- df.original$obs

#Sample sizes
sam.size <- c(3, 10, 30, 100)
# sam.size <- c(10, 30, 60)

df <- data.frame(sample = do.call(rep, list(sam.size, sam.size)),
                 values =  do.call(c, lapply(sam.size, function (x) sample(pop, size = x)))
                )

df.sum <- plyr::ddply(df, c("sample"), summarise,
                mean = mean(values),
                sd = sd(values),
                ci = 1.96*(sd(values)/sqrt(length(values))),
                sem = sd(values)/sqrt(length(values))
               )

df.sum <- reshape2::melt(df.sum, id.vars = c("sample","mean"))

dodge <- position_dodge(width = 4)


m <- seq(-2,2,1)
xValues <- c(H.Means[1]+(m*H.SDs[1]))
# xLabels <- paste(xValues)
sign <- c(rep("-",2), "", rep("+", 2))

xLabels <- paste0(xValues,
                 "\n$\\mu$",
                 sign,
                 abs(m),
                 "$\\sigma$")

xLabels[3] <- paste0(xValues[3],
                 "\n$\\mu$")



ggplot(df.sum, aes(x = sample, y = mean, col = variable, group = variable)) +
  geom_point(position = dodge) +
  scale_x_continuous("Sample Size",
                     labels = sam.size,
                     breaks = sam.size,
                     limits = c(0,105),
                     expand = c(0,0)) +
  scale_y_continuous("Mean and spread",
                     limits = range(xValues),
                     labels = xLabels,
                     breaks = xValues,
                     expand = c(0,0)) +
  scale_color_brewer("Type of error bar",
           palette = "Dark2",
           labels = c("Standard Deviation",
                "95\\% Confidence Interval",
                "SEM")) +
  geom_errorbar(aes(ymin=mean-value, ymax=mean+value), width=.1, position=dodge) +
  geom_hline(yintercept = mean(pop)) +
  geom_hline(yintercept = mean(pop) + sd(pop), linetype="dotted") +
  geom_hline(yintercept = mean(pop) - sd(pop), linetype="dotted") +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
#         legend.direction = "horizontal",
        legend.key = element_blank(),
        axis.ticks = element_line(colour = "black"),
        axis.text = element_text(colour = "black", vjust = 0.75),
        axis.title.y = element_text(colour = "black", size = 9),
        axis.title.x = element_text(colour = "black", size = 9, vjust=-0.2),
        axis.line = element_line(colour = "black"),
#         legend.position = "top",
       legend.position = c(0.75, 0.90),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
# ,
#     plot.margin = unit(c(3, 1, 0.5, 0.5), "lines"),
#     panel.spacing = unit(0.25, "lines")
)


```

```{r dif-types-spread, fig.cap='Comparing different types of error bars.', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/ErrorBars.png")
```
The horizontal lines depict our population mean (solid) and standard deviation (dashed) from a simulated population. Selecting samples of different sizes, i.e. 10, 30 or 60 observations, affects error bars in different ways. Standard deviation is relatively stable across sample sizes, whereas the 95\% confidence interval and the SEM get smaller as sample size increases.

Geoff Cumming and colleagues have a couple of very informative papers on the topic of error bars.^["Error bars in experimental biology". Geoff Cumming, Fiona Fidler, and David L. Vaux. The Journal of Cell Biology, Vol. 177, No. 1, April 9, 2007, pp 7-11.] The authors distinguish between descriptive and inferential error bars.

```{exercise, name="Types of error bars"}
Refer to Figure _dif-types-spread_ and fill in Table \@ref(tab:cat-types). In the _type_ column, indicate if the error bar is descriptive or inferential. In the _description_ column, define, using simple language, what each error bar tells us (no equations or text book definitions!).
```

Table: (\#tab:cat-types) Types of Error bars.

| Error Bar                  | Type | Explanation |
|----------------------------|------|-------------|
| Standard Deviation         |      |             |
| Standard Error of the mean |      |             |
| 95% Confidence Interval    |      |             |

```{exercise, name="Reporting uncertainty"}
What should you think if you see a paper that contains a figure with error bars, without specifying what they represent?
```

To end our discussion about error bars, consider this quote from the GraphPad Prism website^[graphpad.com/guides/prism/6/ statistics/index.htm?statwh entoplotsdvssem.htm]:

> If your goal is to emphasise small and unimportant differences in your data, show your error bars as $SEM$, and hope that your readers think they are SD. If our goal is to cover-up large differences, show the error bars as the SD for the groups, and hope that your readers think they are $SEM$s. 

Some word of advice: if you don't understand the joke, review the differences between SD and SEM.

## 95% Confidence Intervals Exercise

Follow this [link](https://docs.google.com/spreadsheets/d/1YTwFP_ko39FSULcK0H1e4m8Qp7FuopJtlgA4ak6djiY/edit?usp=sharing) to enter you values for the in-class exercise.
