# Two--sample t--test -- comparing a continuous variable across two levels of a categorical variable {#two_sample_t_tests}

```{r echo = FALSE}
martian <- read.delim("data/martian.txt")
H.Means <- tapply(martian$Height, martian$Site, mean)
H.SDs <- plyr::round_any(tapply(martian$Height, martian$Site, sd),1)

Height.mean <- H.Means[1]

Height.sd <- H.SDs[1]
nValue <- length(martian$Height[martian$Site == "Site I"])
Height.SEM <- (Height.sd/sqrt(nValue))

Height.95CI <- c(Height.mean - 1.96*Height.SEM, Height.mean + 1.96*Height.SEM)

distance <- 5 # Set distance away from x-bar for tests

HypoValue <- Height.mean - distance # Set less than
HypoValue2 <- Height.mean + distance # Set greater than.

mu1 <- H.Means[1]
mu2 <- H.Means[2]
sigma1 <- H.SDs[1]
sigma2 <- H.SDs[2]

tValue <- (Height.mean - HypoValue)/Height.SEM
tValue2 <- (H.Means[1] - H.Means[2])/Height.SEM

nValue <- 10
nValue1 <- 10
nValue2 <- 10
denominator <- sqrt(sigma1^2/nValue + sigma2^2/nValue)

NewValues <- martian$Height_UnequalVar[martian$Site == "Site II"]
# c(207.2, 195.7, 198.6, 211.0, 204.7, 227.1, 205.1, 214.7, 212.3, 223.6) # sd = 10.00
# NewValues8 <- c(199.8, 221.6, 208.6, 221.9, 210.3, 200.1, 204.9, 210.1, 217.1, 205.6) # sd = 8.00
muNewValues <- mean(NewValues)
sNewValues <- round(sd(NewValues),2)
nNewValues <- length(NewValues)
# mean(martian$Height[martian$Site == "Site I"])
# t.test(martian$Height[martian$Site == "Site I"], NewValues)
# t.test(martian$Height[martian$Site == "Site I"], NewValues, var.equal = TRUE)


denominator2 <- sqrt(sigma1^2/nValue + sNewValues^2/nNewValues)


```

The first scenario we'll consider is when we want to quantify the effect that the levels of a categorical variable have on the values of a continuous variable. For now, the categorical variable is restricted to have only two levels (for more than two levels you'll have to wait until ANOVA). In particular, we want to see if a change in categories has an effect on the mean of the continuous variable, which is the sort of question you can tackle with a _two--sample t--test_. 

For this we return to our Martian dataset. We'd like to compare height at the two sites, cf. also Table \@ref(tab:heightVSsite).

