## Background

The t-test is a [statistical hypothesis test](https://en.wikipedia.org/wiki/Statistical_hypothesis_testing) that is 
commonly used when the data are normally distributed (follow a normal 
distribution) if the value of the population standard deviation were known. When
the population standard deviation is not known and is replaced by an estimate
based no the data, the test statistic follows a Student's t distribution. 

T-tests are 
handy hypothesis tests in statistics when you want to compare means. You can 
compare a sample mean to a hypothesized or target value using a one-sample 
t-test. You can compare the means of two groups with a two-sample t-test. If 
you have two groups with paired observations (e.g., before and after 
measurements), use the paired t-test.

A t-test looks at the t-statistic, the t-distribution values, and 
the degrees of freedom to determine the statistical significance. 
To conduct a test with three or more means, we would use an analysis 
of variance.

The distriubution that the t-statistic follows was described in a famous 
paper by "Student", a pseudonym for [William Sealy Gosset](https://en.wikipedia.org/wiki/William_Sealy_Gosset), published in 
1908 in Biometrica^[http://seismo.berkeley.edu/~kirchner/eps_120/Odds_n_ends/Students_original_paper.pdf]. 

## The Z-score and probability

Before talking about the t-distribution and t-scores, lets review the Z-score, 
its relation to the normal distribution, and probability.

The Z-score is defined as:

$$Z = \frac{x - \mu}{\sigma}$$
where $\mu$ is a the population mean from which $x$ is drawn and $\sigma$ is
the population standard deviation (taken as known, not estimated from the data).

The probability of observing a $Z$ score of $z$ or greater can be calculated by
$pnorm(z,\mu,\sigma)$.

For example, let's assume that our "population" is known and it truly
has a mean 0 and standard deviation 1. If we have observations drawn from
that population, we can assign a probability of seeing that observation
by random chance _under the assumption that the null hypothesis is **TRUE**_. 

```{r}
zscore = seq(-5,5,1)
```

For each value of zscore, let's calculate the p-value and put the results in a `data.frame`.

```{r}
df = data.frame(
    zscore = zscore,
    pval   = pnorm(zscore, 0, 1)
)
df
```

Why is the p-value of something 5 population standard deviations away from
the mean (zscore=5) nearly 1 in this calculation? What is the default for
`pnorm` with respect to being one-sided or two-sided?

Let's plot the values of probability vs z-score:

```{r}
plot(df$zscore, df$pval, type='b')
```

This plot is the *empirical* cumulative density function (cdf) for our data. How can we use it? If we know the z-score, we can look up the probability of 
observing that value. Since we have constructed our experiment to follow
the standard normal distribution, this cdf also represents the cdf of 
the standard normal distribution.

### Small diversion: two-sided pnorm function

The `pnorm` function returns the "one-sided" probability of having
a value at least as extreme as the observed $x$ and uses the "lower" tail
by default. Let's create a function that computes two-sided p-values.

1. Take the absolute value of x
2. Compute `pnorm` with `lower.tail=FALSE` so we get lower p-values with
larger values of $x$.
3. Since we want to include both tails, we need to multiply the area 
(probability) returned by pnorm by 2.

```{r}
twosidedpnorm = function(x,mu=0,sd=1) {
    2*pnorm(x,mu,sd,lower.tail=FALSE)
}
```

And we can test this to see how likely it is to be 2 or 3 standard deviations
from the mean:

```{r}
twosidedpnorm(2)
twosidedpnorm(3)
```

## The t-distribution

We spent time above working with z-scores and probability. An important
aspect of working with the normal distribution is that we MUST assume
that we know the standard deviation. Remember that the Z-score is defined
as:

$$Z = \frac{x - \mu}{\sigma}$$
The formula for the *population*
standard deviation is:

$$\sigma = \sqrt{\frac{1}{N}\sum_{i=1}^{N}({xi - \mu)^2}}$$
In general, the population standard deviation is taken as "known" as
we did above. 

If we do not but only have a *sample*
from the population, instead of using the Z-score, we use the t-score 
defined as:

$$t = \frac{x - \bar{x}}{s}$$
This looks quite similar to the formula for Z-score, but here we have to
*estimate* the standard deviation, $s$ from the data. The formula for $s$
is:

$$s = \sqrt{\frac{1}{N-1}\sum_{i=1}^{N}({x_{i} - \bar{x})^2}}$$
Since we are estimating the standard deviation from the data, this leads
to extra variability that shows up as "fatter tails" for smaller sample
sizes than for larger sample sizes. We can see this by comparing the
_t-distribution_ for various numbers of degrees of freedom (sample sizes).

We can look at the effect of sample size on the distributions graphically
by looking at the densities for 3, 5, 10, 20 degrees of freedom and
the normal distribution:

```{r}
library(dplyr)
library(ggplot2)
t_values = seq(-6,6,0.01)
df = data.frame(
    value = t_values,
    t_3   = dt(t_values,3),
    t_6   = dt(t_values,6),
    t_10  = dt(t_values,10),
    t_20  = dt(t_values,20),
    Normal= dnorm(t_values)
) |>
    tidyr::gather("Distribution", "density", -value)
ggplot(df, aes(x=value, y=density, color=Distribution)) + 
    geom_line()
```

The `dt` and `dnorm` functions give the density of the distributions
for each point. 


```{r}
df2 = df |> 
    group_by(Distribution) |>
    arrange(value) |> 
    mutate(cdf=cumsum(density))
ggplot(df2, aes(x=value, y=cdf, color=Distribution)) + 
    geom_line()
```

### p-values based on Z vs t

When we have a "sample" of data and want to compute the statistical 
significance of the difference of the mean from the population mean, 
we calculate the standard deviation of the sample means (standard error).

$$z = \frac{x - \mu}{\sigma/\sqrt{n}}$$

Let's look at the relationship between the p-values of Z (from the normal distribution) vs t for a **sample** of data.

```{r}
set.seed(5432)
samp = rnorm(5)
z = sqrt(length(samp)) * mean(samp) #simplifying assumption (sigma=1, mu=0)
```

And the p-value if we assume we know the standard deviation:

```{r}
pnorm(z)
```

```{r}
ts = sqrt(length(samp)) * mean(samp) / sd(samp)
pnorm(ts)
pt(ts,5)
```

### Experiment

1. Simulate a bunch of samples of size `n` from the normal distribution
2. Calculate the p-value distribution for those samples based on the 
normal.
3. Calculate the p-value distribution for those samples based on the 
normal, but with the *estimated* standard deviation.
4. Calculate the p-value distribution for those samples based on the 
t-distribution.


Create a function that draws a sample of size `n` from the standard
normal distribution.

```{r}
zf = function(n) {
    samp = rnorm(n)
    z = sqrt(length(samp)) * mean(samp)
    z
}
```

And give it a try:

```{r}
zf(5)
```

Perform 10000 replicates of our sampling and z-scoring.

```{r}
z10k = replicate(10000,zf(5))
hist(pnorm(z10k))
```

And do the same, but now creating a t-score function. 

```{r}
tf = function(n) {
    samp = rnorm(n)
    ts = sqrt(length(samp)) * mean(samp) / sd(samp)
    ts
}
```

```{r}
t10k = replicate(10000,tf(5))
hist(pnorm(t10k))
```

```{r}
hist(pt(t10k,5))
```

And the qqplot:

```{r}
qqplot(z10k,t10k)
abline(0,1)
```

## t.test

### One-sample

```{r}
x = rnorm(20,1)
# small sample
t.test(x[1:5])
```

Increase sample size:

```{r}
t.test(x[1:20])
```

### two-sample

```{r}
x = rnorm(10,0.5)
y = rnorm(10,-0.5)
t.test(x,y)
```

### from a data.frame

```{r}
df = data.frame(value=c(x,y),group=as.factor(rep(c('g1','g2'),each=10)))
t.test(value ~ group, data=df)
```

### Equivalence to linear model

```{r}
t.test(value ~ group, data=df, var.equal=TRUE)
```

This is *equivalent* to:

```{r}
res = lm(value ~ group, data=df)
summary(res)
```