# Introduction to Hypothesis Testing {#Testing}

We've already come pretty far in our understanding of statistics and have even taken our first steps into the territory of inferential statistics. 

<!-- % You may be wondering how we managed to go on for this long without once mentioning the word $p$--value. Unfortunately, in many disciplines \emph{doing statistics} is eponymous with calculating $p$--values. We hope that we've convinced you by now that there is much more to statistics than counting the number of stars behind the letter $p$!\marginnote{Before we proceed, one further note on hypothesis testing. We placed importance at the beginning of the workshop on having a hypothesis. This may lead you to believe that everything in research leads to hypothesis testing. This is all very well if your data is suited for that, but with small datasets or data that was collected for exploration purposes, it's totally fine -- and even more honest and correct -- to not produce a $p$--value as the end result.} -->

<!-- % However, whether we like it or not, we can't deny that $p$--values play a big role in research, so -->
<!-- Now we're equipped to proceed to the arena of hypothesis testing, which is where the famous $p$--values live in the statistical landscape. -->

## From confidence intervals to hypothesis testing

Confidence intervals and hypothesis testing are close relatives. To see this, let us first remind ourselves how we calculate and interpret confidence intervals. If you have a hard time with the following exercises, refer back to page _sec:conf-int_.

```{r echo = FALSE}

# Data:
martian <- read.delim("data/martian.txt")

# Descriptive stats for both sites
H.Means <- tapply(martian$Height, martian$Site, mean)
H.SDs <- plyr::round_any(tapply(martian$Height, martian$Site, sd),1)
mystery <- 215

# Descriptive stats for site I
Height.mean <- H.Means[1]
Height.sd <- H.SDs[1]
nValue <- length(martian$Height[martian$Site == "Site I"])

# Inferential stats:
Height.SEM <- (Height.sd/sqrt(nValue))
Height.95CI <- c(Height.mean - 1.96*Height.SEM, Height.mean + 1.96*Height.SEM)

# Values for t-tests
distance <- 5 # Set distance away from x-bar for tests
HypoValue <- Height.mean - distance # Set less than
HypoValue2 <- Height.mean + distance # Set greater than.

# Signal:Noise ratio
tValue <- (Height.mean - HypoValue)/Height.SEM
tValue2 <- (Height.mean - HypoValue2)/Height.SEM

# t-tests
# t.test(martian$Height[martian$Site == "Site I"], mu = HypoValue)
# t.test(martian$Height[martian$Site == "Site I"], mu = HypoValue, alternative = "greater")
# t.test(martian$Height[martian$Site == "Site I"], mu = HypoValue, alternative = "less")

# Get p-values:
pValue_2Side <- 2 * (1 - pt(tValue, nValue - 1))
pValue_1SideLarge <- 1 - pt(tValue, nValue - 1)
pValue_1SideSmall <- pt(tValue, nValue - 1)

pValue_2Side <- plyr::round_any(pValue_2Side, 0.001)
pValue_1SideLarge <- plyr::round_any(pValue_1SideLarge, 0.0001, f = floor)
pValue_1SideSmall <- plyr::round_any(pValue_1SideSmall, 0.0001, f = ceiling)

Z_pValue_2Side <- 2 * (1 - pnorm(tValue))
Z_pValue_1SideLarge <- 1 - pnorm(tValue)
Z_pValue_1SideSmall <- pnorm(tValue)

Z_pValue_2Side <- plyr::round_any(Z_pValue_2Side, 0.001)
Z_pValue_1SideLarge <- plyr::round_any(Z_pValue_1SideLarge, 0.0001, f = floor)
Z_pValue_1SideSmall <- plyr::round_any(Z_pValue_1SideSmall, 0.0001, f = ceiling)
```

```{exercise, name="Defining the 95\\% Confidence Interval"}
Imagine that a 95% confidence interval for the height of Martians at Site I is [xxx, xxx]. Write a sentence that correctly describes this finding.
```

```{exercise, name="Calculating from the 95\\% Confidence Interval"}
What is the observed mean? What is the SEM?
```

To see the connection between confidence intervals and hypothesis testing, image that, before we began our exploration to Mars, we did some background research. It seems that Earthlings have always had a notion of what Martians are supposed to look like. Based on this, we formulate the _a priori_ hypothesis (i.e. an hypothesis that is based on a notion we obtained _before_ we collected the sample) that Martians should on average be $`r HypoValue`$ cm tall.


```{exercise, name="Assessing the plausability of a hypothesis"}
Based on our sample, do you believe this? Do our data support this hypothesis? Give reasons as to why or why not.
```

Congratulations! You have just performed a hypothesis test without even noticing! Let's explain:
Recall how we derived the interpretation of the confidence interval above?^[Whenever we talk about the CLT, we implicitly assume that we are in the situation where the sample size $n$ is big enough so that the it holds true!] Our starting point was the CLT, from which we know that 95% of all sample means $\bar{x}$ that we could have gotten from the population will fall into the interval $\mu \pm 1.96 \cdot SEM$. Then we noted that this implies that 95% of all the intervals $\bar{x} \pm 1.96 \cdot SEM$ will in turn contain the true underlying $\mu$, i.e. for the single confidence interval that you get from your sample, the chance that it covers $\mu$ is 95%.

To highlight the relationship between confidence intervals and hypothesis testing, we can use the exact same trick with one further gimmick which relates to the main difference between the two methods, i.e. the presence of an explicit hypothesis! In our example, the hypothesis is that the true mean Martian height is $`r HypoValue`$ cm. An hypothesis of this kind is commonly called _null hypothesis_. We can write $H_0: \mu=`r HypoValue` cm$.

Now suppose that $H_0$ indeed describes the true state of the world and the mean height of the entire Martian population at Site I really _is_ `r HypoValue`cm. In this case we know from the CLT, that the distribution of the sample mean $\bar{X}$ is ^[In hypothesis testing, we pretend that $H_0$ is true, and then look how likely it is that we observe the given data under this premise! In other words, we find out how much our data support the null hypothesis.]

$$\bar{X} \sim N(`r HypoValue`, SEM^2) \quad \mbox{or} \quad \bar{X} \sim N(`r HypoValue`, `r Height.SEM`^2).$$

This implies that 95\% of all sample means, $\bar{x}$, would lie in the interval $`r HypoValue` \pm 1.96 \cdot `r Height.SEM` = `r HypoValue` \pm `r 1.96 * Height.SEM` = (`r c(HypoValue - 1.96 * Height.SEM, HypoValue + 1.96 * Height.SEM)`)$, i.e. the chance of observing an $\bar{x}$ _outside_ this interval is only 5\%. Therefore, if we _do_ observe an $\bar{x}$ outside of this interval, we are inclined to _reject the null hypothesis- with 95\% confidence.

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

mu <- HypoValue
sigma <- Height.SEM
xRangeMin <- plyr::round_any(mu-(3*sigma), 0.5, f = floor) # Set the upper limit of the x-axis
xRangeMax <- plyr::round_any(mu+(3*sigma), 0.5, f = ceiling) # Set the upper limit of the x-axis

x <- seq(xRangeMin, xRangeMax, length=5000) # Define an x series based on inputed mu and xRange
hx <- dnorm(x, mean = mu, sd = sigma) # Define the density of the normal curve at each x position

test.df <- data.frame(xval = x, value = hx)

# Base plot:
q <- ggplot(test.df, aes ( x= xval, y = value)) +
  geom_line() +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )

baseValues <- c(rep(HypoValue, 6), Height.mean, HypoValue)
m <- c(seq(-3,2,1),0,3)

xValues <- plyr::round_any(baseValues + m * Height.SEM,.01)

sign <- c("-","-","-","","+","+","","+")
m2 <- abs(m)
m2[c(4,7)] <- ""
symbols1 <- c(rep("$\\mu_0$",6), "$\\bar{x}$", "$\\mu_0$")

# char <- "SEM"

# These next two symbols just make the SEM symbols in the plot, insead of writing "SEM"
# char <- "$\\sigma / \\sqrt{n}$"
char <- "$\\frac{\\sigma}{\\sqrt{n}}$"
symbols2 <- c(rep(char, 3), "", char, char, "", char)

q_base <- q +
  coord_cartesian(xlim = c(xRangeMin, xRangeMax), ylim = c(0, 0.25)) +
  scale_x_continuous("Martian Height at Site I (cm) under premise of H0",
                     breaks = xValues,
                     limits = c(xRangeMin,xRangeMax),
                     labels = paste0(xValues,"\n",symbols1,sign,m2,symbols2),
                     expand = c(0,0)) +
  #  geom_vline(xintercept = mu, colour = "#C42126")
  # geom_ribbon(data = subset(test.df, test.df$variable == "dist"), aes(x=xval,ymax=value),ymin=0,fill="#F5B319") +
  geom_ribbon(data = subset(test.df, test.df$xval > mu-(1.96*sigma) & test.df$xval < mu+(1.96*sigma)),
              aes(x=xval,ymax=value),ymin=0,fill=MidBlue) +
  # ggtitle(paste0("$",HypoValue,"\\pm 1.96SEM = 95 \\%$ Confidence Interval")) +
#   geom_vline(xintercept = mu + (1.96*sigma), colour = MidRed, linetype = 2) +
#   geom_vline(xintercept = mu - (1.96*sigma), colour = MidRed, linetype = 2) +
#   geom_vline(xintercept = HypoValue, colour = DarkOrange) +
  geom_vline(xintercept = Height.mean, colour = DarkGreen)

# q_base

