---
title: "fruitfly"

output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r libraries}
library(faraway)
library(ggplot2)
library(dplyr)
``` 
This was a designed experiment: 125 fruitflies were divided into 5 groups of 25 each, and each group was assigned a different treatment. (It is conventional to refer to ``treatments'' in DoE, even if it's not a drug.) In this case the treatments were expsoure to female fruitflies. One group was solitary, another with a single virgin female, another with eight virgin females, another with 1 pregnant female and another with 8 pregnant females. Pregnant females will not mate.

The response variable was `longevity`, measured in days. A continuous covariate `thorax` was also recorded: it is the length of the thorax in millimeters, and known to be associated with longevity. The goal of the study was to see if sexual activity has a cost to male fruitflies, in terms of longevity.^[Partridge and Farquhar, 1981, *Nature* **284**, p.580.Available [here](http://www.utstat.utoronto.ca/reid/sta2101f/fruitfly-paper.pdf)]

First, a look at the data:
```{r explore, echo = TRUE}
head(fruitfly);dim(fruitfly)
is.factor(fruitfly$activity)
ggplot(aes(x=thorax, y=longevity), data = fruitfly) + geom_point()  + geom_smooth(method = lm) + facet_wrap(~ activity)
``` 

It is clear that there is a linear relationship between thorax length and longevity, as expected. 

Two questions of interest are whether there is a difference in longevity between activity groups, and whether the linear relationship between thorax and longevity is the same in each group. The latter is perhaps not that important in this example, but if it is the same, the resulting model will be simpler. 

```{r modelling}
lmod <- lm(longevity ~ thorax * activity, data = fruitfly)
summary(lmod)
```

(I've just done what you're not supposed to do, shown all the output in its ugliness.) The notation `a*b` on the right side is shorthand for `a + b + a:b`, where `a:b` is the *interaction* term.

How does `lm` get the output? Try below, and then look at `model.matrix(lmod2)` or `model.matrix(lmod)`. (They are long, you might want just the first 10 rows or so.)

With factor variables, an analysis of variance table is helpful, because you can see the amount of variation explained by the factor by an $F$-test. 

```{r testing} 
anova(lmod)
lmod2 <- lm(longevity ~ thorax + activity, data = fruitfly)
anova(lmod2)
anova(lm(longevity ~ activity + thorax, data = fruitfly))
summary(lmod2)
```

What is the difference between the two anova tables? 
Note that when we omit the interaction term, `activity` shows a large difference between levels; in particular those in the `high` category have a much reduced lifetime, relative to the baseline category. (What is the baseline category?)

The $p$-value for this in the table is $10^{-9}$ which shouldn't be taken very seriously (why not?). A table of sample means in the five groups is more informative (I found this code by googling):

```{r table of means}
group_by(fruitfly, activity) %>% summarize(mean = mean(longevity))
```

To test differences between any two group means we will need to estimate the standard error of the difference; this is not difficult, but a bit tedious. The original paper didn't do this, just reported a "$p < 0.001$" for the high activity group, and "$p < 0.05$" for the low activity group. 

Often with several groups the textbook recommendation is to do all pairwise comparisons of the means, and I'm sure that there is ```R``` code out there somewhere that will quickly do that for you, and offer a choice of adjustments for multiple testing.

For this particular experiment, though, the groups "one" and "many" were described as controls, as was the group "isolated", so the comparisons of interest are fairly limited. The  comparisons of most interest are probably "isolated" against "one" and "many", as these are slightly different types of controls, and then "low" and "high" against the average of the 3 control groups. The effect is so striking here, though, that this detail may not be needed. 