# GAM fitting to rock data

# Downlad the package gam from the appropriate CRAN mirror site, then:
library(gam)
attach(rock)

# First do a linear fit:
rock.lm = lm(log(perm) ~ area + peri + shape)



rock.gam1 = gam(log(perm) ~ s(area) + s(peri) + s(shape))
# Omitting the s() will result in linear terms being fitted.
# An option is lo() for loess fits.

# Compare the two fits:
anova(rock.lm, rock.gam1)

par(mfrow=c(2,2))
plot.gam(rock.gam1, se=T)


# Fit only one nonlinear term:
rock.gam2 = gam(log(perm) ~ area + peri + s(shape))
win.graph()
dev.set(which=3)
plot.gam(rock.gam2, se=T)
anova(rock.lm, rock.gam2, rock.gam1)
## tests model 1 against model 2 against model 3


######## PPR fitting 

shape1 = 100*shape 
area1 = area/400
peri1 = peri/100
c(mean(shape1), mean(area1), mean(peri1))

rock.ppr1 = ppr(log(perm) ~ area1 + peri1 + shape1, nterms = 1, max.terms = 4)
summary(rock.ppr1)
rock.ppr2 = ppr(log(perm) ~ area1 + peri1 + shape1, nterms = 2, max.terms = 4)
summary(rock.ppr2)
rock.ppr3 = ppr(log(perm) ~ area1 + peri1 + shape1, nterms = 3, max.terms = 4)
summary(rock.ppr3)
rock.ppr4 = ppr(log(perm) ~ area1 + peri1 + shape1, nterms = 4, max.terms = 4)
summary(rock.ppr4)


rock.lm2 = lm(log(perm) ~ (area1 + peri1)^2 + I(area1^2) + I(peri1^2))
summary(rock.lm2)

######## Plotting

par(mfrow=c(3,2))
plot(rock.ppr2)
plot(rock.ppr2$fitted, rock.ppr2$resid)
plot(rock.gam2$fitted, rock.gam2$resid)
plot(rock.lm$fitted, rock.lm$resid)
plot(rock.lm2$fitted, rock.lm2$resid)


win.graph()
par(mfrow=c(3,2))
plot(rock.ppr4)
plot(rock.ppr4$fitted, rock.ppr4$resid)