library(faraway) data(nes96) sPID <- nes96$PID levels(sPID) <- c("Democrat","Democrat","Independent","Independent","Independent","Republican","Republican") summary(sPID) 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) nincome <- inca[unclass(nes96$income)] summary(nincome) table(nes96$educ) matplot(prop.table(table(nes96$educ,sPID),1),type="l",xlab="Education",ylab="Proportion",lty=c(1,2,5)) cutinc <- cut(nincome,7) il <- c(8,26,42,58,74,90,107) matplot(il,prop.table(table(cutinc,sPID),1),lty=c(1,2,5),type="l",ylab="Proportion",xlab="Income") cutage <- cut(nes96$age,7) al <- c(24,34,44,54,65,75,85) matplot(al,prop.table(table(cutage,sPID),1),lty=c(1,2,5),type="l",ylab="Proportion",xlab="Age") library(nnet) mmod <- multinom(sPID ~ age + educ + nincome, nes96) mmodi <- step(mmod) mmode <- multinom(sPID ~ age + nincome, nes96) deviance(mmode) - deviance(mmod) pchisq(16.206,mmod$edf-mmode$edf,lower=F) predict(mmodi,data.frame(nincome=il),type="probs") predict(mmodi,data.frame(nincome=il)) summary(mmodi) cc <- c(0,-1.17493,-0.95036) exp(cc)/sum(exp(cc)) predict(mmodi,data.frame(nincome=0),type="probs") (pp <- predict(mmodi,data.frame(nincome=c(0,1)),type="probs")) log(pp[1,1]*pp[2,2]/(pp[1,2]*pp[2,1])) log(pp[1,1]*pp[2,3]/(pp[1,3]*pp[2,1])) sPID[1:4] cm <- diag(3)[unclass(sPID),] cm[1:4,] y <- as.numeric(t(cm)) resp.factor <- gl(944,3) cat.factor <- gl(3,1,3*944,labels=c("D","I","R")) rnincome <- rep(nincome,each=3) head(data.frame(y,resp.factor,cat.factor,rnincome)) nullmod <- glm(y ~ resp.factor + cat.factor, family=poisson) glmod <- glm(y ~ resp.factor + cat.factor + cat.factor:rnincome, family=poisson) deviance(glmod) deviance(mmodi) coef(glmod)[c(1,945:949)] coef(mmodi) 0.016087-0.017665 data(cns) cns cns$CNS <- cns$An+cns$Sp+cns$Other plot(log(CNS/NoCNS) ~ Water, cns, pch=as.character(Work)) binmodw <- glm(cbind(CNS,NoCNS) ~ Water + Work, cns, family=binomial) binmoda <- glm(cbind(CNS,NoCNS) ~ Area + Work, cns, family=binomial) anova(binmodw,binmoda,test="Chi") halfnorm(residuals(binmodw)) summary(binmodw) exp(-0.339058) cmmod <- multinom(cbind(An,Sp,Other) ~ Water + Work, cns) nmod <- step(cmmod) nmod cc <- c(0,0.28963,-0.98083) names(cc) <- c("An","Sp","Other") exp(cc)/sum(exp(cc)) multinom(cbind(NoCNS,An,Sp,Other) ~ Water + Work, cns) library(MASS) pomod <- polr(sPID ~ age + educ + nincome, nes96) c(deviance(pomod),pomod$edf) c(deviance(mmod),mmod$edf) pomodi <- step(pomod) deviance(pomodi)-deviance(pomod) pchisq(11.151,pomod$edf-pomodi$edf,lower=F) pim <- prop.table(table(nincome,sPID),1) logit(pim[,1])-logit(pim[,1]+pim[,2]) summary(pomodi) ilogit(0.209) ilogit(1.292)-ilogit(0.209) predict(pomodi,data.frame(nincome=il,row.names=il),type="probs") 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) opmod <- polr(sPID ~ nincome, method="probit") summary(opmod) dems <- pnorm(0.128-il*0.008182) demind <- pnorm(0.798-il*0.008182) cbind(dems,demind-dems,1-demind) polr(sPID ~ nincome, method="cloglog")