> ############## > # Section 7.1 > ############## > > data(nes96, package="faraway") > party <- nes96$PID > levels(party) <- c("Democrat","Democrat","Independent","Independent", "Independent","Republican","Republican") > inca <- c(1.5,4,6,8,9.5,10.5,11.5,12.5,13.5,14.5,16,18.5,21,23.5, 27.5,32.5,37.5,42.5,47.5,55,67.5,82.5,97.5,115) > income <- inca[unclass(nes96$income)] > rnes96 <- data.frame(party, income, education=nes96$educ, age=nes96$age) > summary(rnes96) party income education age Democrat :380 Min. : 1.50 MS : 13 Min. :19.00 Independent:239 1st Qu.: 23.50 HSdrop: 52 1st Qu.:34.00 Republican :325 Median : 37.50 HS :248 Median :44.00 Mean : 46.58 Coll :187 Mean :47.04 3rd Qu.: 67.50 CCdeg : 90 3rd Qu.:58.00 Max. :115.00 BAdeg :227 Max. :91.00 MAdeg :127 > rnes96[1:10,] party income education age 1 Republican 1.5 HS 36 2 Democrat 1.5 Coll 20 3 Democrat 1.5 BAdeg 24 4 Democrat 1.5 BAdeg 28 5 Democrat 1.5 BAdeg 68 6 Democrat 1.5 Coll 21 7 Democrat 1.5 Coll 77 8 Independent 1.5 Coll 21 9 Independent 1.5 Coll 31 10 Democrat 1.5 HS 39 > > ################## > # Education: Middle-school (MS), High-School Drop (HSdrop), > # High-Schole (HS), College (Coll), College-degress (CCdeg), > # Bachelor Degree (BAdeg), Master Degree (MAdeg). > ################## > > > > library(nnet) > mmod <- multinom(party ~ age + education + income, rnes96) # weights: 30 (18 variable) initial value 1037.090001 iter 10 value 990.568608 iter 20 value 984.319052 final value 984.166272 converged > summary(mmod) Call: multinom(formula = party ~ age + education + income, data = rnes96) Coefficients: (Intercept) age education.L education.Q education.C education^4 education^5 education^6 income Independent -1.197260 0.0001534525 0.06351451 -0.1217038 0.1119542 -0.07657336 0.1360851 0.15427826 0.01623911 Republican -1.642656 0.0081943691 1.19413345 -1.2292869 0.1544575 -0.02827297 -0.1221176 -0.03741389 0.01724679 Std. Errors: (Intercept) age education.L education.Q education.C education^4 education^5 education^6 income Independent 0.3265951 0.005374592 0.4571884 0.4142859 0.3498491 0.2883031 0.2494706 0.2171578 0.003108585 Republican 0.3312877 0.004902668 0.6502670 0.6041924 0.4866432 0.3605620 0.2696036 0.2031859 0.002881745 Residual Deviance: 1968.333 AIC: 2004.333 > > op <- options(contrasts=c("contr.treatment", "contr.treatment")) > mmod <- multinom(party ~ age + education + income, rnes96) # weights: 30 (18 variable) initial value 1037.090001 iter 10 value 990.364722 iter 20 value 984.508641 final value 984.166272 converged > summary(mmod) Call: multinom(formula = party ~ age + education + income, data = rnes96) Coefficients: (Intercept) age educationHSdrop educationHS educationColl educationCCdeg educationBAdeg educationMAdeg income Independent -1.373895 0.0001539014 0.2704482 0.2458744 0.09119446 0.3269554 0.1082654 0.1933497 0.01623914 Republican -3.048576 0.0081945031 0.9876547 1.6915600 1.95336096 1.8835335 1.8708213 1.4539589 0.01724696 Std. Errors: (Intercept) age educationHSdrop educationHS educationColl educationCCdeg educationBAdeg educationMAdeg income Independent 0.766464 0.005374592 0.7460643 0.6992364 0.713322 0.7382877 0.7151263 0.7287344 0.003108590 Republican 1.111205 0.004902674 1.1237993 1.0721857 1.076931 1.0940829 1.0792352 1.0914118 0.002881749 Residual Deviance: 1968.333 AIC: 2004.333 > ############# > # It has two groups of coefficients since the response has three levels. > # Democratic party is the baseline > ############# > > mmodi <- step(mmod) Start: AIC=2004.33 party ~ age + education + income trying - age # weights: 27 (16 variable) initial value 1037.090001 iter 10 value 988.896864 iter 20 value 985.822223 final value 985.812737 converged trying - education # weights: 12 (6 variable) initial value 1037.090001 iter 10 value 992.269502 final value 992.269484 converged trying - income # weights: 27 (16 variable) initial value 1037.090001 iter 10 value 1009.025560 iter 20 value 1006.961593 final value 1006.955275 converged Df AIC - education 6 1996.539 - age 16 2003.625 18 2004.333 - income 16 2045.911 # weights: 12 (6 variable) initial value 1037.090001 iter 10 value 992.269502 final value 992.269484 converged Step: AIC=1996.54 party ~ age + income trying - age # weights: 9 (4 variable) initial value 1037.090001 final value 992.712152 converged trying - income # weights: 9 (4 variable) initial value 1037.090001 final value 1020.425203 converged Df AIC - age 4 1993.424 6 1996.539 - income 4 2048.850 # weights: 9 (4 variable) initial value 1037.090001 final value 992.712152 converged Step: AIC=1993.42 party ~ income trying - income # weights: 6 (2 variable) initial value 1037.090001 final value 1020.636052 converged Df AIC 4 1993.424 - income 2 2045.272 > mmode <- multinom(party ~ age + income, rnes96) # weights: 12 (6 variable) initial value 1037.090001 iter 10 value 992.269502 final value 992.269484 converged > summary(mmode) Call: multinom(formula = party ~ age + income, data = rnes96) Coefficients: (Intercept) age income Independent -1.187181 0.0002260719 0.0161244 Republican -1.154787 0.0041328499 0.0178682 Std. Errors: (Intercept) age income Independent 0.2941042 0.005143979 0.002861903 Republican 0.2736512 0.004700426 0.002666931 Residual Deviance: 1984.539 AIC: 1996.539 > deviance(mmode) - deviance(mmod) [1] 16.20642 > pchisq(16.206,mmod$edf-mmode$edf,lower=F) [1] 0.181982 > > ############## > # Age and Income are continuous variables. > # The comparison between the previous two can tell us the > # significance of education. > ############## > > inclevels <- 0:110 > preds <- data.frame(income=inclevels,predict(mmodi,data.frame(income=inclevels),type="probs")) > preds[1:10,] income Democrat Independent Republican 1 0 0.5898168 0.1821588 0.2280244 2 1 0.5857064 0.1838228 0.2304708 3 2 0.5815839 0.1854890 0.2329271 4 3 0.5774498 0.1871572 0.2353929 5 4 0.5733047 0.1888271 0.2378682 6 5 0.5691492 0.1904984 0.2403524 7 6 0.5649837 0.1921708 0.2428455 8 7 0.5608089 0.1938442 0.2453469 9 8 0.5566253 0.1955183 0.2478565 10 9 0.5524335 0.1971927 0.2503738 > rnes96[1:10,] party income education age 1 Republican 1.5 HS 36 2 Democrat 1.5 Coll 20 3 Democrat 1.5 BAdeg 24 4 Democrat 1.5 BAdeg 28 5 Democrat 1.5 BAdeg 68 6 Democrat 1.5 Coll 21 7 Democrat 1.5 Coll 77 8 Independent 1.5 Coll 21 9 Independent 1.5 Coll 31 10 Democrat 1.5 HS 39 > mmodi Call: multinom(formula = party ~ income, data = rnes96) Coefficients: (Intercept) income Independent -1.1749331 0.01608683 Republican -0.9503591 0.01766457 Residual Deviance: 1985.424 AIC: 1993.424 > > eta1 <- -1.1749331+0.01608683*inclevels[1] > eta2 <- -0.9503591+0.01766457*inclevels[1] > > aa <- c(1,exp(eta1),exp(eta2)) > aa/sum(aa) [1] 0.5898168 0.1821588 0.2280244 > > ############ > # Look at the ways in computing the predicted probabilities. > ############ > > predict(mmodi,data.frame(income=inclevels)) [1] Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat [16] Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat [31] Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat [46] Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Democrat Republican Republican Republican Republican Republican Republican [61] Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican [76] Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican [91] Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican Republican [106] Republican Republican Republican Republican Republican Republican Levels: Democrat Independent Republican > xtabs( ~ predict(mmodi) + rnes96$party) rnes96$party predict(mmodi) Democrat Independent Republican Democrat 284 123 166 Independent 0 0 0 Republican 96 116 159 > (284+0+159)/nrow(rnes96) [1] 0.4692797 > summary(mmodi) Call: multinom(formula = party ~ income, data = rnes96) Coefficients: (Intercept) income Independent -1.1749331 0.01608683 Republican -0.9503591 0.01766457 Std. Errors: (Intercept) income Independent 0.1536103 0.002849738 Republican 0.1416859 0.002652532 Residual Deviance: 1985.424 AIC: 1993.424 > > > cc <- c(0,-1.17493,-0.95036) > exp(cc)/sum(exp(cc)) [1] 0.5898166 0.1821593 0.2280242 > predict(mmodi,data.frame(income=0),type="probs") Democrat Independent Republican 0.5898168 0.1821588 0.2280244 > > ############### > # Those are prediced probabilities when income=0. > ############## > > (pp <- predict(mmodi,data.frame(income=c(0,1)),type="probs")) Democrat Independent Republican 1 0.5898168 0.1821588 0.2280244 2 0.5857064 0.1838228 0.2304708 > log(pp[1,1]*pp[2,2]/(pp[1,2]*pp[2,1])) [1] 0.01608683 > log(pp[1,1]*pp[2,3]/(pp[1,3]*pp[2,1])) [1] 0.01766457 > > ############ > # log ddds ratio > ############ > > ##################### > # More for multinomial models > ##################### > ######################## > # Another Example: Mice Toxicity Study > ######################## > > R1 <- c(15,17,22,38,144) > R2 <- c(1,0,7,59,132) > R3 <- c(281,225,283,202,9) > x <- c(0,62.5,125,250,500) > Y <- cbind(R2,R1,R3) > > freq <- c(R1,R2,R3) > yy <- as.factor(c(rep("Nonlive",5),rep("Malformation",5),rep("Normal",5))) > xx <- c(x,x,x) > yy [1] Nonlive Nonlive Nonlive Nonlive Nonlive Malformation Malformation Malformation Malformation Malformation Normal Normal Normal [14] Normal Normal Levels: Malformation Nonlive Normal > levels(yy) <- c("Nonlive","Malformation","Normal") > > Y R2 R1 R3 [1,] 1 15 281 [2,] 0 17 225 [3,] 7 22 283 [4,] 59 38 202 [5,] 132 144 9 > yy [1] Malformation Malformation Malformation Malformation Malformation Nonlive Nonlive Nonlive Nonlive Nonlive Normal Normal Normal [14] Normal Normal Levels: Nonlive Malformation Normal > freq [1] 15 17 22 38 144 1 0 7 59 132 281 225 283 202 9 > xx [1] 0.0 62.5 125.0 250.0 500.0 0.0 62.5 125.0 250.0 500.0 0.0 62.5 125.0 250.0 500.0 > > ####### > # Multinomial Model > ####### > > g <- multinom(yy~xx,weights=freq) # weights: 9 (4 variable) initial value 1576.508634 iter 10 value 753.227671 final value 753.220332 converged > summary(g) Call: multinom(formula = yy ~ xx, weights = freq) Coefficients: (Intercept) xx Malformation 0.9833348 -0.002100678 Normal 4.9527859 -0.014009587 Std. Errors: (Intercept) xx Malformation 0.2706012 0.0006455004 Normal 0.2493653 0.0007867021 Residual Deviance: 1506.441 AIC: 1514.441 > g$dev [1] 1506.441 > g$edf [1] 4 > > ######## > # Proportional Odds Model > ######## > > library(MASS) > > g <- polr(yy~xx,weights=freq) > summary(g) Re-fitting to get Hessian Call: polr(formula = yy ~ xx, weights = freq) Coefficients: Value Std. Error t value xx -0.0102 0.0004508 -22.64 Intercepts: Value Std. Error t value Nonlive|Malformation -5.0144 0.1949 -25.7322 Malformation|Normal -3.2072 0.1402 -22.8679 Residual Deviance: 1577.477 AIC: 1583.477 > g$dev [1] 1577.477 > g$edf [1] 3 > > ######## > # Saturated Model > ######## > > g.saturated <- multinom(yy~factor(xx),weights=freq) # weights: 18 (10 variable) initial value 1576.508634 iter 10 value 750.025106 iter 20 value 724.633177 iter 30 value 724.577708 iter 40 value 724.469453 final value 724.468130 converged > summary(g.saturated) Call: multinom(formula = yy ~ factor(xx), weights = freq) Coefficients: (Intercept) factor(xx)62.5 factor(xx)125 factor(xx)250 factor(xx)500 Malformation 2.707031 8.790348 -1.561982 -3.146892 -2.620011 Normal 5.637388 8.442839 -1.937959 -4.406626 -8.322952 Std. Errors: (Intercept) factor(xx)62.5 factor(xx)125 factor(xx)250 factor(xx)500 Malformation 1.032330 76.10868 1.119822 1.053076 1.039339 Normal 1.001296 76.10791 1.071900 1.012173 1.058905 Residual Deviance: 1448.936 AIC: 1468.936 > g.saturated$dev [1] 1448.936 > g.saturated$edf [1] 10 > g.saturated$fitted Nonlive Malformation Normal 1 3.370257e-03 0.05050238 0.94612736 2 7.135109e-07 0.07025044 0.92974885 3 2.243822e-02 0.07051429 0.90704749 4 1.973179e-01 0.12709771 0.67558434 5 4.631557e-01 0.50526509 0.03157923 6 3.370257e-03 0.05050238 0.94612736 7 7.135109e-07 0.07025044 0.92974885 8 2.243822e-02 0.07051429 0.90704749 9 1.973179e-01 0.12709771 0.67558434 10 4.631557e-01 0.50526509 0.03157923 11 3.370257e-03 0.05050238 0.94612736 12 7.135109e-07 0.07025044 0.92974885 13 2.243822e-02 0.07051429 0.90704749 14 1.973179e-01 0.12709771 0.67558434 15 4.631557e-01 0.50526509 0.03157923 > g.saturated$fitted[1:5,] Nonlive Malformation Normal 1 3.370257e-03 0.05050238 0.94612736 2 7.135109e-07 0.07025044 0.92974885 3 2.243822e-02 0.07051429 0.90704749 4 1.973179e-01 0.12709771 0.67558434 5 4.631557e-01 0.50526509 0.03157923 > > Y/apply(Y,1,sum) R2 R1 R3 [1,] 0.003367003 0.05050505 0.94612795 [2,] 0.000000000 0.07024793 0.92975207 [3,] 0.022435897 0.07051282 0.90705128 [4,] 0.197324415 0.12709030 0.67558528 [5,] 0.463157895 0.50526316 0.03157895 > > ########## > # Why the deviance of the saturated model is not zero. > # Let us try > ########## > > -2*sum(Y*log(g.saturated$fitted[1:5,])) [1] 1448.936 > > ######### > # The residual deviance is computed based the Bernoulli model > # instead of a multinomial model. > ######## > > g$fitted[1:5,] Nonlive Malformation Normal 1 0.006597497 0.03229814 0.9611044 2 0.012411453 0.05872236 0.9288662 3 0.023229073 0.10334420 0.8734267 4 0.078475024 0.26316261 0.6583624 5 0.521972485 0.34737381 0.1306537 > Y.hat <- g$fitted[1:5,]*apply(Y,1,sum) > index <- Y>0 > G2 <- 2*sum(Y[index]*log(Y[index]/Y.hat[index])) > G2 [1] 128.5409 > g$dev-g.saturated$dev [1] 128.5406 > > ############### > # Section 7.2: Simply introduce the problem. > ############### > > library(MASS) > mlda <- lda(party ~ age + income, rnes96) > mlda Call: lda(party ~ age + income, data = rnes96) Prior probabilities of groups: Democrat Independent Republican 0.4025424 0.2531780 0.3442797 Group means: age income Democrat 47.06579 37.63158 Independent 46.50628 51.65481 Republican 47.41231 53.29846 Coefficients of linear discriminants: LD1 LD2 age 0.005325725 0.0608207675 income 0.033215484 -0.0001139976 Proportion of trace: LD1 LD2 0.9926 0.0074 > > ############# > # Linear Dicriminant Analysis > ############# > > preds <- predict(mlda) > head(preds$posterior) Democrat Independent Republican 1 0.5866620 0.1894525 0.2238855 2 0.5953707 0.1913508 0.2132785 3 0.5932182 0.1908863 0.2158956 4 0.5910492 0.1904151 0.2185357 5 0.5684622 0.1853332 0.2462046 6 0.5948341 0.1912353 0.2139306 > xtabs( ~ predict(mlda)$class + rnes96$party) rnes96$party predict(mlda)$class Democrat Independent Republican Democrat 298 141 184 Independent 0 0 0 Republican 82 98 141 > > ############## > # Section 7.3 > ############## > > data(cns, package="faraway") > cns Area NoCNS An Sp Other Water Work 1 Cardiff 4091 5 9 5 110 NonManual 2 Newport 1515 1 7 0 100 NonManual 3 Swansea 2394 9 5 0 95 NonManual 4 GlamorganE 3163 9 14 3 42 NonManual 5 GlamorganW 1979 5 10 1 39 NonManual 6 GlamorganC 4838 11 12 2 161 NonManual 7 MonmouthV 2362 6 8 4 83 NonManual 8 MonmouthOther 1604 3 6 0 122 NonManual 9 Cardiff 9424 31 33 14 110 Manual 10 Newport 4610 3 15 6 100 Manual 11 Swansea 5526 19 30 4 95 Manual 12 GlamorganE 13217 55 71 19 42 Manual 13 GlamorganW 8195 30 44 10 39 Manual 14 GlamorganC 7803 25 28 12 161 Manual 15 MonmouthV 9962 36 37 13 83 Manual 16 MonmouthOther 3172 8 13 3 122 Manual > cns$CNS <- cns$An+cns$Sp+cns$Other > plot(log(CNS/NoCNS) ~ Water, cns, pch=as.character(Work)) > > ########## > # The new variable CNS is the sum of three variables > ########## > > binmodw <- glm(cbind(CNS,NoCNS) ~ Water + Work, cns, family=binomial) > summary(binmodw) Call: glm(formula = cbind(CNS, NoCNS) ~ Water + Work, family = binomial, data = cns) Deviance Residuals: Min 1Q Median 3Q Max -2.65570 -0.30179 -0.03131 0.57213 1.32998 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.4325803 0.0897889 -49.367 < 2e-16 *** Water -0.0032644 0.0009684 -3.371 0.000749 *** WorkNonManual -0.3390577 0.0970943 -3.492 0.000479 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 41.047 on 15 degrees of freedom Residual deviance: 12.363 on 13 degrees of freedom AIC: 102.49 Number of Fisher Scoring iterations: 4 > > ############ > # Water represented a subspace of area. It caused > # that variable Area and Water were confounded. > ############ > > binmoda <- glm(cbind(CNS,NoCNS) ~ Area + Work, cns, family=binomial) > summary(binmoda) Call: glm(formula = cbind(CNS, NoCNS) ~ Area + Work, family = binomial, data = cns) Deviance Residuals: Min 1Q Median 3Q Max -0.8161 -0.3412 0.0202 0.3201 0.7632 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.84428 0.10445 -46.381 < 2e-16 *** AreaGlamorganC 0.01798 0.14705 0.122 0.902684 AreaGlamorganE 0.34026 0.12795 2.659 0.007831 ** AreaGlamorganW 0.28045 0.14338 1.956 0.050477 . AreaMonmouthOther -0.02721 0.20226 -0.135 0.892988 AreaMonmouthV 0.12710 0.14199 0.895 0.370704 AreaNewport -0.33496 0.20450 -1.638 0.101435 AreaSwansea 0.16427 0.15950 1.030 0.303033 WorkNonManual -0.34580 0.09733 -3.553 0.000381 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 41.0475 on 15 degrees of freedom Residual deviance: 3.0771 on 7 degrees of freedom AIC: 105.2 Number of Fisher Scoring iterations: 4 > > anova(binmodw,binmoda,test="Chi") Analysis of Deviance Table Model 1: cbind(CNS, NoCNS) ~ Water + Work Model 2: cbind(CNS, NoCNS) ~ Area + Work Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 13 12.3628 2 7 3.0771 6 9.2857 0.1581 > > ########### > # The likelihood ratio test > ########### > > library(faraway) > halfnorm(residuals(binmodw)) > sumary(binmodw) Estimate Std. Error z value Pr(>|z|) (Intercept) -4.4325803 0.0897889 -49.3667 < 2.2e-16 Water -0.0032644 0.0009684 -3.3709 0.0007492 WorkNonManual -0.3390577 0.0970943 -3.4920 0.0004793 n = 16 p = 3 Deviance = 12.36281 Null Deviance = 41.04749 (Difference = 28.68468) > exp(-0.339058) [1] 0.7124411 > > ############# > # The above is the odds ratio > ############# > > cmmod <- multinom(cbind(An,Sp,Other) ~ Water + Work, cns) # weights: 12 (6 variable) initial value 762.436928 iter 10 value 685.762336 final value 685.762238 converged > summary(cmmod) Call: multinom(formula = cbind(An, Sp, Other) ~ Water + Work, data = cns) Coefficients: (Intercept) Water WorkNonManual Sp 0.3752018 -0.001297048 0.1157557 Other -1.1225496 0.002182689 -0.2702791 Std. Errors: (Intercept) Water WorkNonManual Sp 0.1900329 0.002032944 0.2086875 Other 0.2795628 0.002896841 0.3247247 Residual Deviance: 1371.524 AIC: 1383.524 > > cmmod$edf [1] 6 > cmmod$fitted.value An Sp Other 1 0.3659780 0.5184562 0.11556582 2 0.3644207 0.5229897 0.11258958 3 0.3636278 0.5252470 0.11112519 4 0.3546765 0.5487745 0.09654896 5 0.3541420 0.5500838 0.09577426 6 0.3732862 0.4949607 0.13175309 7 0.3616868 0.5306386 0.10767453 8 0.3677950 0.5129833 0.11922171 9 0.3737549 0.4715979 0.15464719 10 0.3727047 0.4764122 0.15088314 11 0.3721606 0.4788118 0.14902756 12 0.3656545 0.5039186 0.13042682 13 0.3652479 0.5053207 0.12943146 14 0.3782856 0.4467622 0.17495218 15 0.3708042 0.4845501 0.14464574 16 0.3749470 0.4657954 0.15925760 > > summary(cmmod)$edf [1] 6 > summary(cmmod)$deviance [1] 1371.524 > > ########## > # There are two groups of coefficients in the outout. > # Look at the residual deviance. You need to compare it with the saturated model. > ########## > > mod.sat <- multinom(cbind(An,Sp,Other) ~factor(Water)*factor(Work), cns) # weights: 51 (32 variable) initial value 762.436928 iter 10 value 676.578660 iter 20 value 668.935610 iter 30 value 668.593789 iter 40 value 668.563493 final value 668.563246 converged > summary(mod.sat) Call: multinom(formula = cbind(An, Sp, Other) ~ factor(Water) * factor(Work), data = cns) Coefficients: (Intercept) factor(Water)42 factor(Water)83 factor(Water)95 factor(Water)100 factor(Water)110 factor(Water)122 factor(Water)161 factor(Work)NonManual Sp 0.3829891 -0.12764895 -0.35559238 0.07376707 1.226416 -0.3204569 0.1025095 -0.2696494 0.3102572 Other -1.0986122 0.03570577 0.08003919 -0.45951920 1.791710 0.3037187 0.1176847 0.3646658 -0.5107928 factor(Water)42:factor(Work)NonManual factor(Water)83:factor(Work)NonManual factor(Water)95:factor(Work)NonManual factor(Water)100:factor(Work)NonManual Sp -0.1237784 -0.04994875 -1.354802 0.02627559 Other 0.4750860 1.12393561 -12.916330 -17.58100181 factor(Water)110:factor(Work)NonManual factor(Water)122:factor(Work)NonManual factor(Water)161:factor(Work)NonManual Sp 0.2149784 -0.1026033 -0.3366222 Other 1.3056611 -13.0394844 -0.4600733 Std. Errors: (Intercept) factor(Water)42 factor(Water)83 factor(Water)95 factor(Water)100 factor(Water)110 factor(Water)122 factor(Water)161 factor(Work)NonManual Sp 0.2367712 0.2971983 0.3329645 0.3768627 0.6753145 0.3444155 0.5079182 0.3630085 0.5967157 Other 0.3651480 0.4518260 0.4878873 0.6602732 0.7958172 0.4868452 0.7692143 0.5066217 1.1547161 factor(Water)42:factor(Work)NonManual factor(Water)83:factor(Work)NonManual factor(Water)95:factor(Work)NonManual factor(Water)100:factor(Work)NonManual Sp 0.7555628 0.8381787 8.678397e-01 1.378020e+00 Other 1.3596413 1.3618868 2.422528e-05 2.591136e-07 factor(Water)110:factor(Work)NonManual factor(Water)122:factor(Work)NonManual factor(Water)161:factor(Work)NonManual Sp 0.8542479 1.028587e+00 0.7784765 Other 1.3553794 1.291124e-05 1.4309541 Residual Deviance: 1337.126 AIC: 1401.126 > > mod.sat$edf [1] 32 > > count.obs <- cbind(cns$An,cns$Sp,cns$Other) > rate.obs <- count.obs/apply(count.obs,1,sum) > > mod.sat$fitted.value An Sp Other 1 0.2631620 0.4736827 2.631554e-01 2 0.1249969 0.8750031 3.473291e-09 3 0.6428575 0.3571423 1.995729e-07 4 0.3461564 0.5384582 1.153854e-01 5 0.3124800 0.6250220 6.249806e-02 6 0.4400100 0.4799933 7.999668e-02 7 0.3333273 0.4444467 2.222260e-01 8 0.3333321 0.6666678 1.629503e-07 9 0.3974313 0.4230771 1.794916e-01 10 0.1250041 0.6250002 2.499957e-01 11 0.3584907 0.5660366 7.547272e-02 12 0.3793122 0.4896543 1.310335e-01 13 0.3571434 0.5238087 1.190478e-01 14 0.3846119 0.4307701 1.846179e-01 15 0.4186053 0.4302322 1.511625e-01 16 0.3333391 0.5416710 1.249899e-01 > round(abs(mod.sat$fitted.value-rate.obs),4) An Sp Other 1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 0 0 0 7 0 0 0 8 0 0 0 9 0 0 0 10 0 0 0 11 0 0 0 12 0 0 0 13 0 0 0 14 0 0 0 15 0 0 0 16 0 0 0 > > cmmod$deviance-mod.sat$deviance [1] 34.39798 > mod.sat$edf-cmmod$edf [1] 26 > > 1-pchisq(cmmod$deviance-mod.sat$deviance,mod.sat$edf-cmmod$edf) [1] 0.1252942 > > ############# > # The above is the saturated model. > # The p-value of deviance goodness-of-fit is 0.1252. > ############# > > nmod <- step(cmmod) Start: AIC=1383.52 cbind(An, Sp, Other) ~ Water + Work trying - Water # weights: 9 (4 variable) initial value 762.436928 iter 10 value 686.562074 final value 686.562063 converged trying - Work # weights: 9 (4 variable) initial value 762.436928 final value 686.580556 converged Df AIC - Water 4 1381.124 - Work 4 1381.161 6 1383.524 # weights: 9 (4 variable) initial value 762.436928 iter 10 value 686.562074 final value 686.562063 converged Step: AIC=1381.12 cbind(An, Sp, Other) ~ Work trying - Work # weights: 6 (2 variable) initial value 762.436928 final value 687.227416 converged Df AIC - Work 2 1378.455 4 1381.124 # weights: 6 (2 variable) initial value 762.436928 final value 687.227416 converged Step: AIC=1378.45 cbind(An, Sp, Other) ~ 1 > nmod Call: multinom(formula = cbind(An, Sp, Other) ~ 1, data = cns) Coefficients: (Intercept) Sp 0.2896333 Other -0.9808293 Residual Deviance: 1374.455 AIC: 1378.455 > > ######## > # The above is the intercept only model. > ######## > > cc <- c(0,0.28963,-0.98083) > names(cc) <- c("An","Sp","Other") > exp(cc)/sum(exp(cc)) An Sp Other 0.3688767 0.4927946 0.1383287 > multinom(cbind(NoCNS,An,Sp,Other) ~ Water + Work, cns) # weights: 16 (9 variable) initial value 117209.801938 iter 10 value 27783.394982 iter 20 value 20770.391265 iter 30 value 11403.855681 iter 40 value 5238.502280 iter 50 value 4717.065810 final value 4695.482597 converged Call: multinom(formula = cbind(NoCNS, An, Sp, Other) ~ Water + Work, data = cns) Coefficients: (Intercept) Water WorkNonManual An -5.455142 -0.0029088415 -0.3638609 Sp -5.071037 -0.0043231569 -0.2435803 Other -6.594755 -0.0005127307 -0.6420771 Residual Deviance: 9390.965 AIC: 9408.965 > > ############# > # Section 7.4 > ############# > > x <- seq(-5,5,by=0.1) > xa <- c(21,41,71) > plot(x,dlogis(x),type="l",axes=FALSE,ylab="",xlab="") > segments(x[xa],rep(0,3),x[xa],dlogis(x[xa])) > axis(1,at=x[xa],c(expression(theta[1]),expression(theta[2]),expression(theta[3]))) > > ########### > # Read the interpretation of the plot by yourself. > # It is on the textbook. > ########### > > ########## > # Proportional Odds models. > # you need to pay attention the sign of beta > ########## > > library(MASS) > pomod <- polr(party ~ age + education + income, rnes96) > c(deviance(pomod),pomod$edf) [1] 1984.211 10.000 > c(deviance(mmod),mmod$edf) [1] 1968.333 18.000 > > ########### > # The above provides the likelihood ratio test > ########### > > pomodi <- step(pomod) Start: AIC=2004.21 party ~ age + education + income Df AIC - education 6 2002.8 2004.2 - age 1 2004.4 - income 1 2038.6 Step: AIC=2002.83 party ~ age + income Df AIC - age 1 2001.4 2002.8 - income 1 2047.2 Step: AIC=2001.36 party ~ income Df AIC 2001.4 - income 1 2045.3 > deviance(pomodi)-deviance(pomod) [1] 11.15136 > pchisq(11.151,pomod$edf-pomodi$edf,lower=F) [1] 0.1321668 > > pim <- with(rnes96,prop.table(table(income,party),1)) > logit(pim[,1])-logit(pim[,1]+pim[,2]) 1.5 4 6 8 9.5 10.5 11.5 12.5 13.5 14.5 16 18.5 21 23.5 27.5 -0.9007865 -2.0614230 -0.7576857 -1.0033021 -2.3025851 -0.3083014 -0.7985077 -1.8971200 -1.2527630 -1.1786550 -0.4128452 -0.3542428 -1.5141277 -1.6534548 -0.7467847 32.5 37.5 42.5 47.5 55 67.5 82.5 97.5 115 -0.5225217 -0.9232594 -1.0296194 -0.8219801 -1.4276009 -1.1826099 -0.9867640 -1.4829212 -1.7066017 > > ############ > # How to interpret the last model. > # Note that there is only one group of slopes, although it has two intecepts. > ############ > > summary(pomodi) Re-fitting to get Hessian Call: polr(formula = party ~ income, data = rnes96) Coefficients: Value Std. Error t value income 0.01312 0.001971 6.657 Intercepts: Value Std. Error t value Democrat|Independent 0.2091 0.1123 1.8627 Independent|Republican 1.2916 0.1201 10.7526 Residual Deviance: 1995.363 AIC: 2001.363 > ilogit(0.209) [1] 0.5520606 > ilogit(1.292)-ilogit(0.209) [1] 0.2324249 > > ########### > # The above is the prediction when income is zero. > ########## > > ilogit(0.209-20*0.01312) [1] 0.4866532 > ilogit(1.292-20*0.01312)-ilogit(0.209-20*0.01312) [1] 0.2501852 > > ########## > # The above is the prediction when income is 20 > ########## > > inclevels <- seq(0,100,by=20) > predict(pomodi,data.frame(income=inclevels,row.names=inclevels), type="probs") Democrat Independent Republican 0 0.5520865 0.2323241 0.2155895 20 0.4866801 0.2500729 0.2632470 40 0.4217267 0.2610935 0.3171797 60 0.3593732 0.2641117 0.3765151 80 0.3014331 0.2587658 0.4398011 100 0.2492008 0.2456926 0.5051067 > > ###################### > # Try to understand how to do the prediction. > ###################### > > x <- seq(-4,4,by=0.05) > plot(x,dlogis(x),type="l") > abline(v=c(0.209,1.292)) > abline(v=c(0.209,1.292)-50*0.013120,lty=2) > abline(v=c(0.209,1.292)-100*0.013120,lty=5) > > ############ > # The above is just the derivation of a plot > ############ > > ############## > # The following is just an extension. > # It has theoretical concern. > # Please read the book by yourself. > # I am going to skip that. > ############## > > library(VGAM) > nmod <- vglm(party ~ income, family=cumulative(parallel=FALSE),rnes96) > summary(nmod) Call: vglm(formula = party ~ income, family = cumulative(parallel = FALSE), data = rnes96) Pearson residuals: Min 1Q Median 3Q Max logit(P[Y<=1]) -1.752 -0.7586 -0.3645 0.9717 2.102 logit(P[Y<=2]) -1.664 -1.1919 0.3590 1.1379 1.293 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept):1 0.328868 0.121560 2.705 0.00682 ** (Intercept):2 1.148270 0.128719 8.921 < 2e-16 *** income:1 -0.016186 0.002358 -6.864 6.70e-12 *** income:2 -0.010486 0.002203 -4.761 1.93e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Number of linear predictors: 2 Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2]) Residual deviance: 1987.539 on 1884 degrees of freedom Log-likelihood: -993.7693 on 1884 degrees of freedom Number of iterations: 4 No Hauck-Donner effect found in any of the estimates Exponentiated coefficients: income:1 income:2 0.9839442 0.9895689 > > ################ > # It has two slopes and two intercepts. > # The domian of slopes is tricky. > # After it is solved, the rest is straightforward. > ################ > > pmod <- vglm(party ~ income, family=cumulative(parallel=TRUE), rnes96) > deviance(pmod) - deviance(nmod) [1] 7.824074 > 1-pchisq(7.82,1) [1] 0.00516712 > (1.148-0.329)/(0.0105-0.0162) [1] -143.6842 > > ############# > # The following is proportional probit and cloog models. > # Try to understand them by yourself. > ############# > > opmod <- polr(party ~ income, method="probit", rnes96) > summary(opmod) Re-fitting to get Hessian Call: polr(formula = party ~ income, data = rnes96, method = "probit") Coefficients: Value Std. Error t value income 0.008182 0.001208 6.775 Intercepts: Value Std. Error t value Democrat|Independent 0.1284 0.0694 1.8510 Independent|Republican 0.7976 0.0722 11.0399 Residual Deviance: 1994.892 AIC: 2000.892 > dems <- pnorm(0.128-inclevels*0.008182) > demind <- pnorm(0.798-inclevels*0.008182) > data.frame(Democrat=dems,Independent=demind-dems,Republican=1-demind,row.names=inclevels) Democrat Independent Republican 0 0.5509255 0.2366392 0.2124352 20 0.4857847 0.2512923 0.2629230 40 0.4210219 0.2600578 0.3189204 60 0.3583323 0.2622764 0.3793912 80 0.2992496 0.2577791 0.4429713 100 0.2450342 0.2469077 0.5080581 > polr(party ~ income, method="cloglog") Call: polr(formula = party ~ income, method = "cloglog") Coefficients: income 0.008271135 Intercepts: Democrat|Independent Independent|Republican -0.2866703 0.4564409 Residual Deviance: 2003.96 AIC: 2009.96 > > ########### > > levels(nes96$income) [1] "$3Kminus" "$3K-$5K" "$5K-$7K" "$7K-$9K" "$9K-$10K" "$10K-$11K" "$11K-$12K" "$12K-$13K" "$13K-$14K" "$14K-$15K" "$15K-$17K" "$17K-$20K" "$20K-$22K" [14] "$22K-$25K" "$25K-$30K" "$30K-$35K" "$35K-$40K" "$40K-$45K" "$45K-$50K" "$50K-$60K" "$60K-$75K" "$75K-$90K" "$90K-$105K" "$105Kplus" > inca <- c(1.5,4,6,8,9.5,10.5,11.5,12.5,13.5,14.5,16,18.5,21,23.5, 27.5,32.5,37.5,42.5,47.5,55,67.5,82.5,97.5,115) > nes96$sincome <- inca[unclass(nes96$income)] > lmod <- lm(sincome ~ age + educ + PID, nes96) > sumary(lmod) Estimate Std. Error t value Pr(>|t|) (Intercept) 45.837964 3.353281 13.6696 < 2.2e-16 age -0.049637 0.058214 -0.8527 0.39407 educ.L 39.285404 5.085048 7.7257 2.875e-14 educ.Q 2.855412 4.750839 0.6010 0.54796 educ.C 3.448575 4.008271 0.8604 0.38981 educ^4 2.529881 3.282474 0.7707 0.44107 educ^5 -4.470540 2.799246 -1.5971 0.11059 educ^6 -0.577267 2.341919 -0.2465 0.80536 PID.L 12.710325 2.163432 5.8751 5.882e-09 PID.Q -6.314047 2.861348 -2.2067 0.02758 PID.C 1.030562 2.374037 0.4341 0.66432 PID^4 6.241300 2.929586 2.1304 0.03340 PID^5 -2.517442 2.558453 -0.9840 0.32539 PID^6 -0.792289 3.653186 -0.2169 0.82835 n = 944, p = 14, Residual SE = 27.89922, R-Squared = 0.2 > lmod2 <- lm(sincome ~ unclass(educ) + unclass(PID), nes96) > anova(lmod2,lmod) Analysis of Variance Table Model 1: sincome ~ unclass(educ) + unclass(PID) Model 2: sincome ~ age + educ + PID Res.Df RSS Df Sum of Sq F Pr(>F) 1 941 737108 2 930 723881 11 13227 1.5448 0.1102 > levels(nes96$educ) [1] "MS" "HSdrop" "HS" "Coll" "CCdeg" "BAdeg" "MAdeg" > levels(nes96$PID) [1] "strDem" "weakDem" "indDem" "indind" "indRep" "weakRep" "strRep" > sumary(lmod2) Estimate Std. Error t value Pr(>|t|) (Intercept) 3.51014 3.03844 1.1552 0.2483 unclass(educ) 7.35939 0.57297 12.8443 < 2.2e-16 unclass(PID) 2.46345 0.40308 6.1115 1.444e-09 n = 944, p = 3, Residual SE = 27.98792, R-Squared = 0.19