##---------------------------------------------## ## Script for Part 3: Logit and Probit Models ## ## John Fox ## ## Applied Statistics With R ## ## WU Wien ## ## May/June 2006 ## ##---------------------------------------------## # Logit, Probit, and Linear Probability Models, Chile Data library(car) data(Chile) attach(Chile) yes <- ifelse(vote == "Y", 1, ifelse(vote=="N", 0, NA)) chile.logit <- glm(yes ~ statusquo, family=binomial) # logit summary(chile.logit) chile.probit <- glm(yes ~ statusquo, family=binomial(probit)) # probit summary(chile.probit) chile.linear <- lm(yes ~ statusquo) # linear summary(chile.linear) coef(chile.linear)[2] # compare slope coefficients (coef(chile.logit)/4)[2] coef(chile.probit)[2] (coef(chile.logit)/(pi/sqrt(3)))[2] detach(Chile) # Binary Logit Model, SLID data SLID <- read.table(file.choose(), header=TRUE) attach(SLID) some(SLID) region <- factor(region, levels=c("Atlantic", "Quebec", "Ontario", "Prairies", "BC")) kids0004 <- recode(kids0004, ' 0="No"; 1="Yes" ', as.factor=TRUE) kids0509 <- recode(kids0509, ' 0="No"; 1="Yes" ', as.factor=TRUE) income <- husbandIncome/1000 slid.mod.1 <- glm(working ~ region + (kids0004 + kids0509 + kids1014)* (income + schooling), family=binomial) Anova(slid.mod.1) slid.mod.2 <- glm(working ~ region + kids0004 + kids0509 + income + schooling, family=binomial) summary(slid.mod.2) library(effects) plot(all.effects(slid.mod.2), ask=FALSE) detach(SLID) # Multinomial Logit Model, British Election Panel Study library(nnet) BEPS <- read.table(file.choose()) BEPS$vote <- factor(BEPS$vote, c("Liberal Democrat", "Labour", "Conservative")) some(BEPS) BEPS.mod <- multinom(vote ~ age + men + economic.cond.national + economic.cond.household + Blair + Hague + Kennedy + Europe*political.knowledge, data=BEPS) summary(BEPS.mod, cor=FALSE) # suppress coefficient correlations Anova(BEPS.mod) detach(BEPS) # Proportional Odds Logit Model, World Values Survey library(MASS) WVS <- read.table(file.choose()) WVS$poverty <- ordered(WVS$poverty, levels=c('Too Little', 'About Right', 'Too Much')) some(WVS) WVS.mod.1 <- polr(poverty ~ country*(men + religion + degree + age), data=WVS) Anova(WVS.mod.1) WVS.mod.2 <- polr(poverty ~ men + country*(religion + degree + age), data=WVS) summary(WVS.mod.2) WVS.multinom.2 <- multinom(poverty ~ men + country*(religion + degree + age), data=WVS) # approximate test of proportional-odds assumption deviance(WVS.mod.2) - deviance(WVS.multinom.2) pchisq(deviance(WVS.mod.2) - deviance(WVS.multinom.2), 34 - 18, lower.tail=FALSE) # Binomial Logit Model, American Voter Data closeness <- factor(rep(c('one.sided', 'close'), c(3,3)), levels=c('one.sided', 'close')) preference <- ordered(rep(c('weak', 'medium', 'strong'), 2), levels=c('weak', 'medium', 'strong')) voted <- c(91, 121, 64, 214, 284, 201) did.not.vote <- c(39, 49, 24, 87, 76, 25) data.frame(closeness, preference, logit=log(voted/did.not.vote)) options(contrasts=c("contr.sum", "contr.poly")) mod.campbell <- glm(cbind(voted, did.not.vote) ~ closeness*preference, family=binomial) Anova(mod.campbell) Anova(mod.campbell, type="III")