```

```{r Draw-normal-testing1, fig.cap='Distribution of the sample mean under the null hypothesis.', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/Hypo1.png")
```
i.e. under the premise that the mean height in the population is $\mu_0 = `r HypoValue`$. The observed mean value $\bar{x} = `r Height.mean`$, lies outside the region $\mu_0 \pm 1.96 \cdot SEM$ which provides evidence against $H_0$.



This situation is depicted in Figure _fig:Draw-normal-testing1_ where the green bar shows our observed sample mean $\bar{x} = `r Height.mean`$ and the blue region is the interval $\mu_0 \pm 1.96 \cdot SEM$ that we calculated with our hypothesised value $\mu_0 = `r HypoValue`$ and our $SEM = `r Height.SEM`$. We see that the probability of observing data like ours (i.e. data that yield an observed sample mean $\bar{x} = `r Height.mean`$) under the premise of $H_0$ is not very big. In fact, it's smaller than 5\% as $\bar{x}$ lies outside the blue region. Therefore the data provide evidence against this null hypothesis. We would reject it with 95\% confidence!

```{exercise, name="Rejeting hypotheses"}
Would you reject $H_0: \mu = 196 cm$?
```

The cool thing is that the above reasoning holds for _all_ hypothesised values $\mu_0$ that are outside of the 95\% confidence interval. This is because if $\mu_0$ does not lie inside $\bar{x} \pm 1.96 \cdot SEM$, it's also^[All values $\mu_0$ that are not inside the 95\% confidence interval can be rejected with 95\% confidence.] true that $\bar{x}$ does not lie inside $\mu_0 \pm 1.96 \cdot SEM$! So in our example, even without being much of a mathematician, you can say right away that all hypothesised mean values outside the 95\% CI range $[`r c(Height.mean - 1.96 * Height.SEM, Height.mean + 1.96 * Height.SEM)`]$ are not consistent with our data. All values _within_ the confidence interval, however, remain potential candidates for the true mean Martian height in our population (at site I). Magic! 

Before we move on deeper into the realm of hypothesis testing, let's summarise what we've learned so far about the relationship between confidence intervals and hypothesis testing:

1. A confidence interval contains a range of possible values for the parameter of interest that are consistent with the data, i.e. that are credible possibilities for the true underlying population value of the parameter.

2. In significance tests only one possible value for the parameter, called the hypothesised value, is tested. We determine whether our data are consistent with the proposition that the hypothesised value is the true underlying parameter value.

3. All values outside the confidence interval would be rejected if they were hypothesised in a significance test.^[This is is a general statement! Here you can choose the confidence level as you like. If you are happy with the idea to accidentally reject $H_0$ even if it's true with a chance of $\alpha$\%, then choose a ($1-\alpha$)\% confidence interval. More on $\alpha$ later.] The confidence level determines the certainty with which you'd reject the null hypothesis, e.g. all values outside a 95\% confidence interval can be rejected if we accept that there is a 5\% chance that we make a mistake.



## A general recipe for hypothesis testing {#recipeTesting}

Understanding the previous paragraph is crucial, as it provides a general recipe for hypothesis tesing that ermerges when we disassemble the steps that we took to reach a decision about our null hypothesis:

<dl>
  <dt>1. Pretend that your null hypothesis is true</dt>
  <dd>in our example we pretended that $\mu_0=`r HypoValue`$ cm is the true underlying mean height in the Martian population at site I</dd>
  
  <dt>2. Test distribution</dt>
  <dd>Derive the distribution of your statistic under $H_0$: if $\mu_0=`r HypoValue`$ cm is true, we know that $\bar{X} \sim N(`r HypoValue`, SEM^2)$, i.e. $\bar{X} \sim N(`r HypoValue`, `r Height.SEM`^2)$.</dd>

  <dt>3. Asses</dt>
  <dd>Assess how unlikely it is that you would see your data under the premise of $H_0$: in our example, the probability to obtain a sample with a mean height of $\bar{x} = `r Height.mean`$ cm is less than 5\%.</dd>

  <dt>4. Reach a decision</dt>
  <dd>as we considered 'less then 5\%' above to be small enough, we concluded that there is enough evidence in the data to reject $H_0$.</dd>

</dl>

In principle, that's all there is to is -- case closed! We'll spend a bit more time on this topic though to investigate a few subtleties and introduce some official stats jargon.


## The signal and the noise -- what we do in hypothesis testing {#signalnoise}

At the beginning of our chapter on estimation on page _chap:estimation_, we stated that in hypothesis testing, we address the general question 'In case I see something interesting in my data, is the effect real (i.e. due to some systematic pattern) or is it just due to chance?'. Let's see what we meant by this!

In our example where we observe a mean height of $\bar{x} = `r Height.mean`$ cm, and we want to assess the credibility of the null hypothesis $H_0: \mu = `r HypoValue`$ cm, what do you think is _the effect_ that we talk about?

```{exercise, name="Assessing credibility"}
```

Here, _the effect_ is the fact that the sample mean isn't equal to the hypothesised population mean, i.e. $\bar{x} - \mu_0 \neq 0$. In hypothesis testing we want to know if this effect is real (i.e. if our data and the null hypothesis are far enough apart to count as a systematic pattern) or if the distance $\bar{x} - \mu_0$ we observe could have occurred just by chance despite $\mu_0$ being the true underlying value. 

<!-- % In our example, we therefore ask ourselves whether observing a sample mean of $\bar{x} = \Sexpr{Height.mean}$ cm could just be due to random sampling variation even if the underlying population mean is equal to the hypothesised value $\mu_0 = \Sexpr{HypoValue}$ cm. -->

If you rephrase this question a bit, you are asking whether there is a real signal among the noise we have in the data, where the signal corresponds to the effect we defined above and the noise corresponds to the random sampling variation. 

In our example, the reason why observing $\bar{x} = `r Height.mean`$ cm is so unlikely under $H_0: \mu = `r HypoValue`$ cm is that the difference between the two values (i.e. the signal $\bar{x} - \mu_0$) is quite a bit higher than the $SEM$ which, as we know, can also be seen as the noise we have in estimating $\mu$ using the sample mean, $\bar{x}$. The signal to noise ratio we are looking at here can therefore be written as

$$ \frac{\mbox{signal}}{\mbox{noise}} =  \frac{\bar{x} - \mu_0}{SEM} = \frac{`r Height.mean` - `r HypoValue`}{`r Height.SEM`} = `r tValue`. $$

Does that remind you of something? It should! This is essentially what you saw when we described $Z$--scores! ^[Recall that N(0,1) means the standard Normal distribution, with a $\mu$ of 0 and a $\sigma^2$ of 1. The notation $\stackrel{H_0}{\sim}$ is a quick and lazy way of saying 'distributed under the null hypothesis'. This formula tells us that $(\bar{X} - \mu_0)/SEM$ has a standard Normal distribution if $H_0: \mu = \mu_0$ is true.]
This is really good news because it provides us with a standardised way to assess whether the signal to noise ratio is big enough to reject the null hypothesis! Instead of writing $\bar{X} \sim N(`r HypoValue`, `r Height.SEM`^2)$ under the null hypothesis (cf. point 2 in section _sec:recipeTesting_), we can equivalently state that
$$ \frac{\bar{X} - \mu_0}{SEM} \, \stackrel{H_0}{\sim} \, N(0,1). $$
This means that if we want to know if our observed signal--to--noise ratio, $(`r Height.mean` - `r HypoValue`)/`r Height.SEM` = `r tValue`$ is big, i.e. unlikely to have occurred under $H_0$, we can just have a look where it sits on the standard normal curve! If it's far off from zero, we are inclined to reject $H_0$ and if it's close to the center, the hypothesised value $\mu_0 = `r HypoValue`$ is a plausible value for the true mean height!

This is a neat standardised way to carry out step 3 in section _sec:recipeTesting_! It is depicted in Figure _fig:Draw-normal-testing2_ where our observed $z-$score^[The signal--to--noise ratio in the current context is essentially the same thing as a $Z$--score!!!] of `r tValue` lies outside the interval [-1.96, 1.96] which, as you may recall, is the interval where 95\% of all values fall when the underlying distribution if the standard normal distribution. As before, this means if $\mu_0 = `r HypoValue`$ really was the true population mean, it would be less than 5\% likely to observe $\bar{x}=`r Height.mean`$. However, instead of  wondering if $\bar{x} = `r Height.mean`$ and $\mu_0 = `r HypoValue`$ are far apart given the sampling variation, we can right away say that $z=`r tValue`$ is far off on the standard normal curve! ^[And this is because with $Z$--scores, the limit for a 95\% certainty will always be 1.96.]




```{r echo = FALSE}


mu <- 0
sigma <- 1
xRangeMin <- -3 # Set the upper limit of the x-axis
xRangeMax <- 3 # Set the upper limit of the x-axis

x <- seq(xRangeMin, xRangeMax, length=5000) # Define an x series based on inputed mu and xRange
hx <- dnorm(x, mean = mu, sd = sigma) # Define the density of the normal curve at each x position
hstd <- dnorm(x, mu) # Define the density of the standard normal curve at each x position

test.df <- data.frame(xval = x, dist = hx, std = hstd)
test.df <- reshape2::melt(test.df, id.var = "xval")

# Base plot:
q <- ggplot(subset(test.df, test.df$variable == "std"), aes ( x= xval, y = value)) +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  coord_cartesian(xlim = c(xRangeMin, xRangeMax), ylim = c(0, 0.5)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text.y = element_text(colour = "black"),
        axis.text.x = element_text(colour = "black", hjust = 0.5),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )
tValue <- 2.64
q +
  scale_x_continuous("z-scores", breaks = c(seq(-3,2,1),tValue,3),
                     limits = c(-3,3),
                     labels = c(seq(-3,2,1),tValue,3), expand = c(0,0)) +
  #  geom_vline(xintercept = mu, colour = "#C42126")
#  geom_ribbon(data = subset(test.df, test.df$variable == "dist"), aes(x=xval,ymax=value),ymin=0,fill="#F5B319") +
#  geom_ribbon(data = subset(test.df, test.df$xval > mu-(1.96*sigma) & test.df$xval < mu+(1.96*sigma) & test.df$variable == "dist"),
#              aes(x=xval,ymax=value),ymin=0,fill="#005AA5") +
#  ggtitle("$15.1 \\pm 1.96SEM = 95 \\%$") +
  geom_vline(xintercept = mu + (1.96*sigma), colour = "#D67A64", linetype = 2) +
  geom_vline(xintercept = mu - (1.96*sigma), colour = "#D67A64", linetype = 2) +
  geom_vline(xintercept = mu, colour = DarkOrange) +
  geom_vline(xintercept = tValue, colour = DarkGreen) +
  geom_line()

```

_fig:Draw-normal-testing2_
Distribution of the signal--to--noise ratio under the null hypothesis. The observed signal--to--noise ratio $z= `r tValue`$ lies outside the region [-1.96, 1.96] which provides evidence that the signal is bigger than expected by chance and the null hypothesis should be rejected.


## The test statistic and the null distribution

In the previous paragraph we have, in passing, introduced very important concepts that deserve their own section, namely the test statistic and the null distribution. Both are crucial ingredients in hypothesis testing so let's mention them explicitly!

In our example, the test statistic was $(\bar{X} - \mu_0)/SEM$. Under the null hypothesis, we know that it follows a standard normal distribution which is the null distribution in this case. In more general terms we can write:

>The test statistic of a hypothesis test measures the discrepancy between what you observe and what you expect to see under the null hypothesis. It has this general form^[This is true for any hypothesis test, including all the practical examples that we'll talk about in the next chapter!]

$$\frac{\mbox{(observed value)} - H_0}{\mbox{Scale factor}} \quad = \quad \frac{\mbox{signal}}{\mbox{noise}}.$$

Under $H_0$ it follows a distribution (called the _null distribution_) that assigns a lot of probability to the region around zero, because if $H_0$ is true, what you observe shouldn't be all that different from zero. Therefore we can conclude that:

> The further away from zero the observed value of the test statistic is, the more evidence the data provide against $H_0$!

This last statement is closely related to the definition of a $p$--value, but before we go there, we have one confession to make.

## The problem of not knowing $\sigma$ -- the t-distribution

Do you recall the definition of the $SEM$? It's the standard deviation of the population $\sigma$ divided by the square root of your sample size $n$. If you plug this into the generic formula for the test statistic above, we get $(\bar{X} - \mu_0)/SEM = (\bar{X} - \mu_0)/(\sigma/\sqrt{n})$. The problem is, however, that we _never_ know $\sigma$. So far, we have conveniently ignored this fact, but now it's time to face reality. How can we deal with this?

One obvious solution is to estimate it from the sample! There is some good news and some bad news about this approach. Good news first: if we replace $\sigma$ by its estimator $S$ we get
$$\frac{\bar{X} - \mu_0}{S/\sqrt{n}}  \stackrel{H_0}{\sim} t_{n-1},$$
i.e. a test statistic that looks very similar to the one we talked about, but apparently it follows a different distribution. This distribution is the famous $t$--distribution! Actually, it is very similar to the standard Normal distribution. However, there are two key differences to note. The t-distribution is:

<!-- % \item Always \emph{centered} around $0$, like the standard normal distribution. It's values are normalized in a way similar to Z-scores. -->
- Flatter and has \emph{thicker tails} than the standard normal distribution.
<!-- % This means that for an equivalent percent of observations, a wider x-range is necessary. -->
- Influenced by a parameter that is called \emph{degrees of freedom}, which in this case is $n-1$.

```{r tDistShape, fig.cap = 'The shape of the t-distribution is determined by the degrees of freedom.', echo = FALSE, fig.width = large_width, echo = FALSE}

xRange <- 4 # Set the upper limit of the x-axis
samples <- 10
df <- samples - 1
x <- seq(-1*(xRange*2), (xRange*2), length=5000) # Define an x series based on inputed mu and xRange
hstd <- dnorm(x) # Define the density of the standard normal curve at each x position
tdist1 <- dt(x, df = 19) # Define the density of the standard normal curve at each x position
tdist2 <- dt(x, df = 9) # Define the density of the standard normal curve at each x position
tdist3 <- dt(x, df = 4) # Define the density of the standard normal curve at each x position

#    test.df <- data.frame(value = x, std = hstd, tdist = tdist)
#    test.m <- melt(test.df)

test.df <- data.frame(x, std = hstd, tdist1 = tdist1, tdist2 = tdist2, tdist3 = tdist3)
library(reshape2)
test.m <- melt(test.df, id.vars = c("x"))

# Base plot:
ggplot(test.m, aes ( x = x, y = value, col = variable)) +
  geom_line() +
  scale_x_continuous("x", breaks = seq(-xRange, xRange, 1), labels = seq(-xRange, xRange, 1), expand = c(0,0)) +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  coord_cartesian(xlim = c(-xRange, xRange), ylim = c(0, 0.5)) +
 # ggtitle(paste("The t-distribution for",df,"degrees of freedom\nand the standard normal distribution", sep = " ")) +
  scale_color_manual("Distribution", labels = c("Normal", "t, df = 19", "t, df = 9", "t, df = 4"), values = c("grey80","#F5B319","#005AA5","#00824A")) +

  geom_vline(xintercept = mu + qnorm(0.5 + 0.95/2), colour = "grey80", linetype = 2) +
  geom_vline(xintercept = mu - qnorm(0.5 + 0.95/2), colour = "grey80", linetype = 2) +

  geom_vline(xintercept = mu + qt(0.5 + 0.95/2, df = 19), colour = "#F5B319", linetype = 2) +
  geom_vline(xintercept = mu - qt(0.5 + 0.95/2, df = 19), colour = "#F5B319", linetype = 2) +

  geom_vline(xintercept = mu + qt(0.5 + 0.95/2, df = 9), colour = "#005AA5", linetype = 2) +
  geom_vline(xintercept = mu - qt(0.5 + 0.95/2, df = 9), colour = "#005AA5", linetype = 2) +

  geom_vline(xintercept = mu + qt(0.5 + 0.95/2, df = 4), colour = "#00824A", linetype = 2) +
  geom_vline(xintercept = mu - qt(0.5 + 0.95/2, df = 4), colour = "#00824A", linetype = 2) +


  geom_vline(xintercept = 0, colour = "#C42126") +
  theme(panel.background = element_blank(),
        axis.text.x = element_text(colour = "black", angle=0),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        legend.position = c(0.92, 0.9),
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )

```

The dashed lines mark the interval covering 95\% of the potential X values, centered around $0$. We see that as the df increases, this interval shortens, and the $t$--distribution gets closer and closer to the standard Normal distribution.

The second point means that the shape of the $t$--distribution changes dependent on the sample size. From Figure \@ref(fig:tDistShape), we note that the larger the degrees of freedom (i.e. the larger the sample size) the^[Another way of understanding this is to think about being penalised for having a small sample size. Smaller sample sizes equate to fatter tails and wider 95\% probability ranges.] more the $t$--distribution resembles the standard Normal distribution. This is due to the fact that the extra uncertainty we have from estimating instead of knowing $\sigma$ decreases with a bigger sample size. This extra uncertainty is what is responsible for the fatter tails of the $t$--distribution when we compare it to the standard Normal distribution. Apart from this slight change in distribution, everything we said so far about testing stays the same!

That was the good news, now the bad news: in a strict sense, we can use this nice $t$--distribution only if the variable of interest is Normally distributed in our population. So whenever we are planning to do a $t$--test we need to make sure this underlying assumption is met.^[This is quite important when you have a small sample size. With big sample sizes, we're back to the CLT and we can revert to the Normal distribution! We'll talk more about underlying assumptions when we get to non--parametric tests.]


Let's explore what all the means using the web app: use the t-distribution web app to answer the following questions:


- What are the boundaries of the central interval of the $t$--distribution that encompasses 95\% probability if you have a sample size of $n=$5?:

- What happens to these boundaries if $n$=10?

<!-- % The answer to the last exercise tells us that the probability to observe a $t$--value outside the interval [-2.26, 2.26] is only 5\%. -->

Going back to our Martian height example, we saw that the observed $t$--statistic is $t=(`r Height.mean` - `r HypoValue`)/`r Height.SEM` = `r tValue`$ which is outside the interval [-2.26, 2.26]. This tells us once more that it is very unlikely that our assumption of Martian height being $`r HypoValue`$ cm based on popular culture references was correct. However, we also see that accounting for the extra uncertaincy of not knowing $\sigma$, our result is not quite as clear cut as before, when we compared $(`r Height.mean` - `r HypoValue`)/`r Height.SEM` = `r tValue`$ to the interval [-1.96, 1.96].


Given this new insight, would you now reject $H_0: \mu = 196$ cm? How did you come to this conclusion?

<!-- % **Exercise** write in your own words how and why the sample size influences the shape of the $t$--distribution. \\ -->
<!-- % \noindent \hdashrule[-2.5ex]{10.5cm}{1pt}{1pt} \\ -->
<!-- % \noindent \hdashrule[-2.5ex]{10.5cm}{1pt}{1pt} \\ -->
<!-- % \noindent \hdashrule[-2.5ex]{10.5cm}{1pt}{1pt} \\ -->
<!-- % \\ -->

## The 'p--word'

At this stage we have shown in different ways that our popularly held belief that the true mean Martian height is `r HypoValue` cm, is probably wrong. In particular, we have demonstrated that if this hypothesis were correct, it would be less than 5\% likely to observe 10 Martians that on average have a height of $\bar{x} = `r Height.mean`$. 
While this amount of information might satisfy some people, most researchers are more curious and want to know more: what does less than 5\% mean? Exactly how likely is it then?

The answer to this question can actually be quite easily obtained by looking at the  $t$--distribution again. Just as we can find out that all values outside [-2.26, 2.26] are less than 5\% likely, we could try to find out how likely it is to observe values outside the interval [$`r c(-1*tValue,tValue)`$]. This probability would tell us, how likely it is to observe a sample with the same or more extreme signal to noise ratio, i.e. the same or bigger standardised deviation between the observed mean $\bar{x}$ and the hypothesized value $\mu_0$. This is what the $p$--value is in our case!!

Use the web app to do exactly this.

More generally, the $p$--value is defined as

>The $p$--value is the probability of observing data that show at least as much deviation from the null hypothesis (i.e. that are at least as unlikely to have been observed by chance) as the data we have at hand, assuming that the null hypothesis were true.

This definition explains why small $p$--values result in a rejection of the null hypothesis. It may sound a bit clumpsy though so let's reiterate reiterate things using an example close to home.



## Understanding 'p--values'

Let us suppose that it's late Autumn, so the weather is on the brink of getting very cold. We have a fancy new sweater that we really want to wear, but if it's too cold, we'll have to settle with our ugly old winter jacket. To tackle this serious problem in an appropriate way, we formulate the null hypothesis that it's very cold outside, e.g. $H_0: \mu\leq10$ celsius ^[Note that we use the state that we want to reject (it being cold) as the null hypothesis. This is a convention in hypothesis testing which also has a mathematical reason as we'll see later.]. Since we're scientists, we want to probe this hypothesis by collecting some data. In this case, we look out the window and note what people are wearing. The $p$--value that we then mentally obtain relates what we observe with what we expect under the null hypothesis via the question

> If the null hypothesis, $H_0: \mu\leq10$ celsius, were true, how likely is it to observe what I observe?

Clearly, the likelihood of seeing people in summer clothes when it's $\leq10$ celsius is not very large. Therefore, the more people we see wearing T--shirts, the bigger the deviation between our data and what we expect under $H_0$, the smaller the $p$--value, and the bigger our inclination to reject $H_0$ and leave our ugly jacket at home when we leave the house!

% So imagine, we saw many people in T--shirts. This would not be very consistent with $H_0$, so we would be inclined to reject it and leave our ugly jacket at home when we leave the house. However, and this is an important point, there is still a chance that we were just really unlucky when we collected the data and accidentally the people that we observed just happened to be wearing summer clothes for some reason _despite it being very cold_. Admittedly, that's not very likely, but the probability of accidently rejecting the null hypothesis is not zero! We'll talk about this kind of mistake in a bit. First, let's reiterate what we learned in this section:

<!-- % We have introduced several formal and not so formal definitions of the $p$--value. You can memorize the one that makes most sense to you. The $p$--value: -->

<!-- % \begin{itemize} -->
<!-- % \item is the probability to observe what I observe, or something more extreme, if the null hypothesis were true. -->
<!-- % \item quantifies the evidence in favour of the null hypothesis, i.e. a low $p$--value means very little evidence for the $H_0$. -->
<!-- % \item is the probability to sample values that deviate from what I'd expect under the null hypothesis at least as much as the ones I observed, if the null hypothesis were true. -->
<!-- % \end{itemize} -->

## The significance level, $\alpha$ {#SigAlpha}

If the $p$--value is small, we are inclined to reject the null hypothesis. But how small does it need to be?!
Everyone knows that the usual strategy is to reject the null hypothesis if the $p$--value is less than 0.05. Sticking to this procedure would mean that we accept at most a 5\% chance of accidentally rejecting the $H_0$ despite it being true. However, it is important to point out that 0.05 is not a God--given magic number, and just choosing it out of convention does not make much sense. Instead, we need to make up our own mind as to how important it is to not accidentally reject $H_0$ in our particular study. We therefore need to define a (study--specific) threshold for the $p$--value beyond which we will reject the $H_0$! This threshold is called the _significance level_ and is commonly denoted by $\alpha$.

To understand $\alpha$, let's return to our sweater weather example. Here, we would choose $\alpha$ depending on how much we hate being cold. If we don't mind being cold, then our strategy could be to set a quite liberal (i.e. large) significance level like $\alpha = 0.2$. With this, we don't require too much evidence against $H_0$ before we reject it, a $p$--value $< 0.2$ would be enough. However, if we absolutely detest being cold, then we are not very inclined to reject the $H_0$. In this case, our strategy would be to select a very stringent significance level, e.g. $\alpha = 0.001$, which means we will only reject $H_0$ if we saw lots of people in $T$--shirts and the evidence against it being cold is absolutely compelling. 

% Similarly, if you make a discovery in your research that would make it necessary to re-write the textbooks you'd better be really sure about your results before you submit them to Nature (small $\alpha$). If your results are only supposed to be the starting point for new experiments in your lab you don't need to be quite so certain that they are based on something real (large $\alpha$)!

Similarly, in the Martian height example, the $p$--value threshold we choose beyond which we reject the hypothesis  $H_0: \mu=`r HypoValue`$ cm should depend on how bad an accidental rejection would be. To stick with our slightly made up explanations: if the person who proposed this hypothesis in the first place was a famous and well-respected editor of many important journals in your field, then you want to be really sure before you claim that he's wrong. In this case, you'd choose a stringent significance level $\alpha$ and only reject the null hypothesis for very small $p$-values, i.e. when the evidence against it in your data is very strong.

By the way, rejecting the null hypothesis is what's commonly called 'the test was significant'. The long (and more correct) version of saying this would be to say

> The data showed statistically significant deviation from the null hypothesis of xxx at the yyy \% level.^[Obviously replace xxx with your specific null hypothesis and yyy with your chosen significance level $\alpha$.]


## Summary thus far

We have covered quite a bit of material regarding hypothesis testing already. Before we move on, let's take a little break and summarise what we've learned so far.

Test statistics tell us something about the size of the observed difference between the data and the null hypothesis (i.e. the signal) relative to the uncertainty in the sample (i.e. the noise). The farther away our observations are from $H_0$, the larger (in absolute terms) the value of the test statistic will be, i.e. a large absolute value of the test statistic gives us a hint that the deviance between the data and the null hypothesis is due to some systematic pattern and not just due to unlucky sampling. A standardised way to make a judgement about this is the $p$--value: the larger the value of the test statistic, the more evidence against $H_0$ is contained in the data and the smaller the $p$--value. If it is smaller than a certain threshold $\alpha$, the null hypothesis is rejected based on the data. The significance level, $\alpha$, thus controls the probability that you accidently reject $H_0$, i.e. a false positive result. We'll return to this soon.



## One and two--sided tests

In hypothesis testing, it's not only important to specify the null hypothesis, $H_0$, but also the departure from $H_0$, i.e. the _alternative hypothesis_ $H_1$ which can be thought of as the opposite of $H_0$. Recall that $H_0$ usually describes the state that is _not_ interesting, i.e. the state of boredom. Conversely, $H_1$  matches our _a priori_ desire of what would be interesting to discover in the data. Taken together, $H_0$ and $H_1$ are called the _test problem_, and as we'll see now, it can come in different versions.

In the sweater weather example, temperatures $\leq10$ celsius weren't interesting, since we couldn't wear our sweater. Therefore, the desirable state is the case where the temperature _exceeds_ $10$ celsius. The entire test problem can be thus be formalised as

$$H_0: T \leq  10celsius \quad \quad \mbox{vs.} \quad \quad H_1: T >  10 celsius.$$

In the Martian height example, _any_ departure from $H_0$ is interesting, we have no _a priori_ reason to be more excited about Martians to be taller or shorter than the hypothesised value $`r HypoValue`$ cm. Thus, the test problem can be formalised as

$$ H_0: \mu = `r HypoValue`cm \quad \quad \mbox{vs} \quad \quad H_1: \mu \neq `r HypoValue`cm.$$

These cases demonstrate the difference between one--sided and two--sided test problems. In total, there are three possible scenarios that are summarised in Table _ScenarioDescriptions_.

Table: (\#tab:ScenarioDescriptions) Three scenarios in hypothesis testing and two types of tests: two--sided and one--sided.

| Scenario | $H_0$          | $H_1$           | _a priori_ expectations             | Type       |
|----------|----------------|-----------------|-------------------------------------|------------|
| 1        | $\mu=\mu_0$    | $\mu\neq\mu_0$  | $\bar{x}$ _differs from_ $\mu_0$    | two--sided |
| 2        | $\mu\geq\mu_0$ | $\mu<\mu_0$     | $\bar{x}$ is _less than_ $\mu_0$    | one--sided |
| 3        | $\mu\leq\mu_0$ | $\mu>\mu_0$     | $\bar{x}$ is _greater than_ $\mu_0$ | one--sided |

```{exercise, name="One-sided tests"}
Can you think of a weather problem that is one--sided but in the other direction?
```

```{exercise, name="Two-sided tests"}
Can you think of a two--sided weather problem?
```

```{exercise, name="Scenarios in your research"}
What are one-- and two--sided scenarios in your own research?
```

The distinction between one-- and two--sided test problems is important because it influences the $p$--value. Too see this, let's return to the Martian height example and consider all 3 scenarios listed in Table \@ref(ScenarioDescriptions). The corresponding graphical explanations are shown in Figure \@ref(fig:3ScenCase3).

### Scenario 1: two--sided tests, $H_0: \mu=\mu_0$

We already defined our Martian height situation as being two--sided. That means we are looking at the test problem:

$$H_0: \mu =  `r HypoValue` cm \quad \quad \mbox{vs.} \quad \quad H_1: \mu \neq `r HypoValue` cm.$$

We observed that $\bar{x}>\mu_0$ (i.e. $`r Height.mean`$ cm > $`r HypoValue`$ cm, i.e. look at the right--hand side of the top row of Figure \@ref(fig:3ScenCase3).). Our standardised way of assessing the deviance between the sample and the null hypothesis gives us the test statistic

$$t=(\bar{x}-\mu_0)/SEM = (`r Height.mean`-`r HypoValue`)/`r Height.SEM` = `r tValue`.$$

The dashed red line shows the position of this value on the corresponding null distribution, $t_{`r nValue-1`}$. We find that the area under the curve to the right of `r tValue` (shown in red) is very small. It corresponds to $0.0135$, i.e. if the null hypothesis were true there is only a `r 0.0135*100`\% chance to observe a sample of `r nValue` Martians with a mean height of `r Height.mean` cm or larger (i.e. more extreme). But wait! In the two--sided case, this is only _half_ of the story. The term _more extreme_ refers to deviations from $H_0$ to either side. Therefore, we also need to consider the area under the curve to the left of `r tValue2` which by the way would correspond to an observed $\bar{x}$ of $190$ cm as

$$t=(`r 190`- `r HypoValue`)/ `r Height.SEM` = `r tValue2`.$$
Since our test is two--sided, $\bar{x} = 190$ cm provides _as much evidence against_ $H_0$ as $\bar{x} = `r Height.mean`$ cm. Due to the symmetry of the $t$--distribution we know that the area under the curve to the left of `r tValue2` is also $0.0135$, i.e. if $H_0$ were true, there is only a `r 0.0135*100`\% chance to observe `r nValue` Martians with a mean height of $190$ cm or smaller (i.e. more extreme). Taken together, this means that there is a $2.7$\% chance to observe a sample that provides at least as much evidence against $H_0$ as the one we observed. This is exactly the $p$--value we're interested in in the two--sided scenario. We add up the area under the curve beyond the t-statistic from _both_ directions. As $p < 0.05$, we reject the null hypothesis at the chosen significance level $\alpha = 0.05$.

- Important point number 1: In the two--sided case, exactly the same considerations hold true if we had observed $\bar{x} = 190$ cm which corresponds to the left--hand side of the top row of Figure \@ref(fig:3ScenCase3).

- Important point number 2: Note the position of the dashed blue lines. These the values of the test statistic would result exactly in a $p$--value of that matches the chosen significance level, i.e. $p = 0.05$ in our case. They are called the _critical values_ for $\alpha = 0.05$ because any test statistic beyond these values would yield a significant result, i.e. result a $p$--value less than $0.05$.

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


xRange <- 3 # Set the upper limit of the x-axis
df <- nValue - 1
x <- seq(-1*(xRange*2), (xRange*2), length=5000) # Define an x series based on inputed mu and xRange
tdist1 <- dt(x, df = df) # Define the density of the standard normal curve at each x position

# test.df <- data.frame(x, tdist1 = tdist1)

###### Set case 1:
test.df.case1 <- data.frame(x, tdist1 = tdist1)
test.df.case1$test1 <- "geq"
test.df.case1$test2 <- "leq"

test.df.case1 <- reshape2::melt(test.df.case1, measure.vars = c("test1", "test2"), variable.name = "test")

test.df.case1$ribbon <- 0
test.df.case1$ribbon[test.df.case1$x > tValue & test.df.case1$value == "leq"] <- 1
test.df.case1$ribbon[test.df.case1$x < tValue & test.df.case1$value == "geq"] <- 1

test.df.case1$test <- factor(test.df.case1$value)
test.df.case1$obs <- "greater"

###### Set case 2:
test.df.case2 <- data.frame(x, tdist1 = tdist1)
test.df.case2$test1 <- "leq"
test.df.case2$test2 <- "geq"

test.df.case2 <- reshape2::melt(test.df.case2, measure.vars = c("test1", "test2"), variable.name = "test")

test.df.case2$ribbon <- 0
test.df.case2$ribbon[test.df.case2$x < -1*tValue & test.df.case2$value == "geq"] <- 1
test.df.case2$ribbon[test.df.case2$x > -1*tValue & test.df.case2$value == "leq"] <- 1

test.df.case2$test <- factor(test.df.case2$value)
test.df.case2$obs <- "lesser"

###############
##### case 3:

test.df.case3 <- data.frame(x, tdist1 = tdist1)
test.df.case3$test1 <- "eq"
test.df.case3$test2 <- "neq"

test.df.case3 <- reshape2::melt(test.df.case3, measure.vars = c("test1", "test2"), variable.name = "test")

test.df.case3$ribbon <- 0
test.df.case3$ribbon[test.df.case3$x < -1*tValue & test.df.case3$value == "eq"] <- 1
test.df.case3$ribbon[test.df.case3$x > tValue & test.df.case3$value == "eq"] <- 2

test.df.case3$ribbon[test.df.case3$x < -1*tValue & test.df.case3$value == "neq"] <- 1
test.df.case3$ribbon[test.df.case3$x > tValue & test.df.case3$value == "neq"] <- 2

# test.df.case3$test <- factor(test.df.case3$value)
test.df.case3$test <- "eq"
test.df.case3$test <- factor(test.df.case3$test) # Restablish the factor to drop unused levels.

test.df.case3$obs <- factor(test.df.case3$value)
levels(test.df.case3$obs) <- c("lesser", "greater")

###########################
###########################
# Cumulative plot:

test.df <- rbind(test.df.case3,
                 test.df.case2,
                 test.df.case1)


test.df <- test.df[-4] # remove value column, unnecessary.

test.df$obs <- factor(test.df$obs)
levels(test.df$obs) <- c("Observation:\n$\\bar{x}<\\mu_0$","Observation:\n$\\bar{x}>\\mu_0$")

levels(test.df$test) <- c("Scenario 1:\n$H_0:\\mu=\\mu_0$\n$H_1:\\mu\\neq\\mu_0$",
                          "Scenario 2:\n$H_0:\\mu\\geq\\mu_0$\n$H_1:\\mu<\\mu_0$",
                          "Scenario 3:\n$H_0:\\mu\\leq\\mu_0$\n$H_1:\\mu>\\mu_0$")

test.df.complete <- test.df # save for use later on


base.plot <-
  ggplot(test.df, aes ( x = x, y = tdist1)) +
#     ggtitle("$\\mu\\pm1\\sigma=68.2\\%$") +
    geom_line() +
    scale_x_continuous("t", breaks = seq(-xRange, xRange, 1), labels = seq(-xRange, xRange, 1), expand = c(0,0)) +
    scale_y_continuous("f(x)", expand = c(0,0)) +
    coord_cartesian(xlim = c(-xRange, xRange), ylim = c(0, 0.4)) +
#     geom_vline(xintercept = mu + qt(0.5 + 0.95/2, df = df), colour = DarkRed, linetype = 2) +
#     geom_vline(xintercept = mu - qt(0.5 + 0.95/2, df = df), colour = DarkRed, linetype = 2) +
    theme(panel.background = element_blank(),
          axis.text.x = element_text(colour = "black", angle=0),
          axis.text.y = element_text(colour = "black"),
          axis.ticks.x = element_line(colour = "black"),
          legend.position = c(0.92, 0.9),
          legend.background = element_blank(),
          legend.key = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text.y = element_text(angle = 0),
          strip.background = element_blank(),
          panel.spacing = unit(2, "lines")
    )

drawlines.df <- plyr::ddply(test.df, c("test", "obs"), summarise,
                      x = tValue,
                      tdist1 = 0,
                      test = test[1],
                      value = 0,
                      ribbon = 0,
                      obs = obs[1]
                      )

drawlines.df$x[drawlines.df$obs == "Observation:\n$\\bar{x}<\\mu_0$"] <- -tValue

drawlines.df$p <- c(pValue_2Side,pValue_2Side,
                    pValue_1SideLarge,pValue_1SideSmall,
                    pValue_1SideSmall,pValue_1SideLarge)

drawlines.df$cutoff1 <- c(qt(0.975, df = df), qt(0.975, df = df),
                          qt(0.05, df = df), qt(0.05, df = df),
                          qt(0.95, df = df), qt(0.95, df = df))

drawlines.df$cutoff2 <- c(qt(0.025, df = df),qt(0.025, df = df),
                          NA,NA,
                          NA,NA)

base.plot +
  geom_ribbon(data=subset(test.df, ribbon == 1),
              aes(x=x,ymax=tdist1),ymin=0,fill=MidRed) +
  geom_ribbon(data=subset(test.df, ribbon == 2),
              aes(x=x,ymax=tdist1),ymin=0,fill=MidRed) +
  geom_vline(data = drawlines.df, aes(xintercept = x), colour = DarkRed, linetype = 2) +
  geom_vline(data = drawlines.df, aes(xintercept = cutoff1), colour = DarkBlue, linetype = 2) +
  geom_vline(data = drawlines.df, aes(xintercept = cutoff2), colour = DarkBlue, linetype = 2) +
  geom_text(data = drawlines.df, aes(label = paste("t =",plyr::round_any(x, .01))), x = -2.5, y = 0.38, size = 3, hjust = 0, vjust = 1) +
  geom_text(data = drawlines.df, aes(label = paste("p =",p)), x = -2.5, y = 0.33, size = 3, hjust = 0, vjust = 1) +
  facet_grid(test ~ obs)

```

```{r 3ScenCase3, fig.cap='Two--sided and one--sided tests. ', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/3ScenCase3.png")
```
Three scenarios and two possible outcomes of statistical tests. Each scenario is plotted on a $t_{`r nValue-1`}$ distribution with critical values at $\alpha=0.05$ depicted with dashed blue lines. The $p$--value corresponds to the red area.


### Scenario 2: One--sided tests, $H_0: \mu\geq\mu_0$

If we had some _a priori_ reason why it only would be interesting if the mean Martian height is smaller than what is currently believed, we would have formulated our test problem as

$$H_0: \mu \geq `r HypoValue` cm \quad \quad \mbox{vs.} \quad \quad H_1: \mu < `r HypoValue` cm.$$

This case is presented in the middle row of Figure \@ref(fig:3ScenCase3). Again, the right panel shows our current situation, where $\bar{x}>\mu_0$ (i.e. $`r Height.mean`$ cm > $`r HypoValue`$ cm). The $t$--statistic, the dashed red line, is the same as before: $t=`r tValue`$. However, due to the different nature of the test problem, this situation results in a _very large_ $p$--value compared to the two--sided test: $p = 0.9866$, which is reflected in the red-shaded area under the curve to the left of our $t$-score. This is because our observation, $\bar{x}>\mu_0$, is perfectly in agreement with $H_0: \mu\geq\mu_0$! There is no evidence in the data against the null hypothesis. Unfortunately, this is corresponds to the situation that is not of interest for us.


Conversly, if $\bar{x} = 190$ cm, i.e. $\bar{x}<\mu_0$, the situation is reversed as depicted in the left panel of the middle row of \@ref(fig:3ScenCase3). The $t$--statistic is large and the $p$--value is tiny: $0.0134$, because our data are not in agreement with $H_0$. It would be very unlikely to observe $\bar{x} = 190$ cm in a sample of 10 Martians if $H_0: \mu \geq `r HypoValue`$ cm were true.

^[A convenient way to remember which part of the curve will contribute to the $p$--value is to look at $H_1$. For $H_1: \mu\neq\mu_0$, it's the area beyond the $t$--statistic in _both_ directions. For $H_1: \mu<\mu_0$, it's always that which is _lower than_ the $t$--statistic.] Given that the $p$--value for the previous observation case $\bar{x}>\mu_0$, is $p = 0.9866$, what is the formula for this $p$--value? What is the relationship of this $p$--value to the $p$--value we observed for the two-sided test, i.e. $p=0.027$?

<!-- % You can see how this scenario would only be valid if we really did care about the direction of change. There are lots of subtle differences here, so let's take a look. -->
<!-- %  -->
<!-- % With $H_0: \mu\geq\mu_0$ we only consider the area under the curve to the \emph{left} of the observed $t$--statistic. In the two-sided test, $H_0: \mu=\mu_0$ we took the area under the curve the \emph{left} and to the \emph{right} of the observed $t$--statistic, equidistant to the center, on both sides of the curve.  -->


Before we move on the scenario 3, also take a look at the position of the critical values, i.e. the blue lines. In a one-sided test, there's only one critical value but as before it's positioned on the side of the curve in accordance with $H_1$, i.e. in a region where we'd reject $H_0$. 

```{exercise, name="Position of the critical values"}
Why are the critical values in the one--sided test closer to zero than int he two--sided scenario?
```

<!-- % Not only does the critical value get smaller in one-sided tests, but the $p$--value is halved! This makes sense, since we are only interested in one direction, but you can imagine that this can sometimes decide if a test is "significant" or not. -->
<!-- %  -->
<!-- % Sometimes these considerations doesn't make a big difference - but they may, so you should be aware of them! It's another good reason why we need to have some \emph{a priori} idea about the direction of our test. -->

### Scenario 3: One--sided tests, $H_0: \mu\leq\mu_0$

This situation is exactly the same as in scenario 2, except that everything is reversed. In this situation, we have the test problem

$$ H_0: \mu \leq `r HypoValue` \quad \quad \mbox{vs.} \quad \quad H_1: \mu > `r HypoValue`,$$

which is depicted as scenario 3 in the bottom row of Figure \@ref(fig:3ScenCase3).

Notice that the critical values, the dashed blue lines, and the area under the curve contributing to the $p$--value are on the right side of the curve. This reflects the direction of the test, since only high $t$--statistics will result in a "significant" result.
Our current situation, $\bar{x}>\mu_0$, shown in the right panel, speaks against $H_0$ and results in a small $p$-value. This is _exactly_ the same as the left panel of scenario 2, just on the opposite side!
The converse outcome, $\bar{x}<\mu_0$, shown in the left panel, is in agreement with $H_0$ and results in a very large $p$--value. This is _exactly_ the same as the right panel of scenario 2, just on the opposite side!

## Errors in hypothesis testing

As mentioned above, the significance level $\alpha$ controls the probability to accidently reject $H_0$. This is one of two types of error that can occur in hypothesis testing. The two error types can have different names depending on the text book you look at which can sometimes be a source of confusion, so let's spell them out in detail:

<dl>
  <dt>Type I error</dt>
  <dd>also known as the $\alpha$--error, corresponds to rejecting the null hypothesis $H_0$ although it is true. Type I errors are referred to as _false positives_.^[Note the confusing use of the term 'positive' here. A _positive_ result is we _reject_ the null hypothesis. This is because usually $H_0$ is the state of boredom you want to reject. By the same token, a _negative_ result is _not_ being able to reject the null hypothesis -- boring!] </dd>

  <dt>Type II error</dt>
  <dd>also known as the $\beta$--error, corresponds to not rejecting the null hypothesis $H_0$ although it is false. Type II errors are referred to as _false negatives_.</dd>
  </dl>

<!-- % As we saw, type I errors can be managed by choosing an appropriate significance level for our test, i.e. a threshold for the $p$--value beyond which we reject the null hypothesis. You can probably guess that there is a flip-side to type I errors: -->

The two error types above are flip--sides of each other. Actually, they are two out of four possible outcomes that we can get in hypothesis testing as shown in Figure \@ref(fig:fouroptions).

```{r fouroptions, fig.cap='The four possible outcomes of hypotheses testing -- there are two ways to get it right and two possible pitfalls.', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/GridError.png")
```


To understand the two types of error, we'll return to the sweater weather example. Here, the null hypothesis is $H_0: \mu\leq10$ celsius, i.e. the state we hope to disprove. The null hypothesis can be either true or false, this is the underlying reality that is unknown to us and that we want to find out about by collecting the data. Based on data, we get a test result, i.e. we can either reject the null hypothesis or not. In total, we get the four outcomes shown in Figure \@ref(fig:fouroptionssweater).

```{r fouroptionssweater, fig.cap='The four types of outcomes in the sweater weather example.', echo=FALSE}
knitr::include_graphics("figures/GridSweater.png")
```

The desired options are show in green, and our most favourite is the true positive where it indeed is warm and we leave the ugly coat inside.

Note that we deliberately _never_ say that we _prove_ $H_0$. In hypothesis testing, we either _reject_ or _cannot reject_ the null hypothesis $H_0$, but we never prove it. Why is that? We give the answer to this at the end of the chapter!

<dl>
<dt>A false positive</dt>
<dd>here would be that we accidently reject $H_0$ based on the data, i.e. we think we have seen a sufficient number of people in summer clothes to decide that it's probably not so cold, so we walk outside with a sweater and get a bad surprise because in fact $H_0$ is true and we and up freezing. In a research setting, a false positive is akin to thinking you have a fantastic new result, when you actually don't!</dd>

<dt>A true positive</dt>
<dd>is our preferred outcome, i.e. we reject $H_0$ wear our sweater and it's nice Autumn weather.</dd>

  <dt>A true negative</dt>
  <dd>would be OK as well, it is the situation where we don't find enough evidence against the hypothesis that it's cold and decide to wear a coat when indeed it is cold. We have to wear the ugly coat but at least we are not freezing.</dd>
  
  <dt>A false negative</dt>
  <dd>Here would be the situation where we reject $H_0$ based on our data, wear the ugly coat although actually it's warm outside and we and up sweating all day.</dd>
</dl>


In the sweater weather example, both types of error are not desirable because you end up either freezing or sweating. So how do we deal with the errors?! In hypothesis testing, the error that you actively deal with is the type I error, by choosing your significance level $\alpha$. Adhering to your choice will make sure that the probability of a type I error is not going to be bigger than $\alpha$. However, if you choose a small $\alpha$, you need to be aware that this increases the chance of a type II error! \@ref(pg:typeIandII) So there is a trade--off between them.

```{exercise, name="Types of error bars"}
In the sweater weather example, what would happen if you chose a very stringent (i.e. small) significance level?
```



## Power

There is only one more term that we need to explain before you know the most relevant vocabulary in hypothesis testing, and this is _power_. As you the name suggests, power is something that we want in a test, so let's talk about how we can get lots of it! So far we have mainly talked about type I error, i.e. the error that corresponds to the left--hand side of Figure \@ref(fig:fig:fouroptions) when $H_0$ is true. But what happens on the other side, when $H_0$ is false? Either we incorrectly fail to reject $H_0$ (type II error) or we correctly reject $H_0$. The convention is to denote the probability of a type II error with $\beta$ which directly leads to a probability of $1-\beta$ for a correct rejection of $H_0$ -- this is the famous power of a test:

>_Statistical power_ is the probability of rejecting a false null hypothesis $H_0$. This corresponds to the probability of a _true positive_. In other words, _it is the ability to detect an effect when it is there_.

<!-- % \item[Type I error ($\alpha$)] The probability of rejecting a true null hypothesis $H_0$, i.e. \emph{false positives}. -->
<!-- % -->
<!-- % \item[Type II error ($\beta$)] The probability of failing to reject a false null hypothesis $H_0$, i.e. \emph{false negatives}. -->

So for a given test, how do we determine the probability of the type II error, and thus the power? Is there is procedure to control for it in a way that is similar to setting a significance level which sets an upper limit on the probability of conducting a type I error? Unfortunately the answer is no -- $\beta$ cannot be controlled for directly. It is indirectly determined by three factors: the significant level $\alpha$, the effect size of interest, and the sample size.

In the following, we will demonstrate these dependencies using the Martian height example, and to make it easier to understand we'll just stick to the one--sided test problem
$$ H_0: \mu \leq `r HypoValue`cm \quad \quad \mbox{vs} \quad \quad H_1: \mu > `r HypoValue`cm.$$
The top panel of Figure \@ref(fig:PowerOneSided) shows the null distribution for this test problem.^[This is just a repetition of Scenario 3 from Figure \@ref(fig:3ScenCase3), but we're going to make things easier by simply using the Normal distribution.] It is centered around the hypothesised value. 
The additional panels below correspond to different versions of a situation where the alternative hypothesis is true, namely to $\mu = 197$ cm, $\mu = 200$ cm, and $\mu = 203$ cm, respectively. In these situations, we would like to reject the null hypothesis.


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

# Make one faceted plot for all power options:

Height.SEM <- (Height.sd/sqrt(nValue))

# power.t.test(n = 10,
#              delta = 5,
#              sd = 6,
#              sig.level = .05,
#              type = "one.sample",
#              alternative = "one.sided")



means <- c(195, 197, 200, 203)
labels <- c("Null", "Scenario 1", "Scenario 2", "Scenario 3")

xRange <- 4# Set the upper limit of the x-axis
x <- seq(-1*(xRange*4), (xRange*5)+4, length=10000) + means[1] # Define an x series based on inputed mu and xRange

# Make distributions:
dh <- do.call(dnorm, list(x = x, mean = means, sd = Height.SEM))

# Calculate cut off value for null distributuion
es <- qnorm(0.95,mean = means[1], sd = Height.SEM) # 2   # observed effect (deviation from mean under H0), i.e. alpha

# collect all distributions into a data frame:
df.power <- data.frame(variable = as.vector(t(mapply(rep, labels, each = 10000))),
                       xval = x,
                       value = dh)

# Set colours for ribbons later on, first null:
df.power$fill <- DarkYellow # Normal part 1, yellow
df.power$fill[df.power$variable == "Null" & df.power$xval > es] <- DarkRed # Normal part 2, red

# then the rest of the curves:
df.power <- rbind(df.power[df.power$variable == "Null",],
                  plyr::ddply(df.power[df.power$variable != "Null",], c("variable"), transform,
                        fill = ifelse(xval <= es, DarkBlue, DarkPurple)))

# Make nice labels
alpha <- pnorm(es, means[1], Height.SEM, lower.tail=FALSE)


# Labels for legends:
value0 <- 100-round(alpha*100, 1) #### true negatives
value0 <- paste(value0, "\\% true negatives", sep = "")

value1 <- round(alpha*100, 1) #### Alpha, Type I Error, false positives
value1 <- paste("$\\alpha=$ ", value1, "\\% - Type I Error, false positives", sep = "")

value2 <- paste("$\\beta$, Type II Error, false negatives", sep = "")
value3 <- paste("Power, $1-\\beta$",", true positives", sep = "")

greeks <- c(value0, value1, value2, value3)

levels(df.power$variable) <- c("Null\nDistribution",
                               "Scenario 1:\nEffect size = $2cm$\nLower power",
                               "Scenario 2:\nEffect size = $5cm$",
                               "Scenario 3:\nEffect size = $8cm$\nHigher power")


################## Labels for beta and power:
# names(df.power)
# "variable" "xval"     "value"    "fill"

# Different beta values:
beta <- pnorm(es, means[2:4], Height.SEM)
power <- 1 - beta

# value2 <- round(beta*100, 1) #### Type I Error, false positives
# value2 <- paste("$\\beta=$", value2, "\\% - Type II Error, false negatives", sep = "")

# value3 <- round((1-beta)*100, 1) #### Power, 1-beta, true positives
# value3 <- paste(value3, "\\% Power, $1-\\beta$",", true positives", sep = "")

drawlines.df <- data.frame(variable = levels(df.power$variable)[2:4],
                           beta = beta,
                           power = power)

drawlines.df$xval <- means[2:4]
drawlines.df$value = 0.2
drawlines.df$fill = "#f3b21a"

##############################################

# The four curves, coloured
g <- ggplot(df.power, aes (x = xval, y = value, group = fill, fill = fill, ymax = value)) +
  geom_ribbon(ymin=0, alpha = 0.4) +
  facet_grid(variable ~ .) +

  geom_text(data = drawlines.df, aes(label = paste("$\\beta=$", round(beta*100,2), "\\%", sep = "")), x = 190, y = 0.2, size = 3, hjust = 0, vjust = 1) +
  geom_text(data = drawlines.df, aes(label = paste(round(power*100, 2), "\\% Power", sep = "")), x = 190, y = 0.16, size = 3, hjust = 0, vjust = 1) +

  geom_vline(xintercept = es, colour = DarkRed, linetype = 2) + # Draw effect size line.
  geom_vline(xintercept = means[1], colour = DarkRed) + # Draw mean line.
  geom_vline(data = drawlines.df, aes(xintercept = xval), colour = DarkBlue) +

  scale_x_continuous("x", breaks = seq(189, 205, 1), expand = c(0,0)) +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  scale_fill_identity("",
                      breaks = c(DarkPurple, DarkBlue, DarkRed, DarkYellow),
                      guide = "legend",
                      labels = rev(greeks) ) +

  coord_cartesian(xlim = c(189, 205), ylim = c(0, 0.3)) +

  guides(fill = guide_legend(reverse=TRUE)) +
  theme(panel.background = element_blank(),
        axis.text.x = element_text(colour = "black", angle=0),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        legend.position = c(0.6, 0.95),
        legend.text = element_text(size = 10),
        legend.key.size = unit(0.55, "cm"),
        # legend.direction = "vertical",
        # legend.justification = "left",
        legend.justification = "left",
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.y = element_text(angle = 0),
        strip.background = element_blank()
    )

```



```{r PowerOneSided, fig.cap='Calculating power.', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/PowerScenarios.png")
```
_top panel_: The null distribution is what we expect to observe if the null hypothesis were true, it is therefore centered around the hypothesised $\mu$ (solid red line). The dotted red line marks the critial value, which is determined by $\alpha$. _lower panels_: The distributions depict what we expect to observe in different situations where the alternative hypothesis is true. Note that these situations exhibit a sequentially increasing departure from $H_0$. In each case, the area of the region beyond the critical value corresponds the probability of obtaining a true positive result, i.e. the power of our test under these circumstances.



### Dependence on the significance level, $\alpha$

In Figure _PowerOneSided_, the significance level is chosen as $\alpha = 0.05$. This choice determines the critical value that is displayed by the dashed line. Let's call this $x_{crit}$. This critical value separates the null distribution (top panel) into two parts, yellow and red. The yellow part makes up 95\% of the area under the curve. Conversely, the red part makes up 5\% of the area under the curve. Whether or not we reject the null hypothesis, depends on the observed $\bar{x}$. We reject if $\bar{x} > x_{crit}$ which -- by design of the critical value -- happends in only 5\% of the cases if the null hypothesis is true. 

Now consider the bottom panels which depict situations where $H_1$ is true and we _want_ to reject the null hypothesis. The probability of this is shown in purple in these situations. It corresponds to the power respective power. So what happens to the power if we change $\alpha$? Imagine, we decrease it to say $\alpha = 0.01$. This moves the critical value to the right, and will decrease the area of the purple part in all bottom panels! This is the trade--off between type I and type II error that we previously talked about. If we decrease $\alpha$, we decrease the probability of a type I error because the evidence in the data against the $H_0$ needs to be very large before we reject $H_0$. However, this also means that the probability to _not_ rejecting $H_0$, although it is false, can be quite big. 


### Dependence on the effect size

The effect size is defined as the departure of reality from the null hypothesis. The bottom panels of Figure \@ref(fig:PowerOneSided) therefore show three different effect sizes, namely 2 cm, 5 cm, and 8 cm, respectively. We can see that the power (the purple area under the curves) changes depending on the position of the bottom curve which in turn is governed by the effect size. The bigger the effect size, the farther the distribution of $\bar{X}$ moves away from the null distribution and power increases. In other words, the bigger the true $\mu$, the more likely it is to observe a $\bar{x}>x_{crit}$ in which case we reject $H_0$.

To demonstrate this further, imagine the extreme case where $\mu = 195.0001cm$. Although we are still in the case where we'd like to reject the null hypothesesis, the signal $\mu - \mu_0$ is very small, so even without math we intuitively know that this effect will be hard to detect. Visually, we can also confirm this intuition by imagining the lower panel curves in Figure \@ref(fig:PowerOneSided) moving to the left until they are centered around $\mu = 195.0001cm$. In this case the blue area will grow a lot and with it the probabilty to conduct a type II error. Conversely the power decreases!

This raises the question: if power changes depending on the position of the true $\mu$ (which is unknown to us!!), how can it be calculated? It can't unfortunately, but to rescue the situation, we use a worst--case--scenario where we consider the _minimal effect size_ that we would like to be able to detect and caclulate the _maximally possible_ power. In our example, imagine that we only consider our finding exciting if the true Martian mean height is at least $2$ cm greater than $\mu_0=195$ cm. This means that we would carry out our power calcuations based on Scenario 1 where $\mu = 197$ cm.

### Dependence on the sample size

Let's stay with the Scenario 1 in Figure \@ref(fig:PowerOneSided), where the effect size is 2 cm and power is $`r power[1]`$. What if we are not happy with this value and we would like to have a greater chance to detect a signal of 2 cm with our test? Say we would like to a power of $`r power[2]`$ instead (which is the power of Scenario 2 where the effect size is 5 cm)! This is where sample size comes in. To see this, remember our good old friend the CLT - the $SEM$ decreases as sample size increases.^[You can see this in very large studies, where tiny effects are considered significant, primarily because the sample size is really big.] This ensures that you can estimate $\mu$ with greater precision, and in terms of Figure \@ref(fig:PowerOneSided) it means that the overlap between the top and the bottom curves decreases although the effect size (i.e. the distance between the centers of the curves) stays the same.


In Figure \@ref(fig:PowerOneSided), we calculated the shown power values using a sample size of $n=10$. To achieve the desired power $`r power[2]`$ for detecting an effect size of 2 cm, we'd need to have a sample size of 62! This is shown in  Figure \@ref(fig:PowerHighN). Note how high and narrow the distributions are which reduces the overlap between the two curves. Also notice that the critical value has shifted.


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

# Make one faceted plot for all power options:

newN <- power.t.test(delta = 2,
             sd = 6,
             sig.level = .05,
             power = .839,
             type = "one.sample",
             alternative = "one.sided")
newN <- round(newN$n)

# 65, 0.8514252
# With n = 64, power = 0.8465653
# 63, power = 0.8415619
# 62, power = 0.8364114

Height.SEM <- (Height.sd/sqrt(62))

means <- c(195, 197)
labels <- c("Null", "Scenario 1")

xRange <- 4# Set the upper limit of the x-axis
x <- seq(-1*(xRange*4), (xRange*5)+4, length=10000) + means[1] # Define an x series based on inputed mu and xRange

# Make distributions:
dh <- do.call(dnorm, list(x = x, mean = means, sd = Height.SEM))

# Calculate cut off value for null distributuion
es <- qnorm(0.95,mean = means[1], sd = Height.SEM) # 2   # observed effect (deviation from mean under H0), i.e. alpha

# collect all distributions into a data frame:
df.power <- data.frame(variable = as.vector(t(mapply(rep, labels, each = 10000))),
                       xval = x,
                       value = dh)

# Set colours for ribbons later on, first null:
df.power$fill <- DarkYellow # Normal part 1, yellow
df.power$fill[df.power$variable == "Null" & df.power$xval > es] <- DarkRed # Normal part 2, red

# then the rest of the curves:
df.power <- rbind(df.power[df.power$variable == "Null",],
                  plyr::ddply(df.power[df.power$variable != "Null",], c("variable"), transform,
                        fill = ifelse(xval <= es, DarkBlue, DarkPurple)))

# Make nice labels
alpha <- pnorm(es, means[1], Height.SEM, lower.tail=FALSE)


# Labels for legends:
value0 <- 100-round(alpha*100, 1) #### true negatives
value0 <- paste(value0, "\\% true negatives", sep = "")

value1 <- round(alpha*100, 1) #### Alpha, Type I Error, false positives
value1 <- paste("$\\alpha=$ ", value1, "\\% - Type I Error, false positives", sep = "")

value2 <- paste("$\\beta$, Type II Error, false negatives", sep = "")
value3 <- paste("Power, $1-\\beta$",", true positives", sep = "")

greeks <- c(value0, value1, value2, value3)

levels(df.power$variable) <- c("Null\nDistribution",
                               "Scenario 1:\nEffect size = $2cm$\nn = $62$\nIncreased power")


################## Labels for beta and power:
# names(df.power)
# "variable" "xval"     "value"    "fill"

# Different beta values:
beta <- pnorm(es, means[2], Height.SEM)
power <- 1 - beta

# value2 <- round(beta*100, 1) #### Type I Error, false positives
# value2 <- paste("$\\beta=$", value2, "\\% - Type II Error, false negatives", sep = "")

# value3 <- round((1-beta)*100, 1) #### Power, 1-beta, true positives
# value3 <- paste(value3, "\\% Power, $1-\\beta$",", true positives", sep = "")

drawlines.df <- data.frame(variable = levels(df.power$variable)[2],
                           beta = beta,
                           power = power)

drawlines.df$xval <- means[2]
drawlines.df$value = 0.2
drawlines.df$fill = "#f3b21a"

##############################################

# The four curves, coloured
g <- ggplot(df.power, aes (x = xval, y = value, group = fill, fill = fill, ymax = value)) +
  geom_ribbon(ymin=0, alpha = 0.4) +
  facet_grid(variable ~ .) +

  geom_text(data = drawlines.df, aes(label = paste("$\\beta=$", round(beta*100,2), "\\%", sep = "")), x = 190, y = 0.4, size = 3, hjust = 0, vjust = 1) +
  geom_text(data = drawlines.df, aes(label = paste(round(power*100, 2), "\\% Power", sep = "")), x = 190, y = 0.3, size = 3, hjust = 0, vjust = 1) +

  geom_vline(xintercept = es, colour = "#C42126", linetype = 2) + # Draw effect size line.
  geom_vline(xintercept = means[1], colour = DarkRed) + # Draw mean line.
  geom_vline(data = drawlines.df, aes(xintercept = xval), colour = DarkBlue) +

  scale_x_continuous("x", breaks = seq(189, 205, 1), expand = c(0,0)) +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  scale_fill_identity("",
                      breaks = c(DarkPurple, DarkBlue, DarkRed, DarkYellow),
                      guide = "legend",
                      labels = rev(greeks) ) +

  coord_cartesian(xlim = c(189, 205), ylim = c(0, 0.6)) +

  guides(fill = guide_legend(reverse=TRUE)) +

  theme(panel.background = element_blank(),
        axis.text.x = element_text(colour = "black", angle=0),
        axis.text.y = element_text(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        legend.position = c(0.6, 0.95),
        legend.text = element_text(size = 10),
        legend.key.size = unit(0.55, "cm"),
        # legend.direction = "vertical",
        # legend.justification = "left",
        legend.justification = "left",
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.y = element_text(angle = 0),
        strip.background = element_blank()
    )

```


```{r PowerHighN, fig.cap='Calculating power.', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/PowerLargeN.png")
```
_top_: The null distribution with $n = 62$. _bottom_: The disribution of $\bar{X}$ with an effect size of 2 cm, i.e. $\mu = 197$ cm. Compared to Figure \@ref(fig:PowerOneSided) the curves are very narrow and don't overlap as much which increases the power.

```{exercise, name="Influencers of power"}
Given that the sample size influences power via the $SEM$, can you derive how population variance influences power?
```

In summary, the coloured partition of the curves in Figures \@ref(fig:PowerOneSided) and \@ref(fig:PowerHighN) is dependent on four factors. Figure \@ref(fig:fig:PowerGrid) gives an overview about the effect of these factors on the power!

<!-- \includegraphics[scale = 0.5]{./figures/PowerGridAlpha} -->

```{r PowerGrid, fig.cap='Calculating power.', out.width='100%', fig.align='center', echo=FALSE}
knitr::include_graphics("figures/PowerGrid.png")
```

The power of a statistical test is dependent on four factors: the sample size, the population variance, the significance level, and the effect size under investigation. For a given significance level, the highest power is achieved with a large sample size, a small population variance and large effect sizes. Here, the effect of significance level, effect size and sample size are depicted for a fixed value of the population variance. A qualitative colour scale depicting power under different conditions.



<!-- % \begin{description} -->
<!-- % \item[Effect size] Is the difference between the location of two distributions we are interested in detecting. Here we want to detect any effect that is at least $5cm$. -->
<!-- % \item[Sample size] Directly influences the width of the distributions, since the curves are defined as $N(\Sexpr{HypoValue}, SEM^2)$ and $N(\Sexpr{Height.mean}, SEM^2)$. -->
<!-- % \item[Population variance] like the sample size it determines the width of the distribution. -->
<!-- % \item[Significance level, $\alpha$] Determines the position of our critical value on the null distribution, the dotted red line in figure \ref{PowerOneSided}. -->
<!-- % \end{description} -->



<!-- % Sample size increases power by reducing the standard error of the sampling distribution used in the test. -->

<!-- % \begin{figure*} -->
<!-- % \begin{center} -->
<!-- % \includegraphics[scale = 0.5]{./figures/PowerGrid} -->
<!-- % \end{center} -->
<!-- % \caption[sam var]{The power of a statistical test is dependent on three factors: the sample size, the population variance and the effect size under investigation. The highest power is achieved with a large sample size, small population variance and large effect sizes. A qualitative colour scale depicting power under different conditions.} -->
<!-- % \label{fig:PowerGrid} -->
<!-- % \end{figure*} -->

## Sample size calculation in hypothesis testing

So how can we determine the sample size we need to detect a certain effect size with a certain power? For different tests, the formula comes with different flavours. However, in general, we can use the following recipe:

- Choose the significance level, i.e. the maximum type I error rate that you're willing to accept in your study.
- Choose the minimal effect size that you are interested in.
- Identify the sources of variation in your experiment and make a best-guess of their magnitude, either from previous experiments or from literature.
- Choose the power that you want your experiment to have.
- Calculate the sample size that would allow you to achieve this power, given the first three points.

<!-- %' \begin{figure} -->
<!-- %' <<Draw-alpha, eval = TRUE, echo = FALSE, fig.height = 3.5, fig.width = 7>>= -->
<!-- %' -->
<!-- %' mha <- Height.mean # mean under the alternative hypothesis -->
<!-- %' mho <- HypoValue # mean under the null hypothesis -->
<!-- %' # mha2 <- HypoValue - (Height.mean - HypoValue) # Mean of second direction -->
<!-- %' -->
<!-- %' xRange <- 4# Set the upper limit of the x-axis -->
<!-- %' x <- seq(-1*(xRange*4), (xRange*5)+4, length=5000) + mho # Define an x series based on inputed mu and xRange -->
<!-- %' # x <- seq(180, 200, length=5000) + mho # Define an x series based on inputed mu and xRange -->
<!-- %' -->
<!-- %' # summary(x) -->
<!-- %' dh0 <- dnorm(x, mho, Height.SEM) # Draw normal curve under the null hypothesis, mu = mho, sigma = SEM -->
<!-- %' dh1 <- dnorm(x, mha, Height.SEM) # Draw alternative hypothesis, sigma = Height.sd -->
<!-- %' # dh2 <- dnorm(x, mha2, Height.SEM) # Draw alternative hypothesis, sigma = Height.sd -->
<!-- %' -->
<!-- %' -->
<!-- %' es <- qnorm(0.975,mean = mho, sd = Height.SEM) # 2   # observed effect (deviation from mean under H0), i.e. alpha -->
<!-- %' es2 <- qnorm(0.025,mean = mho, sd = Height.SEM) # 2   # observed effect (deviation from mean under H0), i.e. alpha -->
<!-- %' -->
<!-- %' -->
<!-- %' # dh0 <- dt(x, df = 59) # Draw t-dist, df = 59, mu = 0, sigma = 1 -->
<!-- %' # dh1 <- dt(x, df = 59) # Draw alternative hypothesis, sigma = 1 -->
<!-- %' -->
<!-- %' # library(reshape2) -->
<!-- %' df.power <- data.frame(xval = x, dh0 = dh0, dh1 = dh1) -->
<!-- %' df.power <- reshape2::melt(df.power, id.var = "xval") -->
<!-- %' #names(df.power) -->
<!-- %' -->
<!-- %' df.power$fill <- "#F5B319" # Normal part 1, yellow -->
<!-- %' df.power$fill[df.power$variable == "dh0" & df.power$xval > es] <- "#C42126" # Normal part 2, red -->
<!-- %' df.power$fill[df.power$variable == "dh0" & df.power$xval < es2] <- "#C42126" # Normal part 2, red -->
<!-- %' -->
<!-- %' -->
<!-- %' df.power$fill[df.power$variable == "dh1" & df.power$xval < es] <- "#005AA5" # Ha part 1, blue -->
<!-- %' df.power$fill[df.power$variable == "dh1" & df.power$xval > es] <- "#662D91" # Ha part 2, purple -->
<!-- %' -->
<!-- %' alpha <- pnorm(es, mho, Height.SEM, lower.tail=FALSE) -->
<!-- %' alpha2 <- pnorm(es2, mho, Height.SEM, lower.tail=FALSE) -->
<!-- %' beta <- pnorm(es, mha, Height.SEM) -->
<!-- %' -->
<!-- %' # alpha <- pt(es, df = 0.95, lower.tail=FALSE) -->
<!-- %' # beta <- pt(es, mha, 1) -->
<!-- %' -->
<!-- %' labs.h <- "\\% true negatives" -->
<!-- %' labs.h0 <- "\\% - Type I Error, false positives" -->
<!-- %' labs.h1a <- "\\% - Type II Error, false negatives" -->
<!-- %' labs.h1b <- "\\% Power, " -->
<!-- %' -->
<!-- %' value0 <- 100-round((alpha*2)*100, 1) -->
<!-- %' value0 <- paste(value0, labs.h, sep = "") -->
<!-- %' -->
<!-- %' value1 <- round((alpha*2)*100, 1) -->
<!-- %' value1 <- paste("$\\alpha=$", value1, labs.h0, sep = "") -->
<!-- %' -->
<!-- %' value2 <- round(beta*100, 1) -->
<!-- %' value2 <- paste("$\\beta=$", value2, labs.h1a, sep = "") -->
<!-- %' -->
<!-- %' value3 <- round((1-beta)*100, 1) -->
<!-- %' value3 <- paste(value3, labs.h1b,"$1-\\beta$",", true positives", sep = "") -->
<!-- %' -->
<!-- %' greeks <- c(value1, value2, value3) -->
<!-- %' -->
<!-- %' #q <- ggplot(subset(df.power, df.power$variable == "dh0"), aes ( x= xval, y = value)) + -->
<!-- %' #q <- ggplot(df.power, aes (x= xval, y = value, group = fill, fill = fill, ymax = value)) + -->
<!-- %' -->
<!-- %' ggplot(subset(df.power, df.power$variable == "dh0"), aes (x= xval, y = value, group = fill, fill = fill, ymax = value)) + -->
<!-- %'   scale_x_continuous("Martian height (cm)", breaks = seq(189, 205, 1), expand = c(0,0)) + -->
<!-- %'   scale_y_continuous("f(x)", expand = c(0,0)) + -->
<!-- %'   ggtitle("Null distribution") + -->
<!-- %'   coord_cartesian(xlim = c(189, 205), ylim = c(0, 0.6)) + -->
<!-- %'   theme(panel.background = element_blank(), -->
<!-- %'         plot.background = element_blank(), -->
<!-- %'         legend.background = element_blank(), -->
<!-- %'         legend.key = element_blank(), -->
<!-- %'         axis.text = element_text(colour = "black"), -->
<!-- %'         axis.ticks.y = element_line(colour = "black"), -->
<!-- %'         axis.ticks.x = element_line(colour = "black"), -->
<!-- %'         panel.grid.major = element_blank(), -->
<!-- %'         panel.grid.minor = element_blank(), -->
<!-- %'         legend.position = c(0.05, 0.9), -->
<!-- %'         legend.justification = "left" -->
<!-- %'   ) + -->
<!-- %'   scale_fill_identity("", guide = "legend", -->
<!-- %'                       labels = c(greeks[1], # Label for Normal 2, alpha -->
<!-- %'                                  value0 # Label for Normal 1 -->
<!-- %'                       ) #close c() -->
<!-- %'   ) + # close scale_fill_identity() -->
<!-- %'   geom_ribbon(ymin=0) + # , alpha = 0.4 -->
<!-- %'   guides(fill = guide_legend(reverse=TRUE)) + -->
<!-- %'   geom_vline(xintercept = es, colour = "#C42126", linetype = 2) + # Draw effect size line. -->
<!-- %'   geom_vline(xintercept = mho, colour = "#C42126") # Draw mean line. -->
<!-- %' -->
<!-- %' -->
<!-- %' ggplot(subset(df.power, df.power$variable == "dh1"), aes (x= xval, y = value, group = fill, fill = fill, ymax = value)) + -->
<!-- %'   scale_x_continuous("Martian height (cm)", breaks = seq(189, 205, 1), expand = c(0,0)) + -->
<!-- %'   scale_y_continuous("f(x)", expand = c(0,0)) + -->
<!-- %'   ggtitle("Observed distribution") + -->
<!-- %'   coord_cartesian(xlim = c(189, 205), ylim = c(0, 0.6)) + -->
<!-- %'   theme(panel.background = element_blank(), -->
<!-- %'         plot.background = element_blank(), -->
<!-- %'         legend.background = element_blank(), -->
<!-- %'         legend.key = element_blank(), -->
<!-- %'         axis.text = element_text(colour = "black"), -->
<!-- %'         axis.ticks.y = element_line(colour = "black"), -->
<!-- %'         axis.ticks.x = element_line(colour = "black"), -->
<!-- %'         panel.grid.major = element_blank(), -->
<!-- %'         panel.grid.minor = element_blank(), -->
<!-- %'         legend.position = c(0.05, 0.9), -->
<!-- %'         legend.justification = "left" -->
<!-- %'   ) + -->
<!-- %'   scale_fill_identity("", guide = "legend", -->
<!-- %'                       labels = c(greeks[2:3] # Label for Normal 2, alpha -->
<!-- %'                       ) #close c() -->
<!-- %'   ) + # close scale_fill_identity() -->
<!-- %'   geom_ribbon(ymin=0) +  # , alpha = 0.4 -->
<!-- %'   geom_vline(xintercept = es, colour = "#C42126", linetype = 2) + # Draw effect size line. -->
<!-- %'   geom_vline(xintercept = mha, colour = "#C42126") # Draw mean line. -->
<!-- %' -->
<!-- %' @ -->
<!-- %' \caption[alpha]{Calculating power for a two--sided test.} -->
<!-- %' \label{PowerTwoSided} -->
<!-- %' \end{figure} -->



<!-- % From http://www.psychology.emory.edu/clinical/bliwise/Tutorials/SPOWER/SPOWAPP.HTM -->
<!-- % 1) when we are planning a study and -->
<!-- % 2) when we fail to reject the null hypothesis. -->
<!-- % -->
<!-- % We always want to think about statistical power when we are planning a research study. -->
<!-- % It makes no sense to develop a strong research design, collect and analyze data, and not have enough statistical power to be able to reject the null hypothesis. So, we should always think about how many participants we need in order to have enough statistical power to test the hypothesis of interest. What is "enough" statistical power? Many researchers think that a probability of .80 of rejecting the null hypothesis is good statistical power. How do we figure out how many participants we need? We estimate a likely population effect size from the current literature. We then calculate the sample size needed to detect an effect of that size with the desired level of power. We use the same procedures outlined in our power calculations but estimate "N" instead of estimating "1-b". -->

<!-- % We also need to think about statistical power when we fail to reject the null. Remember that when we reject the null, we can have two possible outcomes. We reject the null because there is no effect in the population; this would be a "correct" result. We also could be making a Type II error. In this case, we would not have enough sufficient statistical power for the test. How do we know which one is operating? We can never know for certain, but if we calculate our effect size, this can give us a clue. If our effect size is very small, then there is a good chance that there is no effect in the population. If our effect size is moderate to large, then we most likely do not have enough statistical power in our study. -->

## Some words of caution

Within this book, this chapter contains the toughest material -- well done for seeing it through! Things will fall into place from now on and will become easier to understand. Before we move on though, we'd like to draw your attention to three more things!

### What makes a Good Hypothesis?

Usually, you don't just collect data for fun, but with a certain question in mind. Good hypotheses, in the context of hypothesis testing, are testable questions that are derived from your research question. Sometimes the research question needs to be broken down so that it can be written (with a formula) in terms of a null hypothesis and a suitable alternative. That's a shame but if you manage, your test results will be meaningful!

### Statistical versus practical significance {#sigandsize}

We saw that the magnitude of a $p$--value depends on several factors, and not all of them tell you anything about the _practical_ significance of your results.^[If the sample size is large, then it is possible to find very tiny deviations from the null hypothesis in your data that may not be of any practical importance.]

<!-- % In fact, the Martian height case is a good example for this phenomenon. Are we really interested in differences between the proposed and measured Martian height in the order of a couple centimeters as in our example? Does it have any practical importance for us? If not, then we shouldn't get excited about low $p$--values. Actually, the average height of our Martians doesn't seem all that different from the the proposed height, even if the $p$--value is small! -->

In a way, statisticians have hijacked the ordinary meaning of the word _significance_, which implies important and practical relevance. Statisticians have given the word _significance_ its own meaning. Some not very statistically inclined people may not realise this and may get excited about a _statistically significant_ finding, when actually:

>The technical term _statistically significant_ relates to the mere existence of an effect, not about the size. To give yourself and the readers of your publications the chance to also assess the practical significance of your results ALWAYS PROVIDE A CONFIDENCE INTERVAL.

### Interpreting non-significant results

What's the right way of interpreting a non-significant $p$--value^[Actually, the term 'a non-significant $p$--value' doesn't make much sense if we are strict! What people mean when they say that is a $p$--value that exceeds the chosen significance level $\alpha$, which in most cases corresponds to $p > 0.05$. As you know now, however, $\alpha$ can in principle vary and should always be included in the sentence when you talk about statistical significance.]? In some papers you can read something like 'we accept the null hypothesis'. If this statement is based on a statistical hypothesis test, it is actually flawed! As we mentioned before, a non-significant $p$--value merely shows that you cannot rule out the null hypothesis. Based on the duality between confidence intervals and hypothesis testing that we covered in section \@ref(sec:ci2test), however, we know that there are many values that we cannot rule out, i.e. _all_ values inside your confidence interval! So we cannot accept just one of them!

Knowing this will also help you understand why the null hypothesis is usually chosen as the state of boredom that you want to reject. It is because you cannot _prove_ anything that you want to see, all you can do is to reject the stuff that is of no interest for you!

There are two lessons to be learned from this:
- Never state that the null hypothesis is true, or that you accept the null hypothesis.
- Never quote a $p$--value without a confidence interval.