---
title: "Intro to the $t$-distribution"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r top_pic, echo=F, message=F}
library(openintro)
data(COL)

#par(
      #mar = c(2, 1, 1, 1),
      #mgp = c(5, 0.6, 0))
plot(c(-3.5, 8),
     c(0, dnorm(0)),
     xlab="",ylab="",
     type = 'n',
     #sub="Open Intro, figure 5.2",
     axes = FALSE)
at <- seq(-4, 10, 2)
axis(1, at)
axis(1, at - 1, rep("", length(at)), tcl = -0.1)
abline(h = 0)

COL. <- fadeColor(COL[1], c("FF", "89", "68", "4C", "33"))
COLt <- fadeColor(COL[1], c("FF", "AA", "85", "60", "45"))
DF   <- c('normal', 8, 4, 2, 1)

X <- seq(-5, 10, 0.02)
Y <- dnorm(X)
lines(X, Y, col = COL.[1])

for (i in 2:5) {
  Y <- dt(X, as.numeric(DF[i]))
  lines(X, Y, col = COL.[i])
}

legend(4.5, 0.4,
       legend = c(DF[1],
       paste('t, df = ', DF[2:5], sep = '')),
       col = COL.,
       text.col = COLt,
       lty = rep(1, 5))
```

### A problem

The following data records the average mercury content in dolphin muscle of 19 Risso’s dolphins from the Taiji area in Japan:

```{r}
d = c(2.57,4.43,2.09,7.68,4.77,2.12,5.13,5.71,5.33,3.31,7.49,4.91,2.58,1.08,6.60,3.91,3.97,6.18,5.90)
```

While it might be important for researchers to assess the threat of mercury 
in the ocean, they do *not* want to go kill more dolphins to get that data.
What can they conclude from this data set, even though it's a bit too small
to use a normal distribution?

## The basics

The normal distribution, as awesome as it is, requires that we work with large sample sizes.

The $t$-distribution is similar but better suited to small sample sizes.

Just as with the normal distribution, there's not just one $t$-distribution but, rather, a *family* of distributions.

Just as there's a formula for the normal distribution, there's a formula for the $t$-distribution. It's a bit more complicated, though.

Like all continuous distributions, we compute probabilities with the $t$-distribution by computing the area under a curve. We do so using either a computer or a table.

The degree to which a $t$-distribution deviates from the normal is determined by
its "degrees of freedom" parameter $df$, which is one less than the sample size.

The mean of the $t$-distribution is zero and its variance is related to the degrees of freedom $df$ by
$$\sigma^2 = \frac{df}{df-2}.$$

Unlike the normal distribution, there's no easy way to translate from a
$t$-distribution with one standard deviation to another standard one. As a result,
it's less common to use tables and more common to use software than it is with the
normal.

Given a particular number of degrees of freedom, however, there is a standard way 
to derive a $t$-score that's analogous to the $z$-score for the normal distribution.
This $t$-score is a crucial thing that you need to know when using tables for
the $t$-distribution.


## The $t$-test for one sample mean

Our main uses of the $t$-distribution will be to construct confidence intervals
and perform hypothesis tests for relatively small data sets. The main tool for 
doing so is the so called $t$-test. The $t$-test is to the $t$-distribution as 
the $p$-test is to the normal distribution. 

The $t$-test a little tricky due to the lack of a "standard" version of the
$t$-distribution. R has a built command, called `t.test` that automates the
procedure, though.
While `t.test` uses the $t$ distribution to deal with small sample sizes, it 
actually works well with samples of all sizes. The reason is that, as the degrees
of freedom grows, the corresponding $t$-distribution approaches the normal 
distribution.

### Example

Here's how we can examine that data set on mercury levels in Risso’s dolphins
using a $t$-test in R:

```{r}
d = c(2.57,4.43,2.09,7.68,4.77,2.12,5.13,5.71,5.33,3.31,7.49,4.91,2.58,1.08,6.60,3.91,3.97,6.18,5.90)
t.test(d, mu=3.6)
```


There's quite a bit going on here that we need to parse. First, there's the command itself: 

    t.test(d, mu=3.6)

I guess it's pretty clear that we're applying `t.test` to our data set `d`. The `mu=3.6` option is used to specify an assumed null-value that we are comparing the data against. Thus, in the output we see the lines:

    ## t = 2.1177, df = 18, p-value = 0.04838
    ## alternative hypothesis: true mean is not equal to 3.6

Note that the value of $3.6$ appears in the statement of the alternative hypothesis; the associated $p$-value for the alternative hypothesis is, evidently, $0.04838$. Thus, we can reject the null hypothesis that the true mean is $3.6$.

The computation of the $p$-value is based on the $t$-score in much the same way that the $p$-value for the normal distribution is based on the $z$-score. In the output above, the $t$-score is indicated the `t = 2.1177` in the output.

Finally, the $95\%$ confidence interval is just what you'd probably think - 
we have a $95\%$ level of confidence that the interval contains the actual mean. 
This confidence interval is constructed in exactly the same way that we've
constructed confidence intervals before - it has the form
$$[\bar{x} - ME, \bar{x} + ME].$$

For this particular problem, we might be interested in a one-sided hypothesis test to check, for example, if the mercury concentration is too high. Thus, we might
run the command

```{r}
t.test(d, mu=3.6, alternative = "greater")
```

Not surprisingly, the $p$-value is exactly half of the previous $p$-value and the
confidence interval is a semi-infinite interval that lies entirely to the right 
of the mean.

### Using the distribution itself

Direct use of the `t.test` command is the easiest way to apply a $t$-test and
often the way to use the $t$-distribution. However, we should be able to gain a
deeper understanding of the $t$-distribution and it's relationship to the normal
distribution by using the fundamental distribution functions:

 * `dt` - the density function
 * `pt` - the cumulative probability function
 * `qt` - the inverse of `pt` known as the quantile function

These are all analogous to the functions the corresponding functions `dnorm`,
`pnorm`, and `qnorm` that we've used when exploring the normal distribution.

In fact, every distribution in R has functions that following this naming scheme.
 
### `dt`

This function defines the density for the $t$-distribution. Thus, we can plot it
like so:

```{r}
f = function(x){dt(x,5)}
plot(f, -3,3)
```

### `pt`

The *cumulative probability distribution function* for the $t$-distribution 

Given a random variable $X$ with a $t$-distribution with `df` degrees of freedom and a real number $x$, `pt(x,df)` computes the area under the graph from $-\infty$ to $x$. In the picture below, $df=18$ and $x=-1$.

For example, the following command computes the probability that a random variable $X$ generated by the $t$-distribution with $18$ degrees of freedom is less than $-1$:

```{r}
pt(-1,18)
```

This is illustrated by the following picture:

```{r cummulative_pic, echo=F}
dx = 0.01
x <- seq(-6,6, by=dx)
y <- dt(x, 18)
plot(x,y, type='l', ylim=c(0,0.6), yaxs='i', xlim=c(-5,5), xlab="", ylab="", yaxt='n', lwd=3)
x <- seq(-6,-1, by=dx)
y <- dt(x, 18)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='lightgray')
title(main = 't-density function with df=18 \n and the area representing P(X < -1)')
text(-3,0.22, 'p=0.165')
arrows(-3,0.2, -1.5,0.05)
```


### `qt` 
 
The *quantiles* for the $t$-distribution
 
This is an inverse to `pt`. Thus, given a probability $p$ and some number $df$ of degrees of freedom, `qt(p,df)` computes the real number $x$ such that `pt(x,df)` produces that probability $p$. Symbolically,
 $$\texttt{pt}(\texttt{qt}(p,df),df) = p.$$
Referring again to the picture above, we saw that $\texttt{pt}(-1,18)=0.165$; thus, we should have $\texttt{qt}(0.165,18)=-1$ (or maybe just close due to round off):
 
```{r}
qt(0.165,18)
```

The interpretation is that just over $16.5\%$ of the population lies below $x=-1$ in this $t$-distribution.

Note that the `qt` command is quite handy for finding $t^*$ multiplier values, as 
we'll see below.

### The $t$-score and the $p$-value

The $t$-score for a $t$-distribution is computed in much the same way that a $z$-score is computed for a normal distribution. That is, given a sample of size $n$, with mean $\bar{x}$, and standard deviation $s$ the $t$-score of the observation $x$ is
$$T = \frac{x-\bar{x}}{s/\sqrt{n}}.$$
In our dolphin example above, the assumed mean is $3.6$ and the observed mean is $4.513684$. We can compute the standard error easily enough:

```{r}
n = length(d)
se = sd(d)/sqrt(n)
se
```

Thus, the $t$-score is
```{r}
t_score = (mean(d)-3.6)/se
t_score
```
In agreement with the first line from our application of `t.test`.

The $p$-value is then computed from the $t$-score in the same way that the normal $p$-value is computed from the $z$-score. In this case, we have an observation of the mean that is larger than the assumed mean. Thus, the probability that the mean might be this large or larger, *assuming* the null value is
```{r}
1 - pt(t_score, 18)
```
Of course, the alternative hypothesis, as stated by the output of `t.test` is that the true mean is *unequal* to 3.6; thus, we should double this number to obtain the result in that output.

### Computing confidence intervals 

Computing a confidence interval by hand is fairly easy, if you have the right 
information. The general form is

$$[\bar{x} - ME, \bar{x} + ME],$$
where the margin of error $ME$ has the form
$$ME = t^* \times SE.$$


#### Getting $t^*$

We get $t^*$ from the $t$-distribution in the same way that we get $z^*$ from the
normal distribution. Thus, to find a
$95\%$ confidence interval, we examine the following picture.:

```{r confidence, echo=F}
dx = 0.01
x <- seq(-4,4, by=dx)
y <- dt(x, 18)
plot(x,y, type='l', ylim=c(0,0.6), yaxs='i', xlim=c(-4,4), xlab="", ylab="", yaxt='n', lwd=3)

