# R-code for Snoring and Heart Disease 3.2 - 3.2.4
#                                                  by Yu-Ling Tseng

snoring <- read.table('snoringdata.txt')
names(snoring) <- c('snoring','disease', 'total')
snoring.x <- snoring$snoring       ## use scores 0, 2, 4, 5 for snoring levels 
snoring.y <- cbind(ndisease=snoring$disease,
                   nnormal=snoring$total-snoring$disease)
snoring.lprob <- glm(snoring.y ~ snoring.x,
                   family=binomial(link='identity'))
snoring.glm <- glm(snoring.y ~ snoring.x, family=binomial)
snoring.glm.probit <- glm(snoring.y ~ snoring.x,
                   family=binomial(link='probit'))

summary(snoring.lprob)      # c.f. p. 69
summary(snoring.glm)        # c.f. p. 71
summary(snoring.glm.probit) # c.f. p. 72

# Table 3.1 p. 69

tb31<-data.frame("snoring"=snoring.x, "H-Disease"=snoring$disease, 
"prop.Yes"=snoring$disease/snoring$total, 
"Linear fit"=predict(snoring.lprob, type='response'), 
"logit fit"=predict(snoring.glm ,type='response'),
"probit fit"=predict(snoring.glm.probit, type='response'))
tb31

# Figure 3.1  p. 70

plot(snoring$snoring, snoring$disease/snoring$total, ylim=c(0, .2))
legend(locator(1),c("Linear prob.","Logistic","Probit"),lty=c(3,2,4),lwd=c(2,2,2),col=c(1,2,4))
xx <- seq(0,5, length=100)
lines(xx,
      predict(snoring.lprob, newdata=data.frame(snoring.x=xx),
              type='response'),
      col=1, lty=3, lwd=2)
lines(xx,
      predict(snoring.glm, newdata=data.frame(snoring.x=xx), type='response'),
      col=2, lty=2,lwd=2)
lines(xx,
      predict(snoring.glm.probit, newdata=data.frame(snoring.x=xx),
              type='response'),
      col=4, lty=4, lwd=2)