---
title: "t-Tests"
output: html_document
---

```{r setup, include=TRUE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(Sleuth3)
require(ggplot2)
require(tigerstats)
require(xtable)

finches <- case0201
```

## The two-sample t-Test

We can now compute the pooled SD estimate of the difference between the pre-draught and post-draught samples

```{r pooledvar}
X <- subset(finches, Year == 1976)$Depth
Y <- subset(finches, Year == 1978)$Depth

pooled_sd <- sqrt((sum((X-mean(X))^2) + sum((Y-mean(Y))^2))/(length(X)+length(Y)-2))
```

Let's compute some p-values for the Null Hypothesis that $\mu_X = \mu_Y$. We have that

$$\frac{\bar{X}-\bar{Y}}{s_p\sqrt{1/N_X+1/N_Y}}$$

is t-distributed with $N_X+N_Y-2$ d.f.

```{r pval}
Nx = length(X)
Ny = length(Y)
est = mean(X) - mean(Y)
pooled_sd = sqrt((sum((X-mean(X))^2) + sum((Y-mean(Y))^2))/(length(X)+length(Y)-2))
diff_sd = pooled_sd*sqrt(1/Nx + 1/Ny)
df = Nx + Ny - 2

```

Let's compute the two-sided and one-sided p-values. One-sided:

```{r pvals_manual_onesided}
df = Nx + Ny - 2
pt(est/diff_sd, df = df)
```

Two-sided:
```{r pvals_manual_twosided}
2*pt(est/diff_sd, df = df)
```

This was easy: in this case (and in most cases we'll see), the one-sided p-value is just twice the one-sided p-value.


It's possible to use R's `t-test` function to compute the same values:

```{r pvals_auto}
t.test(X, Y, alternative="two.sided", var.equal=TRUE)
t.test(X, Y, alternative="less", var.equal=TRUE)
```


Before, we were assuming equal variances. We didn't have to. `t-test` uses an approximation when the two variances are not assumed to be equal.

```{r pvals_auto_diff_vars}
t.test(X, Y, alternative="two.sided", var.equal=FALSE)
t.test(X, Y, alternative="less", var.equal=FALSE)
```




## Robustness of the t-test to deviations from assumptions
### Long tails

Generally, thanks to the CLT, the t-tests are pretty robust to deviations from assumptions. However, long-tailed distrbutions  will present problems.

Here's an example of a long-tailed distrbution. In this case, it looks basically like a normal distribution, except occasionally it has outliers. 

```{r longtail}
ind1 <- rbinom(100, 1, .02)
X1 <- (1-ind1)*rnorm(100, sd=5) + ind1*rnorm(100, sd=100)
hist(X1)
```
  
Let's simulate the t-test for this kind of case

```{r longtailsimul}
p = c()
for (i in 1:1000){
  ind1 <- rbinom(100, 1, .01)
  X1 <- (1-ind1)*rnorm(100, sd=5) + ind1*rnorm(100, sd=100)
  
  ind2 <- rbinom(100, 1, .01)
  X2 <- (1-ind2)*rnorm(100, sd=5) + ind2*rnorm(100, sd=100)
  
  p <- c(p, t.test(X1, X2, alternative="two.sided", var.equal=FALSE)$p.value)
}

mean(p<.05)

```
As you can see, we rejected the null hypothesis less often than expected -- not a problem if we're only worried about Type I errors, but a big problem if we're worried about Type II errors! 


###Non-independence

Let's now make our samples non-independent (What do you think the effect would be?)

```{r non-independence}
p = c()
for (i in 1:1000){
  X1 <- c(rep(rnorm(1), 20 ), rep(rnorm(1), 20 ), rep(rnorm(1), 20 ), rep(rnorm(1), 20 ), rep(rnorm(1), 20 ))
  X2 <- c(rep(rnorm(1), 20 ), rep(rnorm(1), 20 ), rep(rnorm(1), 20 ), rep(rnorm(1), 20 ), rep(rnorm(1), 20 ))
  
  p = c(p, t.test(X1, X2, alternative="two.sided", var.equal=FALSE)$p.value)
}

mean(p<.05)
```

A disaster! Of course, the data looks anything but normal if you plot it:

```{r histnonindep}
hist(X1)
```


### Multiplicative effects

Rainfall was measured on days when cumulus clouds were seeded and unseeded with silver iodide.

```{r showclouds, message=FALSE}
clouds <- case0301
ggplot(clouds, aes(x=Rainfall)) + 
  geom_histogram(data=subset(clouds, Treatment == "Unseeded"), fill = "red", alpha = 0.2) + 
  geom_histogram(data=subset(clouds, Treatment == "Seeded"), fill = "blue", alpha = 0.2)
```

There's clearly a difference there, but the disibutions are very skewed. Can we do better if we apply the log transformation to the data?


```{r showclouds2, message=FALSE}
clouds[, "lograin"] <- log(1+clouds$Rainfall)
ggplot(clouds, aes(x=lograin)) + 
  geom_histogram(data=subset(clouds, Treatment == "Unseeded"), fill = "red", alpha = 0.2) + 
  geom_histogram(data=subset(clouds, Treatment == "Seeded"), fill = "blue", alpha = 0.2)
```

Looks better, but in this case a density plot would be better

```{r clouddensity}
densityplot(~lograin, groups = Treatment, auto.key = TRUE, data = clouds)

```

Close enough!

```{r ttests}
t.test(lograin~Treatment, alternative="two.sided", data=clouds)
```

Interpretation: seeding clouds increases rainfall by a factor of somewhere between 1.3 and 7.4.


Let's try the wrong t-test, just for fun

```{r ttests_bad}
t.test(X1, X2, alternative="two.sided", var.equal=FALSE)
```

The p-value in this case is meaningless, but let's point and laugh: it's larger than 0.05.

## The Sampling Distribution of the Sample Variance

Let's now return to guessing at an appropriate standard deviation. 
```{r estsd}
upper = (89*2-2)/(pooled_sd^2*qchisq(.05, df=2*89-2))
lower = (89*2-2)/(pooled_sd^2*qchisq(.95, df=2*89-2))
c(lower, upper)
```

So a guess of 1.0 is completely appropriate.
