### 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