###   Soc 6708
###   Advanced Statistics

### Robert Andersen, University of Toronto

###  Wednesday, May 23, 2012

#############################
#   Logit and Probit Models
#       Part I
#############################

library(car)
library(MASS)

library(foreign)
India<-read.spss('c:/data/india99.sav', , 
        use.value.labels=TRUE, to.data.frame=TRUE)
summary(India)
attach(India)

#Fitting a OLS model for vote for the BJP in India
#Using the GLM function
bjp.lm<-glm(BJPVOTE~LIBERAL, family=gaussian)
summary(bjp.lm)
summary(LIBERAL)
library(xtable)
print(xtable(bjp.lm))
detach(India)

#A model that gives wrong predictions
Convote<-read.table('c:/teaching/icpsr/data/convote.dat', header=T)
summary(Convote)
attach(Convote)
Convote.lm<-glm(CONVOTE~LRSCALE, family=gaussian)
summary(Convote.lm)
print(xtable(Convote.lm))
detach(Convote)


#Fitting a logit model for vote for the BJP in India
#use the glm function, sepcifying the binomial family
attach(India)
bjp.logit<-glm(BJPVOTE~CASTE, family=binomial)
summary(bjp.logit)
print(xtable(bjp.logit))


#Adding variables and interactions
bjp.logit2<-glm(BJPVOTE~MEN+CASTE+LIBERAL+CASTE*LIBERAL, family=binomial)
summary(bjp.logit2)
print(xtable(bjp.logit2))

#loading car library to perform analysis of deviance
#which is done with the Anova function,
library(car)
Anova(bjp.logit2)


#we can obtain the R^2 analogue by
#1-(G^2 for the model/G^2 for null model)
bjp.null<-update(bjp.logit2, .~1)


#Function for Pseudo R2
pseudoR2<-function(model) {
     1-(deviance(model)
     /model$null.deviance)}
pseudoR2(bjp.logit2)

#Effect displays
#The car library functions for effects displays
#rely on the lattice package, which by
#default gives a gray background to the graphs
#I prefer white so specify that as follows:
library(lattice)
trellis.device
lset(col.whitebg())
#car allows me to make two types of
#effect displays--one that gives 
#a separate graph for each effect, and 
#includes a confidence envelope
library(effects)
plot(all.effects(bjp.logit2))
plot(effect('CASTE*LIBERAL', bjp.logit2))

#The other plots all the effects 
#on a single graph but without 
#confidence envelopes
      
plot(effect('CASTE*LIBERAL', bjp.logit2), 
      multiline=TRUE)
#rescale.axis=FALSE evenly spaces the probilities
#on the vertical axis
plot(effect('CASTE*LIBERAL', bjp.logit2), 
      multiline=TRUE, rescale.axis=FALSE)

plot(all.effects(bjp.logit2), multiline=TRUE)



#Probit model

bjp.probit<-glm(BJPVOTE~MEN+CASTE+LIBERAL+CASTE*LIBERAL, 
     family=binomial(link = "probit"))
summary(bjp.probit)
Anova(bjp.probit)

plot(all.effects(bjp.probit), multiline=TRUE,
    rescale.axis=FALSE)


print(xtable(bjp.logit2))


#Comparing coefficients from logit 
#and probit models

probitmod<-glm(BJPVOTE~LIBERAL, 
     family=binomial(link = "probit"))
logitmod<-glm(BJPVOTE~LIBERAL, 
     family=binomial)
probitmod; logitmod
logitmod$coefficients/probitmod$coefficients