pdf(file = "gam-plots.pdf")
 
 # stealing stuff from Chapter 5 of
 #     Simon Wood (2006)
 #     Generalized Additive Models: An Introduction with R
 #     Chapman & Hall/CRC Texts, Boca Raton, FL.
 # uses the R package mgcv which is one of the "recommended packages that
 # comes with every R installation (no need for install.packages or for
 # using the "packages" menu in the GUI's -- it is already installed).
 
 # note: in this package have multiple ways to specify smooths in formulas
 # s(x), s(x, y), te(x, y), ti(x, y), see help for s, te, ti, t2, smooth.terms
 
 # Cherry trees, see help(trees)
 
 # using the default, thin-plate splines
 
 library(mgcv)
 data(trees)
 out1 <- gam(Volume ~ s(Height) + s(Girth),
     family = Gamma(link = log), data = trees)
 print(out1)
 summary(out1)
 plot(out1, residuals = TRUE)
 
 # using cubic regression splines (number of knots fixed ???)
 
 out2 <- gam(Volume ~ s(Height, bs = "cr") + s(Girth, bs = "cr"),
     family = Gamma(link = log), data = trees)
 summary(out2)
 plot(out2, residuals = TRUE)

 predict(out1) - predict(out2)
 # not a lot of difference

 # using 2-D thin plate splines

 out3 <- try(gam(Volume ~ s(Height, Girth),
     family = Gamma(link = log), data = trees))
 # humpf, guess we have to follow the book
 out3 <- try(gam(Volume ~ s(Height, Girth, k = 25),
     family = Gamma(link = log), data = trees))
 # it works, but why do we need k = 25, what is it doing?
 summary(out3)
 plot(out3, residuals = TRUE)

 # is this any better ?
 anova(out1, out3)
 # apparently this is nonsense ?????????
 anova(out3)
 # no use
 # read p. 195 of Wood book !!!!!!!!!!!!!!!!!!!!!!!
 # would have to bootstrap to get P-value, but how ?????
 # just getting a valid P-value is an open research question

 # still following book
 out4 <- try(gam(Volume ~ te(Height, Girth, k = 5),
     family = Gamma(link = log), data = trees))
 summary(out4)
 plot(out4, residuals = TRUE)

 # still following book, new data set, brain image
 library(gamair)
 data(brain)
 brain <- brain[brain$medFPQ > 5e-3, ]
 bout0 <- gam(medFPQ ~ s(X, Y, k = 100), data = brain)
 summary(bout0)
 gam.check(bout0)
 # looks terrible so go with book.
 bout1 <- gam(medFPQ ~ s(X, Y, k = 100), data = brain,
     family = Gamma(link = log))
 summary(bout1)
 gam.check(bout1)

 vis.gam(bout1, plot.type = "contour")