heart <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/621/data/heart.txt")
require(locfit)

## Local likelihood vs. linear vs. kernel density
d0 <- with(heart, density(sbp[chd==0],bw="ucv"))
d1 <- with(heart, density(sbp[chd==1],bw="ucv"))
x <- seq(100,220,len=101)
f0 <- approx(d0$x,d0$y,xout=x)
f1 <- approx(d1$x,d1$y,xout=x)
p1 <- mean(heart$chd)
y1 <- p1*f1$y/(p1*f1$y+(1-p1)*f0$y)
fit.glm <- glm(chd~sbp, data=heart, family="binomial")
y2 <- predict(fit.glm, data.frame(sbp=x), type="response")
fit <- locfit(chd~lp(sbp), data=heart, family="binomial")
y3 <- predict(fit, data.frame(sbp=x), type="response")
Y <- cbind(y1, y2, y3)
col <- c("#FF4E3740","#00B50040","#008DFFFF")
matplot(x, Y, lty=1, lwd=3, col=col, type="l", xlab="SBP", ylab="Pr(CHD)", ylim=range(Y))
legend("topleft", c("Kernel\ndensity", "Logistic\nregression", "Local\nlikelihood"), col=col, lwd=3, ncol=3, bty="n")

## AIC
## SBP
AICp <- aicplot(chd~sbp, data=heart, family="binomial", alpha=seq(.2,1,.05), deg=1)
with(AICp, {
  plot(alpha, values, pch=19, xlab="Bandwidth", ylab="AIC", las=1, cex=0.7, type="o")
  ind <- seq(1, length(alpha), len=5)
  axis(3, at=alpha[ind], labels=formatC(df[ind], digits=1, format="f"))
})
mtext("Effective degrees of freedom", line=2.5)
## Age
AICp <- aicplot(chd~age, data=heart, family="binomial", alpha=seq(.2,1,.05), deg=1)
with(AICp, {
  plot(alpha, values, pch=19, xlab="Bandwidth", ylab="AIC", cex=0.7, type="o", las=1)
  ind <- seq(1, length(alpha), len=5)
  axis(3, at=alpha[ind], labels=formatC(df[ind], digits=1, format="f"))
})
mtext("Effective degrees of freedom", line=2.5)

## Confidence bands
p <- function(x){exp(x)/(1+exp(x))}
par(mfrow=c(1,2),mar=c(5,5,1,1))
ci <- scb(chd~lp(sbp), data=heart, family="binomial", simul=FALSE, ev=lfgrid(51))
plot(rep(ci$xev,2),c(ci$upper,ci$lower),type="n",xlab="Systolic blood pressure",ylab=expression(hat(theta)(x)),ylim=c(-2.2,3.5))
with(ci,polygon(c(xev,rev(xev)),c(lower,rev(upper)),col="gray80",border=FALSE))
lines(ci$xev,ci$coef)
plot(rep(ci$xev,2),p(c(ci$upper,ci$lower)),type="n",ylim=c(0,1),xlab="Systolic blood pressure",ylab=expression(hat(pi)(x)))
with(ci,polygon(c(xev,rev(xev)),p(c(lower,rev(upper))),col="gray80",border=FALSE))
lines(ci$xev,p(ci$coef))

ci <- scb(chd~lp(sbp), data=heart, family="binomial", simul=TRUE)
plot(rep(ci$xev,2),c(ci$upper,ci$lower),type="n",xlab="Systolic blood pressure",ylab=expression(hat(theta)(x)),ylim=c(-2.2,3.5))
with(ci,polygon(c(xev,rev(xev)),c(lower,rev(upper)),col="gray80",border=FALSE))
lines(ci$xev,ci$coef)
plot(rep(ci$xev,2),p(c(ci$upper,ci$lower)),type="n",ylim=c(0,1),xlab="Systolic blood pressure",ylab=expression(hat(pi)(x)))
with(ci,polygon(c(xev,rev(xev)),p(c(lower,rev(upper))),col="gray80",border=FALSE))
lines(ci$xev,p(ci$coef))

ci <- scb(chd~lp(age, nn=0.55, deg=1), data=heart, family="binomial", simul=FALSE, ev=lfgrid(51))
plot(rep(ci$xev,2),c(ci$upper,ci$lower),type="n",xlab="Systolic blood pressure",ylab=expression(hat(theta)(x)), ylim=c(-5,1))
with(ci,polygon(c(xev,rev(xev)),c(lower,rev(upper)),col="gray80",border=FALSE))
lines(ci$xev,ci$coef)
plot(rep(ci$xev,2),p(c(ci$upper,ci$lower)),type="n",ylim=c(0,1),xlab="Systolic blood pressure",ylab=expression(hat(pi)(x)))
with(ci,polygon(c(xev,rev(xev)),p(c(lower,rev(upper))),col="gray80",border=FALSE))
lines(ci$xev,p(ci$coef))

ci <- scb(chd~lp(age, nn=0.55, deg=1), data=heart, family="binomial", simul=TRUE)
plot(rep(ci$xev,2),c(ci$upper,ci$lower),type="n",xlab="Systolic blood pressure",ylab=expression(hat(theta)(x)), ylim=c(-5,1))
with(ci,polygon(c(xev,rev(xev)),c(lower,rev(upper)),col="gray80",border=FALSE))
lines(ci$xev,ci$coef)
plot(rep(ci$xev,2),p(c(ci$upper,ci$lower)),type="n",ylim=c(0,1),xlab="Systolic blood pressure",ylab=expression(hat(pi)(x)))
with(ci,polygon(c(xev,rev(xev)),p(c(lower,rev(upper))),col="gray80",border=FALSE))
lines(ci$xev,p(ci$coef))

require(gam)
fit2 <- gam(chd~lo(sbp), data=heart, family="binomial")
fit1 <- gam(chd~sbp, data=heart, family="binomial")
fit0 <- gam(chd~1, data=heart, family="binomial")
anova(fit0, fit1, fit2)

fit2 <- gam(chd~lo(age), data=heart, family="binomial")
fit1 <- gam(chd~age, data=heart, family="binomial")
fit0 <- gam(chd~1, data=heart, family="binomial")
anova(fit0, fit1, fit2)