#logistic regression
#Generalized linear models, McCullagh and Nelder
#section 4.6

lizard<-c(20, 13, 8, 6, 34, 31, 17, 12, 2, 0, 3, 0, 11, 5, 15, 1,
8, 8, 4, 0, 69, 55, 60, 21, 1, 0, 1, 0, 20, 4, 32, 5,
4, 12, 5, 1, 18, 13, 8, 4, 4, 0, 3, 1, 10, 3, 8, 4)
lizard<-matrix(lizard,byrow=F,ncol=6)
> lizard
     [,1] [,2] [,3] [,4] [,5] [,6] 
[1,]   20    2    8    1    4    4
[2,]   13    0    8    0   12    0
[3,]    8    3    4    1    5    3
[4,]    6    0    0    0    1    1
[5,]   34   11   69   20   18   10
[6,]   31    5   55    4   13    3
[7,]   17   15   60   32    8    8
[8,]   12    1   21    5    4    4

y<-c(lizard[,1],lizard[,3],lizard[,5])
m<-y+c(lizard[,2],lizard[,4],lizard[,6])
	
d<-rep(c(0,0,1,1),6)  #diameter of perch
> d
 [1] 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1

h<-rep(c(0,1),12)  #height of perch
> h
 [1] 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1

s<-rep(c(rep(0,4),rep(1,4)),3)  #sunny (0) or shade(1)
> s
 [1] 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1

t<-factor(rep(c(1,2,3),rep(8,3)))  #time of day: early, mid-day, late
> t
 [1] 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3

fit1<-glm(cbind(y,m-y)~d+h+s+t,family=binomial)
> summary(fit1)
Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.66015  -0.37800   0.04488   0.62644   1.48717  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.9447     0.3415   5.695 1.23e-08 ***
d            -0.7626     0.2113  -3.610 0.000306 ***
h             1.1300     0.2571   4.395 1.11e-05 ***
s            -0.8473     0.3224  -2.628 0.008585 ** 
t2            0.2271     0.2502   0.908 0.363984    
t3           -0.7368     0.2990  -2.464 0.013730 *    # see table 4.8 of GLM
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 70.102  on 22  degrees of freedom
Residual deviance: 14.205  on 17  degrees of freedom
AIC: 83.029

Number of Fisher Scoring iterations: 4

# Estimated dispersion
14.205/17  # 0.8355882

## check for interactions
fit2<-glm(cbind(y,m-y)~d+h+s+t+h:s,family=binomial)
summary(fit2)   ## see table 4.7 in GLM

plot(predict(fit1), resid(fit1,type='pearson'),
 main="Residual plot for the logistic model")
abline(0,0)