y = c(24,35,21,30)
n = c(1379,638,213,254)
snore = c(0,2,4,5)

snoring<-data.frame(snore,y,n)
M.logit<-glm(cbind(y,n-y)~snore,family=binomial(link=logit),data=snoring)
M.logit
coef(M.logit)
confint(M.logit)
summary(M.logit)
fitted(M.logit)


M.probit<-glm(cbind(y,n-y)~snore,family=binomial(link=probit),data=snoring)
M.probit
fitted(M.probit)


M.lm<-glm(y/n~snore,weight=n,data=snoring)
M.lm
fitted(M.lm)

beta0<-coef(M.logit)[1]
beta1<-coef(M.logit)[2]
plot(snore,y/n,type="p",pch=16)
abline(coef(M.lm),lty=2,col=2)
curve(exp(beta0+beta1*x)/(1+exp(beta0+beta1*x)),0,5,add=T,col=4)
beta01<-coef(M.probit)[1]
beta11<-coef(M.probit)[2]
curve(pnorm(beta01+beta11*x),0,5,add=T,col=8)
legend("topleft",legend=c(NA,"lm","logit","probit"),lty=c(NA,2,1,1),col=c(NA,2,4,8),bty="n")

fitted(M.logit)
fitted(M.probit)
fitted(M.lm)

#####################Snoring Example #####################
snoring <-matrix(c(24,1355,35,603,21,192,30,224), ncol=2, byrow=TRUE)
dimnames(snoring) <-list(snore=c("never","sometimes","often","always"),heartdisease=c("yes","No"))
scores.a <- c(0,2,4,5) 
scores.b <- c(0,2,4,6) 
scores.c <- 0:3 
scores.d <- 1:4 
snoring.lg.a <- glm( snoring ~ scores.a, family=binomial() ) 
summary(snoring.lg.a)
snoring.lg.b <- glm( snoring ~ scores.b, family=binomial() ) 
summary(snoring.lg.b)
snoring.lg.c <- glm( snoring ~ scores.c, family=binomial() ) 
summary(snoring.lg.c)
snoring.lg.d <- glm( snoring ~ scores.d, family=binomial ()) 
summary(snoring.lg.d)
coef(snoring.lg.a)
coef(snoring.lg.b)
coef(snoring.lg.c )
coef(snoring.lg.d) 
predict(snoring.lg.a, type="response" )
predict(snoring.lg.b, type="response" )
predict(snoring.lg.c, type="response" )
predict(snoring.lg.d, type="response" )
snoring.logit <-glm( snoring ~ scores.a, family=binomial(link="logit"))
summary(snoring.logit)
predict(snoring.logit, type="response")


#################################################################
## glm 
snore.f<-factor(1:4,labels=c("never","sometime","almost","all"))
snore.f

M<-glm(cbind(y,n-y)~snore.f,family=binomial(link=logit))
M
summary(M)
exp(coef(M))
exp(confint(M))


anova(M,test="Chisq")
anova(M,test="LRT")
anova(M.logit,M,test="Chisq")
anova(M.probit,M,test="Chisq")
anova(M.logit,M.probit,test="Rao")