# Multinomial logit model example library(nnet) BEPS <- read.table("BEPS.txt") BEPS$vote <- factor(BEPS$vote, c("Liberal Democrat", "Labour", "Conservative")) multinom.mod <- multinom(vote ~ age + men + economic.cond.national + economic.cond.household + Blair + Hague + Kennedy + Europe*political.knowledge, data=BEPS) predictors <- data.frame(expand.grid(list( age=mean(BEPS$age), men=.5, economic.cond.national=mean(BEPS$economic.cond.national), economic.cond.household=mean(BEPS$economic.cond.household), Blair=mean(BEPS$Blair), Hague=mean(BEPS$Hague), Kennedy=mean(BEPS$Kennedy), Europe=seq(1:11), political.knowledge=0:3))) effects.multinom <- polytomousEffects(multinom.mod, predictors) attach(effects.multinom) par(mfrow=c(3,4), mar=c(5,5,4,2) + 0.1, cex.main=2, font.lab=par("font.main"), cex.axis=1.5, font.axis=1, cex.lab=par("cex.main")) for (party in c("Labour", "Conservative", "Liberal Democrat")){ for (knowledge in 0:3){ plot(c(1,11), c(0,1), type="n", xlab=if (party == "Liberal Democrat") "Attitude Towards Europe" else "", ylab=party, main=if (party == "Labour") paste("Knowledge =", knowledge) else "") for (prefix in c("p.", "L.", "U.")){ lines(1:11, get(paste(prefix, make.names(party), sep="") )[political.knowledge == knowledge], lty=if (prefix == "p.") 1 else 2, lwd=if (prefix == "p.") 3 else 1, col="red") } } } ### Proportional-odds model example library(MASS) WVS <- read.table("WVS.txt") WVS$poverty <- ordered(WVS$poverty, levels=c('Too Little', 'About Right', 'Too Much')) WVS$country <- factor(WVS$country, c('Sweden', 'Norway', 'Australia', 'USA')) polr.mod <- polr(poverty ~ men + religion + degree + country*poly(age,3), data=WVS) predictors <- data.frame(expand.grid(list( men=.5, religion=mean(WVS$religion), degree=mean(WVS$degree), age=18:87, country=c('Sweden', 'Norway', 'Australia', 'USA')))) effects.polr <- polytomousEffects(polr.mod, predictors) attach(effects.polr) par(mfrow=c(3,4), mar=c(5,5,4,2) + 0.1, cex.main=2, font.lab=par("font.main"), cex.axis=1.5, font.axis=1, cex.lab=par("cex.main")) for (response in c("Too Little", "About Right", "Too Much")){ for (ctry in c("Sweden", "Norway", "Australia", "USA")){ plot(c(18, 87), c(0,1), type="n", xlab=if(response == "Too Much") "Age" else "", ylab=if(ctry == "Sweden") response else "", main=if(response == "Too Little") ctry else "") for (prefix in c("p.", "L.", "U.")){ lines(18:87, get(paste(prefix, make.names(response), sep="") )[country == ctry], lty=if (prefix == "p.") 1 else 2, lwd=if (prefix == "p.") 3 else 1, col="red") } } }