---
title: "Anova"
output: html_document
---

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

fish = case0602
```

First, let's plot the data

```{r exploratory_analysis}
bwplot(Pair~Percentage, data=fish)
densityplot(~Percentage, groups = Pair, auto.key = TRUE, data=fish)
```

We might be wondering whether the length matters, or whether there is a lot of variation between the different male pairs. Let's plot the data:

```{r plotlength}
plot(fish$Length, fish$Percentage)
```

It certainly doesn't look like there is a *linear* trend there... Recalling STA302, let's try to see anyway:

```{r linreg}
summary(lm(Percentage~Length, data=fish))

```

It doesn't look like there is evidence that the length matters for the proportion of time spent with the yellow-swordtailed male. (Which isn't to say that there isn't any nonlinear relationship; but it doesn't really look like that either.)


Now, let's compute the average percentage of time spent with the yellow-swordtailed male. (We are running two regressions, forcing one of them to have 0 intercept)


```{r runreg}
summary(lm(Percentage~Pair, data=fish))
summary(lm(Percentage~0+Pair, data=fish))
```

Note that really we are probably wondering about the differences between the different pairs because the males are matched by length. We can convert the lengths into factors, and run the regression again in the same way.

```{r runreg_factor}
summary(lm(Percentage~0+factor(Length), data=fish))
```

Now, the lengths were treated like *categorical* variables (i.e., like indicator variables.)




Let's compute a couple of values to make sure that we have things right

```{r confirm_reg}
coefs <- summary(lm(Percentage~Pair, data=fish))$coefficients[, "Estimate"]
c(mean(subset(fish, Pair=="Pair1")$Percentage), coefs["(Intercept)"])
c(mean(subset(fish, Pair=="Pair2")$Percentage)-
  mean(subset(fish, Pair=="Pair1")$Percentage), coefs["PairPair2"])
c(mean(subset(fish, Pair=="Pair3")$Percentage)-
  mean(subset(fish, Pair=="Pair1")$Percentage), coefs["PairPair3"])
```


## F-test

Let's go back to the boxplot
```{r boxplot}
bwplot(Pair~Percentage, data=fish)
```

Do the variances look equal here? I'd say they look a bit *fishy*, especially with Pair 3. I'm kind of wary of running an F-test here. On the bright side, the density plot indicates that the distributions are mostly normal.

Let's look at some diagnostic plots
```{r qqplot}
fit <- lm(Percentage~Pair, data=fish)
plot(fit, 1)
plot(fit, 2)
```

There are definitely outliers there -- maybe it's love!

Still, let's forge on -- it didn't stop one of the textbooks from trying

```{r anova}
anova(lm(Percentage~Pair, data=fish))
```

The F-test (which you may or may not trust) doesn't say that the pairs are different.

## Pairwise t-test

Let's run the unadjusted pairwise t-test:

```{r pairwise_t}
with(fish, pairwise.t.test(Percentage, Pair, p.adj = "none"))
```

No difference!

Let's try the Bonferroni correction (of course there will be no difference either...)

```{r pairwise_bonf}
with(fish, pairwise.t.test(Percentage, Pair, p.adj = "bonf"))

```

Finally, let's use Tukey's HSD:

```{r tukeyhsd}
TukeyHSD(aov(lm(Percentage ~ Pair, data= fish)), conf.level = 0.95)
```

## Confidence Intervals

Maybe we're interested in whether any of the pairs are significantly different from the mean. We can't just look at each one individually (though we can use Tukey's HSD). If we want to use Bonferroni, we do go

```{r bonf_confint}
fit <- lm(Percentage~Pair, data=fish)
confint(fit, level= 1-0.05/nlevels(fish$Pair))
```

## Finally, what about the yellow swordtails?

Let's check for normality one last time, and run the test

```{r swordtails}
hist(fish$Percentage)
t.test(fish$Percentage, mu=.5)
```
