#===========================================================
# Illustraton of different link-functions for the beetle data
#============================================================

# Read grouped beetle data from the web:
data="http://www.stat.ufl.edu/~aa/glm/data/Beetles2.dat"
beetle = read.table(data, header = T)

# Plot the proportion deads a function of dose:
plot(beetle$logdose,beetle$dead/beetle$n,pch=16,col='blue',
     ylim=c(0,1),xlab="log(dose)",ylab="Proportion dead")


#=================================
# Logistic regression (logit-link)
#=================================

# Fit logistic regression with logdose as covariate:
fit.logistic=glm(cbind(dead,n-dead)~logdose,family=binomial,data=beetle)
summary(fit.logistic)

# Plot fitted line:
x.grid = seq(1.65, 1.95, 0.01)
pi.hat.logistic = predict(fit.logistic, newdata=data.frame(logdose=x.grid), type="response")
lines(x.grid,pi.hat.logistic,lwd=2,col='red')

#=========================
# Probit model(probit-link)
#=========================

#Fit probit model:
fit.probit=glm(cbind(dead,n-dead)~logdose,family=binomial(link=probit),data=beetle)
summary(fit.probit)
#Estimates almost equal to the ones for logistic divided by 1.8 (ref lecture):
fit.logistic$coef/1.8


# Plot fitted line:
pi.hat.probit = predict(fit.probit, newdata=data.frame(logdose=x.grid), type="response")
lines(x.grid,pi.hat.probit,lwd=2,col='black')

#===========================
# Complementary log-log link
#===========================

#Fit model with complementary log-log link:
fit.cloglog=glm(cbind(dead,n-dead)~logdose,family=binomial(link=cloglog),data=beetle)
summary(fit.cloglog)

# Plot fitted line:
pi.hat.cloglog = predict(fit.cloglog, newdata=data.frame(logdose=x.grid), type="response")
lines(x.grid,pi.hat.cloglog,lwd=2,col='green')