# The data file extinction.txt contains measurements on
# breeding pairs of land-bird species collected from 16 islands are
# Britain over the course of several decades.  For each species, the
## data set contains an average time to extinction on those islands
# where it appeared (this is actually the reciprocal of the average of
# 1/T_i, where T_i is the length of time species i remained on
# the island, and 1/T_i is set to 0 if the species did not become
# extinct on the island.  This mean is sometimes called the geometric
# mean). Other variables are the average number of nesting pairs (the
# average, over all islands where the bird appeared, of the number of
# nesting pairs per year), the size of the species (either large or
# small), and the migratory status of the species (migrant or
# resident). The data are from Pimm, S.L., Jones, H.L., Diamond, J.
# (1988) On the risk of extinction, American Naturalist,
# 132, 757-85. It is expected that species with large numbers
# of nesting pairs will tend to remain longer before becoming extinct.
# Of interest is whether, after accounting for the number of nesting
# pairs, size or migratory status is related to extinction time.  There may also be
# some interest in whether the effect of size differs depending on the
# number of nesting pairs.


Data <- read.csv(file="C:/Documents and Settings/steeleb/My Documents/Math 543/extinction.csv",header=T,sep=",")

head(Data,5)
y <-Data$TIME

attach(Data)
n <- length(y)
plot(Data$PAIRS,y,ylab="Extinction time",xlab="N pairs",pch=16,col=2)

library(lattice)
xyplot(y~PAIRS|SIZE,groups=STATUS,pch=16,col=1:2)
legend("topright",legend=levels(STATUS),pch=16,col=1:2)

xyplot(log(y)~PAIRS|SIZE,groups=STATUS,pch=16,col=1:2,ylab="log-Extinction time",xlab="N pairs")
legend("topright",legend=levels(STATUS),pch=16,col=1:2)


panl <- rep("L and R",n)
panl <- replace(panl,(SIZE=="S")&(STATUS=="R"),"S and R")
panl <- replace(panl,(SIZE=="L")&(STATUS=="M"),"L and M")
panl <- replace(panl,(SIZE=="S")&(STATUS=="M"),"S and M")
xyplot(log(y)~Data$PAIRS|panl,pch=16,
	panel=function(x,y){
	panel.lmline(x,y,type="l")
	panel.xyplot(x,y,col=1,pch=16)})


levels(SIZE)
contrasts(SIZE) <- contr.treatment(n=2, base = 2, contrasts = TRUE)
levels(STATUS)
contrasts(STATUS) <- contr.treatment(n=2, base = 1, contrasts = TRUE)

y <- log(Data$TIME)
lm.obj <- lm(y~PAIRS+SIZE+STATUS)
model.matrix(lm.obj)[1:10,]
head(Data,10)

anova(lm.obj)
summary(lm.obj)
confint(lm.obj)


################## order matters - a little #################

anova(lm(y~PAIRS+STATUS+SIZE))
anova(lm(y~STATUS+SIZE+PAIRS))
anova(lm(y~PAIRS+STATUS*SIZE))


Residuals <- lm.obj$residuals
qqnorm(Residuals)
abline(a=0,b=1)
Fitted <- lm.obj$fitted
plot(Fitted,Residuals)
identify(Fitted,Residuals,labels=SPECIES)

xyplot(Fitted~PAIRS,groups=panl,col=1:4,pch=16)
legend("bottomright",pch=16,col=1:4,legend=levels(as.factor(panl)))

xyplot(Residuals~Fitted,groups=panl,col=1:4,pch=16)
legend("bottomright",pch=16,col=1:4,legend=levels(as.factor(panl)))