# loglin_ex_corr.r

# Chargement des fonctions et packages

library(MASS)

##############
# Exercice 3.1

load("Titanic.rda")
Ti <- Titanic
summary(Ti)

# Création des deux formats de données nécessaires
# Table pour utiliser la fonction loglm
T <- xtabs(~Ti$classe+Ti$age+Ti$sex+Ti$living)
T

# Data Frame pour utiliser la fonction glm
DF <- as.data.frame(T)
DF

# Modèles avec toutes les interactions jusqu'à un certain ordre
M_O3 <- glm(Freq~(Ti.classe+Ti.age+Ti.sex+Ti.living)^3,data=DF,family=poisson)
summary(M_O3)

M_O2 <- glm(Freq~(Ti.classe+Ti.age+Ti.sex+Ti.living)^2,data=DF,family=poisson)
summary(M_O2)
1-pchisq(116.59,13)
# p-valeur = 0 -> modèle rejeté

# 1) Procédure de sélection backward à partir du modèle M_O3 en se basant sur la déviance
anova(M_O3)
1-pchisq(4.02,3)
1-pchisq(1.69,1)
# -> on supprime Ti.classe:Ti.age:Ti.sex

M_O3a <- glm(Freq~(Ti.classe+Ti.age+Ti.sex+Ti.living)^3
			-Ti.classe:Ti.age:Ti.sex,data=DF,family=poisson)
summary(M_O3a)
1-pchisq(9.7834,6)
anova(M_O3a)
1-pchisq(29.21,3)
1-pchisq(12.17,1)
# Aucune de ces deux interactions ne peut être supprimée.
# Le modèle final est le modèle M_O3a.

# 2) Idée: se baser sur AIC.
M_AIC <- step(M_O3,direction="backward")
# -> on arrive à un second modèle final contenant aussi 3 interactions d'ordre 3

# 3) Idée: se baser sur BIC.
M_BIC <- step(M_O3,direction="backward",k=log(2201))
# -> on arrive à un modèle final contenant 3 interactions d'ordre 3
# ce modèle est le même que celui obtenu précédemment en se basant sur la déviance


##############
# Exercice 3.2

load("Syndicalisation.rda")
Sy <- Syndicalisation
summary(Sy)

# Création des deux formats de données nécessaires
# Table pour utiliser la fonction loglm
T <- xtabs(~Sy$sexe+Sy$syndiqué+Sy$travail+Sy$état_civil+Sy$expérience)
T

# Data Frame pour utiliser la fonction glm
DF <- as.data.frame(T)
DF

# Modèle saturé d'ordre 5
M_O5 <- glm(Freq ~ (Sy.sexe+Sy.syndiqué+Sy.travail+Sy.état_civil+Sy.expérience)^5,
			data=DF,family=poisson)

# 1) Meilleur modèle en terme de BIC
M_BIC <- step(M_O5,direction="backward",k=log(534))

# 2) Meilleur modèle en terme de AIC
M_AIC <- step(M_O5,direction="backward")

# 3) Meilleur modèle en terme de déviance

# Modèle d'ordre 4
M_O4 <- glm(Freq ~ (Sy.sexe+Sy.syndiqué+Sy.travail+Sy.état_civil+Sy.expérience)^4,
			data=DF,family=poisson)
summary(M_O4)
1-pchisq(17.021,10)

# Modèle d'ordre 3
M_O3 <- glm(Freq ~ (Sy.sexe+Sy.syndiqué+Sy.travail+Sy.état_civil+Sy.expérience)^3,
			data=DF,family=poisson)
summary(M_O3)
1-pchisq(49.314,47)

# Modèle d'ordre 2
M_O2 <- glm(Freq ~ (Sy.sexe+Sy.syndiqué+Sy.travail+Sy.état_civil+Sy.expérience)^2,
			data=DF,family=poisson)
summary(M_O2)
1-pchisq(107.96,99)

# Modèle d'ordre 1 (modèle d'indépendance)
M_O1 <- glm(Freq ~ (Sy.sexe+Sy.syndiqué+Sy.travail+Sy.état_civil+Sy.expérience),
			data=DF,family=poisson)
summary(M_O1)
1-pchisq(331.23,133)

# Procédure de sélection backward manuelle à partir de M_O3
anova(M_O3)
# Etape 1
M_O3_1 <- glm(Freq ~ (Sy.sexe+Sy.syndiqué+Sy.travail+Sy.état_civil+Sy.expérience)^3
			-Sy.sexe:Sy.syndiqué:Sy.expérience,
			data=DF,family=poisson)
anova(M_O3_1)

# Etape 2
M_O3_2 <- glm(Freq ~ (Sy.sexe+Sy.syndiqué+Sy.travail+Sy.état_civil+Sy.expérience)^3
			-Sy.sexe:Sy.syndiqué:Sy.expérience-Sy.syndiqué:Sy.état_civil:Sy.expérience,
			data=DF,family=poisson)
anova(M_O3_2)
summary(M_O3_2)
1-pchisq(50.551,51)

# Etc ...


##############
# Exercice 3.3

load("PoidsBebe.rda")
PB <- PoidsBebe
summary(PB)

# Création des deux formats de données nécessaires
# Table pour utiliser la fonction loglm
T <- xtabs(~PB$LOW+PB$SMOKE+PB$HT+PB$UI)

# Data Frame pour utiliser la fonction glm
DF <- as.data.frame(T)
DF

# Modèle avec interactions jusqu'à l'ordre 4
M_O4 <- glm(Freq ~ (PB.LOW+PB.SMOKE+PB.HT+PB.UI)^4,data=DF,family=poisson)
summary(M_O4)

# Modèle avec interactions jusqu'à l'ordre 3
M_O3 <- glm(Freq ~ (PB.LOW+PB.SMOKE+PB.HT+PB.UI)^3,data=DF,family=poisson)
summary(M_O3)

# Modèle avec interactions jusqu'à l'ordre 2
M_O2 <- glm(Freq ~ (PB.LOW+PB.SMOKE+PB.HT+PB.UI)^2,data=DF,family=poisson)
summary(M_O2)
1-pchisq(0.6587,5)

# Modèle avec interactions jusqu'à l'ordre 1
M_O1 <- glm(Freq ~ (PB.LOW+PB.SMOKE+PB.HT+PB.UI),data=DF,family=poisson)
summary(M_O1)
1-pchisq(20.110,11)

# --> Modèle le plus simple acceptable: M_O2

# Meilleur modèle sur le principe du BIC
M_BIC <- step(M_O4,direction="backward",k=log(189))

# --> Il y a 2 "meilleurs modèles" possibles:
# 1) Le modèle d'indépendance (BIC minimal)
# 2) Le modèle 
# Freq ~ PB.LOW + PB.SMOKE + PB.HT + PB.UI + PB.LOW:PB.HT + PB.LOW:PB.UI + PB.HT:PB.UI
# Meilleur BIC avant que les améliorations deviennent marginales (< 2)


# Sélection en termes de déviance en prenant M_O2 comme point de départ
anova(M_O2)

M_O2a <- glm(Freq ~ (PB.LOW+PB.SMOKE+PB.HT+PB.UI)^2
		-PB.SMOKE:PB.HT,data=DF,family=poisson)
anova(M_O2a)

M_O2b <- glm(Freq ~ (PB.LOW+PB.SMOKE+PB.HT+PB.UI)^2
		-PB.SMOKE:PB.HT-PB.SMOKE:PB.UI,data=DF,family=poisson)
anova(M_O2b)

# --> En terme de déviance, c'est ce modèle qui est le meilleur (différent du modèle M_BIC)