ex.data<-data.frame(expand.grid(
personality=factor(c("A","B"),levels=c("A","B")),
cholesterol=factor(c("Normal","High"),levels=c("Normal","High")),
dbp=factor(c("Normal","High"),levels=c("Normal","High"))),
count=c(716,819,207,186,79,67,25,22))


library(MASS)
#short hand writing
# 1*2 is equivalent to write 1+2+1*2
# 1:2 only include the interaction between 1 and 2
# 1*2*3  is equivalent to write 1+ 2+ 3+1*2 + 1*3 + 2*3 + 1*2*3
# (1+2+3)^2 is equivalent to write 1+2+3+1*2+1*3+2*3

#example 3.2.1
fitp.c.d<-loglm(count~personality+cholesterol+dbp,data=ex.data,param=T,fit=T) #[1][2][3] 
fitp.c.d
fitp.c.d$fit

#example 3.2.3
fitpc.pd<-loglm(count~personality*cholesterol+personality*dbp,data=ex.data,param=T,fit=T) #[12][13] 
fitpcpd.array<-fitted(fitpc.pd)
odds.ratio<-function(x) x[1,1]*x[2,2]/(x[2,1]*x[1,2])
apply(fitpcpd.array,1,odds.ratio) #given level of personality
apply(fitpcpd.array,2,odds.ratio) #given level of cholesterol
apply(fitpcpd.array,3,odds.ratio) #given level of dbp

#example 3.4.2
fitpcd<-loglm(count~dbp*cholesterol*personality,data=ex.data,param=T,fit=T) # pcd
fitpc.pd.cd<-update(fitpcd, .~. - dbp:cholesterol:personality) # pc, pd, cd
fitpc.pd<-update(fitpc.pd.cd, .~. - dbp:cholesterol) #pc, pd
fitpc.cd<-update(fitpc.pd.cd, .~. - dbp:personality) #pc, cd
fitpd.cd<-update(fitpc.pd.cd, .~. - cholesterol:personality) #pd, cd
fitp.cd<-update(fitpc.pd.cd, .~. - cholesterol:personality -dbp:personality)#p, cd
fitc.pd<-update(fitpc.pd.cd, .~. - cholesterol:personality -dbp:cholesterol)#c,pd
fitd.pc<-update(fitpc.pd.cd, .~. - dbp:personality -dbp:cholesterol)#c,pd
fitp.c.d<-update(fitd.pc, .~. - cholesterol:personality) #p, c, d

#example 3.6.1, model selection, 
#do forwards selection or backwards selection or both;
#start from a model of our choice;
#use AIC (the default) or BIC
upper<-formula(~dbp*cholesterol+cholesterol*personality+dbp*personality)
model.aic = step(fitp.c.d, scope=list(lower= ~., upper= ~.^2))
model.aic = step(fitp.c.d, scope=list(lower= ~., upper=upper))
model.aic = step(fitpcd, direction="backward")