knitr::opts_chunk$set(echo = TRUE)
library(faraway)
library(SMPracticals)
## Loading required package: ellipse
## 
## Attaching package: 'ellipse'
## The following object is masked from 'package:graphics':
## 
##     pairs
## 
## Attaching package: 'SMPracticals'
## The following objects are masked from 'package:faraway':
## 
##     beetle, bliss, seeds
data(troutegg)
ftable(xtabs(cbind(survive,total) ~ location + period, data = troutegg))
##                  survive total
## location period               
## 1        4            89    94
##          7            94    98
##          8            77    86
##          11          141   155
## 2        4           106   108
##          7            91   106
##          8            87    96
##          11          104   122
## 3        4           119   123
##          7           100   130
##          8            88   119
##          11           91   125
## 4        4           104   104
##          7            80    97
##          8            67    99
##          11          111   132
## 5        4            49    93
##          7            11   113
##          8            18    88
##          11            0   138
bmod <- glm(cbind(survive,total-survive) ~ 
              location + period, 
            family = binomial, 
            data = troutegg)

summary(bmod)
## 
## Call:
## glm(formula = cbind(survive, total - survive) ~ location + period, 
##     family = binomial, data = troutegg)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.8305  -0.3650  -0.0303   0.6191   3.2434  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.6358     0.2813  16.479  < 2e-16 ***
## location2    -0.4168     0.2461  -1.694   0.0903 .  
## location3    -1.2421     0.2194  -5.660 1.51e-08 ***
## location4    -0.9509     0.2288  -4.157 3.23e-05 ***
## location5    -4.6138     0.2502 -18.439  < 2e-16 ***
## period7      -2.1702     0.2384  -9.103  < 2e-16 ***
## period8      -2.3256     0.2429  -9.573  < 2e-16 ***
## period11     -2.4500     0.2341 -10.466  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1021.469  on 19  degrees of freedom
## Residual deviance:   64.495  on 12  degrees of freedom
## AIC: 157.03
## 
## Number of Fisher Scoring iterations: 5
plot.glm.diag(bmod)

# plot(bmod)

bmod2 <- glm(cbind(survive,total-survive) ~ 
               location + period,
             family = quasibinomial,
             data = troutegg)

summary(bmod2)
## 
## Call:
## glm(formula = cbind(survive, total - survive) ~ location + period, 
##     family = quasibinomial, data = troutegg)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.8305  -0.3650  -0.0303   0.6191   3.2434  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.6358     0.6495   7.138 1.18e-05 ***
## location2    -0.4168     0.5682  -0.734 0.477315    
## location3    -1.2421     0.5066  -2.452 0.030501 *  
## location4    -0.9509     0.5281  -1.800 0.096970 .  
## location5    -4.6138     0.5777  -7.987 3.82e-06 ***
## period7      -2.1702     0.5504  -3.943 0.001953 ** 
## period8      -2.3256     0.5609  -4.146 0.001356 ** 
## period11     -2.4500     0.5405  -4.533 0.000686 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 5.330358)
## 
##     Null deviance: 1021.469  on 19  degrees of freedom
## Residual deviance:   64.495  on 12  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
plot.glm.diag(bmod2) ## see SM Section 10.2.3

The summary tables from either model are not that helpful. We can get an idea of the relative importance of location and period using analysis of the deviance. I am following the text (2nd edition p.59).

phihat <- sum(residuals(bmod, type="pearson")^2/12) #12 is the residual degrees of freedom
drop1(bmod, scale=phihat, test = "F")  #F-test
## Warning in drop1.glm(bmod, scale = phihat, test = "F"): F test assumes
## 'quasibinomial' family
## Single term deletions
## 
## Model:
## cbind(survive, total - survive) ~ location + period
## 
## scale:  5.330322 
## 
##          Df Deviance    AIC F value    Pr(>F)    
## <none>         64.50 157.03                      
## location  4   913.56 308.32  39.494 8.142e-07 ***
## period    3   228.57 181.81  10.176  0.001288 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally the plots of the can be helpful to see what’s going on.

elogits <- with(troutegg, log((survive+0.5)/(total-survive+0.5)))
with(troutegg, interaction.plot(period,location,elogits))

Location 5 appears to have lower survival rates than the other 4 locations, and survival rates decrease somewhat over time.