#######################################################################
#######################################################################
#######################################################################
######                                                           ######
######             3rd VALUE training workshop, 2014             ######
######                                                           ######
######         Preparatory material for Rglimclim training:      ######
######                     self-study exercises                  ######
######                                                           ######
######              Analysis of Galapagos Islands data           ######
######                                                           ######
######  This script is a "clean" version of the file             ######
######  R_SelfStudy.r: it is designed to be run using            ######
######  source(), in which case it writes all results to the     ######
######  file Galapagos_Results.txt. See the header to file       ######
######  R_SelfStudy.r for more details.                          ######
######                                                           ######
#######################################################################
#######################################################################
#######################################################################
##
##	Divert output to file 
##
sink("Galapagos_Results.txt")
##
##	First step: read the data and produce summaries of each variable
##
tortoise.data <- read.table("galapagos.dat",header=TRUE)
cat("\nSUMMARY OF DATA SET:\n")
cat("====================\n")
print(summary(tortoise.data))
##
##	Plot the data ...
##
plot(tortoise.data)
##
##	... and identify the outlying data point
##
big.island <- (tortoise.data$Area > 3000)
cat("\nOUTLYING DATA POINT:\n")
cat("====================\n")
print(tortoise.data[big.island,])
##
##	Plot Endemics vs Elevation (publication-quality), and output to PDF
##	
plot(tortoise.data$Elevation,tortoise.data$Endemics,
     xlab="Elevation (m)",ylab="No. of species",
     main="Variation of endemic species numbers with\nisland elevation",
     pch=15,col="blue")
box(lwd=2)
dev.copy(pdf,"endemics.pdf",width=6,height=6)
dev.off()
#
#	Fit a linear regression model and output results
#
endemics.model <- lm(Endemics ~ Elevation, data=tortoise.data)
cat("\nREGRESSION OF ENDEMICS UPON ELEVATION:\n")
cat("======================================\n")
print(summary(endemics.model))
#
#	Produce the previous plot again, with the fitted regression 
#	line superimposed
#
plot(tortoise.data$Elevation,tortoise.data$Endemics,
     xlab="Elevation (m)",ylab="No. of species",
     main="Variation of endemic species numbers with\nisland elevation",
     pch=15,col="blue")
box(lwd=2)
abline(endemics.model, col="red", lty=2, lwd=2)
#
#	And plot some diagnostics (R produces 4 plots by default, 
#	so make space for them)
#
par(mfrow=c(2,2))
plot(endemics.model)
##
##	Return output to Console, and reset graphics window for 1 plot 
##	per page
##
sink(); par(mfrow=c(1,1))
