# Loglin.r

# Chargement des fonctions et packages

library(descr)
library(MASS)

# Exemple pour une table de contingence à deux dimensions: Infections
load("Infections.rda")

summary(Infections)
crosstab(Infections$ocean,Infections$r_infect,plot=FALSE)

crosstab(Infections$ocean,Infections$r_infect,expected=TRUE,plot=FALSE)

# Transformation du format des données
Itable <- xtabs(~Infections$ocean+Infections$r_infect)
Itable
Idt <- as.data.frame(Itable)
Idt

# Modèle d'indépendance
Indep1 <- glm(Freq~relevel(Infections.ocean,2)+relevel(Infections.r_infect,3),data=Idt,family=poisson)
summary(Indep1)

Indep2 <- glm(Freq~Infections.ocean+Infections.r_infect,data=Idt,family=poisson)
summary(Indep2)

Indep3 <- loglm(~Infections$ocean+Infections$r_infect,data=Itable)
Indep3

# Modèle saturé

Sat <- glm(Freq~relevel(Infections.ocean,2)*relevel(Infections.r_infect,3),data=Idt,family=poisson)
summary(Sat)


# Exemple pour une table de contingence à trois dimensions: Infections
table(Infections$ocean,Infections$r_infect,Infections$sexe)

Itable3 <- xtabs(~Infections$sexe+Infections$ocean+Infections$r_infect)
Idt3 <- as.data.frame(Itable3)

# Modèle d'association homogène
I3_ass_hom <- glm(Freq~Infections.sexe*Infections.ocean+Infections.sexe*Infections.r_infect+Infections.ocean*Infections.r_infect,data=Idt3,family=poisson)
summary(I3_ass_hom)
anova(I3_ass_hom)

I3_ass_hom_2 <- loglm(~Infections$sexe*Infections$ocean+Infections$sexe*Infections$r_infect+Infections$ocean*Infections$r_infect,data=Itable3)
I3_ass_hom_2

# Modèle d'indépendance conditionnelle
I3_ind_cond <- glm(Freq~Infections.sexe*Infections.ocean+Infections.ocean*Infections.r_infect ,data=Idt3,family=poisson)
summary(I3_ind_cond)

I3_ind_cond_2 <- loglm(~Infections$sexe*Infections$ocean+Infections$ocean*Infections$r_infect,data=Itable3)
I3_ind_cond_2

anova(I3_ind_cond_2,I3_ass_hom_2)


##################################################################
# Procédure de sélection backward basée sur le test de la déviance

# Modèle d'association homogène (toutes les interactions d'ordre 2)
I3_ass_hom_2

# Modèle d'indépendance
I3_ind_2 <- loglm(~Infections$sexe+Infections$ocean+Infections$r_infect,data=Itable3)
I3_ind_2

# Composantes du modèle d'association homogène
anova(I3_ass_hom)

1-pchisq(0.109,1)
1-pchisq(0.239,2)

# Modèle d'indépendance conditionnelle entre sexe et r_infect:
anova(I3_ind_cond)

1-pchisq(0.109,1)
1-pchisq(8.698,2)

# Modèle avec une seule interaction d'ordre 2 entre r_infect et ocean
I3_final <- glm(Freq~Infections.ocean*Infections.r_infect+Infections.sexe,data=Idt3,family=poisson)
anova(I3_final)

1-pchisq(8.698,2)


# Procédure backward partant du modèle saturé et fondée sur AIC 
I3_Sat <- glm(Freq~Infections.ocean*Infections.r_infect*Infections.sexe,data=Idt3,family=poisson)
I3_AIC <- step(I3_Sat,direction="backward")
summary(I3_AIC)


# Procédure backward partant du modèle saturé et fondée sur BIC 
n <- sum(Idt3$Freq)
I3_BIC <- step(I3_Sat,direction="backward",k=log(n))
summary(I3_BIC)


###############################################
# Exemple GPEM
load("GPEM.rda")

summary(GPEM)

# Transformation du format des données
T_GPEM <- xtabs(~GPEM$GENDER+GPEM$PMS+GPEM$EMS+GPEM$M.STATUS)
T_GPEM
DF_GPEM <- as.data.frame(T_GPEM)
DF_GPEM

# Modèle saturé
GPEM_Sat <- loglm(~(GPEM$GENDER+GPEM$PMS+GPEM$EMS+GPEM$M.STATUS)^4,data=T_GPEM)
GPEM_Sat

# Modèle jusqu'aux interactions d'ordre 3
GPEM_O3 <- loglm(~(GPEM$GENDER+GPEM$PMS+GPEM$EMS+GPEM$M.STATUS)^3,data=T_GPEM)
GPEM_O3

