# Bliss example page 155
library(faraway)
data(bliss)
modl=glm(cbind(dead,alive)~conc,family=binomial,bliss)
summary(modl)$coef


y=bliss$dead/30
mu=y
eta=logit(mu)
z=eta+(y-mu)/(mu*(1-mu))
w=30*mu*(1-mu)
lmod=lm(z~conc,weights=w,bliss)
coef(lmod)

for(k in 1:5){
eta=lmod$fit
mu=ilogit(eta)
z=eta+(y-mu)/(mu*(1-mu))
w=30*mu*(1-mu)
lmod=lm(z~conc,weights=w,bliss)
cat(k,coef(lmod),'\n')
}



summary(lmod)

xm=model.matrix(lmod)
wm=diag(w)
sqrt(diag(solve(t(xm) %*% wm %*% xm)))

summary(lmod)$coef[,2]/summary(lmod)$sigma

# tests of fit

sumary(modl)

1-pchisq(deviance(modl),df.residual(modl))

anova(modl,test='Chi')

modl2=glm(cbind(dead,alive)~conc+I(conc^2),family=binomial,bliss)
anova(modl,modl2,test='Chi')

anova(modl2,test='Chi')

# diagnostics

residuals(modl)

residuals(modl,'pearson')

residuals(modl,'response')

bliss$dead/30-fitted(modl)

residuals(modl,'working')

modl$residuals

residuals(lmod)

# leverage and influence

influence(modl)$hat

rstudent(modl)

influence(modl)$coef

cooks.distance(modl)

# new dataset

data(gala)
gala=gala[,-2]
modp=glm(Species~.,family=poisson,gala)

# plot residuals three ways (p. 162)

par(mfrow=c(1,3))
plot(residuals(modp)~predict(modp,type='response'),
xlab=expression(hat(mu)),ylab='Deviance Residuals')

plot(residuals(modp)~predict(modp,type='link'),
xlab=expression(hat(eta)),ylab='Deviance Residuals')

plot(residuals(modp,type='response')~predict(modp,type='link'),
xlab=expression(hat(mu)),ylab='Response Residuals')

# plot species v. area three ways (p. 164)

plot(Species~Area,gala)

plot(Species~log(Area),gala)

mu=predict(modp,type='response')
z=predict(modp)+(gala$Species-mu)/mu
plot(z~log(Area),gala,ylab='Linearized Response')


# also plot linearized response against Area (not good)

plot(z~Area,gala,ylab='Linearized Response')


# take logs in all independent variables and compare with earlier deviance 
modpl=glm(Species~log(Area)+log(Elevation)+log(Nearest)+log(Scruz+0.1)+log(Adjacent),
family=poisson,gala)
c(deviance(modp),deviance(modpl))

# modpl is much better as assessed by deviance

# partial residual plot
par(mfrow=c(1,1))

mu=predict(modpl,type='response')
u=(gala$Species-mu)/mu+coef(modpl)[2]*log(gala$Area)
plot(u~log(Area),gala,ylab='Partial Residual')
abline(0,coef(modpl)[2])



# diagnostic for link function
z=predict(modpl)+(gala$Species-mu)/mu
plot(z~predict(modpl),xlab='Linear Predictor',
ylab='Linearized Response')

# unusual points - page 166
halfnorm(rstudent(modpl))
# a.k.a. "jackknife residuals"

# leverage plot
gali=influence(modpl)
halfnorm(gali$hat)

halfnorm(cooks.distance(modpl))


modplr=glm(Species~log(Area)+log(Elevation)+log(Nearest)+log(Scruz+0.1)+
log(Adjacent),family=poisson,gala,subset=-25)
cbind(coef(modpl),coef(modplr))

# final model
modpla=glm(Species~log(Area)+log(Adjacent),family=poisson,gala)
dp=sum(residuals(modpla,type='pearson')^2)/modpla$df.res
sumary(modpla,dispersion=dp)