> library(SMPracticals) > data(jacamar) > jacamar[1:5,] species colour N S E 1 Ab Unpainted 0 0 14 2 Pa Unpainted 6 1 0 3 Di Unpainted 1 0 2 4 Pl Unpainted 4 1 5 5 Cf Unpainted 0 0 0 > dim(jacamar) [1] 48 5 > fit1 = glm(cbind(E, N+S) ~ colour + species, family = binomial, data = jacamar) # response is a two-column vector of 'successes' and 'failures'; see Details section of ?glm > summary(fit1) Call: glm(formula = cbind(E, N + S) ~ colour + species, family = binomial, data = jacamar) Deviance Residuals: Min 1Q Median 3Q Max -2.3912 -0.6250 0.0000 0.7751 3.3933 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.6759 0.4590 1.473 0.140837 colourBrown -1.5973 0.6139 -2.602 0.009267 ** colourYellow -1.3930 0.5823 -2.392 0.016740 * colourBlue -2.1680 0.7141 -3.036 0.002398 ** colourGreen -2.3071 0.6837 -3.374 0.000740 *** colourRed -3.2633 0.7896 -4.133 3.58e-05 *** colourOrange -3.1498 0.7260 -4.339 1.43e-05 *** colourBlack -4.4381 0.9500 -4.672 2.99e-06 *** speciesPa -2.0666 0.8504 -2.430 0.015094 * speciesDi 0.1284 0.5961 0.215 0.829390 speciesPl 0.9094 0.4619 1.969 0.048975 * speciesCf 2.8103 0.8046 3.493 0.000478 *** speciesSs 2.3337 0.6776 3.444 0.000573 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 134.235 on 43 degrees of freedom Residual deviance: 67.282 on 31 degrees of freedom AIC: 133.19 Number of Fisher Scoring iterations: 5 > anova(fit1) # reproduces part of Table 10.3 Analysis of Deviance Table Model: binomial, link: logit Response: cbind(E, N + S) Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 43 134.235 colour 7 25.775 36 108.460 species 5 41.178 31 67.282 > fit2 = glm(cbind(E, N+S) ~ species + colour, family = binomial, data = jacamar) # the other order > anova(fit2) # gives us the rest of Table 10.3 (with full data) Analysis of Deviance Table Model: binomial, link: logit Response: cbind(E, N + S) Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 43 134.235 species 5 19.641 38 114.595 colour 7 47.313 31 67.282 > plot.glm.diag(fit1) # a program special to SMPracticals > plot(fit2$fit, glm.diag(fit2)$rd) # > identify(fit2$fit, glm.diag(fit2)$rd) # points 1 and 47 > fit3 = glm(cbind(E, N+S) ~ species + colour, family = binomial, data = jacamar[-c(3,47),]) # this doesn't work; check the output for speciesCf > summary(fit3) > which(jacamar$species == "Cf") # as per note on p.485 [1] 5 11 17 23 29 35 41 47 > jacamar.small = jacamar[-c(1,5,11,17,23,29,35,41,47),] > dim(jacamar.small) [1] 39 5 > fit3 = glm(cbind(E,N+S) ~ colour +species, family = binomial, data = jacamar.small) > anova(fit3) # now we get the rest of Table 10.3 Df Deviance Resid. Df Resid. Dev NULL 35 73.681 colour 7 10.477 28 63.205 species 4 35.183 24 28.022 > fit4 = glm(cbind(E,N+S) ~ species +colour, family = binomial, data = jacamar.small) > anova(fit4) # and the other order Df Deviance Resid. Df Resid. Dev NULL 35 73.681 species 4 27.635 31 46.047 colour 7 18.025 24 28.022 > pchisq(18.025, 7) # is colour significant? [1] 0.9881423 > 1-.Last.value [1] 0.01185772 > summary(fit4) > fit5 = glm(cbind(E,N+S) ~ species +colour - 1, family = binomial, data = jacamar.small) # notice the -1; output now agrees with Table 10.4 > summary(fit5) Coefficients: Estimate Std. Error z value Pr(>|z|) speciesAb -1.9894 0.7910 -2.515 0.01190 * speciesPa -2.2187 0.8455 -2.624 0.00869 ** speciesDi -0.5597 0.6697 -0.836 0.40333 speciesPl 0.1622 0.5358 0.303 0.76204 speciesSs 1.5019 0.7763 1.935 0.05303 . colourBrown 0.1588 0.7306 0.217 0.82792 colourYellow 0.3347 0.6824 0.490 0.62380 colourBlue -0.5349 0.8084 -0.662 0.50814 colourGreen -0.8330 0.7453 -1.118 0.26369 colourRed -1.9257 0.8827 -2.182 0.02914 * colourOrange -1.9385 0.8475 -2.287 0.02218 * colourBlack -1.2552 0.8609 -1.458 0.14484 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 113.681 on 36 degrees of freedom Residual deviance: 28.022 on 24 degrees of freedom AIC: 91.928 Number of Fisher Scoring iterations: 5 > residuals(fit4, type="pearson") # these are the unstandardized Pearson residuals; used for computing Pearson's statistic 2 3 4 6 7 -0.872479346 1.091144192 -0.256805318 0.471918599 0.566642011 8 9 10 12 13 -0.618373942 0.285355994 -0.455980801 0.435895791 -0.521000637 14 15 16 18 19 1.457857007 -1.367207694 0.504915054 0.691434568 -0.693297640 20 21 22 24 25 0.000000000 1.728592971 0.111353430 -0.710689786 1.684695506 > sum(.Last.value^2) [1] 25.58358 # agrees with text p. 485