#  Here is something more on logistic regression
#  Beetle mortality data on page 127
y=c(6,13,18,28,52,53,61,60)
n=c(59,60,62,56,63,59,62,60)
x=c(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)
n_y=n-y
beetle.mat=cbind(y,n_y)
# Figure 7.3
plot(x,y/n,xlab="Dose",ylab="Proportion Killed")
res.glm1=glm(beetle.mat~x,family=binomial(link="logit"))
fit_y=n*fitted.values(res.glm1)
D=2*sum( y *log(y/fit_y) + (n-y)*log((n-y)/(n-fit_y)) )
# Carefull with log((n-y))
p1=y *log(y/fit_y) 
p2=(n-y)*log((n-y)/(n-fit_y))
D=2*sum(p1)+2*sum(p2[-8])
#  Checking for null deviance
fit_y_0=n*(sum(y)/sum(n))
p1=y *log(y/fit_y_0)
p2=(n-y)*log((n-y)/(n-fit_y_0))
D=2*sum(p1)+2*sum(p2[-8])
# 
res.glm2=glm(beetle.mat~x,family=binomial(link="probit"))
res.glm3=glm(beetle.mat~x,family=binomial(link="cloglog"))
fit_y1=fit_y
fit_y2=n*fitted.values(res.glm2)
fit_y3=n*fitted.values(res.glm3)
result=cbind(y,fit_y1,fit_y2,fit_y3)
round(result,digits=2)
lines(x,fitted.values(res.glm2),col=2)
lines(x,fitted.values(res.glm3),col=3)
lines(x,fitted.values(res.glm1),col=1)
#
res.glm1$deviance
res.glm2$deviance
res.glm3$deviance
#
res.glm1$aic
res.glm2$aic
res.glm3$aic
#
r=rstandard(res.glm3)
m=influence.measures(res.glm3)
d=cooks.distance(res.glm3)
plot(d)
# Pearson's chi-square
pis=fitted.values(res.glm3)
X=sum ( ((y-n*pis)^2)/(n*pis*(1-pis)) )
1-pchisq(X,6)
# Computing the C statistic
# can be done from AIC
p=2
Lhat=(res.glm1$aic-2*p)/(-2)
res.glm0=glm(beetle.mat~1,family=binomial(link="logit"))
p=1
Ltilde=(res.glm0$aic-2*p)/(-2)
C=-2*(Ltilde-Lhat)
1-pchisq(C,1)
# For Rsq Ltilde needs a little modification 
dif=(Ltilde-Lhat)
Ltilde=Ltilde-sum(log(choose(n,y)))
Rsq=dif/Ltilde