rm(list = ls())
cush <- read.table(
url("http://stat.unm.edu/~fletcher/LLM/DATA/TAB4-11.DAT"),
#"C:\\E-drive\\Books\\ANREG2\\newdata\\TAB21-11.DAT",
#"C:\\E-drive\\Books\\LOGLIN3\\DATA\\TAB4-11.DAT",
		  sep="",col.names=c("Syn","Tetra","Preg"))
attach(cush)
cush

#Create a 3 x 21 table of 0-1 entries,
#each row has 1's for a different type of syndrome
j=rep(seq(1,21),3)
i=c(rep(1,21),rep(2,21),rep(3,21))
Tet=c(Tetra,Tetra,Tetra)
Pre=c(Preg,Preg,Preg)
y=c(Syn,Syn,Syn)
y[1:6]=1
y[7:21]=0
y[22:27]=0
y[28:37]=1
y[38:58]=0
y[59:63]=1
datal=c(y,i,j,Tet,Pre)
datl=matrix(datal,63,5,dimnames =
list(NULL,c("y", "i", "j","Tet","Pre")))
datl

#Fit the log-linear model for logistic discrimination.
i=factor(i)
j=factor(j)
lp=log(Pre)
lt=log(Tet)
ld <- glm(y ~ i + j + i:lt +i:lp ,family = poisson)
ldp=summary(ld)
ldp
anova(ld)

# Table 4.12
q=ld$fit
# Divide by sample sizes
p1=ld$fit[1:21]/6
p2=ld$fit[22:42]/10
p3=ld$fit[43:63]/5
# Produce table
estprob = c(Syn,p1,p2,p3)
EstProb=matrix(estprob,21,4,dimnames =
list(NULL,c("Group", "A", "B","C")))
EstProb

# Table 4.13 Proportional prior probabilities.
post = c(Syn,ld$fit)
PropProb=matrix(post,21,4,dimnames =
list(NULL,c("Group", "A", "B","C")))
PropProb

# Table 4.13 Equal prior probabilities.
p=p1+p2+p3
pp1=p1/p
pp2=p2/p
pp3=p3/p
post = c(Syn,pp1,pp2,pp3)
EqProb=matrix(post,21,4,dimnames=
list(NULL,c("Group","A","B","C")))
EqProb