#Exercise 3 (Logit and probit link)
beetle = read.table("http://www.uio.no/studier/emner/matnat/math/STK3100/h14/data/beetle.dat", header=F)
colnames(beetle) = c("dose", "tot", "dead")
attach(beetle)

beetle

#a: look at two glm models, and plot the fitted curve
glm1 <- glm(cbind(dead,tot-dead)~dose,family=binomial)           # I(x^2)
glm2 <- glm(cbind(dead,tot-dead)~dose+I(dose^2),family=binomial) # x + I(x^2)

plot(dose, glm1$fitted.values, type="l", col="blue", ylab="fitted values")
lines(dose, glm2$fitted.values, col="red")
points(dose, dead/tot)
legend("topleft",legend=c("x_i","x_i + x_i^2", "the data"),col=c("blue","red", "black"),lty=1)

summary(glm2) # H0: beta_I(dose^2) = 0
Zstat =  156.41/sqrt(3348.128)                      #\hat{beta2}/hat{sd_beta2}
p_value = pnorm(Zstat, lower.tail = FALSE)*2        #p-value 0.00687 < 0.05 , signifikant



#b: look at the covariance and correlation between the beta estimates
summary(glm2)$cov.scaled 
summary(glm2,correlation=T)$correlation

glm3 <- glm(cbind(dead,tot-dead)~I(dose^2),family=binomial)     #I(x^2)
summary(glm3)

#since x and x^2 is highly correlated, it is enough to include only one of the variables.



#c: look at the probit link model, compare the likelihood with the logit link model:
glm4 <-  glm(cbind(dead,tot-dead)~dose,family=binomial(link=probit))

plot(dose, glm1$fitted.values, type="l", col="blue")
lines(dose, glm4$fitted.values, col="red")
points(dose, dead/tot)
legend("topleft",legend=c("logit","probit", "observed data"),col=c("blue","red", "black"),lty=1)

logLik(glm1)
logLik(glm4)