library(HSAUR)

# VARIOUS DATASETS FOR BINARY DATA (LOGISTIC REGRESSION)
?plasma;                 
?womensrole
?polyps
#
#############
# PLASMA DATA
#############
#
data("plasma",package="HSAUR");           head(plasma)
layout(matrix(1:2,ncol=2))
  cdplot(ESR~fibrinogen,data=plasma)
  cdplot(ESR~globulin,data=plasma)
layout(matrix(1:1,ncol=1))
#
plot(plasma$fibrinogen,plasma$ESR)
layout(matrix(1:2,ncol=2))
  plot(plasma$ESR,plasma$fibrinogen,main="Fib")
  plot(plasma$ESR,plasma$globulin,main="Glo")
layout(matrix(1:1,ncol=1))
#
plasma1 = glm( ESR ~ fibrinogen, data=plasma, family=binomial())
summary(plasma1)
confint(plasma1,parm="fibrinogen")			    # Based on profile likelihood (Venables and Ripley's MASS book, ?8.4, pp. 220-221)
confint.default(plasma1,parm="fibrinogen")	# Based on approx normality
exp(coef(plasma1)["fibrinogen"])
exp(confint.default(plasma1,parm="fibrinogen"))
#
plasma2 = glm(ESR~fibrinogen+globulin, data=plasma, family=binomial())
summary(plasma2)
confint.default(plasma2,parm="globulin")
#
# RESIDUAL DEVIANCE
plasma1["deviance"] 
plasma2["deviance"]
anova(plasma1,plasma2,test="Chisq")
#
prob=predict(plasma1,type="response")
prob
#
# BUBBLE PLOT
plot(globulin~fibrinogen,data=plasma,xlim=c(2,6),ylim=c(25,50),pch=".")
symbols(plasma$fibrinogen,plasma$globulin,circles=prob,add=TRUE)
#
###################
# WOMENSROLE DATA
###################
# (aggregated data)
#
data("womensrole",package="HSAUR")
head(womensrole)
#
women1 = glm( cbind(agree,disagree) ~ sex + education, data=womensrole, family=binomial())
summary(women1)
#
role1=predict(women1,type="response")
role1
#
confint.default(women1)
exp(coef(women1))
exp(confint.default(women1))
#
myplot = function(role.fitted){
	f = womensrole$sex == "Female"
	plot(womensrole$education,role.fitted,type="n",ylab="Prob of ageeing", xlab="Education",ylim=c(0,1))
	lines(womensrole$educatio[!f],role.fitted[!f],lty=1)
	lines(womensrole$educatio[f],role.fitted[f],lty=2)
	lgtxt=c("Fitted (Males)", "Fitted (Females)")
	legend("topright",lgtxt,lty=1:2,bty="n")
	y = womensrole$agree/(womensrole$agree+womensrole$disagree)
	text(womensrole$education,y,ifelse(f,"F","M"),family="HerseySerif",cex=1.25)
	}
#
women2 = glm( cbind(agree,disagree) ~ sex * education, data=womensrole, family=binomial())
summary(women2)
role2=predict(women2,type="response")
role2
#
myplot(role1)
myplot(role2)
#
res=residuals(women2,type="deviance")
plot(role2,res,xlab="Fitted Values",ylab="Residuals")
abline(h=0,lty=2)
#
# Probit
#
women3 = glm( cbind(agree,disagree) ~ sex * education, data=womensrole, family=binomial(link=probit))
summary(women3)
role3=predict(women3,type="response")
role3
#
# Complementary log-log
#
women4 = glm( cbind(agree,disagree) ~ sex * education, data=womensrole, family=binomial(link=cloglog))
summary(women4)
role4=predict(women4,type="response")
role4
#
f = womensrole$sex == "Female"
plot(womensrole$education,role2,type="n",ylab="Prob of ageeing", xlab="Education",ylim=c(0,1))
lines(womensrole$educatio[!f],role2[!f],lty=1)
lines(womensrole$educatio[!f],role3[!f],lty=1,col=2)
lines(womensrole$educatio[!f],role4[!f],lty=1,col=4)
legend("topright",legend=c("logit","probit","cloglog"),lty=1,col=c(1,2,4))