########## R script: galapagos.Rs ##########

# For performing GLM analysis of the
# data on species diversity on the Galapagos
# Islands.

# Last changed: 31 MAR 2008

# Define wait() function to permit pausing.

wait <- function() 
{
   cat("Hit Enter to continue\n")
   ans <- readline()
}

# Set flag for saving graphics in PDF files.

createPDF <- FALSE

# Read in data from file.

galapagos <- read.table("galapagos.txt",header=TRUE)

# Fit a Poisson regression model.

fitPois <- glm(Species~.,data=galapagos,family=poisson)
print(summary(fitPois))
wait()

# Obtain vectors useful for diagnosis.

y <- galapagos$Species
etaHat <- predict(fitPois,type="link")
muHat <-  predict(fitPois,type="response")
pearsonResids <- residuals(fitPois,type="pearson")
devianceResids <- residuals(fitPois,type="deviance")

# Do some basic residual plots

if (createPDF) pdf("galapagos1.pdf",width=8)
par(mfrow=c(2,2))
qqnorm(pearsonResids,main="Pearson residuals")
abline(0,1)
qqnorm(devianceResids,main="Deviance residuals")
abline(0,1)
plot(etaHat,pearsonResids)
abline(0,0)
plot(etaHat,devianceResids)
abline(0,0)
wait()

# Do a graphical check of the mean=variance assumption.

if (createPDF) pdf("galapagos2.pdf",width=8)
par(mfrow=c(1,1))
plot(log(muHat),log((y-muHat)^2))
abline(0,1)
wait()

# Fit an Poisson quasi-likelihood model.

fitQuasPois <- glm(Species~.,data=galapagos,family=quasipoisson)
print(summary(fitQuasPois))
wait()

# Do some basic residual plots

etaHat <- predict(fitQuasPois,type="link")
muHat <-  predict(fitQuasPois,type="response")
pearsonResids <- residuals(fitQuasPois,type="pearson")
devianceResids <- residuals(fitQuasPois,type="deviance")
if (createPDF) pdf("galapagos3.pdf",width=8)
par(mfrow=c(1,2))
plot(etaHat,pearsonResids)
abline(0,0)
plot(etaHat,devianceResids)
abline(0,0)
wait()

# Check that the mean/variance relationship
# using the over-dispersion parameter.

if (createPDF) pdf("galapagos4.pdf",width=8)
par(mfrow=c(1,1))
phiHat <- sqrt(summary(fitQuasPois)$dispersion)
plot(log(phiHat*muHat),log((y-muHat)^2))
abline(0,1)
wait()

# Do partial residual plots.

partialResids <- residuals(fitQuasPois,type="partial")
if (createPDF) pdf("galapagos5.pdf",width=8)

par(mfrow=c(2,3))
for (i in 1:ncol(partialResids))
{
   xvar <- galapagos[,i+1]
   yvar <- partialResids[,i]

   plot(xvar,yvar,main=dimnames(partialResids)[[2]][i],
        xlab=dimnames(partialResids)[[2]][i],
        ylab="partial residuals")

   fitSS <- smooth.spline(xvar,yvar,df=3)
   xvarg <- seq(min(xvar),max(xvar),length=101)

   lines(xvarg,predict(fitSS,xvarg)$y,col="red")
}   

# Do plots, but restriction of the range of "Area" and "Adjacent".

if (createPDF) pdf("galapagos6.pdf",width=8)
par(mfrow=c(2,3))
for (i in 1:ncol(partialResids))
{
   xvar <- galapagos[,i+1]
   yvar <- partialResids[,i]

   xlimVal <- range(xvar)
   if ((i==1)|(i==5)) xlimVal <- c(0,1000)

   plot(xvar,yvar,main=dimnames(partialResids)[[2]][i],
        xlab=dimnames(partialResids)[[2]][i],
        ylab="partial residuals",xlim=xlimVal)

   fitSS <- smooth.spline(xvar,yvar,df=3)
   xvarg <- seq(min(xvar),max(xvar),length=101)

   lines(xvarg,predict(fitSS,xvarg)$y,col="red")
}   

########## End of galapagos.Rs ##########