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.