###############  Handout #5 for ST434/534  ###############
##########################################################
################Logit models#####################
##########################################################
#################Example Abortion Data######################
library(MASS)

#write abortion data from Table 3.1 
race.labs<-factor(c("white","nonwhite"),levels=c("white","nonwhite"))
sex.labs<-factor(c("male","female"),levels=c("male","female"))
opinion.labs<-factor(c("yes","no","und"),levels=c("yes","no","und"))
age.labs<-factor(c("18-25","26-35","36-45","46-55","56-65","66+"),levels=
c("18-25","26-35","36-45","46-55","56-65","66+"))

temp<-c(96,24,140,21,44,5,43,4,1,2,1,1,138,18,171,25,64,7,65,6,2,1,4,2,
117,16,152,20,56,7,58,5,6,3,9,1,75,12,101,17,48,6,51,5,5,4,9,1,72,6,102,14,49,8,
58,5,6,3,10,1,83,4,111,13,60,10,67,5,8,4,16,1)

ex.data<-data.frame(expand.grid(race=race.labs,sex=sex.labs,opinion=opinion.labs,age=age.labs),count=temp)
#fit model [RSA][RSO][OA]
fit1<-loglm(count~race*sex*age+race*sex*opinion+opinion*age,data=ex.data,param=T,fit=T) 
#another way
fit2 <- glm(count~race*sex*age+race*sex*opinion+opinion*age,data=ex.data, family=poisson)
fit1.array<-fitted(fit1) 
#generate table 4.8
aa<-ftable(fit1.array, col.vars = "age") 
aa
aa[1,]/aa[2,] #row 1 of table 4.9
aa[4,]/aa[5,] #row 2 of table 4.9
aa[7,]/aa[8,] #row 3 of table 4.9
aa[10,]/aa[11,] #row 4 of table 4.9

#Table4.9
temp2<-c(aa[1,]/aa[2,],aa[4,]/aa[5,],aa[7,]/aa[8,],aa[10,]/aa[11,])
tabletemp<-array(temp2,dim=c(6,2,2))
dimnames(tabletemp)<-list(age=c("18-25","26-35","36-45","46-55","56-65","66+"),
sex=c("male","female"),race=c("white","nonwhite"))
table4.9<-ftable(tabletemp,row.vars=c("race","sex"),col.vars = "age")
table4.9

#now delete the category undecided
d<-which(ex.data$opinion=='und')
newdata<-ex.data[-d,]
newdata

#fit model [RSA][RSO][OA] using data with undecided deleted
fit3<-loglm(count~race*sex*age+race*sex*opinion+opinion*age,data=newdata,param=T,fit=T)
fit3.array<-fitted(fit3) 
ftable(fit3.array, col.vars = "age") #table 4.10
anova(fit3)
fit3$deviance #9.104
#another way
fit4 <- glm(count~race*sex*age+race*sex*opinion+opinion*age,data=newdata, family=poisson)


#use logit model directly, I create the dataset "abortion.txt"
ex.data2<-read.table(file="W:/teaching/stat434534/datasets/abortion.txt",header=T,sep="",col.names=c("ID","race","sex","age","yes","no"))

#model 5
fit5<-glm(yes/(yes+no)~race*sex+age,data=ex.data2,family=binomial,weight=yes+no) 
#note, the left side of glm function is yes/(yes+no), this the probability of support
fit5

#recode age to incorporate a linear trend
ex.data2$k<-as.numeric(ex.data2$age)
#model 6
fit6<-glm(yes/(yes+no)~race*sex+k,data=ex.data2,family=binomial,weight=yes+no)
fit6
anova(fit6,fit5)


#recode g as (f,e) to incorporate the fact that males have similar odds of support
for(i in 1:24){
if (ex.data2$race[i]=='white'& ex.data2$sex[i]=='male') {ex.data2$f[i]<-1} 
if (ex.data2$race[i]=='white'& ex.data2$sex[i]=="female") {ex.data2$f[i]<-2} 
 if (ex.data2$race[i]=='nonwhite'& ex.data2$sex[i]=="male") {ex.data2$f[i]<-1}
 if (ex.data2$race[i]=='nonwhite'& ex.data2$sex[i]=="female") {ex.data2$f[i]<-3} 
 }

for(i in 1:24){
if (ex.data2$race[i]=='white'& ex.data2$sex[i]=='male') {ex.data2$e[i]<-1} 
if (ex.data2$race[i]=='white'& ex.data2$sex[i]=="female") {ex.data2$e[i]<-1} 
 if (ex.data2$race[i]=='nonwhite'& ex.data2$sex[i]=="male") {ex.data2$e[i]<-2}
 if (ex.data2$race[i]=='nonwhite'& ex.data2$sex[i]=="female") {ex.data2$e[i]<-1} 
 }


#model 7 

fit7<-glm(yes/(yes+no)~as.factor(f)+age,data=ex.data2,family=binomial,weight=yes+no)
fit7

#model 8
fit8<-glm(yes/(yes+no)~as.factor(f)+k,data=ex.data2,family=binomial,weight=yes+no)
fit8
summary(fit8)
fit8.array<-fitted(fit8) 
 
tableprob<-array(fit8.array,dim=c(4,6))
dimnames(tableprob)<-list(racesex=c("wm","nm","wf","nf"),age=c("18-25","26-35","36-45","46-55","56-65","66+"))
tableprob #probability of support table 

tablem8<-tableprob/(1-tableprob)
#no.l5, white female, 46-55, fitted value=0.677.
dimnames(tablem8)<-list(racesex=c("wm","nm","wf","nf"),age=c("18-25","26-35","36-45","46-55","56-65","66+"))
tablem8 #odds ratio table in page 158