require(gam)
bmd <- read.delim("../../data/bmd.txt")
heart <- read.delim("../../data/heart.txt")

## Basic usage:
fit <- gam(spnbmd~lo(age),data=bmd)
plot(fit)
fit <- gam(spnbmd~lo(age,span=.8,degree=2),data=bmd)
plot(fit)

## CHD
fit <- gam(chd~lo(sbp)+lo(tobacco)+lo(ldl)+lo(obesity)+lo(adiposity)+lo(typea)+lo(alcohol)+lo(age)+famhist,data=heart,family="binomial")
summary(fit)
par(mfrow=c(3,3))
plot(fit,se=TRUE)
plot(fit,se=TRUE,ask=TRUE)
fit <- gam(chd~lo(sbp)+lo(tobacco)+lo(ldl)+lo(obesity)+lo(adiposity)+lo(typea)+lo(alcohol)+lo(age)+famhist,data=heart,family="binomial")
## Overall effect: Tobacco
fit0 <- gam(chd~lo(sbp)+lo(ldl)+lo(obesity)+lo(adiposity)+lo(typea)+lo(alcohol)+lo(age)+famhist,data=heart,family="binomial")
anova(fit0,fit)

## Slide 9
gcv <- function(fit)
  {
    n <- length(fit$y)
    df <- 1+fit$df.null-fit$df.residual
    fit$deviance/(n*(1-df/n)^2)
  }
lsv <- function(fit)
  {
    n <- length(fit$y)
    df <- 1+fit$df.null-fit$df.residual
    crossprod(fit$y-fit$fitted)/(n*(1-df/n)^2)
  }
Span <- seq(.1,1,by=.02)
CV <- array(NA,dim=c(length(Span),3),dimnames=list(Span,c("GCV","AIC","LSV")))
for (i in 1:length(Span))
  {
    fit <- gam(chd~lo(tobacco,span=Span[i],degree=2),data=heart,family=binomial)
    CV[i,] <- c(gcv(fit),fit$aic,lsv(fit))
  }
for (j in 2:3) CV[,j] <- (CV[,j] - mean(CV[,j]))*sd(CV[,1])/sd(CV[,j]) + mean(CV[,1])
col <- c("black","blue","red")
matplot(Span,CV,type="l",lty=1,lwd=3,col=col)
for (j in 1:3) abline(v=Span[which.min(CV[,j])],col=col[j])
legend("top",col=col,legend=c("GCV (Deviance)","AIC","GCV (Least squares)"),lty=1,lwd=3,bg="white")

## Slide 18 (Extension)
m <- bmd[bmd$gender=="male",]
f <- bmd[bmd$gender=="female",]
Span <- seq(.1,1,by=.02)
Degree <- 1:2
AIC <- array(NA,dim=c(length(Span),2,2),dimnames=list(Span,paste("Deg",Degree,sep=""),c("M","F")))
Data <- list(m,f)
for (i in 1:length(Span))
  {
    for (j in 1:length(Degree))
      {
        for (k in 1:2)
          {
            fit <- gam(spnbmd~lo(age,span=Span[i],degree=Degree[j]),data=Data[[k]])
            AIC[i,j,k] <- fit$aic
          }
      }
  }
par(mfrow=c(2,2))
for (j in 1:2){for (k in 1:2) plot(Span,AIC[,j,k],type="l",main=paste(dimnames(AIC)[[3]][k],", Deg=",j,sep=""),ylab="AIC")}
Span[apply(AIC[,,1],2,which.min)]
Span[apply(AIC[,,2],2,which.min)]
fit.m <- gam(spnbmd~lo(age,span=.52),data=m)
fit.f <- gam(spnbmd~lo(age,span=.32),data=f)
par(mfrow=c(1,2))
plot(fit.m,se=TRUE)
plot(fit.f,se=TRUE)

## Slide 19
fit <- gam(chd~lo(sbp)+lo(tobacco)+famhist,data=heart,family="binomial")
summary(fit)
fit0 <- gam(chd~lo(sbp)+famhist,data=heart,family="binomial")
anova(fit0,fit)
fit0 <- gam(chd~lo(sbp)+lo(tobacco),data=heart,family="binomial")
anova(fit0,fit)

## Slide 21-22
require(locfit)
pband <- function(fit,xx,col.p="gray80",col.l="black",add=FALSE,...)
  {
    n <- fit$mi["n"]
    rdf <- n-2*fit$dp["df1"]+fit$dp["df2"]
    hat <- predict(fit,newdata=xx,se.fit=TRUE)
    lwr <- hat$fit - qt(.975,rdf)*hat$se
    upr <- hat$fit + qt(.975,rdf)*hat$se
    if (!add) plot(NULL,xlim=range(xx),...)
    polygon(c(xx,rev(xx)),c(lwr,rev(upr)),border=FALSE,col=col.p)
    lines(xx,hat$fit,lwd=2,col=col.l)
  }
sband <- function(fit,xx,col.p="gray80",col.l="black",add=FALSE,...)
  {
    s <- scb(fit,ev=lfgrid(70))
    if (!add) plot(NULL,xlim=range(xx),...)
    with(s,polygon(c(xev,rev(xev)),c(lower,rev(upper)),border=FALSE,col=col.p))
    with(s,lines(xev,coef,lwd=2,col=col.l))
  }
fit.m <- locfit(spnbmd~age,m)
fit.f <- locfit(spnbmd~age,f)
xx <- seq(min(bmd$age),max(bmd$age),len=101)
pband(fit.m,xx,ylim=c(-.04,.12),col.p=rgb(0,0,1,alpha=.5),col.l="blue")
pband(fit.f,xx,add=TRUE,col.p=rgb(1,0,0,alpha=.5),col.l="red")
sband(fit.m,xx,ylim=c(-.04,.12),col.p=rgb(0,0,1,alpha=.5),col.l="blue")
sband(fit.f,xx,add=TRUE,col.p=rgb(1,0,0,alpha=.5),col.l="red")

## Slide 24
par(mfrow=c(1,2))
fit <- gam(chd~lo(typea,age),data=heart,family=binomial)
summary(fit)
plot(fit,ticktype="detailed",theta=-30,col=rgb(.50,.42,.96),shade=.5)
fit <- gam(chd~lo(tobacco,ldl),data=heart,family=binomial)
summary(fit)
plot(fit,ticktype="detailed",theta=-30,col=rgb(.50,.42,.96),shade=.5)

## Slide 25
fit1 <- gam(chd~lo(sbp)+lo(tobacco)+lo(ldl)+lo(obesity)+lo(adiposity)+lo(typea)+lo(alcohol)+lo(age)+famhist,data=heart,family="binomial")
fit2 <- gam(chd~lo(sbp,tobacco,ldl,obesity,adiposity,typea,alcohol,age)+famhist,data=heart,family="binomial")