library(car) Chile$vote1 <- Chile$vote Chile$vote1[Chile$vote1 == "A"] <- NA Chile$vote1 <- ordered(Chile$vote1) Chile$yes <- factor(with(Chile, ifelse(vote=="Y", "yes", ifelse(vote=="N", "no", NA)))) Chile$education <- factor(Chile$education, levels=c("P", "S", "PS")) modb <- glm(yes ~ region + population + sex + age + education + income, family=binomial, data=Chile) Anova(modb) modb2 <- glm(yes ~ statusquo + region + population + sex + age + education + income, family=binomial, data=Chile) Anova(modb2) library(effects) plot(allEffects(modb2), ask=FALSE) library(MASS) modp <- polr(vote1 ~ region + sex + age + education + log(income) + log(population), data=Chile, Hess=TRUE) summary(modp) Anova(modp) modp2 <- polr(vote1 ~ statusquo + region + sex + age + education + log(income) + log(population), data=Chile, Hess=TRUE) summary(modp2) Anova(modp2) plot(effect("statusquo", modp2), style="stacked") eff <- allEffects(modp) plot(eff, ask=FALSE) plot(allEffects(modp), style="stacked", ask=FALSE) library(nnet) modm <- multinom(vote1 ~ region + sex + age + education + log(income) + log(population), data=Chile) summary(modm) Anova(modm) plot(allEffects(modm), style="stacked", ask=FALSE) logLik(modp) logLik(modm) AIC(modp, modm) AIC(modp, modm, k=log(2260)) modm2 <- multinom(vote1 ~ statusquo + region + sex + age + education + income + population, data=Chile) Anova(modm2) AIC(modm2, modp2, k=log(2260)) plot(allEffects(modm2), ask=FALSE, style="stacked") # diagnostics influencePlot(modb) dfbs <- dfbetas(modb) nms <- colnames(dfbs) par(mfrow=c(3,4)) for (i in 1:11) plot(dfbs[,i], ylab=nms[i]) par(mfrow=c(2,2)) for (var in c("population", "age", "income")) cr.plot(modb, var) table(Chile$population) table(Chile$income) mod.nonlin <- glm(yes ~ region + as.factor(population) + sex + age + education + as.factor(income), family=binomial, data=Chile) mod.poplin <- glm(yes ~ region + population + sex + age + education + as.factor(income), family=binomial, data=Chile) anova(mod.poplin, mod.nonlin, test="Chisq") mod.inclin <- glm(yes ~ region + as.factor(population) + sex + age + education + income, family=binomial, data=Chile) anova(mod.inclin, mod.nonlin, test="Chisq")