#====================================================
# Illustration of deviance for logistic regression 
# with grouped and ungrouped 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)

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


# Read ungrouped beetle data:
data="http://www.stat.ufl.edu/~aa/glm/data/Beetles.dat"
beetle.u = read.table(data, header = T)

# Fit logistic regression model based on ungrouped data:
fit.u=glm(y~x,family=binomial,data=beetle.u)
# Summary for ungrouped data:
summary(fit.u)

# Note that the deviances for ungrouped data differ from 
# the model with grouped data, but differences between 
# deviances remain the same ($null.deviance gives the deviance
#for the null modell with only an intercept):
c(fit$null.deviance,fit$deviance,fit$null.deviance-fit$deviance)
c(fit.u$null.deviance,fit.u$deviance,fit.u$null.deviance-fit.u$deviance)

# Likelihood ratio test for grouped and ungrouped data:
fit.null=glm(cbind(dead,n-dead)~1,family=binomial,data=beetle)
anova(fit.null,fit,test="LRT")
fit.u.null=glm(y~1,family=binomial,data=beetle.u)
anova(fit.u.null,fit.u,test="LRT")

#We did a likelihood ratio test using deviance 
#for the low birthweight data in "week37-commadns-testing.R",
#for M1: Model with lwt,smoke,race,ht (but no age)
#and M2: Model with lwt,smoke,race,ht and age

# Read data:
data="http://www.uio.no/studier/emner/matnat/math/STK3100/data/lowbirthweight.txt"
lowbirthweight=read.table(data,header=T)

# Fit logistic regressions M2 and M1:
fit.all=glm(low~age+lwt+factor(smoke)+race+factor(ht),family=binomial,data=lowbirthweight)
fit.noage=glm(low~lwt+smoke+factor(race)+factor(ht),family=binomial,data=lowbirthweight)

# Likelihood ratio test using deviance (details about deviance later)
anova(fit.noage,fit.all,test="LRT")

#Could fex compare more M1 to models, fex all models 
#with one less covariate than M1
drop1(fit.noage,test="LRT")

#This we could have done for M2 as well 
#(including the LRT test for M1 and M2 above)
drop1(fit.all,test="LRT")

#Can compare all types of nested models (can differ more than just one covariate) 
# for example M1 and a model M3 with age and ftv 
#(number of physician visits during the first trimester ):
fit.m3=glm(low~age+lwt+factor(smoke)+race+factor(ht)+ftv,family=binomial,data=lowbirthweight)
anova(fit.noage,fit.m3,test="LRT")


#===================================================
# Illustration of residuals with grouped beetle data
#===================================================

# Plot proportions dead and fitted model
plot(beetle$logdose,beetle$dead/beetle$n,pch=16,col='blue',
     ylim=c(0,1),xlab="log(dose)",ylab="Proportion dead")
x.grid = seq(1.65, 1.95, 0.01)
pi.hat = predict(fit, newdata=data.frame(logdose=x.grid), type="response")
lines(x.grid,pi.hat,lwd=2,col='red')

# Deviance residuals for grouped data:
residuals(fit)

# Pearson residuals for grouped data
residuals(fit,"pearson")

# Plot residuals for grouped data versus logdose:
plot(beetle$logdose,residuals(fit),pch=16,col='blue',
     ylim=c(-2,2),xlab="log(dose)",ylab="Residuals")
points(beetle$logdose,residuals(fit,"pearson"),pch=17, col='red')
legend("bottomright",c("deviance","pearson"),pch=c(16,17),col=c("blue","red"))

# To see the tendency, we add a local regression plot
lines(lowess(beetle$logdose,residuals(fit)),col="blue",lwd=2)
lines(lowess(beetle$logdose,residuals(fit,"pearson")),col="red",lwd=2)

# Residual plots for grouped data from the car package 
# (package must be installed)
library(car)
residualPlots(fit)


#====================================================
# Illustration of residuals with low birthweight data
#====================================================

# Plot residuals versus linear predictor:
plot(predict(fit.noage),residuals(fit.noage),pch=16,col='blue',
     ylim=c(-3,3),xlab="Linear predictor",ylab="Residuals")
points(predict(fit.noage),residuals(fit.noage,"pearson"),pch=17, col='red')
legend("bottomright",c("deviance","pearson"),pch=c(16,17),col=c("blue","red"))


# Plot residuals versus mother's weight:
plot(lowbirthweight$lwt,residuals(fit.noage),pch=16,col='blue',
     ylim=c(-3,3),xlab="Mother's weight",ylab="Residuals")
lines(lowess(lowbirthweight$lwt,residuals(fit.noage)),col="blue",lwd=2)

# Residual plots from the car package 
library(car)
residualPlots(fit.noage)