> #------------------------------------------------------------------- > # 7. LOGISTIC REGRESSION > #------------------------------------------------------------------- > > #------------------------------------------------------------------- > # KYPHOSIS > #------------------------------------------------------------------- > m<-matrix(scan("/glim/datasets/kyphosis"),81,5,T) Read 405 items > y<-m[,2] > age<-m[,3] > cbind(y,age) y age [1,] 0 71 [2,] 0 158 [3,] 1 128 [4,] 0 2 [78,] 0 26 [79,] 0 120 [80,] 1 42 [81,] 0 36 > par(mfrow=c(2,2)) > plot(age,y,yaxs="i",ylim=c(0,1),main="Figure 7.1 Kyphosis") > summary(glm(y~1,family=binomial)) Call: glm(formula = y ~ 1, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -0.6864 -0.6864 -0.6864 -0.6864 1.7671 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.3257 0.2727 -4.862 1.16e-06 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 83.234 on 80 degrees of freedom Residual deviance: 83.234 on 80 degrees of freedom AIC: 85.234 Number of Fisher Scoring iterations: 3 > glm1<-glm(y~age,family=binomial) > summary(glm1) Call: glm(formula = y ~ age, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -0.9023 -0.7397 -0.6028 -0.5521 1.9449 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.809352 0.530332 -3.412 0.000646 *** age 0.005442 0.004820 1.129 0.258896 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 83.234 on 80 degrees of freedom Residual deviance: 81.932 on 79 degrees of freedom AIC: 85.932 Number of Fisher Scoring iterations: 3 > i<-order(age) > i [1] 5 6 14 16 37 4 54 57 29 27 26 52 75 25 70 66 13 69 39 21 78 20 31 42 81 [26] 8 80 59 38 10 7 44 51 1 56 41 17 19 63 11 73 40 36 23 45 28 60 22 34 9 [51] 62 64 65 58 79 49 32 72 3 33 61 24 48 47 46 53 35 55 43 12 30 77 2 71 68 [76] 15 18 50 76 67 74 > lines(age[i],fitted(glm1)[i]) > a<-seq(-500,1000,100) > eta<- -1.809076886 + 0.005439848 *a > phat<-exp(eta)/(1+exp(eta)) > plot(a,phat,yaxs="i",ylim=c(0,1),type="l",main="Figure 7.2 Kyphosis enlarged") > points(age,y) > plot(age,y,yaxs="i",ylim=c(0,1),main="Figure 7.3 Quadratic in age") > age2<-age^2 > glm2<-glm(y~age+age2,family=binomial) > summary(glm2) Call: glm(formula = y ~ age + age2, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -1.0079 -0.8412 -0.4155 -0.2210 2.3919 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.7699500 1.1398759 -3.307 0.000942 *** age 0.0700279 0.0267689 2.616 0.008896 ** age2 -0.0003651 0.0001468 -2.487 0.012865 * --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 83.234 on 80 degrees of freedom Residual deviance: 72.739 on 78 degrees of freedom AIC: 78.739 Number of Fisher Scoring iterations: 4 > lines(age[i],fitted(glm2)[i]) > > #------------------------------------------------------------------- > # 245-T > #------------------------------------------------------------------- > misc<-c(10,17,18,11,16,24,20,17,9,15,16,15) > totf<-c(222,221,188,183,197,218,208,219,198,216,244,218) > herb<-c(0,0,0,0,1454,3280,788,0,304,560,0,0) > p<-misc/totf > cbind(month.name,herb,misc,totf,p) month.name herb misc totf p [1,] "January" "0" "10" "222" "0.045045045045045" [2,] "February" "0" "17" "221" "0.076923076923077" [3,] "March" "0" "18" "188" "0.0957446808510638" [4,] "April" "0" "11" "183" "0.0601092896174863" [5,] "May" "1454" "16" "197" "0.0812182741116751" [6,] "June" "3280" "24" "218" "0.110091743119266" [7,] "July" "788" "20" "208" "0.0961538461538462" [8,] "August" "0" "17" "219" "0.0776255707762557" [9,] "September" "304" "9" "198" "0.0454545454545455" [10,] "October" "560" "15" "216" "0.0694444444444444" [11,] "November" "0" "16" "244" "0.0655737704918033" [12,] "December" "0" "15" "218" "0.0688073394495413" > plot(herb,p,main="Figure 7.4 245-T") > summary(glm(p~1,family=binomial,weight=totf)) Call: glm(formula = p ~ 1, family = binomial, weights = totf) Deviance Residuals: Min 1Q Median 3Q Max -1.78153 -0.58321 -0.06067 0.54583 1.89358 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.5232 0.0758 -33.29 <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: 13.216 on 11 degrees of freedom Residual deviance: 13.216 on 11 degrees of freedom AIC: 1341.4 Number of Fisher Scoring iterations: 5 > glm<-glm(p~herb,family=binomial,weight=totf) > summary(glm) Call: glm(formula = p ~ herb, family = binomial, weights = totf) Deviance Residuals: Min 1Q Median 3Q Max -1.49318 -0.29192 -0.07846 0.53689 1.43971 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.621e+00 8.943e-02 -29.303 <2e-16 *** herb 1.629e-04 7.074e-05 2.303 0.0213 * --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 13.216 on 11 degrees of freedom Residual deviance: 8.310 on 10 degrees of freedom AIC: 1338.5 Number of Fisher Scoring iterations: 5 > cbind(p,fitted(glm),resid(glm)) p 1 0.04504505 0.06782132 -1.43211122 2 0.07692308 0.06782132 0.52756510 3 0.09574468 0.06782132 1.43971084 4 0.06010929 0.06782132 -0.42258532 5 0.08121827 0.08442103 -0.16263388 6 0.11009174 0.11044338 -0.01657168 7 0.09615385 0.07640292 1.03404181 8 0.07762557 0.06782132 0.56488101 9 0.04545455 0.07102066 -1.49318090 10 0.06944444 0.07382237 -0.24836440 11 0.06557377 0.06782132 -0.14035281 12 0.06880734 0.06782132 0.05777124 > i<-order(herb) > lines(herb[i],fitted(glm)[i]) > h<-seq(-5000,35000,1000) > eta<- -2.6206479509+ 0.0001629353 *h > phat<-exp(eta)/(1+exp(eta)) > plot(h,phat,yaxs="i",ylim=c(0,1),type="l",main="Figure 7.5 245-T enlarged") > points(herb,p) > > #------------------------------------------------------------------- > # CLINICS (AGAIN) > #------------------------------------------------------------------- > died<-c(3,4,17,2) > surv<-c(176,293,197,23) > totc<-died+surv > p<-died/totc > clin<-factor(c(1,1,2,2)) > care<-factor(c(1,2,1,2)) > cbind(p,clin,care) p clin care [1,] 0.01675978 1 1 [2,] 0.01346801 1 2 [3,] 0.07943925 2 1 [4,] 0.08000000 2 2 > summary(glm(p~1,family=binomial,weight=totc)) Call: glm(formula = p ~ 1, family = binomial, weights = totc) Deviance Residuals: 1 2 3 4 -1.563 -2.411 2.924 1.011 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.2771 0.1996 -16.41 <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: 17.828 on 3 degrees of freedom Residual deviance: 17.828 on 3 degrees of freedom AIC: 225.38 Number of Fisher Scoring iterations: 5 > summary(glm(p~care,family=binomial,weight=totc)) Call: glm(formula = p ~ care, family = binomial, weights = totc) Deviance Residuals: 1 2 3 4 -2.4024 -0.6923 1.7628 1.6905 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.9258 0.2295 -12.747 <2e-16 *** care2 -1.0381 0.4717 -2.201 0.0277 * --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 17.828 on 3 degrees of freedom Residual deviance: 12.216 on 2 degrees of freedom AIC: 221.77 Number of Fisher Scoring iterations: 6 > summary(glm(p~clin,family=binomial,weight=totc)) Call: glm(formula = p ~ clin, family = binomial, weights = totc) Deviance Residuals: 1 2 3 4 0.223333 -0.179765 -0.003172 0.009271 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.2047 0.3806 -11.048 < 2e-16 *** clin2 1.7555 0.4495 3.906 9.4e-05 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 17.82840 on 3 degrees of freedom Residual deviance: 0.08229 on 2 degrees of freedom AIC: 209.63 Number of Fisher Scoring iterations: 6 > summary(glm(p~care+clin,family=binomial,weight=totc)) Call: glm(formula = p ~ care + clin, family = binomial, weights = totc) Deviance Residuals: 1 2 3 4 0.11107 -0.09263 -0.04706 0.14186 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.1372 0.5075 -8.152 3.57e-16 *** care2 -0.1104 0.5609 -0.197 0.84398 clin2 1.6991 0.5305 3.203 0.00136 ** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 17.828399 on 3 degrees of freedom Residual deviance: 0.043256 on 1 degrees of freedom AIC: 211.60 Number of Fisher Scoring iterations: 6 > glm<-glm(p~clin,family=binomial,weight=totc) > cbind(p*totc,fitted(glm)*totc,resid(glm)) [,1] [,2] [,3] 1 3 2.632354 0.223332781 2 4 4.367649 -0.179765476 3 17 17.012552 -0.003172306 4 2 1.987448 0.009271422 > > #------------------------------------------------------------------- > # 245-T (AGAIN) > #------------------------------------------------------------------- > freq<-c(misc,totf-misc) > herbb<-c(herb,herb) > month<-gl(12,1,24) > surv<-gl(2,12,24) > cbind(month,surv,freq,herbb) month surv freq herbb [1,] 1 1 10 0 [2,] 2 1 17 0 [3,] 3 1 18 0 [4,] 4 1 11 0 [5,] 5 1 16 1454 [6,] 6 1 24 3280 [7,] 7 1 20 788 [8,] 8 1 17 0 [9,] 9 1 9 304 [10,] 10 1 15 560 [11,] 11 1 16 0 [12,] 12 1 15 0 [13,] 1 2 212 0 [14,] 2 2 204 0 [15,] 3 2 170 0 [16,] 4 2 172 0 [17,] 5 2 181 1454 [18,] 6 2 194 3280 [19,] 7 2 188 788 [20,] 8 2 202 0 [21,] 9 2 189 304 [22,] 10 2 201 560 [23,] 11 2 228 0 [24,] 12 2 203 0 > summary(glm(freq~month+surv+surv:herbb,family=poisson)) Call: glm(formula = freq ~ month + surv + surv:herbb, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -1.446215 -0.176585 -0.004789 0.168447 1.383445 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.712e+00 1.070e-01 25.340 <2e-16 *** month2 -4.515e-03 9.502e-02 -0.048 0.9621 month3 -1.662e-01 9.911e-02 -1.677 0.0935 . month4 -1.932e-01 9.984e-02 -1.935 0.0530 . month5 -1.374e-01 9.821e-02 -1.400 0.1617 month6 -6.498e-02 9.828e-02 -0.661 0.5085 month7 -7.439e-02 9.658e-02 -0.770 0.4412 month8 -1.361e-02 9.524e-02 -0.143 0.8864 month9 -1.178e-01 9.776e-02 -1.206 0.2280 month10 -3.386e-02 9.561e-02 -0.354 0.7232 month11 9.449e-02 9.275e-02 1.019 0.3083 month12 -1.818e-02 9.535e-02 -0.191 0.8488 surv2 2.621e+00 8.942e-02 29.307 <2e-16 *** surv1.herbb 1.629e-04 7.074e-05 2.303 0.0213 * --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 2199.03 on 23 degrees of freedom Residual deviance: 8.31 on 10 degrees of freedom AIC: 176.41 Number of Fisher Scoring iterations: 3 > > #------------------------------------------------------------------- > # BIRDS > #------------------------------------------------------------------- > cal<-c( 36.9, 41.1, 35.8, 40.6, 34.6, 38.9, 34.3, 37.9, + 32.8, 37.0, 31.7, 36.1, 31.0, 36.3, 29.8, 34.2, + 29.1, 33.4, 28.2, 32.8, 27.4, 32.0, 27.8, 31.9, + 25.5, 30.7, 24.9, 29.5, 23.7, 28.5, 23.1, 27.7) > spec<-gl(2,1,32) > temp<-rep((0:15)*2,rep(2,16)) > plot(temp,cal,type="n",main="Figure 7.6 Birds") > points(temp[spec==1],cal[spec==1],pch="A") > points(temp[spec==2],cal[spec==2],pch="B") > glm<-glm(cal~spec-1+temp:spec) > summary(glm) Call: glm(formula = cal ~ spec - 1 + temp:spec) Deviance Residuals: Min 1Q Median 3Q Max -0.5243 -0.2305 -0.1347 0.1599 1.1821 Coefficients: Estimate Std. Error t value Pr(>|t|) spec1 36.57941 0.19739 185.32 <2e-16 *** spec2 40.83897 0.19739 206.90 <2e-16 *** spec1.temp -0.45279 0.01121 -40.39 <2e-16 *** spec2.temp -0.43676 0.01121 -38.96 <2e-16 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 0.1709338) Null deviance: 33549.9000 on 32 degrees of freedom Residual deviance: 4.7861 on 28 degrees of freedom AIC: 40.012 Number of Fisher Scoring iterations: 2 > abline(coef(glm)[c(1,3)]) > abline(coef(glm)[c(2,4)]) > glm<-glm(cal~temp*spec) > summary(glm) Call: glm(formula = cal ~ temp * spec) Deviance Residuals: Min 1Q Median 3Q Max -0.5243 -0.2305 -0.1347 0.1599 1.1821 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 36.57941 0.19739 185.315 < 2e-16 *** temp -0.45279 0.01121 -40.388 < 2e-16 *** spec2 4.25956 0.27915 15.259 4.25e-15 *** temp.spec2 0.01603 0.01585 1.011 0.321 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 0.1709338) Null deviance: 705.0550 on 31 degrees of freedom Residual deviance: 4.7861 on 28 degrees of freedom AIC: 40.012 Number of Fisher Scoring iterations: 2 > plot(temp,resid(glm),type="n",main="Figure 7.7 Birds residuals") > points(temp[spec==1],resid(glm)[spec==1],pch="A") > points(temp[spec==2],resid(glm)[spec==2],pch="B") > Temp<-factor(temp) > glm<-glm(cal~Temp+temp*spec) > summary(glm) Call: glm(formula = cal ~ Temp + temp * spec) Deviance Residuals: Min 1Q Median 3Q Max -4.240e-01 -5.191e-02 1.776e-15 5.191e-02 4.240e-01 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 36.87022 0.22215 165.972 < 2e-16 *** Temp2 -0.81603 0.28372 -2.876 0.0122 * Temp4 -2.28206 0.28434 -8.026 1.32e-06 *** Temp6 -2.94809 0.28538 -10.330 6.24e-08 *** Temp8 -4.16412 0.28682 -14.518 7.83e-10 *** Temp10 -5.18015 0.28867 -17.945 4.65e-11 *** Temp12 -5.44618 0.29092 -18.721 2.63e-11 *** Temp14 -7.11221 0.29355 -24.229 7.87e-13 *** Temp16 -7.87824 0.29655 -26.566 2.23e-13 *** Temp18 -8.64426 0.29992 -28.822 7.25e-14 *** Temp20 -9.46029 0.30364 -31.156 2.47e-14 *** Temp22 -9.32632 0.30770 -30.310 3.62e-14 *** Temp24 -11.09235 0.31209 -35.543 4.00e-15 *** Temp26 -12.00838 0.31679 -37.907 1.64e-15 *** Temp28 -13.12441 0.32178 -40.786 5.93e-16 *** Temp30 -13.84044 0.32707 -42.317 3.55e-16 *** spec2 4.25956 0.19142 22.252 2.52e-12 *** temp.spec2 0.01603 0.01087 1.474 0.1625 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 0.0803771) Null deviance: 705.0550 on 31 degrees of freedom Residual deviance: 1.1253 on 14 degrees of freedom AIC: 21.686 Number of Fisher Scoring iterations: 2 > diff<-cal[spec==2]-cal[spec==1] > tem<-(0:15)*2 > plot(tem,diff,main="Figure 7.8 Birds differences") > summary(glm(diff~tem)) Call: glm(formula = diff ~ tem) Deviance Residuals: Min 1Q Median 3Q Max -0.75574 -0.15228 -0.04162 0.06184 0.84809 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.25956 0.19142 22.252 2.52e-12 *** tem 0.01603 0.01087 1.474 0.163 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 0.1607542) Null deviance: 2.6000 on 15 degrees of freedom Residual deviance: 2.2506 on 14 degrees of freedom AIC: 20.023 Number of Fisher Scoring iterations: 2 > > #------------------------------------------------------------------- > # ESKIMO > #------------------------------------------------------------------- > m<-matrix(scan("/glim/datasets/eskimo"),ncol=2,byrow=T) Read 72 items > pres<-m[,1] > abs<-m[,2] > tot<-pres+abs > p<-pres/tot > age<-rep(1:6,6) > plot(age,p,yaxs="i",ylim=c(0,1),main="Figure 7.9 Eskimo") > glm<-glm(p~age,family=binomial,weight=tot) > summary(glm) Call: glm(formula = p ~ age, family = binomial, weights = tot) Deviance Residuals: Min 1Q Median 3Q Max -2.4374 -0.9261 -0.3488 0.9059 2.3244 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.71593 0.24203 -11.22 <2e-16 *** age 0.77734 0.07518 10.34 <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: 203.568 on 35 degrees of freedom Residual deviance: 59.187 on 34 degrees of freedom AIC: 579.51 Number of Fisher Scoring iterations: 4 > lines(age[1:6],fitted(glm)[1:6]) > a2<-age^2 > glm<-glm(p~age+a2,family=binomial,weight=tot) > summary(glm) Call: glm(formula = p ~ age + a2, family = binomial, weights = tot) Deviance Residuals: Min 1Q Median 3Q Max -2.43900 -0.94281 0.02440 0.86376 2.52769 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.49915 0.48987 -7.143 9.13e-13 *** age 1.36435 0.31721 4.301 1.70e-05 *** a2 -0.08972 0.04612 -1.945 0.0517 . --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 203.57 on 35 degrees of freedom Residual deviance: 55.45 on 33 degrees of freedom AIC: 577.77 Number of Fisher Scoring iterations: 4 > Age<-factor(age) > glm<-glm(p~age+a2+Age,family=binomial,weight=tot) > summary(glm) Call: glm(formula = p ~ age + a2 + Age, family = binomial, weights = tot) Deviance Residuals: Min 1Q Median 3Q Max -2.24696 -0.81913 -0.02167 0.87265 2.68329 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.02919 0.75412 -4.017 5.9e-05 *** age 0.88428 0.76241 1.160 0.246 a2 -0.02224 0.11498 -0.193 0.847 Age2 0.19690 0.48593 0.405 0.685 Age3 0.25226 0.62426 0.404 0.686 Age4 0.66887 0.60791 1.100 0.271 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 203.568 on 35 degrees of freedom Residual deviance: 53.128 on 30 degrees of freedom AIC: 581.45 Number of Fisher Scoring iterations: 4 > sex<-gl(2,6,36) > pop<-gl(3,12,36) > glm(p~age+sex,family=binomial,weight=tot) Call: glm(formula = p ~ age + sex, family = binomial, weights = tot) Coefficients: (Intercept) age sex2 -2.6425 0.7774 -0.1588 Degrees of Freedom: 35 Total (i.e. Null); 33 Residual Null Deviance: 203.6 Residual Deviance: 58.59 AIC: 580.9 > glm(p~age+pop,family=binomial,weight=tot) Call: glm(formula = p ~ age + pop, family = binomial, weights = tot) Coefficients: (Intercept) age pop2 pop3 -2.65203 0.80128 -0.07083 -0.60304 Degrees of Freedom: 35 Total (i.e. Null); 32 Residual Null Deviance: 203.6 Residual Deviance: 54.18 AIC: 578.5 > glm(p~age+sex+pop,family=binomial,weight=tot) Call: glm(formula = p ~ age + sex + pop, family = binomial, weights = tot) Coefficients: (Intercept) age sex2 pop2 pop3 -2.57853 0.80118 -0.15784 -0.07261 -0.60339 Degrees of Freedom: 35 Total (i.e. Null); 31 Residual Null Deviance: 203.6 Residual Deviance: 53.59 AIC: 579.9 > glm(p~age+pop+age:pop,family=binomial,weight=tot) Call: glm(formula = p ~ age + pop + age:pop, family = binomial, weights = tot) Coefficients: (Intercept) age pop2 pop3 age.pop2 age.pop3 -3.22631 1.00694 0.21001 1.47078 -0.09977 -0.64809 Degrees of Freedom: 35 Total (i.e. Null); 30 Residual Null Deviance: 203.6 Residual Deviance: 40.46 AIC: 568.8 > glm(p~age+sex+pop+age:pop,family=binomial,weight=tot) Call: glm(formula = p ~ age + sex + pop + age:pop, family = binomial, weights = tot) Coefficients: (Intercept) age sex2 pop2 pop3 age.pop2 -3.1467 1.0101 -0.1909 0.2225 1.4934 -0.1050 age.pop3 -0.6557 Degrees of Freedom: 35 Total (i.e. Null); 29 Residual Null Deviance: 203.6 Residual Deviance: 39.63 AIC: 570 > glm(p~age+sex+pop+age:sex+age:pop,family=binomial,weight=tot) Call: glm(formula = p ~ age + sex + pop + age:sex + age:pop, family = binomial, weights = tot) Coefficients: (Intercept) age sex2 pop2 pop3 age.sex2 -3.155945 1.013427 -0.170388 0.223186 1.493409 -0.007047 age.pop2 age.pop3 -0.105333 -0.655869 Degrees of Freedom: 35 Total (i.e. Null); 28 Residual Null Deviance: 203.6 Residual Deviance: 39.63 AIC: 572 > glm(p~age+sex+pop+age:sex+age:pop+sex:pop,family=binomial,weight=tot) Call: glm(formula = p ~ age + sex + pop + age:sex + age:pop + sex:pop, family = binomial, weights = tot) Coefficients: (Intercept) age sex2 pop2 pop3 age.sex2 -3.23214 0.99725 0.01333 0.41556 1.80258 0.02030 age.pop2 age.pop3 sex2.pop2 sex2.pop3 -0.09919 -0.65592 -0.46004 -0.68089 Degrees of Freedom: 35 Total (i.e. Null); 26 Residual Null Deviance: 203.6 Residual Deviance: 37.65 AIC: 574 > glm(p~age+sex+pop+age:sex+age:pop+sex:pop+age:sex:pop,family=binomial,weight=tot) Call: glm(formula = p ~ age + sex + pop + age:sex + age:pop + sex:pop + age:sex:pop, family = binomial, weights = tot) Coefficients: (Intercept) age sex2 pop2 pop3 -3.7976 1.2001 1.0861 1.5032 2.6494 age.sex2 age.pop2 age.pop3 sex2.pop2 sex2.pop3 -0.3641 -0.4932 -0.9442 -3.1738 -2.4610 age.sex2.pop2 age.sex2.pop3 0.9716 0.5977 Degrees of Freedom: 35 Total (i.e. Null); 24 Residual Null Deviance: 203.6 Residual Deviance: 32.14 AIC: 572.5 > glm<-glm(p~age+pop+age:pop,family=binomial,weight=tot) > summary(glm) Call: glm(formula = p ~ age + pop + age:pop, family = binomial, weights = tot) Deviance Residuals: Min 1Q Median 3Q Max -2.0107 -0.8015 -0.1724 0.9120 1.9330 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.22631 0.35695 -9.039 < 2e-16 *** age 1.00694 0.11747 8.572 < 2e-16 *** pop2 0.21001 0.65185 0.322 0.747318 pop3 1.47078 0.59911 2.455 0.014090 * age.pop2 -0.09977 0.21429 -0.466 0.641510 age.pop3 -0.64809 0.17621 -3.678 0.000235 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 203.568 on 35 degrees of freedom Residual deviance: 40.463 on 30 degrees of freedom AIC: 568.79 Number of Fisher Scoring iterations: 4 > glm<-glm(p~pop-1+age:pop,family=binomial,weight=tot) > summary(glm) Call: glm(formula = p ~ pop - 1 + age:pop, family = binomial, weights = tot) Deviance Residuals: Min 1Q Median 3Q Max -2.0107 -0.8015 -0.1724 0.9120 1.9330 Coefficients: Estimate Std. Error z value Pr(>|z|) pop1 -3.2263 0.3570 -9.039 < 2e-16 *** pop2 -3.0163 0.5454 -5.530 3.20e-08 *** pop3 -1.7555 0.4812 -3.649 0.000264 *** pop1.age 1.0069 0.1175 8.572 < 2e-16 *** pop2.age 0.9072 0.1792 5.062 4.16e-07 *** pop3.age 0.3588 0.1313 2.732 0.006291 ** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 233.661 on 36 degrees of freedom Residual deviance: 40.463 on 30 degrees of freedom AIC: 568.79 Number of Fisher Scoring iterations: 4 > plot(age,p,yaxs="i",ylim=c(0,1),type="n",main="Figure 7.10 Eskimo") > points(age[pop==1],p[pop==1],pch="1") > points(age[pop==2],p[pop==2],pch="2") > points(age[pop==3],p[pop==3],pch="3") > lines(age[1:6],fitted(glm)[1:6],lty=2) > lines(age[13:18],fitted(glm)[13:18],lty=3) > lines(age[25:30],fitted(glm)[25:30],lty=4) >