## ------------------------------------------------------------------------
#dat <- read.table(file = "http://www.simonqueenborough.info/R/data/sparrows.txt", header = TRUE)
dat <- read.table('~/Dropbox/MyLab/website/R/data/sparrows.txt', header = TRUE, sep = '\t')

## ---- sparrow-strip, fig.cap = "**Fig.** sparrow heads."-----------------
stripchart(dat$Head~dat$Sex, vertical = TRUE, method = 'jitter', col = c('darkgrey', 'red'), pch = 21)
points(x = 1:2, y = tapply(dat$Head, dat$Sex, mean), pch = 20, cex = 2)
summary(aov(Head ~ Sex, data = dat))

## ---- sparrow-strip2, fig.cap = "**Fig.** Including Tarsus"--------------
par(mfrow = c(1,2))

stripchart(dat$Head~dat$Sex, vertical = TRUE, method = 'jitter', col = c('darkgrey', 'red'), pch = 21)
points(x = 1:2, y = tapply(dat$Head, dat$Sex, mean), pch = 20, cex = 2)

plot(Head ~ Tarsus, data = dat[dat$Sex == 'Female', ], pch = 21)
points(Head ~ jitter(Tarsus), data = dat[dat$Sex == 'Male',] , col = 'red', pch = 21)
legend('topleft', pch = c(21, 21), col = c(1,2), legend = c('Female', 'Male'))

## ------------------------------------------------------------------------
m0 <- aov(Head ~ Sex + Tarsus, data = dat)
summary(m0)

## ------------------------------------------------------------------------
m1 <- aov(Head ~ Sex * Tarsus, data = dat)
summary(m1)

## ------------------------------------------------------------------------
summary(m1)

## ------------------------------------------------------------------------
summary.lm(m1)

## ------------------------------------------------------------------------
summary.lm(m1)$call

## ------------------------------------------------------------------------
summary(summary.lm(m1)$residuals)

## ------------------------------------------------------------------------
summary.lm(m1)$coef

## ---- fig.cap = "**Fig.** Head size in sparrows"-------------------------
F_intercept <- coef(m1)['(Intercept)']
F_slope     <- coef(m1)['Tarsus']
F_line       <- c(F_intercept, F_slope)

M_intercept <- F_intercept + coef(m1)['SexMale']
M_slope     <- F_slope + coef(m1)['SexMale:Tarsus']
M_line      <- c(M_intercept, M_slope)

par(lwd = 2, las = 1)
plot(Head ~ Tarsus, data = dat[dat$Sex == 'Female', ], pch = 21)
points(Head ~ jitter(Tarsus), data = dat[dat$Sex == 'Male',] , col = 'red', pch = 21)
legend('topleft', pch = c(21, 21), col = c(1,2), legend = c('Female', 'Male'))
abline(F_line, col = 'black')
abline(M_line, col = 'red')