q = qt(0.025, 18)
x <- seq(-4,q, by=dx)
y <- dt(x, 18)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='#eeeeee')
x <- seq(-q,4, by=dx)
y <- dt(x, 18)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='#eeeeee')

title(main = 'A 95% confidence interval')
```

Note that the gray portions must contain $5\%$ of the area under the curve. As a result, the endpoints of the interval can be computed as follows:

```{r}
c(qt(0.025, 18), qt(0.975,18))
```

Of course, you should adjust the 18 to whatever value of $df$ you need.

Note that we can also find that $t^*=2.1$ in the $t$-table on our webpage, 
where we see something that looks like so:

<style>
table {
  text-align: left;
  line-height: 40px;
  border-collapse: separate;
  border-spacing: 0;
  border: 2px solid #ed1c40;
  width: 500px;
  margin: 50px auto;
  border-radius: .25rem;
}

thead tr:first-child {
  background: #ed1c40;
  color: #fff;
  border: none;
}

th:first-child,
td:first-child {
  padding: 0 15px 0 20px;
}

thead tr:last-child th {
  border-bottom: 3px solid #ddd;
}

tbody tr:hover {
  background-color: rgba(237, 28, 64, .1);
  cursor: default;
}

tbody tr:last-child td {
  border: none;
}

tbody td {
  border-bottom: 1px solid #ddd;
}

td:last-child {
  text-align: right;
  padding-right: 10px;
}

.button {
  color: #aaa;
  cursor: pointer;
  vertical-align: middle;
}

.edit:hover {
  color: #0a79df;
}

.delete:hover {
  color: #dc2a2a;
}
</style>

<table>
<tbody>
<tr class="odd">
<td align="right">one tail</td>
<td align="right"> 0.100</td>
<td align="right"> 0.050</td>
<td align="right"> 0.025</td>
<td align="right"> 0.010</td>
<td align="right"> 0.005</td>
</tr>
<tr class="even">
<td align="right">two tails</td>
<td align="right">0.200</td>
<td align="right">0.100</td>
<td align="right">0.050</td>
<td align="right">0.020</td>
<td align="right">0.010</td>
</tr>
<tr class="odd">
<td align="right"><span>df</span> 1</td>
<td align="right"><span>3.08</span></td>
<td align="right"><span>6.31</span></td>
<td align="right"><span>12.71</span></td>
<td align="right"><span>31.82</span></td>
<td align="right"><span>63.66</span></td>
</tr>
<tr class="even">
<td align="right">2</td>
<td align="right"><span>1.89</span></td>
<td align="right"><span>2.92</span></td>
<td align="right"><span>4.30</span></td>
<td align="right"><span>6.96</span></td>
<td align="right"><span>9.92</span></td>
</tr>
<tr class="odd">
<td align="right">...</td>
<td align="right"><span>...</span></td>
<td align="right"><span>...</span></td>
<td align="right"><span>...</span></td>
<td align="right"><span>...</span></td>
<td align="right"><span>...</span></td>
</tr>
<tr class="even">
<td align="right">18</td>
<td align="right"><span>1.33</span></td>
<td align="right"><span>1.73</span></td>
<td align="right"><span>2.10</span></td>
<td align="right"><span>2.55</span></td>
<td align="right"><span>2.88</span></td>
</tr>
</tbody>
</table>

The entries in this table are called critical $t^*$ values. The columns indicate several common choices for confidence level and are alternately labeled either one-sided or two. The rows correspond to degrees of freedom. Now, look in the row where $df=18$ and where the two-sided test is equal to 0.05. we see that $t^*=2.1$. 

#### Our confidence interval

With $t^*$ in hand, we can compute the margin of error:

```{r}
m = mean(d)
me = 2.1*sd(d)/sqrt(19)
c(m-me, m+me)
```

Again, in agreement of the results of `t.test`.

## Another example

### Chips Ahoy!

In 1998, Nabisco announced its 1000 chip challenge, claiming that every bag of Chips
Ahoy cookies contained at least 1000 chips. A group of students counted the chips in
16 bags of cookies and found the following data:

```{r}
chips = c(1219,1214,1087,1200,1419,1121,1325,1345,1244,1258,1356,1132,1191,1270,1295,1135)
length(chips)
```

 * Create a $99\%$ confidence interval for the number of chips in the bag.
 * What is the $p$-value that the average number of chips is larger than 1000?

#### Solutions

There are a couple of approaches that we should understand.

**Solution 1**: The easiest approach is to just use the `t.test`.

For the first part a simple call to `t.test` with the `conf.level` parameter set to $0.99$ should do it.

```{r}
t.test(chips, conf.level=0.99)
```

I guess the confidence interval is $[1168.73, 1307.64]$. The $p$-value that we
see doesn't really help because we haven't set the null value or alternative
hypothesis. We can do this like so:

```{r}
t.test(chips, conf.level=0.99, alternative = "greater", mu=1000)
```

It looks like Nabisco is safe!


**Solution 2**: In order to do the HW, you'll need to
understand this from a slightly different perspective. From that perspective, all
you really need to know is the mean and standard deviation of the data, rather than 
the whole data set. Of course, we can compute this with R:

```{r}
c(mean(chips), sd(chips))
```

In the 2nd HW problem, that's *all* you're given. As before, the 99\% confidence 
interval should have the form
$$[1238.188 - ME, 1238.188 + ME].$$
The margin of error $ME$ has the form
$$ME = t^* \times SE,$$
where $SE = \sigma/\sqrt{n} = 94.282/\sqrt{16} \approx 23.5705$ and $t^*$ is chosen
from the $t$-distribution with 15 degrees of freedom so that it's corresponding 
two sided $p$-value is $0.01$. We can do this in R as follows:

```{r}
qt(0.005,15)
```
We can also look this up in a table to see that the value is about $2.95$ putting
this all together, we find that our confidence interval is
$$[1238.188 - 2.95\times23.5705, 1238.188 + 2.95\times23.5705] \approx [1168.655,1307.721].$$

To compute the $p$-value for the alternative hypothesis that the number of the
chips in the bag is greater than 1000, we first compute the $t$-statistic:

$$\frac{x-\bar{x}}{s/\sqrt(n)} = \frac{1238.188 - 1000}{94.282/\sqrt{16}} \approx 10.10534.$$

Finally, note that `1-pt(10.10534, 15)` return the following computation suggests
that we must reject the null hypothesis that the number of chips is *only* 1000.

```{r}
1-pt(10.10534,15)
```
