R : Copyright 1999, The R Development Core Team Version 0.90.1 (December 15, 1999) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type "?license" or "?licence" for distribution details. R is a collaborative project with many contributors. Type "?contributors" for a list. Type "demo()" for some demos, "help()" for on-line help, or "help.start()" for a HTML browser interface to help. Type "q()" to quit R. > > #------------------------------------------------------------------- > # 1.1 SIMPLE REGRESSION > #------------------------------------------------------------------- > x<-c(6,6.3,6.5,6.8,7,7.1,7.5,7.5,7.6) > y<-c(39,58,49,53,80,86,115,124,104) > par(mfrow=c(2,2)) > plot(x,y,main="Figure 1.1 Simple regression") > glm.out<-glm(y~x) > glm.out Call: glm(formula = y ~ x) Coefficients: (Intercept) x -270.35 50.42 Degrees of Freedom: 8 Total (i.e. Null); 7 Residual Null Deviance: 7612 Residual Deviance: 1014 AIC: 74.06 > summary(glm.out) Call: glm(formula = y ~ x) Deviance Residuals: Min 1Q Median 3Q Max -19.504 -8.378 -1.630 7.202 16.202 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -270.35 51.86 -5.213 0.001236 ** x 50.42 7.47 6.750 0.000265 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 144.8236) Null deviance: 7612.0 on 8 degrees of freedom Residual deviance: 1013.8 on 7 degrees of freedom AIC: 74.059 Number of Fisher Scoring iterations: 2 > lines(x,fitted(glm.out)) > glm(y~1) Call: glm(formula = y ~ 1) Coefficients: (Intercept) 78.67 Degrees of Freedom: 8 Total (i.e. Null); 8 Residual Null Deviance: 7612 Residual Deviance: 7612 AIC: 90.2 > r<-(7612-1013.765)/7612 > r [1] 0.8668202 > tstat<-50.41952/7.469724 > tstat [1] 6.74985 > tstat^2/(7+tstat^2) [1] 0.8668201 > > > #------------------------------------------------------------------- > # 1.2 MAYER'S ESTIMATE OF THE LIBRATION OF THE MOON > #------------------------------------------------------------------- > m<-matrix(scan("/glim/datasets/moon.dat"),byrow=T,ncol=5) Read 135 items > m [,1] [,2] [,3] [,4] [,5] [1,] 13 10 0.8836 -0.4682 1 [2,] 13 8 0.9996 -0.0282 1 [3,] 13 12 0.9899 0.1421 1 [4,] 14 15 0.2221 0.9750 3 [5,] 14 42 0.0006 1.0000 3 [6,] 13 1 0.9308 -0.3654 1 [7,] 14 31 0.0602 0.9982 3 [8,] 14 57 -0.1570 0.9876 2 [9,] 13 5 0.9097 -0.4152 1 [10,] 13 2 1.0000 0.0055 1 [11,] 13 12 0.9689 0.2476 1 [12,] 13 11 0.8878 0.4602 1 [13,] 13 34 0.7549 0.6558 3 [14,] 13 53 0.5755 0.8178 3 [15,] 13 58 0.3608 0.9326 3 [16,] 14 14 0.1302 0.9915 3 [17,] 14 56 -0.1068 0.9943 3 [18,] 14 47 -0.3363 0.9418 2 [19,] 15 56 -0.8560 0.5170 2 [20,] 13 29 0.8002 0.5997 3 [21,] 15 55 -0.9952 -0.0982 2 [22,] 15 39 -0.8409 0.5412 2 [23,] 16 9 -0.9429 0.3330 2 [24,] 16 22 -0.9768 0.2141 2 [25,] 15 38 -0.6262 -0.7797 2 [26,] 14 54 -0.4091 -0.9125 2 [27,] 13 7 0.9284 -0.3716 1 > y<- -m[,1]-m[,2]/60 > y [1] -13.16667 -13.13333 -13.20000 -14.25000 -14.70000 [6] -13.01667 -14.51667 -14.95000 -13.08333 -13.03333 [11] -13.20000 -13.18333 -13.56667 -13.88333 -13.96667 [16] -14.23333 -14.93333 -14.78333 -15.93333 -13.48333 [21] -15.91667 -15.65000 -16.15000 -16.36667 -15.63333 [26] -14.90000 -13.11667 > alpha<-m[,3] > asin<-m[,4] > grp<-m[,5] > glm.out<-glm(y~alpha+asin) > glm.out Call: glm(formula = y ~ alpha + asin) Coefficients: (Intercept) alpha asin -14.55825 1.50580 -0.07192 Degrees of Freedom: 26 Total (i.e. Null); 24 Residual Null Deviance: 32.14 Residual Deviance: 0.5689 AIC: -19.6 > summary(glm.out) Call: glm(formula = y ~ alpha + asin) Deviance Residuals: Min 1Q Median 3Q Max -0.32216 -0.09244 0.01687 0.09447 0.34905 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -14.55825 0.03540 -411.269 <2e-16 alpha 1.50580 0.04173 36.082 <2e-16 asin -0.07192 0.05083 -1.415 0.17 (Intercept) *** alpha *** asin --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 0.02370312) Null deviance: 32.13944 on 26 degrees of freedom Residual deviance: 0.56887 on 24 degrees of freedom AIC: -19.595 Number of Fisher Scoring iterations: 2 > x<-cbind(rep(1,27),alpha,asin) > x alpha asin [1,] 1 0.8836 -0.4682 [2,] 1 0.9996 -0.0282 [3,] 1 0.9899 0.1421 [4,] 1 0.2221 0.9750 [5,] 1 0.0006 1.0000 [6,] 1 0.9308 -0.3654 [7,] 1 0.0602 0.9982 [8,] 1 -0.1570 0.9876 [9,] 1 0.9097 -0.4152 [10,] 1 1.0000 0.0055 [11,] 1 0.9689 0.2476 [12,] 1 0.8878 0.4602 [13,] 1 0.7549 0.6558 [14,] 1 0.5755 0.8178 [15,] 1 0.3608 0.9326 [16,] 1 0.1302 0.9915 [17,] 1 -0.1068 0.9943 [18,] 1 -0.3363 0.9418 [19,] 1 -0.8560 0.5170 [20,] 1 0.8002 0.5997 [21,] 1 -0.9952 -0.0982 [22,] 1 -0.8409 0.5412 [23,] 1 -0.9429 0.3330 [24,] 1 -0.9768 0.2141 [25,] 1 -0.6262 -0.7797 [26,] 1 -0.4091 -0.9125 [27,] 1 0.9284 -0.3716 > beta.ls<-solve(t(x) %*% x, t(x) %*% y) > beta.ls [,1] [1,] -14.55824555 [2,] 1.50579501 [3,] -0.07192118 > g<- cbind(grp==1,grp==2,grp==3)*1 > g [,1] [,2] [,3] [1,] 1 0 0 [2,] 1 0 0 [3,] 1 0 0 [4,] 0 0 1 [5,] 0 0 1 [6,] 1 0 0 [7,] 0 0 1 [8,] 0 1 0 [9,] 1 0 0 [10,] 1 0 0 [11,] 1 0 0 [12,] 1 0 0 [13,] 0 0 1 [14,] 0 0 1 [15,] 0 0 1 [16,] 0 0 1 [17,] 0 0 1 [18,] 0 1 0 [19,] 0 1 0 [20,] 0 0 1 [21,] 0 1 0 [22,] 0 1 0 [23,] 0 1 0 [24,] 0 1 0 [25,] 0 1 0 [26,] 0 1 0 [27,] 1 0 0 > beta.m<-solve(t(g) %*% x, t(g) %*% y) > beta.m [,1] [1,] -14.54719339 [2,] 1.49580457 [3,] -0.09961272 > sum((y-x %*% beta.m)^2) [1] 0.5766745 > summary(glm.out)$cov.unscaled * summary(glm.out)$dispersion (Intercept) alpha asin (Intercept) 0.0012530451 -0.0004108272 -0.000898483 alpha -0.0004108272 0.0017416292 0.000236933 asin -0.0008984831 0.0002369330 0.002583829 > a<-solve(t(g) %*% x, t(g)) > a %*% t(a) * summary(glm.out)$dispersion [,1] [,2] [,3] [1,] 0.0016040221 -0.0005550267 -0.0018779456 [2,] -0.0005550267 0.0020089745 0.0005190049 [3,] -0.0018779456 0.0005190049 0.0053867812 > plot(asin,alpha,type="n",main="Figure 1.2 Mayer's groups") > for(i in 1:3) text(asin[grp==i],alpha[grp==i],i) > > > #------------------------------------------------------------------- > # 1.3 MULTIPLE REGRESSION > #------------------------------------------------------------------- > m<-matrix(scan("/glim/datasets/peru"),byrow=T,ncol=10) Read 390 items > m [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 21 1 71.0 1629 8.0 7.0 12.7 88 170 76 [2,] 22 6 56.5 1569 3.3 5.0 8.0 64 120 60 [3,] 24 5 56.0 1561 3.3 1.3 4.3 68 125 75 [4,] 24 1 61.0 1619 3.7 3.0 4.3 52 148 120 [5,] 25 1 65.0 1566 9.0 12.7 20.7 72 140 78 [6,] 27 19 62.0 1639 3.0 3.3 5.7 72 106 72 [7,] 28 5 53.0 1494 7.3 4.7 8.0 64 120 76 [8,] 28 25 53.0 1568 3.7 4.3 0.0 80 108 62 [9,] 31 6 65.0 1540 10.3 9.0 10.0 76 124 70 [10,] 32 13 57.0 1530 5.7 4.0 6.0 60 134 64 [11,] 33 13 66.5 1622 6.0 5.7 8.3 68 116 76 [12,] 33 10 59.1 1486 6.7 5.3 10.3 72 114 74 [13,] 34 15 64.0 1578 3.3 5.3 7.0 88 130 80 [14,] 35 18 69.5 1645 9.3 5.0 7.0 60 118 68 [15,] 35 2 64.0 1648 3.0 3.7 6.7 60 138 78 [16,] 36 12 56.5 1521 3.3 5.0 11.7 72 134 86 [17,] 36 15 57.0 1547 3.0 3.0 6.0 84 120 70 [18,] 37 16 55.0 1505 4.3 5.0 7.0 64 120 76 [19,] 37 17 57.0 1473 6.0 5.3 11.7 72 114 80 [20,] 38 10 58.0 1538 8.7 6.0 13.0 64 124 64 [21,] 38 18 59.5 1513 5.3 4.0 7.7 80 114 66 [22,] 38 11 61.0 1653 4.0 3.3 4.0 76 136 78 [23,] 38 11 57.0 1566 3.0 3.0 3.0 60 126 72 [24,] 39 21 57.5 1580 4.0 3.0 5.0 64 124 62 [25,] 39 24 74.0 1647 7.3 6.3 15.7 64 128 84 [26,] 39 14 72.0 1620 6.3 7.7 13.3 68 134 92 [27,] 41 25 62.5 1637 6.0 5.3 8.0 76 112 80 [28,] 41 32 68.0 1528 10.0 5.0 11.3 60 128 82 [29,] 41 5 63.4 1647 5.3 4.3 13.7 76 134 92 [30,] 42 12 68.0 1605 11.0 7.0 10.7 88 128 90 [31,] 43 25 69.0 1625 5.0 3.0 6.0 72 140 72 [32,] 43 26 73.0 1615 12.0 4.0 5.7 68 138 74 [33,] 43 10 64.0 1640 5.7 3.0 7.0 60 118 66 [34,] 44 19 65.0 1610 8.0 6.7 7.7 74 110 70 [35,] 44 18 71.0 1572 3.0 4.7 4.3 72 142 84 [36,] 45 10 60.2 1534 3.0 3.0 3.3 56 134 70 [37,] 47 1 55.0 1536 3.0 3.0 4.0 64 116 54 [38,] 50 43 70.0 1630 4.0 6.0 11.7 72 132 90 [39,] 54 40 87.0 1542 11.3 11.7 11.3 92 152 88 > age<-m[,1] > mig<-m[,2] > wt<-m[,3] > ht<-m[,4] > chin<-m[,5] > arm<-m[,6] > calf<-m[,7] > pulse<-m[,8] > sys<-m[,9] > dias<-m[,10] > plot(mig,sys,main="Figure 1.3 Peru, effect of mig") > summary(glm(sys~mig)) Call: glm(formula = sys ~ mig) Deviance Residuals: Min 1Q Median 3Q Max -20.927 -9.495 -1.836 7.153 41.028 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 129.0855 3.7851 34.103 <2e-16 *** mig -0.1136 0.2127 -0.534 0.596 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 175.1744) Null deviance: 6531.4 on 38 degrees of freedom Residual deviance: 6481.5 on 37 degrees of freedom AIC: 316.09 Number of Fisher Scoring iterations: 2 > pt(-0.5341708,37) [1] 0.2982093 > f<-(6531.436-6481.452 )/175.1744 > c(f,1-pf(f,1,37)) [1] 0.2853385 0.5964185 > -0.5341708^2 [1] -0.2853384 > 0.2982093*2 [1] 0.5964186 > plot(wt,sys,main="Figure 1.4 Peru, effect of wt") > summary(glm(sys~wt)) Call: glm(formula = sys ~ wt) Deviance Residuals: Min 1Q Median 3Q Max -20.2943 -8.4910 0.4457 6.6617 35.0399 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 66.5969 16.4639 4.045 0.000255 *** wt 0.9629 0.2591 3.716 0.000665 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 128.5421) Null deviance: 6531.4 on 38 degrees of freedom Residual deviance: 4756.1 on 37 degrees of freedom AIC: 304.02 Number of Fisher Scoring iterations: 2 > summary(glm(sys~wt+mig)) Call: glm(formula = sys ~ wt + mig) Deviance Residuals: Min 1Q Median 3Q Max -17.469 -7.878 1.076 6.292 24.113 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 50.3191 15.8184 3.181 0.00302 ** wt 1.3541 0.2672 5.067 1.22e-05 *** mig -0.5718 0.1879 -3.043 0.00436 ** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 105.0877) Null deviance: 6531.4 on 38 degrees of freedom Residual deviance: 3783.2 on 36 degrees of freedom AIC: 297.09 Number of Fisher Scoring iterations: 2 > (4756.056-3783.157)/105.0877 [1] 9.257972 > -3.042691^2 [1] -9.257969 > pt(-3.042691,36) [1] 0.002179931 > f<-(6531.436-3783.157)/2/105.0877 > c(f,1-pf(f,2,36)) [1] 1.307612e+01 5.385457e-05 > plot(wt,mig,main="Figure 1.5 mig vs wt") > rmig<-reisd(glm(mig~wt)) > rsys<-resid(glm(sys~wt)) > plot(rmig,rsys,main="Figure 1.6 sys vs mig, removing wt") > summary(glm(rsys~rmig-1)) Call: glm(formula = rsys ~ rmig - 1) Deviance Residuals: Min 1Q Median 3Q Max -17.469 -7.878 1.076 6.292 24.113 Coefficients: Estimate Std. Error t value Pr(>|t|) rmig -0.5718 0.1829 -3.126 0.00339 ** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 99.55677) Null deviance: 4756.1 on 39 degrees of freedom Residual deviance: 3783.2 on 38 degrees of freedom AIC: 293.09 Number of Fisher Scoring iterations: 2 > 3783.157/36 [1] 105.0877 > plot(mig,sys,type="n",main="Figure 1.7 Point label = wt/10") > wf<-floor(wt/10); wf [1] 7 5 5 6 6 6 5 5 6 5 6 5 6 6 6 5 5 5 5 5 5 6 5 5 7 7 6 6 6 [30] 6 6 7 6 6 7 6 5 7 8 > for(i in 5:8) text(mig[wf==i],sys[wf==i],i) > summary(glm(sys~ht+mig)) Call: glm(formula = sys ~ ht + mig) Deviance Residuals: Min 1Q Median 3Q Max -24.224 -8.421 -1.743 7.137 37.910 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 40.36194 63.57169 0.635 0.530 ht 0.05639 0.04034 1.398 0.171 mig -0.13500 0.21058 -0.641 0.526 (Dispersion parameter for gaussian family taken to be 170.7686) Null deviance: 6531.4 on 38 degrees of freedom Residual deviance: 6147.7 on 36 degrees of freedom AIC: 316.03 Number of Fisher Scoring iterations: 2 > plot(ht,mig,main="Figure 1.8 mig vs ht") > summary(glm(sys~ht+mig+wt+chin+arm+calf+age)) Call: glm(formula = sys ~ ht + mig + wt + chin + arm + calf + age) Deviance Residuals: Min 1Q Median 3Q Max -20.1440 -6.8827 0.9546 5.9280 21.0118 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 129.38413 56.56907 2.287 0.0292 * ht -0.06731 0.04255 -1.582 0.1238 mig -0.56787 0.22040 -2.577 0.0150 * wt 2.11856 0.46216 4.584 7.05e-05 *** chin -1.35274 0.86637 -1.561 0.1286 arm -0.98328 1.33833 -0.735 0.4680 calf 0.23114 0.61969 0.373 0.7117 age -0.27770 0.28560 -0.972 0.3384 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for gaussian family taken to be 105.4732) Null deviance: 6531.4 on 38 degrees of freedom Residual deviance: 3269.7 on 31 degrees of freedom AIC: 301.40 Number of Fisher Scoring iterations: 2 > f<-(3783.157-3269.668)/5/105.4732 > c(f,1-pf(f,5,31)) [1] 0.9736862 0.4491722 >