# Modèle jusqu'aux interactions d'ordre 2
GPEM_O2 <- loglm(~(GPEM$GENDER+GPEM$PMS+GPEM$EMS+GPEM$M.STATUS)^2,data=T_GPEM)
GPEM_O2

# Modèle d'indépendance
GPEM_indep <- loglm(~GPEM$GENDER+GPEM$PMS+GPEM$EMS+GPEM$M.STATUS,data=T_GPEM)
GPEM_indep

# Sélection backward en partant du modèle GPEM_O3
GPEM_O3_glm <- glm(Freq~(GPEM.GENDER+GPEM.PMS+GPEM.EMS+GPEM.M.STATUS)^3,data=DF_GPEM,family=poisson)
anova(GPEM_O3_glm)

GPEM_O3_1 <- glm(Freq~(GPEM.GENDER+GPEM.PMS+GPEM.EMS+GPEM.M.STATUS)^3-GPEM.GENDER:GPEM.PMS:GPEM.EMS,data=DF_GPEM,family=poisson)
summary(GPEM_O3_1)
1-pchisq(GPEM_O3_1$deviance,GPEM_O3_1$df.residual)
anova(GPEM_O3_1)

GPEM_O3_2 <- glm(Freq~(GPEM.GENDER+GPEM.PMS+GPEM.EMS+GPEM.M.STATUS)^3-GPEM.GENDER:GPEM.PMS:GPEM.EMS-GPEM.GENDER:GPEM.PMS:GPEM.M.STATUS,data=DF_GPEM,family=poisson)
summary(GPEM_O3_2)
1-pchisq(GPEM_O3_2$deviance,GPEM_O3_2$df.residual)
anova(GPEM_O3_2)

GPEM_O3_3 <- glm(Freq~(GPEM.GENDER+GPEM.PMS+GPEM.EMS+GPEM.M.STATUS)^3-GPEM.GENDER:GPEM.PMS:GPEM.EMS-GPEM.GENDER:GPEM.PMS:GPEM.M.STATUS-GPEM.GENDER:GPEM.EMS:GPEM.M.STATUS,data=DF_GPEM,family=poisson)
summary(GPEM_O3_3)
1-pchisq(GPEM_O3_3$deviance,GPEM_O3_3$df.residual)
anova(GPEM_O3_3)

GPEM_O3_4 <- glm(Freq~(GPEM.GENDER+GPEM.PMS+GPEM.EMS+GPEM.M.STATUS)^3-GPEM.GENDER:GPEM.PMS:GPEM.EMS-GPEM.GENDER:GPEM.PMS:GPEM.M.STATUS-GPEM.GENDER:GPEM.EMS:GPEM.M.STATUS-GPEM.GENDER:GPEM.M.STATUS,data=DF_GPEM,family=poisson)
summary(GPEM_O3_4)
1-pchisq(GPEM_O3_4$deviance,GPEM_O3_4$df.residual)
anova(GPEM_O3_4)

GPEM_O3_5 <- glm(Freq~(GPEM.GENDER+GPEM.PMS+GPEM.EMS+GPEM.M.STATUS)^3-GPEM.GENDER:GPEM.PMS:GPEM.EMS-GPEM.GENDER:GPEM.PMS:GPEM.M.STATUS-GPEM.GENDER:GPEM.EMS:GPEM.M.STATUS-GPEM.GENDER:GPEM.M.STATUS-GPEM.GENDER:GPEM.EMS,data=DF_GPEM,family=poisson)
summary(GPEM_O3_5)
1-pchisq(GPEM_O3_5$deviance,GPEM_O3_5$df.residual)
anova(GPEM_O3_5)

GPEM_O3_6 <- glm(Freq~(GPEM.GENDER+GPEM.PMS+GPEM.EMS+GPEM.M.STATUS)^3-GPEM.GENDER:GPEM.PMS:GPEM.EMS-GPEM.GENDER:GPEM.PMS:GPEM.M.STATUS-GPEM.GENDER:GPEM.EMS:GPEM.M.STATUS-GPEM.GENDER:GPEM.M.STATUS-GPEM.GENDER:GPEM.EMS-GPEM.PMS:GPEM.EMS:GPEM.M.STATUS,data=DF_GPEM,family=poisson)
summary(GPEM_O3_6)
1-pchisq(GPEM_O3_6$deviance,GPEM_O3_6$df.residual)
anova(GPEM_O3_6)


# Procédure backward partant du modèle GPEM_O3 et fondée sur BIC 
n <- sum(GPEM_O3_glm$weights)
GPEM_BIC <- step(GPEM_O3_glm,direction="backward",k=log(n))
summary(GPEM_BIC)