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

## Slide 5
densityplot(~sbp, Data, group=chd, auto.key=list(columns=2,text=c("No CHD","CHD")), xlab="Systolic blood pressure", xlim=c(80,250), cex=0.5, lwd=5)

## Slide 6
d0 <- with(Data, density(sbp[chd==0], bw="ucv"))
d1 <- with(Data, 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(Data$chd)

plot(x, p1*f1$y/(p1*f1$y+(1-p1)*f0$y), type="l", ylim=c(0,1), xlab="Systolic Blood Pressure", ylab="Pr(CHD)", lwd=3)
abline(h=.5,lty=2)
fit <- glm(chd~sbp, Data, family="binomial")
lines(x, predict(fit,data.frame(sbp=x), type="response"), col="blue", lwd=3)
legend("topleft", legend=c("Kernel","Logistic"), lty=1, col=c("black","blue"), lwd=3)
