require(mgcv)

## Heart
Data <- read.delim("../../data/heart.txt")
fit <- gam(chd~s(sbp)+s(tobacco)+s(ldl)+s(obesity)+s(adiposity)+s(typea)+s(alcohol)+s(age)+famhist,Data,family="binomial")
summ <- summary(fit)
exp(summ$p.table[2,1]+c(-1,0,1)*qnorm(.975)*summ$p.table[2,2])
sum(fit$edf)
pTerms <- c("famhist")
sTerms <- c("sbp","tobacco","ldl","obesity","adiposity","typea","alcohol","age")

## Slide 5-6
for (i in 1:8)
  {
    plot(fit,select=i,shade=TRUE,ylim=c(-5,5),col="slateblue",lwd=2,rug=FALSE)
    rug(jitter(Data[,sTerms[i]]),col="#6A5ACDBF")
    form0 <- as.formula(paste("chd~",paste(c(pTerms),collapse="+"),"+",paste("s(",sTerms[-i],",sp=",fit$sp[-i],")",sep="",collapse="+")))
    fit0 <- gam(form0,Data,family="binomial")
    aov0 <- anova(fit0,fit,test="Chisq")
    form1 <- as.formula(paste("chd~",paste(c(pTerms,sTerms[i]),collapse="+"),"+",paste("s(",sTerms[-i],",sp=",fit$sp[-i],")",sep="",collapse="+")))
    fit1 <- gam(form1,Data,family="binomial")
    aov1 <- anova(fit1,fit,test="Chisq")
    p <- aov1[2,5]
    mtext(paste("Overall effect: p =",formatP(aov0[2,5])),cex=0.8,line=2)
    if (aov1[2,3] > .1) mtext(paste("Linearity: p =",formatP(p)),cex=0.8,line=1)
  }

## Slide 16
fit <- gam(chd~s(age,ldl,k=9,fx=TRUE),Data,family="binomial")
vis.gam(fit,type="response",plot.type="contour",color="topo",n.grid=100,lty=0,main="Thin-plate spline",contour.col="black")

## Slide 17
par(mfrow=c(2,2))
fit <- gam(chd~age+ldl,Data,family="binomial")
vis.gam(fit,type="response",plot.type="contour",color="topo",n.grid=99,lty=0,main="GLM")
fit <- gam(chd~age*ldl,Data,family="binomial")
vis.gam(fit,type="response",plot.type="contour",color="topo",n.grid=99,lty=0,main="GLM w/ interaction")
fit <- gam(chd~s(age,k=5,fx=TRUE)+s(ldl,k=5,fx=TRUE),Data,family="binomial")
vis.gam(fit,type="response",plot.type="contour",color="topo",n.grid=99,lty=0,main="GAM")
fit <- gam(chd~s(age,ldl,k=9,fx=TRUE),Data,family="binomial")
vis.gam(fit,type="response",plot.type="contour",color="topo",n.grid=99,lty=0,main="Thin-plate")