# Generalized Additive Models (gam) demonstration
windows(record=T)
#
# Generate x values
n=250
p=3
x=((1:n)/n)^p
hist(x)
#
# Generate y.mean values
b=(1:n)/n
d=(x-.6)^2
y.mean=b+d
plot(x,y.mean,main="Expected Value [Y|x]",type="l",col="black",lwd=3,ylab="")
rug(x)
par(new=T)
plot(x,b,axes=F,ylim=range(y.mean,b,d),type="l",col="green",lwd=3,ylab="",xlab="")
par(new=T)
plot(x,d,axes=F,ylim=range(y.mean,b,d),type="l",col="red",lwd=3,ylab="",xlab="")
legend("topleft",legend=c("y.mean = b+d","b","d"),
        lwd=c(5,3,3),col=c("black","green","red"))
# Generate y values as noise on top of y.mean
sig=.15
set.seed(1)
y=y.mean+rnorm(n,sd=sig)
plot(x,y,main="x-y Scatterplot")
rug(x)
#
# Fit a gam with a spline
library(mgcv)                          
g0=gam(y~b+s(x))
plot(g0,main="g0")
summary(g0)
#
# Stiffen the spline
g2=gam(y~b+s(x,k=2))
plot(g2,main="g2")
summary(g2)
g4=gam(y~b+s(x,k=4))
plot(g4,main="g4")
summary(g4)
g6=gam(y~b+s(x,k=6))
plot(g6,main="g6")
summary(g6)
#
# Choose spline with k=2 and plot the fitted values
plot(x,y.mean,type="l",lwd=5,col="black")
rug(x)
par(new=T)
plot(x,g2$fitted.values,ylim=range(y.mean),type="l",col="red",lwd=2,ylab="",
      main="Fitted (g2) vs E[Y|x]")
legend("topleft",legend=c("Mean","Fitted"),,lwd=c(5,3),col=c("black","red"))