> data <- read.table("beetle.dat", header=T)

> data
   dose number killed 
1 1.691     59      6
2 1.724     60     13
3 1.755     62     18
4 1.784     56     28
5 1.811     63     52
6 1.837     59     53
7 1.861     62     61
8 1.884     60     60

> attach(data)
> binom.dat <- matrix(append(killed,number-killed),ncol=2) 
> fit.logit <- glm(binom.dat ~ dose, family=binomial(link=logit))
> fit.logit

Coefficients:
 (Intercept)     dose 
   -60.74002 34.28587

Degrees of Freedom: 8 Total; 6 Residual
Residual Deviance: 11.11558 

> fit.cloglog <- glm(binom.dat ~ dose, family=binomial(link=cloglog))
   # Fits the complementary log-log model
> summary(fit.cloglog)  #  much better fit than logit

                Value Std. Error   t value 
(Intercept) -39.52250   3.232269 -12.22748
       dose  22.01488   1.795086  12.26397

    Null Deviance: 284.2024 on 7 degrees of freedom

Residual Deviance: 3.514334 on 6 degrees of freedom

> predict <- fitted(fit.cloglog)  #  predicted proportions
> pearson <- (killed - number*predict)/sqrt(number*predict*(1-predict))
      #  Pearson residuals 

> matrix(append(append(dose,killed/number),append(predict,pearson)),ncol=4)
      [,1]      [,2]       [,3]       [,4] 
[1,] 1.691 0.1016949 0.09582053  0.1532963
[2,] 1.724 0.2166667 0.18802470  0.5678056
[3,] 1.755 0.2903226 0.33777042 -0.7899455
[4,] 1.784 0.5000000 0.54177561 -0.6274340
[5,] 1.811 0.8253968 0.75684020  1.2684452
[6,] 1.837 0.8983051 0.91843617 -0.5649628
[7,] 1.861 0.9838710 0.98575233 -0.1250008
[8,] 1.884 1.0000000 0.99913569  0.2278237





> binom2.dat <- matrix(append(number-killed,killed),ncol=2) # reverse S, F

> fit.loglog <- glm(binom2.dat ~ dose, family=binomial(link=cloglog))
               #  this is the log-log model for original data
> summary(fit.loglog)  # poorer fit than comp. log-log model

                Value Std. Error   t value 
(Intercept)  37.66401   2.947725  12.77731
       dose -21.58474   1.679123 -12.85477

    Null Deviance: 284.2024 on 7 degrees of freedom

Residual Deviance: 27.57268 on 6 degrees of freedom

> fit.logit2 <- glm(binom2.dat ~ dose, family=binomial(link=logit))
   # Now we fit logit model again, reversing S and F

> summary(fit.logit2)  #  same fit, estimates change sign 

                Value Std. Error   t value 
(Intercept)  60.74002   5.174975  11.73726
       dose -34.28587   2.909317 -11.78485

    Null Deviance: 284.2024 on 7 degrees of freedom

Residual Deviance: 11.11558 on 6 degrees of freedom

> q()