Table: (\#tab:heightVSsite) The form of the dataset that we consider when we want to compare the Martian height across sites.

| Observation | Site   | Height (cm) |
|-------------|--------|-------------|
| 1           | 1      | 192.5       |
| ...      | ... | ...      |
| 10          | 1      | 201.8       |
| 11          | 2      | 205.5       |
| ...      | ... | ...      |
| 20          | 2      | 206.8       |

The sample size at both sites is the same, i.e. $n_1 = n_2 = 10$. Moreover, $\bar{x}=`r mu1`$ and $s=`r sigma1`$ at site I, and $\bar{x}=`r mu2`$ and $s=`r sigma2`$ at site II. This is plotted in Figure \@ref(fig:AltVis). 

```{r AltVis, fig.cap="Dot plots depicting every data point. The mean height and standard deviation for each site is drawn in black.", fig.width = large_width, echo = FALSE}

ggplot(martian, aes (x = Site, y = Height, col = Site)) +
  geom_jitter(position = position_jitter(width = .1), alpha = 0.6) +
  stat_summary(fun.y = "mean", geom = "point", col = "black") +
  stat_summary(fun.data = "mean_sdl", geom = "errorbar", fun.args = list(mult = 1), width = 0.05, col = "black") +
  scale_x_discrete("") +
  scale_y_continuous("",
                     limits = c(185,225),
                     breaks = seq(185, 225, 5),
                     labels = c("", "190", "", "200", "", "210", "", "220", ""),
                     expand = c(0,0)) +
  scale_colour_manual(values = c(MidBlue, MidRed), guide = FALSE) +
  theme(panel.background = element_blank(),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.line.x = element_line(colour = "white"),
        axis.line = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        panel.spacing=unit(1.5, "lines")
  )

```

Plotting the data is actually the very first thing you need to do when you want to start any sort of analysis! When you afterwards plan to carry out a $t$--test, you also need to check the underlying assumption of Normally-distributed data!
<!-- ^[Recall that we can use a Q--Q plot (Figure \@ref(fig:NewQQPlots)) to check whether our data are reasonably Normally-distributed. If your data are not Normal, you can use alternatives to the $t$--test, e.g. the 'Mann--Whitney U test' (also know as the 'Wilcoxon rank--sum test').] -->

Here, we assume that the Martian height at both sites actually is Normally distributed, i.e. $$ X \sim N(\mu_1, \sigma_1) \quad \quad Y \sim N(\mu_2, \sigma_2),$$ where $X$ and $Y$ denote the random variables 'Height at site I' and 'Height at site II', respectively.

We are interested in whether the mean Martian height is different between the sites.

```{exercise, name="Defining statistical hypotheses"}
Write down the null hypothesis $H_0$ and the alternative hypothesis $H_1$ we wish to test given we want to perform a two--sided test! Hint: use the two (unknown) site mean values $\mu_1$ and $\mu_2$.
```



Now that we know what the null and alternative hypotheses are, let's think about the signal--to--noise ratio we are interested in here. Remember in the previous chapter we had the (one--sample) situation where our hypotheses had the form
$$ H_0: \mu = `r HypoValue` \quad \quad \mbox{versus}  \quad \quad H_1: \mu \neq `r HypoValue`.$$
In that situation, we looked at

$$\frac{\mbox{signal}}{\mbox{noise}} = \frac{\mbox{(observed value - hypothesised value)} }{\mbox{scale factor}} = \frac{\bar{X} - \mu_0}{S/\sqrt{n}}$$

which was our $t$--statistic that, under the premise of a true $H_0$, followed a $t$--distribution with $n-1$ degrees of freedom. To test whether the signal was big enough to reject $H_0$ in our sample, we just needed to plug in our sample estimates $\bar{x}$ and $s$, and see where the resulting value $t = (\bar{x} - \mu_0)/(s/\sqrt{n})$ sits on the $t_{n-1}$--distribution.

Here, the situation is slightly more complicated because we have two Normal distributions and four parameters that we somehow must account for. One trick is to rewrite our hypotheses in a way that looks more similar to the one in the one--sample situation:
$$ H_0: \mu_1 - \mu_2 = 0 \quad \quad \mbox{versus}  \quad \quad H_1: \mu_1 - \mu_2 \neq 0.$$

This formulation emphasises that it is the _difference_ between the means $\mu_1 - \mu_2$ we are interested in here. With this, it is easy to see that the ingredients for our test statistic are

**observed value**: We cannot observe the difference between the two population means directly, but like in the one--sample situation we can estimate it using the sample means. Hence we look at $\bar{x}_1 - \bar{x}_2$ which is an estimate for $\mu_1 - \mu_2$.

**hypothesised value**: The hypothesised value is zero (which again corresponds to the state of boredom).

**scale factor**: In the one--sample situation the scale factor was the $SEM$, i.e. the standard deviation of the sample mean $sd(\bar{X})$ which measures the precision that is associated with estimating the population mean by the sample mean. Here, we need the standard deviation of $\bar{X}_1 - \bar{X}_2$ which is the estimate we are interested in now. Let's call this $sd(\bar{X}_1 - \bar{X}_2)$ for now.

From the above, we can now derive the form of the test statistic we are interested in here:

$$\frac{\mbox{signal}}{\mbox{noise}}  =  \frac{\text{(observed value - hypothesised value)} }{\text{scale factor}}  =  \frac{(\bar{X}_1 - \bar{X}_2) - 0}{sd(\bar{X}_1 - \bar{X}_2)}  = \frac{\bar{X}_1 - \bar{X}_2}{sd(\bar{X}_1 - \bar{X}_2)}.$$

Easy! However, we are not quite out of the woods yet, we still don't know the form of $sd(\bar{X}_1 - \bar{X}_2)$ and we also don't know the null distribution of the test statistic.^[Given we're still talking about $t$--tests, it's not hard to guess that the null distribution that we're after will be a $t$--distribution. The real question therefore is 'what's the correct degrees of freedom'?]

Actually, there are different versions of $sd(\bar{X}_1 - \bar{X}_2)$ depending on the nature of your data, and these versions also determine which $t$--distribution you should use. We'll go through the three most important versions in this workshop.

```{R echo = FALSE}
tValue2 <- (mu1 - mu2) /
sqrt((sigma1^2/nValue) + (sigma2^2/nValue))
```

## Version 1 -- unequal variances test

```{r echo = FALSE}

mars.t.1 <- t.test(martian$Height ~ martian$Site)
mars.t.2 <- t.test(martian$Height_UnequalVar8 ~ martian$Site)
mars.t.3 <- t.test(martian$Height_UnequalVar ~ martian$Site)
mars.t <- mars.t.1
```

The most general form of the two--sample $t$--test is called _Welch's variant for unequal variances_.^[Good old Mr. Welch wrote the corresponding paper in 1948, not that long ago when you think about it. Stats really is a very young discipline! Welch, B. L. (1947). "The generalization of Student's problem when several different population variances are involved". Biometrika 34.] The 'unequal variances' refers to the population parameters $\sigma_1^2$ and $\sigma_2^2$. You don't know them of course, but if the corresponding estimates $s_1^2$ and $s_2^2$ obtained from the sample are quite different, it's unlikely that $\sigma_1^2 = \sigma_2^2$. In this case, $sd(\bar{X}_1 - \bar{X}_2)$ has the form

$$sd(\bar{X_1} - \bar{X_2}) = \sqrt{{\frac{S_1^2}{n_1}} + {\frac{S_2^2}{n_2}}}.$$

We won't get into the reasons for this, but the form of this formula shouldn't seem too absurd when you remember that in the one--sample situation the scale factor had the form $SEM = \sqrt{S^2/n}$. With the the above formula plugged in, the entire $t$--statistic looks like

$$ t = \frac{\bar{X}_1 - \bar{X}_2}{\sqrt{{\frac{S_1^2}{n_1}} + {\frac{S_2^2}{n_2}}}} $$

If we plug the values of our Martian height example in the equation, we find the $t$--statistic

$$t = \frac{(`r mu1`-`r mu2`)}
{\sqrt{(\frac{`r sigma1`^2}{`r nValue1`})+(\frac{`r sigma2`^2}{`r nValue2`})}}
= \frac{`r mu1-mu2`}{`r denominator`}
= `r (mu1-mu2)/denominator`.$$

Under the null hypothesis $H_0: \mu_1 - \mu_2 = 0$, it follows a $t$--distribution. In the case of assuming unequal variance, the degrees of freedom have quite a complicated form

$$\mathrm{d.f.} = \frac{(s_1^2/n_1 + s_2^2/n_2)^2}{(s_1^2/n_1)^2/(n_1-1) + (s_2^2/n_2)^2/(n_2-1)}.$$

For our example, there are `r as.character(round(mars.t.1$parameter,2))` degrees of freedom. Remember, the $p$--value is equal to the area under the null distribution for all values that are _further_ from zero, in either direction than the $t$--value we observed. This is illustrated in Figure \@ref(twoSamDist). The $t$--statistic is really far from 0, so the corresponding $p$--value is very low.


```{r twoSamDist, fig.width = large_width, echo = FALSE}

xRange <- 4 # Set the upper limit of the x-axis
#samples <- 10
df <- mars.t$parameter
x <- seq(-1*(xRange*2), (xRange*2), length=5000) # Define an x series based on inputed mu and xRange
tdist1 <- dt(x, df) # Define the density of the standard normal curve at each x position


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

# Base plot:
base33 <-
  ggplot(test.df, aes ( x = x, y = tdist1)) +
  geom_line(col = MidYellow) +
  scale_x_continuous("x", breaks = seq(-xRange, xRange, 2), labels = seq(-xRange, xRange, 2), expand = c(0,0)) +
  geom_ribbon(data=subset(test.df, x < -abs(mars.t$statistic)),
              aes(x=x,ymax=tdist1),ymin=0,fill=MidRed) +
  geom_ribbon(data=subset(test.df, x > abs(mars.t$statistic)),
              aes(x=x,ymax=tdist1),ymin=0,fill=MidRed) +
  geom_vline(xintercept = mars.t$statistic, colour = "black") +
  geom_vline(xintercept = 0, colour = MidRed) +
  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()
  )

base33 + coord_cartesian(xlim = c(-xRange, xRange), ylim = c(0, 0.4)) + scale_y_continuous("f(x)", breaks = seq(0, 0.5, 0.1), expand = c(0,0))

```
The distribution of the $t$--statistic under the null hypothesis $H_0: \mu_1 - \mu_2 = 0$ is the $t_{`r mars.t$parameter`}$--distribution. The position of our observed value for this statistic is shown by the red dashed line. To calculate the $p$--value, the area under the curve is summed for both tails of the distribution because we decided to use a two--sided test.

In principle this is it -- Welch's variant covered!!  However, recall that $s_1 = s_2 = `r sigma1`$ cm for our example, so in reality it actually seems not unlikely that $\sigma_1 = \sigma_2$ in which case we can use a version of the the two--sample $t$--test that is build on this assumption. So to show the impact of different standard devations in the two groups, consider two other scenarios, where all the parameters are the same as before except now, $s_2 = 8$ cm (scenario 2) and $s_2 = `r sNewValues`$ cm (scenario 2), cf. Table \@ref(tab:WelchScenarios).

Table: (\#tab:WelchScenarios) Welch's variant of the two--sample $t$-test for comparing groups having unequal variances. Three scenarios are considered: Two groups with a consistent difference in $\bar{x}$ but with progressively larger differences in $s$.

| Parameter              | Scenario 1                                    | Scenario 2                                    | Scenario 3                                    |
|------------------------|-----------------------------------------------|-----------------------------------------------|-----------------------------------------------|
| Group 1 ($\bar{x}, s$) | (200, 6)                                      | (200, 6)                                      | (200, 6)                                      | 
| Group 2 ($\bar{x}, s$) | (210, 6)                                      | (210, 8)                                      | (210, 10)                                     | 
| t-statistic            | `r mars.t.1$statistic`                        | `r mars.t.2$statistic`                        | `r mars.t.3$statistic`                        | 
| Degrees of Freedom     | `r as.character(round(mars.t.1$parameter,2))` | `r as.character(round(mars.t.2$parameter,2))` | `r as.character(round(mars.t.3$parameter,2))` | 
| $p$-value              | `r as.character(round(mars.t.1$p.value,3))`   | `r as.character(round(mars.t.2$p.value,3))`   |`r as.character(round(mars.t.3$p.value,3))`    | 

These three scenarios are plotted in Figure \@ref(fig:TtestScenarios), and the results of Welch's two-sample $t$--tests for each scenario are summarised in Table \@ref(tab:WelchScenarios). 
 
```{r TtestScenarios, fig.cap = "Mean with standard deviation.", fig.width = large_width, echo = FALSE}

df.Height <- martian[2:5]
df.Height <- melt(df.Height)
levels(df.Height$variable) <- c("Scenario 1\n(Original Data)", "Scenario 2\n", "Scenario 3\n")

ggplot(df.Height, aes (x = Site, y = value, col = Site)) +
  geom_jitter(position = position_jitter(width = .1), alpha = 0.6, size = 3) +
  stat_summary(fun.y = "mean", geom = "point", col = "black") +
  stat_summary(fun.data = "mean_sdl", geom = "errorbar", fun.args = list(mult = 1), width = 0.05, col = "black") +
  scale_x_discrete("") +
  scale_y_continuous("",
                     limits = c(185,230),
                     breaks = seq(185, 230, 5),
                     labels = c("", "190", "", "200", "", "210", "", "220", "", "230"),
                     expand = c(0,0)) +
  scale_colour_manual(values = c(MidBlue, MidRed), guide = FALSE) +
  facet_grid(. ~ variable) +
  theme(panel.background = element_blank(),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.line.x = element_line(colour = "white"),
        axis.line = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        panel.spacing=unit(1.5, "lines")
  )

```

To consider the effects of increasing $\sigma_2$, so that $\sigma_2 \neq \sigma_1$, we compare two simulations with our original scenario, where $\sigma_2 = \sigma_1$.

What can we deduce from this table? As $s_2$ increases, keeping all other parameters constant, two things happen. _Both_ the degrees of freedom _and_ the test statistic _decrease_. Both of these changes have the consequence of increasing the $p$-value. This result makes sense, since the large variance in one group precludes a clear distinction between the two groups. A simple way to understand this is that the error bars have a greater overlap, as seen in Figure \@ref(fig:TtestScenarios).

Notice that we can still get a very small $p$--value with overlapping error bars! It is a common misconception that overlapping errors correspond to high, i.e. _non-significant_ $p$--values.


```{exercise, name="Defining the p-value"}
How to you interpret the $p$--values above? What's your conclusion regarding the null hypothesis?
```

## Version 2 -- equal variances

```{r echo = FALSE}

mars.t.1 <- t.test(martian$Height ~ martian$Site, var.equal = TRUE)
mars.t.2 <- t.test(martian$Height_UnequalVar8 ~ martian$Site, var.equal = TRUE)
mars.t.3 <- t.test(martian$Height_UnequalVar ~ martian$Site, var.equal = TRUE)

```


Mr. Welch's version of the two--sample $t$--test didn't have any restriction regarding the two population variances $\sigma_1^2$ and $\sigma_2^2$, and essentially this resulted in a scale factor $sd(\bar{X}_1 - \bar{X}_2)$, i.e. the noise in the denominator of our signal--to--noise ratio, where each of the two variances were estimated separately. However, if you can assume that the two groups have the some population variance, the situation is a bit different, because in that case your underlying assumption is

$$ X \sim N(\mu_1, \sigma^2) \quad \quad Y \sim N(\mu_2, \sigma^2).$$

So there is only _one_ variance that you need to estimate! And because it's the same for both groups, you can use a combined form to estimate the common variance as

$$S^2_{pooled} = \frac{(n_1-1)S_{1}^2+(n_2-1)S_{2}^2}{n_1+n_2-2}.$$

With this, the common scale factor has the form 

$$sd(\bar{X}_1 - \bar{X}_2) = S_{pooled} \cdot \sqrt{\frac{1}{n_1}+\frac{1}{n_2}}, $$ 

and the test statistic becomes

$$t = \frac{\bar{X}_1 - \bar{X}_2 }{\sqrt{\frac{(n_1-1)S_{1}^2+(n_2-1)S_{2}^2}{n_1+n_2-2}}\cdot \sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}.$$

Under the null hypothesis it follows a $t$--distribution with $n_1 + n_2 - 2$ degrees of freedom. In our example that would be $`r nValue`+`r nValue`-2=`r nValue+nValue-2`$. Table \@ref(tab:NotWelchScenarios) summarises the results for all three scenarios. We see that the main difference to Table \@ref(tab:WelchScenarios) is that the degrees of freedom remains constant. Also in each scenario, it is at least as large is the one in Table \@ref(tab:WelchScenarios) where it decreases as the difference in variance between the groups increases.


Table: (\#tab:NotWelchScenarios) The two--sample $t$-test for comparing groups assuming equal variances. Three scenarios are considered: Two groups with a consistent difference in $\bar{x}$ but with progressively larger differences in $s$..

| Parameter              | Scenario 1                                    | Scenario 2                                    | Scenario 3                                    |
|------------------------|-----------------------------------------------|-----------------------------------------------|-----------------------------------------------|
| Group 1 ($\bar{x}, s$) | (200, 6)                                      | (200, 6)                                      | (200, 6)                                      | 
| Group 2 ($\bar{x}, s$) | (210, 6)                                      | (210, 8)                                      | (210, 10)                                     | 
| t-statistic            | `r mars.t.1$statistic`                        | `r mars.t.2$statistic`                        | `r mars.t.3$statistic`                        | 
| Degrees of Freedom     | `r as.character(round(mars.t.1$parameter,2))` | `r as.character(round(mars.t.2$parameter,2))` | `r as.character(round(mars.t.3$parameter,2))` | 
| $p$-value              | `r as.character(round(mars.t.1$p.value,3))`   | `r as.character(round(mars.t.2$p.value,3))`   |`r as.character(round(mars.t.3$p.value,3))`    | 



   <!-- $t$-statistic & -->
   <!-- \Sexpr{mars.t.1$statistic} & -->
   <!-- \Sexpr{mars.t.2$statistic} & -->
   <!-- \Sexpr{mars.t.3$statistic} \\ -->

   <!-- Degrees of Freedom & -->
   <!-- \Sexpr{as.character(round(mars.t.1$parameter,2))} & -->
   <!-- \Sexpr{as.character(round(mars.t.2$parameter,2))} & -->
   <!-- \Sexpr{as.character(round(mars.t.3$parameter,2))} \\ -->

   <!-- $p$-value & -->
   <!-- \Sexpr{as.character(round(mars.t.1$p.value,3))} & -->
   <!-- \Sexpr{as.character(round(mars.t.2$p.value,3))} & -->
   <!-- \Sexpr{as.character(round(mars.t.3$p.value,3))} \\ -->

   <!-- \bottomrule -->

So what's the effect of assuming a constant variance in both groups? Recall that higher degrees of freedom means that the $t$--distribution more closely resembles the standard Normal distribution. This in turn results in lighter tails^[The expression 'light tail' in the context of a statistical distribution refers to the distributions that have high values around zero and fall off quickly when you move away from zero. For $t$--distributions, the tail gets lighter and lighter if you increase the degrees of freedom. Use the $t$--distribution app to convince yourself about this.]. Like in the previous chapter, these lighter tails have to do with the smaller amount of uncertainty that is inherent in the test procedure. Here, the explanation for this is that we need to estimate only one variance parameter when we assume a common variance for both groups. In general, this will result in smaller $p$--values for the same signal--to--noise ratio. 

<!-- **maybe** a pic with the same value of the test statistic sitting on different $t$--dists with different dfs, and the corresponding area under the curves, two--tailed -->

<!-- % When we assume equal variance, we achieve a lower $p$-value in each scenario. This can be understood in terms of the $t$-statistic and the degrees of freedom for each test. In both versions of the test, as the difference in $s$ increased, the $t$-statistic decreased, resutling in lower $p$-values. However, in Welch's version, the degrees of freedom \emph{also} decreased, meaning that $p$-values were further reduced. Here, where a difference in variance doesn't play a role in calculating the degrees of freedom, it is a constant in all three scenarios. -->

This simulation may lead you to conclude that you should always assume equal variance, since that will result in a lower $p$-value. However, a word of warning: if the variances _are equal_, your procedure will be more precise, BUT if they _are unequal_, you spuriously remove uncertainty, which will result in an inflated type I error rate!

An obvious question is how can you determine if you can assume equal variances. For example, the $F$--test is designed for just that.^[We saw this already when we discussed testing for Normality on page \@ref(sec:qqplots).] However, testing the assumptions with yet another test that also has assumptions is not really optimal. Where do you stop testing assumptions once you are in that loop of testing test assumptions by testing?! The recommendation we give is to eyeball your data to judge whether the assumption of equal variances in your two groups seems OK. If it's unclear, use the Welch variant of the two--sample t--test as the conservative option.

In this simulation, both the sample sizes and the sample means were held constant in the three scenarios. If these values differ, we begin to see further differences between the two methods. Use the two-sample t-test web app to explore the scenarios presented in this section. Also explore the effect of varying sample sizes.

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

A1 <- c(29.97, 29.99, 29.99, 30.01, 30.02, 30.11)
A2 <- c(29.72, 29.89, 29.93, 29.98, 29.98, 30.02)

mean(A2) - mean(A1)
var(A2) - var(A1)

t.test(A1, A2, var.equal = T)
t.test(A1, A2, var.equal = F)

df <- data.frame(A1, A2)
library(reshape2)
library(ggplot2)
df <- melt(df)
ggplot(df, aes(x=variable, y = value )) +
  geom_point(position = position_jitter(0.05), size = 5, alpha = 0.6) +
  stat_summary(fun.y = "mean", geom = "point") +
  stat_summary(fun.data = "mean_sdl", geom = "errorbar", width = 0.1, fun.args = list(mult = 1))
sd(A1)
sd(A2)
t.test(A1, A2, var.equal = T)
t.test(A1, A2, var.equal = F)

t.test(A1*10, A2*10, var.equal = T)
t.test(A1*10, A2*10, var.equal = F)
```

## Version 3 -- paired data

The last type of $t$--test to consider is when the data are paired. This means that we measured the same continuous variable twice in _the same_ individual, for example under two different conditions. In our Martian exaple, the continuous variable we will use is the time to run 100 m. We have already measured this value on Mars, so after returning to Earth, we decided to let our Martians run again. Now we'd like to compare the time on Mars to that on Earth. Therefore, our dataset has the form shown in Table \@ref(tab:pairedspeed).

Table: (\#tab:pairedspeed) The form of the Martian dataset when investigating time to run 100 m on Mars versus Earth, for the same individual.

| Martian | Mars (s) | Earth (s) |
|----|-----|-----|
| 1 | 23.25 | 18.98 |
| 2 | 27.44 | 19.19 |
| 3 | 22.12 | 13.39 |
| 4 | 16.34 | 10.93 |
| ... | ... | ... |
| 9 | 21.32 | 18.18 |
| 10 | 20.84 | 16.17 |

Here, we will restrict our analysis to the `r nValue` Martians from site I. As we measures each Martian twice, we are dealing with `r nValue*2` data-points. Importantly, however, they are _not_ independent. Instead, they come in pairs, and we can assume that these pairs are independent, since they come from different Martians. This data situation is is shown graphically in Figure \@ref(fig:MartianPairedlines). As each measurement on Earth has a corresponding measurement on Mars, it makes sense to connect values within each pair. In a situation like this, you have _paired_ data and should run a _paired $t$--test_ where we are primarily interested in the change in the measured value between the two states _of the same individual_.^[This is actually a much nicer situation because this is more direcly what we were interested in. We don't have to worry about potential confounding effects, e.g. if by chance the group on Mars contained the Martians that would have been fast on Earth as well!] 

What does that mean? Essentially, we are more interested in the _slope_ of the line than the actual values. Before, we were interested in the yellow and red dots in Figure \@ref(fig:MartianPairedlines), now we are interested in the black lines. Once we understand that we have paired data, we realise that we are actually dealing with a one--sample situation again, because we have a single sample of slopes! This is shown in Figure _PairedSeries_ where the difference between paired values is plotted.

<!-- % **Exercise**From table \ref{tab:pairedspeed}, what are the slopes for Martians 1 and 2? -->

Just like in the one--sample $t$--test, the underlying assumption is that the slopes are Normally-distributed, and the null hypothesis is that the mean of the slopes is zero! The test statistic has the same form as that for the one--sample $t$--test we saw in chapter \@ref(chap:onesamplepractice), and its null distribution again is $t_{n-1}$. The results of our paired $t$--test (or a one-sample $t$--test on the slopes, if you like) are summarised in the following table.

(ref:mPair) A scatterplot of the data `cars` using **base** R graphics. 

Paired data for Martian 100 m trials. _Left_: The black lines connect the speed on Mars and that on Earth that were measured for each of our 10 Martians from site I. _Right_: Differences between speed on Mars and speed on Earth for the 10 Martians. The values correspond to the slopes on the black lines in the previous figure.

```{r MartianPairedlines, fig.cap='(ref:mPair)', fig.width = large_width, echo = FALSE}


martian %>% 
  select(c(Site, Speed.Mars, Speed.Earth)) %>% 
  gather(key, value, -Site) %>% 
  filter(Site == "Site I") %>% 
  mutate(key = factor(key, labels = c("Mars", "Earth")),
         obs = rep(1:(length(value)/2), 2),
         jitt = jitter(c(rep(1, (length(value)/2)),rep(2, (length(value)/2))), .2)) -> martian.m

paired1 <- ggplot(martian.m, aes(x = jitt, y = value, group = obs, col = key)) +
  coord_fixed(0.5, xlim = c(0.5,2.5)) +
  geom_line(col = "black", alpha = 0.2) +
  geom_point(alpha = 0.6, size = 3) +
  scale_y_continuous(name = "Time to run 100m (s)",
                     limits = c(16, 28),
                     # labels = seq(9.5, 14.5, 1.5),
                     # breaks = seq(9.5, 14.5, 1.5),
                     expand = c(0,0)) +
  scale_x_continuous(name = "Location",
                     limits = c(0.8,2.2),
                     labels = c("Mars", "Earth"),
                     breaks = 1:2,
                     expand = c(0,0)) +
  scale_colour_manual("location", values = c(DarkYellow, DarkRed), labels = c("Mars", "Earth"), guide = FALSE) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.line.y = element_line(colour = "black"),
        axis.line.x = element_blank(),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black")
  )


# Paired:
martian.m %>% 
  select(-Site, -jitt) %>%
  spread(key, value) %>% 
  mutate(diff = Mars - Earth,
         group = 1) -> martian.m

paired2 <- ggplot(martian.m, aes(x = group, y = diff)) +
  coord_fixed(0.25, xlim = c(0.5,1.5)) +
  geom_point(col = DarkYellow, alpha = 0.6, size = 4, position = position_jitter(0.3, seed = 136)) +
  stat_summary(fun.y = "mean", geom = "point", col = "black") +
  stat_summary(fun.data = "mean_sdl", geom = "errorbar", fun.args = list(mult = 1), width = 0.05, col = "black") +
  scale_y_continuous("Difference in time (Earch - Mars)", limits = c(-8, 4), breaks = seq(-10, 4, 2), expand = c(0,0)) +
    theme(panel.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.ticks.y = element_line(colour = "black"),
          axis.line.y = element_line(colour = "black"),
        axis.text.y = element_text(colour = "black"),
          axis.line.x = element_line(color = "white"),
        axis.ticks.x = element_line(color = "white"),
        axis.title.x = element_text(color = "white"),
        axis.text.x = element_text(color = "white")
        )

cowplot::plot_grid(paired1, paired2, ncol = 2, align = "v", axis = "b")

```

```{r fig.width = large_width, echo = FALSE}

mars.t <- t.test(martian[martian$Site == "Site I",]$Ratio)
# mars.t <- t.test(martian.m$diff)
mars.t$alternative <- "Two-sided"

df.res <- data.frame(Test = mars.t$method,
                     "Test mean, group 1" = mars.t$estimate[1],
                     "Null difference in means" = mars.t$null.value,
                     "t statistic" = mars.t$statistic,
                     "Direction of test" = mars.t$alternative,
                     "Degrees of freedom" = mars.t$parameter,
                     "p--value" = mars.t$p.value,
                     "95 confidence interval" =
                       paste("(",round(mars.t$conf.int[1], 2),", ", round(mars.t$conf.int[2], 2),")", sep = "")
)


df.res <- data.frame(t(df.res), row.names = c("Test", "Sample mean", #"Test mean, group 2",
                                              "Hypothesised value", "t-statistic",
                                              "Direction of test", "Degrees of Freedom",
                                              "p-value", paste("95",I("%")," Confidence Interval", sep ="")))

# uncomment this:
# tli.table <- xtable(df.res)
# pairedT_p_value <- mars.t$p.value
```

```{r echo = FALSE}
df.res %>%
  knitr::kable(col.names = c(""))
# print (tli.table, include.colnames=FALSE, floating=FALSE,hline.after=c(0,nrow(tli.table)))
```

