#AM 020 GAM USING SPLINES
#LOADING ZUUR ET AL'S ISIT EXAMPLE DATA
#EXTRACTING DATA FOR STATIONS 8 & 13 ONLY
#IF PREFERRED THIS MIGHT BE DONE INSTEAD IN MS EXCEL
setwd("c:/DATA/Models")
ISIT=read.table("ISIT.txt",header=T)
ISIT
S8=ISIT$Sources[ISIT$Station==8]
D8=ISIT$SampleDepth[ISIT$Station==8]
S13=ISIT$Sources[ISIT$Station==13]
D13=ISIT$SampleDepth[ISIT$Station==13]
So=c(S8, S13); 
De=c(D8, D13)
ID=rep(c(8, 13), c(length(S8), length(S13)))
mi=max(min(D8), min(D13))
ma=min(max(D8), max(D13))
I1=De > mi & De < ma
ISIT813=data.frame(So=So[I1],De=De[I1],ID=ID[I1])
ISIT813
#NEW DATASET ISIT813 WRITTEN TO DISK
write.table(ISIT813,file="ISIT813.txt")


##IMPORTANT - RESTART R AT THIS POINT
#THIS IS AS IF LOADING THE REVISED DATASET ISIT813 FROM SCRATCH

#READING ISIT813 SUBSAMPLE DATA TABLE FROM DISK
setwd("c:/DATA/Models")
I=read.table("ISIT813.txt",header=T)
I

#PLOTS FROM ISIT813 DATASET	
op=par(mfrow = c(1, 2))
plot(I$De[I$ID==8], I$So[I$ID==8], pch = 16, xlab = "Depth",
    ylab = "Sources", col = 1, main = "Station 8",
    xlim = c(500, 3000), ylim = c(0, 40))
	
plot(I$De[I$ID==13], I$So[I$ID==13], pch = 16, xlab = "Depth",
    ylab = "Sources", col = 1, main = "Station 13",
    xlim = c(500, 3000), ylim = c(0, 40))
par(op)
	
#GAM MODEL FOR ISIT813 SINGLE CURVE WITH OFFSET ONLY	
library(mgcv)
attach(I)
fID=factor(ID)
M4=gam(So~s(De)+fID)
M4	
anova(M4)
summary(M4)

plot(M4)
gam.check(M4)

#GAM MODEL FOR ISIT813 WITH INTERACTION
M5=gam(So~s(De)+s(De,by=as.numeric(ID==13)))
#NOTE: CODE HERE IS MODIFIED ACCORDING TO ERRATA AT www.highstat.com
#NUMERICAL VALUES REPORTED CORRESPOND TO R OUTPUT FROM THE ONLINE
#R SCRIPT, BUT NOT TO THE VALUES REPORTED ON P.60
M5
anova(M5)
summary(M5)

op=par(mfrow = c(1, 2))
plot(M5)
par(op)

gam.check(M5)

#F-TEST COMPARING NESTED MODELS 
anova(M4,M5,test='F')
AIC(M4)
AIC(M5)

#TWO OTHER INTERACTION MODELS
M6a=gam(So~s(De,by = as.numeric(ID== 8))+s(De,by = as.numeric(ID== 13))-1)
anova(M6a)
summary(M6a)

M6b=gam(So~s(De,by = ID) + factor(ID))
anova(M6b)
summary(M6b)

anova(M5,M6a,M6b)
anova(M5,M6a)
anova(M5,M6b)
anova(M6a,M6b)	
AIC(M5)
AIC(M6a)
AIC(M6b)