##---------------------------------------------------------------------##
## Script for Soc 740  Lecture 9                                       ##
##    John Fox                                                         ##
##    Winter 2014                                                      ##
##---------------------------------------------------------------------##


# Chilean plebiscite data

library(car)

?Chile

Chile$yes <- with(Chile, ifelse(vote == "Y", 1, ifelse(vote=="N", 0, NA)))
Chile <- na.omit(Chile[,c("yes", "statusquo")])
table(Chile$yes)


with (Chile,
    plot(statusquo, jitter(yes, factor=1/4), xlab="Support for the Status Quo", 
    ylab="Voting Intention"))
abline(h=c(0,1), lty=2)

mod.lin <- lm(yes ~ statusquo, data=Chile)
summary(mod.lin)

abline(mod.lin, lwd=2) # linear probability model by OLS

mod <- glm(yes ~ statusquo, family=binomial, data=Chile) # logit model by ML
summary(mod)

sq <- with(Chile, seq(min(statusquo), max(statusquo), length=1000))
fit <- predict(mod, newdata=data.frame(statusquo=sq), type="response")
lines(sq, fit, lwd=2)

  # lowess fit using non-robust option (iter=0)
with(Chile, lines(lowess(statusquo, yes, f=0.25, iter=0), lwd=2, lty=2))

# logistic regression for SLID data

SLID <- read.table(
"http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/SLID-women.txt", 
    header=TRUE) 

SLID$region <- factor(SLID$region,       # re-order levels
    levels=c("Atlantic", "Quebec", "Ontario", "Prairies", "BC"))

mod.slid.1 <- glm(working ~ region + kids0004 + kids0509 + kids1014 
    + familyIncome + education, family=binomial, data=SLID)
    
Anova(mod.slid.1)

mod.slid.2 <- glm(working ~ region + kids0004 + kids0509 + familyIncome 
    + education, family=binomial, data=SLID)
summary(mod.slid.2)

confint(mod.slid.2) # LR-based intervals

confint.default(mod.slid.2) # Wald-based intervals

    # Effect plots for the logit model
    
library(effects)

quantile(SLID$education, c(.1, .9))
quantile(SLID$familyIncome, c(.1, .9))

eff <- allEffects(mod.slid.2, xlevels=list(education=10:17, 
    familyIncome=seq(11.6394, 44.3394, length=10)))
plot(eff, ylim=logit(c(0.65, 0.92))) 
# note: forces common y-axis to compare effects
#  and constrains range of income and education

plot(allEffects(mod.slid.2))  # default plot