library(MASS)
library(mlogit)



satisliving<-data.frame(expand.grid(contact = c("Low","High"),  satisfaction = c("Low","Medium","High"), 
                                    housing =c("Tower blocks","Apartments","Atrium houses","Terraced houses"), influence = c("Low","Medium","High")), counts =  c(21,14,21,19,28,37,61,78,23,46,17,43,13,20,9,23,10,20,18,57,6,23,7,13,34,17,22,23,36,40,43,48,35,45,40,86,8,10,8,22,12,24,15,31,13,21,13,13,10,3,11,5,36,23,26,15,18,25,54,62,6,7,7,10,9,21,7,5,5,6,11,13))

satisliving<-satisliving[order(satisliving$contact,satisliving$housing, satisliving$influence, satisliving$satisfaction),]



# Fitting poisson Generalized Linear Model
fit.glm <- glm(counts ~ contact*satisfaction*housing*influence,
                family="poisson",data=satisliving)


# Choose a model by AIC in a Stepwise Algorithm
fit.glm.step <- stepAIC(fit.glm)
summary(fit.glm.step)


#probability table
satisliving.names <- lapply(satisliving[, -5], levels)
x <- predict(fit.glm.step, type = "response")
satis.pm <- matrix(x, ncol = 3, byrow = T,dimnames = list(NULL, satisliving.names[[2]]))
satis.pr <- satis.pm/drop(satis.pm%*% rep(1, 3))
glm.prob.table <- cbind(expand.grid(satisliving.names[-2]), prob = round(satis.pr, 2))
glm.prob.table