beetle <- read.table("c:/course/GLM/GLM_my-course-note/SAS_R-code/beetle.dat",header=T)
mortality <- beetle$y
n <- beetle$n
phat <- beetle$y/beetle$n

plot(beetle$dose,phat)
plot(beetle$dose, log(phat/(1-phat)))

# logistic regression model

beetle.reg <- glm(mortality/n ~ beetle$dose, family=binomial(link="logit"),weight=n)
summary(beetle.reg)

coef(beetle.reg)
vcov(beetle.reg)

fit.p <- c(fitted.values(beetle.reg))
fit.y <- n*fit.p

plot(beetle$dose,phat)
points(beetle$dose,fit.p,pch="f")

# Let's look at residuals
# The residuals() function gives the non-standardized deviance and Pearson residuals.
# The rstandard() function gives the standardized residuals.

residuals(beetle.reg)
residuals(beetle.reg,type="pearson")

plot(beetle$dose,rstandard(beetle.reg),ylab="standardized deviance residuals")


# cloglog model
beetle.log.reg <- glm(mortality/n ~ beetle$dose, family=binomial(link="cloglog"),weight=n)
summary(beetle.log.reg)

fit.p <- c(fitted.values(beetle.log.reg))
fit.y <- beetle$n*fit.p

plot(beetle$dose,phat)
points(beetle$dose,fit.p,pch="f")

residuals(beetle.log.reg)
residuals(beetle.log.reg,type="pearson")

plot(beetle$dose,rstandard(beetle.log.reg),ylab="standardized deviance")